Merge branch 'hdivdiv_bdelems' into 'master'

Hdivdiv quad elements See merge request !254
NGSolve · Oct 20, 2017 · c638b13 · c638b13
2 parents 3be8385 + 24861ab
commit c638b13
Show file tree

Hide file tree

Showing 2 changed files with 424 additions and 19 deletions.
diff --git a/comp/hdivdivfespace.cpp b/comp/hdivdivfespace.cpp
@@ -525,19 +525,39 @@ namespace ngcomp
             ndof += first_facet_dof[f+1] - first_facet_dof[f];            
         }
         break;
+      case ET_QUAD:
+        //ndof += 2*(oi[0]+2)*(oi[0]+1) +1;
+        ndof += (oi[0]+1+HDivDivFE<ET_QUAD>::incsg)*(oi [0]+1+HDivDivFE<ET_QUAD>::incsg)
+          + (oi[0]+2)*(oi[0])*2
+          + 2*(oi[0]+1+HDivDivFE<ET_QUAD>::incsugv) +1;
+        if(discontinuous)
+        {
+          for (auto f : ma->GetElFacets(ei))
+            ndof += first_facet_dof[f+1] - first_facet_dof[f];            
+        }
+        break;
       case ET_PRISM:
         ndof += 3*(oi[0]+1+incrorder_xx1)*(oi[0]+incrorder_xx1)*(oi[2]+1+incrorder_xx2)/2 + 
           (oi[0]+1+incrorder_zz1)*(oi[0]+2+incrorder_zz1)*(oi[2]-1+incrorder_zz2)/2 + 
           (oi[0]+1)*(oi[0]+2)*(oi[2]+1)/2*2;
         if(discontinuous)
         {
-          /*
+          for (auto f : ma->GetElFacets(ei))
-          auto fnums = ma->GetElFacets(ei);
+            ndof += first_facet_dof[f+1] - first_facet_dof[f];            
-          for(int ii=0; ii<fnums.Size(); ii++)
+        }
-          {
+        break;
-            ndof += first_facet_dof[fnums[ii]+1] - first_facet_dof[fnums[ii]];
+      case ET_HEX:
-          }
+        ndof += 3*(oi[0]+2)*(oi[0])*(oi[0]+2) + 3*(oi[0]+1)*(oi[0]+2)*(oi[0]+1);
-          */
+        if(discontinuous)
+        {
+          for (auto f : ma->GetElFacets(ei))
+            ndof += first_facet_dof[f+1] - first_facet_dof[f];            
+        }
+        break;
+      case ET_TET:
+        ndof += (oi[0]+1)*(oi[0]+2)*(oi[0]+1);
+        if(discontinuous)
+        {
           for (auto f : ma->GetElFacets(ei))
             ndof += first_facet_dof[f+1] - first_facet_dof[f];            
         }
@@ -650,6 +670,17 @@ namespace ngcomp
       fe->ComputeNDof();
       return *fe;
     }
+    case ET_QUAD:
+    {
+      auto fe = new (alloc) HDivDivFE<ET_QUAD> (order,plus);
+      fe->SetVertexNumbers (ngel.Vertices());
+      int ii = 0;
+      for(auto f : ngel.Facets())
+        fe->SetOrderFacet(ii++,order_facet[f]);
+      fe->SetOrderInner(order_inner[ei.Nr()]);
+      fe->ComputeNDof();
+      return *fe;
+    }
     case ET_PRISM:
     {
       auto fe = new (alloc) HDivDivFE<ET_PRISM> (order,plus);
@@ -661,6 +692,28 @@ namespace ngcomp
       fe->ComputeNDof();
       return *fe;
     }
+    case ET_HEX:
+    {
+      auto fe = new (alloc) HDivDivFE<ET_HEX> (order,plus);
+      fe->SetVertexNumbers (ngel.vertices);
+      int ii = 0;
+      for(auto f : ngel.Facets())
+        fe->SetOrderFacet(ii++,order_facet[f]);
+      fe->SetOrderInner(order_inner[ei.Nr()]);
+      fe->ComputeNDof();
+      return *fe;
+    }
+    case ET_TET:
+    {
+      auto fe = new (alloc) HDivDivFE<ET_TET> (order,plus);
+      fe->SetVertexNumbers (ngel.vertices);
+      int ii = 0;
+      for(auto f : ngel.Facets())
+        fe->SetOrderFacet(ii++,order_facet[f]);
+      fe->SetOrderInner(order_inner[ei.Nr()]);
+      fe->ComputeNDof();
+      return *fe;
+    }
     default:
       throw Exception(string("HDivDivFESpace::GetFE: element-type ") +
         ToString(ngel.GetType()) + " not supported");