|  | @@ -2470,14 +2470,14 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |        };
 | 
	
		
			
				|  |  |        compute_norm_area_elem(S.GetElemList(), normal, area_elem);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -      S.quadrature_BS  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -      S.quadrature_dBS .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -      S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU   , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | -      S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU   , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | -      S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | -      S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF  , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | -      s.quadrature_dUxD .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -      s.quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | +      S.quadrature_BS   .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | +      S.quadrature_dBS  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | +      S.quadrature_FxU  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU   , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | +      S.quadrature_DxU  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU   , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | +      S.quadrature_FxdU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | +      S.quadrature_dUxF .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF  , order_singular, order_direct, -1.0, comm);
 | 
	
		
			
				|  |  | +      S.quadrature_dUxD .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD  , order_singular, order_direct, -1.0, comm,  0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | +      S.quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |        auto compute_B0_deprecated = [&S](const Real alpha) { // alpha/|r| \hat{\theta}
 | 
	
		
			
				|  |  |          const Vector<ElemBasis> X = S.GetElemList().ElemVector();
 | 
	
	
		
			
				|  | @@ -2949,9 +2949,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |            // dg_dnu1 = (2 B) \cdot (n \cdnot \nabla) \nabla G[sigma]
 | 
	
		
			
				|  |  |            Vector<ElemBasis> d2Gsigma;
 | 
	
		
			
				|  |  | -          Quadrature<Real> quadrature_Fxd2U;
 | 
	
		
			
				|  |  | -          quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -          quadrature_Fxd2U.Eval(d2Gsigma, S.GetElemList(), sigma, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  | +          S.quadrature_Fxd2U.Eval(d2Gsigma, S.GetElemList(), sigma, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  |            for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |              for (Long j = 0; j < Nnodes; j++) {
 | 
	
		
			
				|  |  |                dg_dnu1[i][j] = 0;
 | 
	
	
		
			
				|  | @@ -2990,9 +2988,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |            // dg_dnu3 = (sigma (\nabla D)^T [2 B]
 | 
	
		
			
				|  |  |            Vector<ElemBasis> nablaDtv;
 | 
	
		
			
				|  |  | -          Quadrature<Real> quadrature_dUxD;
 | 
	
		
			
				|  |  | -          quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -          quadrature_dUxD.Eval(nablaDtv, S.GetElemList(), v, S.Laplace_dUxD);
 | 
	
		
			
				|  |  | +          S.quadrature_dUxD.Eval(nablaDtv, S.GetElemList(), v, S.Laplace_dUxD);
 | 
	
		
			
				|  |  |            for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |              for (Long j = 0; j < Nnodes; j++) {
 | 
	
		
			
				|  |  |                dg_dnu3[i][j] = 0;
 | 
	
	
		
			
				|  | @@ -3212,9 +3208,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |            // dAdnu2 = (u n) \cdot (n \cdnot \nabla) \nabla G[v]
 | 
	
		
			
				|  |  |            Vector<ElemBasis> d2G_v;
 | 
	
		
			
				|  |  | -          Quadrature<Real> quadrature_Fxd2U;
 | 
	
		
			
				|  |  | -          quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -          quadrature_Fxd2U.Eval(d2G_v, S.GetElemList(), v, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  | +          S.quadrature_Fxd2U.Eval(d2G_v, S.GetElemList(), v, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  |            for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |              for (Long j = 0; j < Nnodes; j++) {
 | 
	
		
			
				|  |  |                dAdnu2[i][j] = 0;
 | 
	
	
		
			
				|  | @@ -3234,9 +3228,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |            // dAdnu3 = (v n \cdot \nabla D[u]
 | 
	
		
			
				|  |  |            Vector<ElemBasis> nablaDt_u_n;
 | 
	
		
			
				|  |  | -          Quadrature<Real> quadrature_dUxD;
 | 
	
		
			
				|  |  | -          quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -          quadrature_dUxD.Eval(nablaDt_u_n, S.GetElemList(), u_n, S.Laplace_dUxD);
 | 
	
		
			
				|  |  | +          S.quadrature_dUxD.Eval(nablaDt_u_n, S.GetElemList(), u_n, S.Laplace_dUxD);
 | 
	
		
			
				|  |  |            for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |              for (Long j = 0; j < Nnodes; j++) {
 | 
	
		
			
				|  |  |                dAdnu3[i][j] = 0;
 | 
	
	
		
			
				|  | @@ -3620,9 +3612,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |              Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |              Vector<ElemBasis> gradB;
 | 
	
		
			
				|  |  | -            Quadrature<Real> quadrature_Fxd2U;
 | 
	
		
			
				|  |  | -            quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -            quadrature_Fxd2U.Eval(gradB, S.GetElemList(), sigma, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  | +            S.quadrature_Fxd2U.Eval(gradB, S.GetElemList(), sigma, S.Laplace_Fxd2U);
 | 
	
		
			
				|  |  |              Vector<ElemBasis> dphi_dnu(Nelem);
 | 
	
		
			
				|  |  |              for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |                for (Long j = 0; j < Nnodes; j++) {
 | 
	
	
		
			
				|  | @@ -3752,9 +3742,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |              Vector<ElemBasis> dphi_dnu(Nelem);
 | 
	
		
			
				|  |  |              Vector<ElemBasis> nablaDt_nxGv;
 | 
	
		
			
				|  |  | -            Quadrature<Real> quadrature_dUxD;
 | 
	
		
			
				|  |  | -            quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
 | 
	
		
			
				|  |  | -            quadrature_dUxD.Eval(nablaDt_nxGv, S.GetElemList(), nxGv, S.Laplace_dUxD);
 | 
	
		
			
				|  |  | +            S.quadrature_dUxD.Eval(nablaDt_nxGv, S.GetElemList(), nxGv, S.Laplace_dUxD);
 | 
	
		
			
				|  |  |              for (Long i = 0; i < Nelem; i++) {
 | 
	
		
			
				|  |  |                for (Long j = 0; j < Nnodes; j++) {
 | 
	
		
			
				|  |  |                  dphi_dnu[i][j] = 0;
 |