Merge branch 'hcurl_dual_3d' into 'master'

Hcurl dual 3d See merge request jschoeberl/ngsolve!311
NGSolve · Mar 8, 2018 · cc0972260c84f0b9527145a8fee99f39235919aa · cc09722
2 parents d07789c + 62dc596
commit cc0972260c84f0b9527145a8fee99f39235919aa
Unified Split

Showing with 202 additions and 27 deletions.

+202 −27 fem/hcurlhofe_impl.hpp
diff --git a/fem/hcurlhofe_impl.hpp b/fem/hcurlhofe_impl.hpp
@@ -1222,12 +1222,12 @@ namespace ngfem
        for (int i = 0; i < 3; i++)
          ii += order_edge[i];

-        // now come the inner ...
+        /*// now come the inner ...
        Vec<2,AutoDiff<2,T>> adp(mip);
-        /*
-        AutoDiff<2> adx(x,0);
-        AutoDiff<2> ady(y,1);
-        */
+        
+        //AutoDiff<2> adx(x,0);
+        //AutoDiff<2> ady(y,1);
+        
        AutoDiff<2,T> adx = adp(0);
        AutoDiff<2,T> ady = adp(1);
        
@@ -1244,6 +1244,23 @@ namespace ngfem
        for (int i = 1; i <= p; i++)
          for (int j = 1; j <= p-i; j++)
            shape[ii++] = Vec<2,T> (THDiv2Shape<2,T> (uDv_minus_vDu (adpol1[i], adpol2[j])));
+	
+	*/
+	
+	auto xphys = mip.GetPoint()(0);
+        auto yphys = mip.GetPoint()(1);
+        DubinerBasis3::Eval(order-2, x, y,
+                            SBLambda([&] (size_t nr, auto val)
+                                     {
+				       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<2,T> (val, 0);
+                                       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<2,T> (val*x, val*y);
+                                     }));
+	LegendrePolynomial::Eval(order-2,x,
+				 SBLambda([&] (size_t nr, auto val)
+					  {
+					    shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<2,T>(0,val);
+					  }));
+	
      }
  }

@@ -1293,29 +1310,187 @@ namespace ngfem
          ii += order_edge[i];
      }
    if (ip.VB() == BND)
-      { // inner shapes
-        // now come the inner ...
-        Vec<3,AutoDiff<3,T>> adp(mip);
-
-        AutoDiff<3,T> adx = adp(0);
-        AutoDiff<3,T> ady = adp(1);
-        AutoDiff<3,T> adz = adp(1);
-        
-        AutoDiff<3,T> l3 = 1-adx-ady-adz;
+      {
+	T x2, y2;
+	for (int f = 0; f < 4; f++)
+	  {
+	    int p = order_face[f][0];
+	    if (f == facetnr)
+	      {
+		//INT<4> fav = ET_T::GetFaceSort (f, vnums);
+		switch(facetnr)
+		  {
+		  case 0:
+		    x2 = z;//lam[fav[2]];
+		    y2 = y;//lam[fav[1]];
+		    break;
+		  case 1:
+		    x2 = x;//lam[fav[0]];
+		    y2 = z;//lam[fav[2]];
+		    break;
+		  case 2:
+		    x2 = x;//lam[fav[0]];
+		    y2 = y;//lam[fav[1]];
+		    break;
+		  case 3:
+		    x2 = x;//lam[fav[0]];
+		    y2 = y;//lam[fav[1]];
+		  default:
+		    break;
+		  }
+		Vec<2,AutoDiff<2,T>> adp;
+		adp[0] = AutoDiff<2,T>(x2,0);
+		adp[1] = AutoDiff<2,T>(y2,1);
+		AutoDiff<2,T> adx = adp(0);
+		AutoDiff<2,T> ady = adp(1);
+		
+		AutoDiff<2,T> l2 = 1-adx-ady;

-        ArrayMem<AutoDiff<3,T>, 20> adpol1(order+1), adpol2(order+1);
-        LegendrePolynomial::EvalScaled(order, adx-l3, adx+l3, adpol1);
-        LegendrePolynomial::Eval(order, 2*ady-1, adpol2);
-        int p = order_face[0][0];
-        /*
-        for (int i = 0; i < p; i++)
-          for (int j = 0; j < p-i; j++)
-            if (i > 0 || j > 0)
-              shape.Row(ii++) = Vec<3> (THDiv2Shape<3> (Du (adpol1[i]*adpol2[j])));
-        for (int i = 1; i <= p; i++)
-          for (int j = 1; j <= p-i; j++)
-            shape.Row(ii++) = Vec<3> (THDiv2Shape<3> (uDv_minus_vDu (adpol1[i], adpol2[j])));
-        */
+		ArrayMem<AutoDiff<2,T>, 20> adpol1(order+1), adpol2(order+1);
+		LegendrePolynomial::EvalScaled(order, adx-l2, adx+l2, adpol1);
+		LegendrePolynomial::Eval(order, 2*ady-1, adpol2);
+		
+		for (int i = 0; i < p; i++)
+		  for (int j = 0; j < p-i; j++)
+		    if (i > 0 || j > 0)
+		      {
+			auto vec1 = Vec<2,T> (THDiv2Shape<2,T> (Du (adpol1[i]*adpol2[j])));
+			switch(facetnr)
+			  {
+			  case 0:
+			    shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(0.0,vec1(1),vec1(0));
+			    break;
+			  case 1:
+			    shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec1(0),0.0,vec1(1));
+			    break;
+			  case 2:
+			    shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec1(0),vec1(1),0.0);
+			    break;
+			  case 3:
+			    shape[ii++] = 1/sqrt(3.0)*1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec1(0),vec1(1),-vec1(0)-vec1(1));
+			    break;
+			  default:
+			    break;
+			  }
+		      }
+		for (int i = 1; i <= p; i++)
+		  for (int j = 1; j <= p-i; j++)
+		    {
+		      auto vec2 = Vec<2,T> (THDiv2Shape<2,T> (uDv_minus_vDu (adpol1[i], adpol2[j])));
+		      switch(facetnr)
+			{
+			case 0:
+			  shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(0.0,vec2(1),vec2(0));
+			  break;
+			case 1:
+			  shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec2(0),0.0,vec2(1));
+			  break;
+			case 2:
+			  shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec2(0),vec2(1),0.0);
+			  break;
+			case 3:
+			  shape[ii++] = 1/sqrt(3.0)*1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(vec2(0),vec2(1),-vec2(0)-vec2(1));
+			  break;
+			default:
+			  break;
+			}
+		    }
+		/*
+		DubinerBasis3::Eval(order-2, x2, y2,
+				    SBLambda([&] (size_t nr, auto val)
+					     {
+					       switch(facetnr)
+						 {
+						 case 0:
+						   shape[ii++] = Vec<3,T> (0, 0, val);
+						   shape[ii++] = Vec<3,T> (0, val*y2, val*x2);
+						   break;
+						 case 1:
+						   shape[ii++] = Vec<3,T> (val, 0, 0);
+						   shape[ii++] = Vec<3,T> (val*x2, 0, val*y2);
+						   break;
+						 case 2:
+						   shape[ii++] = Vec<3,T> (val, 0, 0);
+						   shape[ii++] = Vec<3,T> (val*x2, val*y2, 0);
+						   break;
+						 case 3:
+						   shape[ii++] = 1/sqrt(3)*Vec<3,T> (val, 0, -val);
+						   shape[ii++] = 1/sqrt(3)*Vec<3,T> (val*x2, val*y2, -val*x2-val*y2);
+						   
+						   break;
+						 }
+					     }));
+		LegendrePolynomial::Eval(order-2,x2,
+					 SBLambda([&] (size_t nr, auto val)
+						  {
+						    switch(facetnr)
+						      {
+						      case 0:
+							shape[ii++] = Vec<3,T>(0, val, 0);
+							break;
+						      case 1:
+							shape[ii++] = Vec<3,T>(0, 0, val);
+							break;
+						      case 2:
+							shape[ii++] = Vec<3,T>(0, val, 0);
+							break;
+						      case 3:
+							shape[ii++] = 1/sqrt(3)*Vec<3,T>(0, val, -val);
+							break;
+						      }
+						  }));
+		*/
+	      }
+	    else
+	      ii += (p+1)*(p-1);
+	  }
+      }
+    else
+      {
+        for (int i = 0; i < 4; i++)
+	  {
+	    int p = order_face[i][0];
+	    ii += (p+1)*(p-1);
+	  }
+      }
+    if (ip.VB() == VOL)
+      {
+	auto xphys = mip.GetPoint()(0);
+        auto yphys = mip.GetPoint()(1);
+	auto zphys = mip.GetPoint()(2);
+	
+	LegendrePolynomial leg;
+	JacobiPolynomialAlpha jac1(1);    
+	leg.EvalScaled1Assign 
+	  (order-3, lam[2]-lam[3], lam[2]+lam[3],
+	   SBLambda ([&](size_t k, T polz) LAMBDA_INLINE
+		     {
+		       // JacobiPolynomialAlpha jac(2*k+1);
+		       JacobiPolynomialAlpha jac2(2*k+2);
+		       
+		       jac1.EvalScaledMult1Assign
+			 (order-3-k, lam[1]-lam[2]-lam[3], 1-lam[0], polz, 
+			  SBLambda ([&] (size_t j, T polsy) LAMBDA_INLINE
+				    {
+				      // JacobiPolynomialAlpha jac(2*(j+k)+2);
+				      jac2.EvalMult(order-3 - k - j, 2 * lam[0] - 1, polsy, 
+						    SBLambda([&](size_t j, T val) LAMBDA_INLINE
+							     {
+							       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(val*x, val*y, val*zphys);
+							       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(val, 0, 0);
+							       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T>(0, val, 0);
+							     }));
+				      jac2.IncAlpha2();
+				    }));
+		       jac1.IncAlpha2();
+		     }));
+
+
+        DubinerBasis3::Eval(order-3, x, y,
+                            SBLambda([&] (size_t nr, auto val)
+                                     {
+				       shape[ii++] = 1/mip.GetMeasure()*mip.GetJacobian()*Vec<3,T> (0, 0, val);
+                                     }));
      }
  }