Merge branch 'h1_dual_shape_3d' into 'master'

H1 dual shape 3d See merge request jschoeberl/ngsolve!314
NGSolve · Mar 14, 2018 · 77cd7a1f95e5288b804e3e65eb311a022336d036 · 77cd7a1
2 parents 2ca1df1 + 00febef
commit 77cd7a1f95e5288b804e3e65eb311a022336d036
Unified Split

Showing with 97 additions and 36 deletions.

+4 −2 comp/h1hofespace.cpp

+84 −26 fem/h1hofe_impl.hpp

+1 −1 fem/scalarfe.cpp

+1 −1 fem/scalarfe.hpp

+3 −2 fem/symbolicintegrator.cpp

+2 −2 fem/tscalarfe.hpp

+2 −2 fem/tscalarfe_impl.hpp
diff --git a/comp/h1hofespace.cpp b/comp/h1hofespace.cpp
@@ -85,7 +85,7 @@ namespace ngcomp
    static void GenerateMatrix (const AFEL & fel, const MIP & mip,
                                MAT & mat, LocalHeap & lh)
    {
-      static_cast<const ScalarFiniteElement<D>&>(fel).CalcDualShape (mip.IP(), mat.Row(0));
+      static_cast<const ScalarFiniteElement<D>&>(fel).CalcDualShape (mip, mat.Row(0));
    }
    template <typename AFEL, typename MIP, typename MAT,
              typename std::enable_if<!std::is_convertible<MAT,SliceMatrix<double,ColMajor>>::value, int>::type = 0>
@@ -224,7 +224,9 @@ namespace ngcomp
        break;
      case 3:
        additional_evaluators.Set ("hesse", make_shared<T_DifferentialOperator<DiffOpHesse<3>>> ());
-	additional_evaluators.Set ("hesseboundary", make_shared<T_DifferentialOperator<DiffOpHesseBoundary<3>>> ());break;
+	additional_evaluators.Set ("hesseboundary", make_shared<T_DifferentialOperator<DiffOpHesseBoundary<3>>> ());
+	additional_evaluators.Set ("dual", make_shared<T_DifferentialOperator<DiffOpDual<3>>> ());
+	break;
      default:
        ;
      }
diff --git a/fem/h1hofe_impl.hpp b/fem/h1hofe_impl.hpp
@@ -35,7 +35,7 @@ namespace ngfem
    template<typename Tx, typename TFA>  
      INLINE void T_CalcShape (TIP<DIM,Tx> ip, TFA & shape) const;
    
-    void CalcDualShape2 (const IntegrationPoint & ip, SliceVector<> shape) const
+    void CalcDualShape2 (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
    { throw Exception ("dual shape not implemented, H1Ho"); }
    
  };
@@ -107,18 +107,20 @@ namespace ngfem
  }


-  template<> inline void H1HighOrderFE_Shape<ET_TRIG> ::
-  CalcDualShape2 (const IntegrationPoint & ip, SliceVector<> shape) const
+  template<>
+  inline void H1HighOrderFE_Shape<ET_TRIG> ::CalcDualShape2 (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
  {
+    auto & ip = mip.IP();
+    shape = 0.0;
    double lam[3] = { ip(0), ip(1), 1-ip(0)-ip(1) };
+    size_t ii = 3;

    if (ip.VB() == BBND)
-      for (size_t i = 0; i < 3; i++)
-        shape[i] = (i == ip.FacetNr()) ? 1 : 0;
-    else
-      for (size_t i = 0; i < 3; i++)
-        shape[i] = 0;
-    size_t ii = 3;
+      {
+	for (size_t i = 0; i < 3; i++)
+	  shape[i] = (i == ip.FacetNr()) ? 1 : 0;
+      }
+    

    // edge-based shapes
    for (int i = 0; i < N_EDGE; i++)
@@ -130,29 +132,17 @@ namespace ngfem
              EdgeOrthoPol::
                EvalScaledMult (order_edge[i]-2, 
                                lam[e[1]]-lam[e[0]], lam[e[0]]+lam[e[1]], 
-                                lam[e[0]]*lam[e[1]], shape+ii);
-            }
-          else
-            {
-              for (size_t j = 0; j < order_edge[i]-1; j++)
-                shape[ii+j] = 0;
+                                1.0/mip.GetMeasure()*lam[e[0]]*lam[e[1]], shape+ii);
            }
          ii += order_edge[i]-1;
 	}

-
    // inner shapes
-    if (order_face[0][0] >= 3)
+    if (ip.VB() == VOL && order_face[0][0] >= 3)
      {
-        if (ip.VB() == VOL)
-          {
-            INT<4> f = GetVertexOrientedFace (0);
-            DubinerBasis3::Eval (order_face[0][0]-3, 
-                                 lam[f[0]], lam[f[1]], shape+ii);
-          }
-        else
-          for ( ; ii < ndof; ii++)
-            shape[ii] = 0;
+	INT<4> f = GetVertexOrientedFace (0);
+	DubinerBasis3::EvalMult(order_face[0][0]-3, 
+				lam[f[0]], lam[f[1]],1.0/mip.GetMeasure(), shape+ii);
      }
  }

@@ -279,6 +269,74 @@ namespace ngfem
  }


+  template<>
+  inline void H1HighOrderFE_Shape<ET_TET> ::
+  CalcDualShape2 (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
+  {
+    auto & ip = mip.IP();
+    shape = 0.0;
+    double lam[4] = { ip(0), ip(1), ip(2), 1-ip(0)-ip(1)-ip(2) };
+    size_t ii = 4;
+
+    if (ip.VB() == BBBND)
+      {
+	for (size_t i = 0; i < 4; i++)
+	  shape[i] = (i == ip.FacetNr()) ? 1 : 0;
+      }
+    // edge-based shapes
+    for (int i = 0; i < N_EDGE; i++)
+      {
+	if (order_edge[i] >= 2 && ip.FacetNr() == i && ip.VB() == BBND)
+	  {
+	    INT<2> e = GetVertexOrientedEdge(i);
+	    EdgeOrthoPol::
+	      EvalScaledMult (order_edge[i]-2, 
+			      lam[e[1]]-lam[e[0]], lam[e[0]]+lam[e[1]], 
+			      1.0/mip.GetMeasure()*lam[e[0]]*lam[e[1]], shape+ii);
+	  }
+	
+	ii += order_edge[i]-1;
+      }
+    // face shapes
+    for (int i = 0; i < N_FACE; i++)
+      {
+	if (order_face[i][0] >= 3 && ip.FacetNr() == i && ip.VB() == BND)
+	  {
+	    INT<4> f = GetVertexOrientedFace (0);
+	    DubinerBasis3::EvalMult (order_face[0][0]-3, 
+				     lam[f[0]], lam[f[1]], 1.0/mip.GetMeasure(), shape+ii);
+	  }
+	ii += (order_face[i][0]-2)*(order_face[i][0]-1)/2;
+      }
+    //inner shapes
+    if (ip.VB() == VOL && order_cell[0][0] >= 4)
+      {
+	int p = order_cell[0][0] - 4;
+	LegendrePolynomial leg;
+	JacobiPolynomialAlpha jac1(1);    
+	leg.EvalScaled1Assign 
+	  (p, lam[2]-lam[3], lam[2]+lam[3],
+	   SBLambda ([&](size_t k, double polz) LAMBDA_INLINE
+		     {
+		       // JacobiPolynomialAlpha jac(2*k+1);
+		       JacobiPolynomialAlpha jac2(2*k+2);
+		       
+		       jac1.EvalScaledMult1Assign
+			 (p-k, lam[1]-lam[2]-lam[3], 1-lam[0], polz, 
+			  SBLambda ([&] (size_t j, double polsy) LAMBDA_INLINE
+				    {
+				      // JacobiPolynomialAlpha jac(2*(j+k)+2);
+				      jac2.EvalMult(p - k - j, 2 * lam[0] - 1, polsy, 
+						    SBLambda([&](size_t j, double val) LAMBDA_INLINE
+							     {
+							       shape[ii++] = 1.0/mip.GetMeasure()*val;
+							     }));
+				      jac2.IncAlpha2();
+				    }));
+		       jac1.IncAlpha2();
+		     }));
+      }
+  }


  /* *********************** Prism  **********************/
diff --git a/fem/scalarfe.cpp b/fem/scalarfe.cpp
@@ -284,7 +284,7 @@ namespace ngfem


  void BaseScalarFiniteElement ::
-  CalcDualShape (const IntegrationPoint & ip, SliceVector<> shape) const
+  CalcDualShape (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
  {
    throw Exception (string("CalcDualShape not overloaded for element ") + typeid(*this).name());
  }
diff --git a/fem/scalarfe.hpp b/fem/scalarfe.hpp
@@ -98,7 +98,7 @@ namespace ngfem
                                                 BareSliceVector<> coefs) const;


-    NGS_DLL_HEADER virtual void CalcDualShape (const IntegrationPoint & ip, SliceVector<> shape) const;
+    NGS_DLL_HEADER virtual void CalcDualShape (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const;
  };

  /**
diff --git a/fem/symbolicintegrator.cpp b/fem/symbolicintegrator.cpp
@@ -1639,8 +1639,9 @@ namespace ngfem

                    // tb.Start();
                    BaseMappedIntegrationRule & bmir = mir.Range(i, i+bs, lh);
+
                    proxy1->Evaluator()->CalcMatrix(fel_trial, bmir, Trans(bbmat1), lh);
-                    
+
                    if (!samediffop)
                      proxy2->Evaluator()->CalcMatrix(fel_test, bmir, Trans(bbmat2), lh);
                    // tb.Stop();
@@ -1794,7 +1795,7 @@ namespace ngfem
          ngfem::ELEMENT_TYPE etfacet = transform.FacetType (k);
          const IntegrationRule& ir_facet = GetIntegrationRule(etfacet, fel_trial.Order()+fel_test.Order());
          IntegrationRule & ir_facet_vol = transform(k, ir_facet, lh);
-          cout << "ir_facet_vol = " << ir_facet_vol << ", facetnr = " << ir_facet_vol[0].FacetNr() << endl;
+          
          BaseMappedIntegrationRule & mir = trafo(ir_facet_vol, lh);
          
          ProxyUserData ud;
diff --git a/fem/tscalarfe.hpp b/fem/tscalarfe.hpp
@@ -124,7 +124,7 @@ namespace ngfem
    // NGS_DLL_HEADER virtual void GetPolOrders (FlatArray<PolOrder<DIM> > orders) const;

    HD NGS_DLL_HEADER 
-    virtual void CalcDualShape (const IntegrationPoint & ip, SliceVector<> shape) const override;
+    virtual void CalcDualShape (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const override;
    
  protected:
    /*
@@ -141,7 +141,7 @@ namespace ngfem
      static_cast<const FEL*> (this) -> T_CalcShape (ip, shape);
    }

-    void CalcDualShape2 (const IntegrationPoint & ip, SliceVector<> shape) const
+    void CalcDualShape2 (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
    {
      throw Exception (string("dual shape not implemented for element ")+typeid(*this).name()); 
    }
diff --git a/fem/tscalarfe_impl.hpp b/fem/tscalarfe_impl.hpp
@@ -959,9 +959,9 @@ namespace ngfem

  template <class FEL, ELEMENT_TYPE ET, class BASE>
  void T_ScalarFiniteElement<FEL,ET,BASE> :: 
-  CalcDualShape (const IntegrationPoint & ip, SliceVector<> shape) const
+  CalcDualShape (const BaseMappedIntegrationPoint & mip, SliceVector<> shape) const
  {
-    static_cast<const FEL*>(this) -> CalcDualShape2 (ip, shape);
+    static_cast<const FEL*>(this) -> CalcDualShape2 (mip, shape);
  }