.

include/sctl/boundary_quadrature.hpp

@@ -2518,364 +2518,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
       SetupQuadrature(S.quadrature_dUxD , S, S.Laplace_dUxD  , order_singular, order_direct, -1.0, comm,  0.01 * pow<-2,Real>(ORDER));
       SetupQuadrature(S.quadrature_Fxd2U, S, S.Laplace_Fxd2U , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 
-      {
-      auto compute_A = [&S,&Jt,&Jp] (const Vector<Real>& x) {
-        Vector<ElemBasis> normal, area_elem;
-        compute_norm_area_elem(S, normal, area_elem);
-
-        auto compute_B0 = [&S, &Jp](const Real alpha) {
-          Vector<ElemBasis> B0;
-          S.quadrature_BS.Eval(B0, S.GetElemList(), Jp, S.BiotSavart);
-          return B0*alpha;
-        };
-        auto compute_half_n_plus_dG = [&S, &normal](const Vector<ElemBasis>& sigma) { // B = n sigma/2 + dG[sigma]
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> B;
-          S.quadrature_FxdU.Eval(B, S.GetElemList(), sigma, S.Laplace_FxdU);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              for (Long k = 0; k < COORD_DIM; k++) {
-                B[i*COORD_DIM+k][j] -= 0.5*sigma[i][j]*normal[i*COORD_DIM+k][j];
-              }
-            }
-          }
-          return B;
-        };
-
-        auto compute_A11 = [&S,&normal,&compute_half_n_plus_dG](Vector<Real>& B_dot_n, const Vector<Real>& sigma) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          B_dot_n.ReInit(Nelem * Nnodes);
-          Vector<ElemBasis> sigma_(Nelem);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              sigma_[i][j] = sigma[i*Nnodes+j];
-            }
-          }
-
-          Vector<ElemBasis> B_dot_n_ = compute_dot_prod(normal, compute_half_n_plus_dG(sigma_));
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              B_dot_n[i*Nnodes+j] = B_dot_n_[i][j];
-            }
-          }
-        };
-        auto compute_A12 = [&S,&normal,&compute_B0](Vector<Real>& B_dot_n, const Real alpha) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          B_dot_n.ReInit(Nelem * Nnodes);
-          Vector<ElemBasis> B_dot_n_ = compute_dot_prod(normal, compute_B0(alpha));
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              B_dot_n[i*Nnodes+j] = B_dot_n_[i][j];
-            }
-          }
-        };
-        auto compute_A21 = [&S,&normal,&compute_half_n_plus_dG](const Vector<Real>& sigma) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> sigma_(Nelem);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              sigma_[i][j] = sigma[i*Nnodes+j];
-            }
-          }
-
-          Vector<ElemBasis> B = compute_half_n_plus_dG(sigma_);
-
-          Vector<ElemBasis> J(Nelem * COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) { // Set J
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> b, n;
-              b(0) = B[i*COORD_DIM+0][j];
-              b(1) = B[i*COORD_DIM+1][j];
-              b(2) = B[i*COORD_DIM+2][j];
-
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
-
-              J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
-              J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
-              J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
-            }
-          }
-
-          Vector<ElemBasis> A;
-          S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
-
-          Real pol_circ = 0;
-          { // compute pol_circ
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            const auto& quad_wts = ElemBasis::QuadWts();
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                pol_circ += A[i*COORD_DIM+0][j] * dX[i*COORD_DIM*2+1][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+1][j] * dX[i*COORD_DIM*2+3][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+2][j] * dX[i*COORD_DIM*2+5][j] * quad_wts[j] / S.NtNp_[0];
-              }
-            }
-          }
-          return pol_circ;
-        };
-        auto compute_A22 = [&S,&compute_B0,&normal](const Real alpha) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> B = compute_B0(alpha);
-
-          Vector<ElemBasis> J(Nelem * COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) { // Set J
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> b, n;
-              b(0) = B[i*COORD_DIM+0][j];
-              b(1) = B[i*COORD_DIM+1][j];
-              b(2) = B[i*COORD_DIM+2][j];
-
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
-
-              J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
-              J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
-              J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
-            }
-          }
-
-          Vector<ElemBasis> A;
-          S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
-
-          Real pol_circ = 0;
-          { // compute pol_circ
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            const auto& quad_wts = ElemBasis::QuadWts();
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                pol_circ += A[i*COORD_DIM+0][j] * dX[i*COORD_DIM*2+1][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+1][j] * dX[i*COORD_DIM*2+3][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+2][j] * dX[i*COORD_DIM*2+5][j] * quad_wts[j] / S.NtNp_[0];
-              }
-            }
-          }
-          return pol_circ;
-        };
-
-        const Vector<Real> sigma(x.Dim()-1,(Iterator<Real>)x.begin(),false);
-        const Real& alpha = x[x.Dim()-1];
-
-        Vector<Real> Ax;
-        Ax.ReInit(x.Dim());
-        Vector<Real> Bdotn(x.Dim()-1,Ax.begin(),false);
-        Real& flux = Ax[x.Dim()-1];
-
-        Vector<Real> Adotn_0, Adotn_1;
-        compute_A11(Adotn_0, sigma);
-        compute_A12(Adotn_1, alpha);
-        Bdotn = Adotn_0 + Adotn_1;
-
-        flux = compute_A21(sigma) + compute_A22(alpha);
-        return Ax;
-      };
-      auto compute_invA = [&S,&comm,&compute_A] (Vector<ElemBasis>& sigma, Real& alpha, Real flux) {
-        typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_A](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
-          (*Ax) = compute_A(x);
-        };
-        const Long Nelem = S.NElem();
-        const Long Nnodes = ElemBasis::Size();
-
-        Vector<Real> rhs_(Nelem * Nnodes + 1);
-        rhs_ = 0;
-        rhs_[Nelem * Nnodes] = flux;
-
-        Vector<Real> x_(Nelem * Nnodes + 1);
-        x_ = 0;
-        ParallelSolver<Real> linear_solver(comm, true);
-        linear_solver(&x_, BIOp, rhs_, 1e-8, 50);
-
-        sigma.ReInit(Nelem);
-        for (Long i = 0; i < Nelem; i++) {
-          for (Long j = 0; j < Nnodes; j++) {
-            sigma[i][j] = x_[i*Nnodes+j];
-          }
-        }
-        alpha = x_[Nelem * Nnodes];
-      };
-
-      auto compute_Aadj = [&S,&Jt,&Jp] (const Vector<Real>& x) {
-        Vector<ElemBasis> normal, area_elem;
-        compute_norm_area_elem(S, normal, area_elem);
-
-        auto compute_B0 = [&S, &Jp](const Real alpha) {
-          Vector<ElemBasis> B0;
-          S.quadrature_BS.Eval(B0, S.GetElemList(), Jp, S.BiotSavart);
-          return B0*alpha;
-        };
-
-        auto compute_A12 = [&S,&normal,&compute_B0](Vector<Real>& B_dot_n, const Real alpha) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          B_dot_n.ReInit(Nelem * Nnodes);
-          Vector<ElemBasis> B_dot_n_ = compute_dot_prod(normal, compute_B0(alpha));
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              B_dot_n[i*Nnodes+j] = B_dot_n_[i][j];
-            }
-          }
-        };
-        auto compute_A22 = [&S,&compute_B0,&normal](const Real alpha) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> B = compute_B0(alpha);
-
-          Vector<ElemBasis> J(Nelem * COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) { // Set J
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> b, n;
-              b(0) = B[i*COORD_DIM+0][j];
-              b(1) = B[i*COORD_DIM+1][j];
-              b(2) = B[i*COORD_DIM+2][j];
-
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
-
-              J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
-              J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
-              J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
-            }
-          }
-
-          Vector<ElemBasis> A;
-          S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
-
-          Real pol_circ = 0;
-          { // compute pol_circ
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            const auto& quad_wts = ElemBasis::QuadWts();
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                pol_circ += A[i*COORD_DIM+0][j] * dX[i*COORD_DIM*2+1][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+1][j] * dX[i*COORD_DIM*2+3][j] * quad_wts[j] / S.NtNp_[0];
-                pol_circ += A[i*COORD_DIM+2][j] * dX[i*COORD_DIM*2+5][j] * quad_wts[j] / S.NtNp_[0];
-              }
-            }
-          }
-          return pol_circ;
-        };
-        auto compute_A11adj = [&S,&normal](Vector<Real>& U, const Vector<Real>& sigma) { // A11adj  = I/2 + D
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> sigma_n(Nelem*COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              sigma_n[i*COORD_DIM+0][j] = sigma[i*Nnodes+j] * normal[i*COORD_DIM+0][j];
-              sigma_n[i*COORD_DIM+1][j] = sigma[i*Nnodes+j] * normal[i*COORD_DIM+1][j];
-              sigma_n[i*COORD_DIM+2][j] = sigma[i*Nnodes+j] * normal[i*COORD_DIM+2][j];
-            }
-          }
-          S.quadrature_dUxF.Eval(U, S.GetElemList(), sigma_n*(-1), S.Laplace_dUxF);
-          U = sigma*(-0.5) + U;
-        };
-        auto compute_A12adj = [&S,&compute_A12,&area_elem](const Vector<Real>& sigma_) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<Real> A12_sigma_;
-          compute_A12(A12_sigma_, 1);
-
-          Vector<ElemBasis> A12_sigma(Nelem), sigma(Nelem);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              sigma[i][j] = sigma_[i*Nnodes+j];
-              A12_sigma[i][j] = A12_sigma_[i*Nnodes+j];
-            }
-          }
-          return compute_inner_prod(area_elem, A12_sigma, sigma);
-        };
-        auto compute_A21adj = [&S,&area_elem,&normal](Vector<Real>& A21adj_flux, Real flux) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> density(Nelem * COORD_DIM);
-          { // Set density
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
-                density[i*COORD_DIM+0][j] = dX[i*COORD_DIM*2+1][j] * s;
-                density[i*COORD_DIM+1][j] = dX[i*COORD_DIM*2+3][j] * s;
-                density[i*COORD_DIM+2][j] = dX[i*COORD_DIM*2+5][j] * s;
-              }
-            }
-          }
-
-          Vector<ElemBasis> Gdensity;
-          S.quadrature_FxU.Eval(Gdensity, S.GetElemList(), density, S.Laplace_FxU);
-
-          Vector<ElemBasis> nxGdensity(Nelem * COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) { // Set nxGdensity
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> Gdensity_, n;
-              Gdensity_(0) = Gdensity[i*COORD_DIM+0][j];
-              Gdensity_(1) = Gdensity[i*COORD_DIM+1][j];
-              Gdensity_(2) = Gdensity[i*COORD_DIM+2][j];
-
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
-
-              nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1);
-              nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2);
-              nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0);
-            }
-          }
-          S.quadrature_dUxF.Eval(A21adj_flux, S.GetElemList(), nxGdensity, S.Laplace_dUxF);
-          A21adj_flux *= flux;
-        };
-        auto compute_A22adj = [&compute_A22] (const Real alpha) {
-          return compute_A22(alpha);
-        };
-
-        const Vector<Real> sigma(x.Dim()-1,(Iterator<Real>)x.begin(),false);
-        const Real& alpha = x[x.Dim()-1];
-
-        Vector<Real> Ax;
-        Ax.ReInit(x.Dim());
-        Vector<Real> Bdotn(x.Dim()-1,Ax.begin(),false);
-        Real& flux = Ax[x.Dim()-1];
-
-        Vector<Real> Adotn_0, Adotn_1;
-        compute_A11adj(Adotn_0, sigma);
-        compute_A21adj(Adotn_1, alpha);
-        Bdotn = Adotn_0 + Adotn_1;
-
-        flux = compute_A12adj(sigma) + compute_A22adj(alpha);
-        return Ax;
-      };
-      auto compute_invAadj = [&S,&comm,&compute_Aadj] (Vector<Real>& b) {
-        typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_Aadj](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
-          (*Ax) = compute_Aadj(x);
-        };
-        const Long Nelem = S.NElem();
-        const Long Nnodes = ElemBasis::Size();
-
-        Vector<Real> x(b.Dim());
-        x = 0;
-        ParallelSolver<Real> linear_solver(comm, true);
-        linear_solver(&x, BIOp, b, 1e-8, 50);
-        return x;
-      };
-      }
-
       Vector<ElemBasis> Bt0, Bp0;
       Vector<ElemBasis> Jt, Jp;
       { // Set Bt0, Bp0