Merge branch 'hex_vectorfacet' into 'master'

Hex vectorfacet See merge request jschoeberl/ngsolve!332
NGSolve · Mar 22, 2018 · e4202a5 · e4202a5
2 parents b2deb9c + 2eefbff
commit e4202a5
Show file tree

Hide file tree

Showing 2 changed files with 49 additions and 3 deletions.
diff --git a/comp/vectorfacetfespace.cpp b/comp/vectorfacetfespace.cpp
@@ -276,7 +276,8 @@ namespace ngcomp
 		  {
 		  case ET_TET: ndof += 4*(order+1)*2; break;
 		  case ET_PRISM: ndof += 2 * (2*(order+1)+3*(2*order+1)); break;
-		  default: throw Exception (string("VectorFacetFESpace: Element type not implemented"));
+		  case ET_HEX: ndof += 6*(2*order+1)*2; break;
+		  default: throw Exception (string("VectorFacetFESpace: Element type not implemented"));		  
 		  }
 	      }
 	    first_inner_dof[ne] = ndof;

diff --git a/fem/vectorfacetfe.cpp b/fem/vectorfacetfe.cpp
@@ -580,9 +580,54 @@ namespace ngfem
 
   template<>
   void VectorFacetVolumeFE<ET_HEX> ::
-  CalcShape ( const IntegrationPoint & ip, int facet, SliceMatrix<> shape ) const
+  CalcShape ( const IntegrationPoint & ip, int fanr, SliceMatrix<> shape ) const
   {
-    throw Exception("VectorFacetVolumeHex::CalcShape not implemented!");
+    AutoDiff<3> x(ip(0), 0), y(ip(1),1), z(ip(2),2);
+
+    // vertex numbering on HEX:
+    // { 0, 0, 0 },
+    // { 1, 0, 0 },
+    // { 1, 1, 0 },
+    // { 0, 1, 0 },
+    // { 0, 0, 1 },
+    // { 1, 0, 1 },
+    // { 1, 1, 1 },
+    // { 0, 1, 1 }
+
+    AutoDiff<3> mux[8]  = { 1-x,   x,   x, 1-x, 1-x,   x,   x, 1-x };  
+    AutoDiff<3> muy[8]  = { 1-y, 1-y,   y,   y, 1-y, 1-y,   y,   y }; 
+    AutoDiff<3> muz[8]  = { 1-z, 1-z, 1-z, 1-z,   z,   z,   z,   z };    
+
+    AutoDiff<3> sigma[8];
+    for (int i = 0; i < 8; i++) sigma[i] = mux[i] + muy[i] + muz[i];
+
+    shape = 0.0;
+    {
+      int p = facet_order[fanr][0];
+
+      INT<4> f = ET_trait<ET_HEX>::GetFaceSort (fanr, vnums);     
+      AutoDiff<3> xi  = sigma[f[0]] - sigma[f[1]]; 
+      AutoDiff<3> zeta = sigma[f[0]] - sigma[f[3]];   
+
+      ArrayMem< double, 10> polx(p+1), poly(p+1);
+      LegendrePolynomial (p, xi.Value(), polx);
+      LegendrePolynomial (p, zeta.Value(), poly);   
+
+      int ii = first_facet_dof[fanr];
+      for (int i = 0; i <= p; i++)
+        for (int j = 0; j <= p; j++)
+        {
+          double val = polx[i] * poly[j];
+          shape(ii,0) = val * xi.DValue(0);
+          shape(ii,1) = val * xi.DValue(1);
+          shape(ii,2) = val * xi.DValue(2);
+          ii++;
+          shape(ii,0) = val * zeta.DValue(0);
+          shape(ii,1) = val * zeta.DValue(1);
+          shape(ii,2) = val * zeta.DValue(2);
+          ii++;
+        } 
+    }
   }
 
   template<>