.

include/sctl/boundary_quadrature.hpp

@@ -3509,8 +3509,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
         Vector<Long> NtNp;
         NtNp.PushBack(20);
         NtNp.PushBack(4);
-        NtNp.PushBack(20);
-        NtNp.PushBack(4);
+        //NtNp.PushBack(20);
+        //NtNp.PushBack(4);
         S = Stellarator<Real,ORDER>(NtNp);
       }
       if (S.Nsurf() == 1) flux_pol = 0.0;
@@ -3519,7 +3519,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       { // Set pressure
         Vector<ElemBasis> normal, area_elem;
         compute_norm_area_elem(S, normal, area_elem);
-        pressure = area_elem;
+        pressure = area_elem*0;
       }
 
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -3771,185 +3771,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
       }
     }
 
-    static void test_() {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
-      Comm comm = Comm::World();
-      Profile::Enable(true);
-
-      Real flux_tor = 1.0, flux_pol = 1.0;
-      Stellarator<Real,ORDER> S;
-      { // Init S
-        Vector<Long> NtNp;
-        NtNp.PushBack(20);
-        NtNp.PushBack(4);
-        NtNp.PushBack(20);
-        NtNp.PushBack(4);
-        S = Stellarator<Real,ORDER>(NtNp);
-      }
-
-      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-      SetupQuadrature(S.quadrature_BS  , S, S.BiotSavart  , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
-      SetupQuadrature(S.quadrature_FxU , S, S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
-      SetupQuadrature(S.quadrature_FxdU, S, S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
-
-      Vector<ElemBasis> Jt, Jp, Bt0, Bp0;
-      compute_harmonic_vector_potentials(Jt, Jp, S);
-      EvalQuadrature(Bt0, S.quadrature_BS, S, Jp, S.BiotSavart);
-      EvalQuadrature(Bp0, S.quadrature_BS, S, Jt, S.BiotSavart);
-
-      auto compute_B = [&S,&Bt0,&Bp0] (const Vector<ElemBasis>& sigma, Real alpha, Real beta) {
-        const Long Nelem = S.NElem();
-        Vector<ElemBasis> B(S.NElem() * COORD_DIM);
-        if (sigma.Dim()) {
-          const Long Nnodes = ElemBasis::Size();
-          Vector<ElemBasis> normal, area_elem;
-          compute_norm_area_elem(S, normal, area_elem);
-          EvalQuadrature(B, S.quadrature_FxdU, S, 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];
-              }
-            }
-          }
-        } else {
-          B = 0;
-        }
-        if (alpha != 0) B += Bt0*alpha;
-        if (beta  != 0) B += Bp0*beta;
-        return B;
-      };
-      auto compute_A = [&S,compute_B] (const Vector<Real>& x) {
-        const Long Nelem = S.NElem();
-        const Long Nnodes = ElemBasis::Size();
-        SCTL_ASSERT(x.Dim() == Nelem*Nnodes+2);
-
-        Vector<ElemBasis> normal, area_elem;
-        compute_norm_area_elem(S, normal, area_elem);
-
-        Vector<ElemBasis> sigma(Nelem);
-        for (Long i = 0; i < Nelem; i++) {
-          for (Long j = 0; j < Nnodes; j++) {
-            sigma[i][j] = x[i*Nnodes+j];
-          }
-        }
-        Real alpha = x[Nelem*Nnodes + 0];
-        Real beta  = x[Nelem*Nnodes + 1];
-
-        Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
-        Vector<ElemBasis> BdotN = compute_dot_prod(B, normal);
-
-        Real flux_tor, flux_pol;
-        { // Compute flux_tor, flux_pol
-          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;
-          EvalQuadrature(A, S.quadrature_FxU, S, J, S.Laplace_FxU);
-
-          Vector<Real> circ_pol(S.Nsurf()), circ_tor(S.Nsurf());
-          { // compute circ
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            const auto& quad_wts = ElemBasis::QuadWts();
-            for (Long k = 0; k < S.Nsurf(); k++) {
-              circ_pol[k] = 0;
-              circ_tor[k] = 0;
-              Long Ndsp = S.ElemDsp(k);
-              for (Long i = 0; i < S.NTor(k)*S.NPol(k); i++) {
-                for (Long j = 0; j < Nnodes; j++) {
-                  circ_pol[k] += A[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+1][j] * quad_wts[j] / S.NTor(k);
-                  circ_pol[k] += A[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+3][j] * quad_wts[j] / S.NTor(k);
-                  circ_pol[k] += A[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+5][j] * quad_wts[j] / S.NTor(k);
-
-                  circ_tor[k] += A[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+0][j] * quad_wts[j] / S.NPol(k);
-                  circ_tor[k] += A[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+2][j] * quad_wts[j] / S.NPol(k);
-                  circ_tor[k] += A[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+4][j] * quad_wts[j] / S.NPol(k);
-                }
-              }
-            }
-          }
-          if (S.Nsurf() == 1) {
-            flux_tor = circ_pol[0];
-            flux_pol = 0;
-          } else if (S.Nsurf() == 2) {
-            flux_tor = circ_pol[1] - circ_pol[0];
-            flux_pol = circ_tor[0] - circ_tor[1];
-          } else {
-            SCTL_ASSERT(false);
-          }
-        }
-
-        Vector<Real> Ax(Nelem*Nnodes+2);
-        for (Long i = 0; i < Nelem; i++) {
-          for (Long j = 0; j < Nnodes; j++) {
-            Ax[i*Nnodes+j] = BdotN[i][j];
-          }
-        }
-        Ax[Nelem*Nnodes + 0] = flux_tor;
-        Ax[Nelem*Nnodes + 1] = flux_pol;
-        return Ax;
-      };
-      auto compute_invA = [&S,&comm,&compute_A] (Vector<ElemBasis>& sigma, Real& alpha, Real& beta, Real flux_tor, Real flux_pol) {
-        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 + 2);
-        rhs_ = 0;
-        rhs_[Nelem * Nnodes + 0] = flux_tor;
-        rhs_[Nelem * Nnodes + 1] = flux_pol;
-
-        Vector<Real> x_(Nelem * Nnodes + 2);
-        x_ = 0;
-        ParallelSolver<Real> linear_solver(comm, true);
-        linear_solver(&x_, BIOp, rhs_, 1e-8, 100);
-
-        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 + 0];
-        beta  = x_[Nelem * Nnodes + 1];
-      };
-
-      Real alpha = 0, beta = 0;
-      Vector<ElemBasis> sigma;
-      compute_invA(sigma, alpha, beta, flux_tor, flux_pol);
-      Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
-      { // Write VTU
-        VTUData vtu;
-        vtu.AddElems(S.GetElemList(), sigma, ORDER);
-        vtu.WriteVTK("sigma", comm);
-      }
-      { // Write VTU
-        VTUData vtu;
-        vtu.AddElems(S.GetElemList(), B, ORDER);
-        vtu.WriteVTK("B", comm);
-      }
-    }
-
   private:
 
     static void tmp() {