Merge branch 'hdivdiv_calcmappedshape' into 'master'

Hdivdiv calcmappedshape See merge request !236
NGSolve · Sep 20, 2017 · dc1159af93c19e70c6308e8c315810bfaca64db1 · dc1159a
2 parents 0704e78 + 3d60d11
commit dc1159af93c19e70c6308e8c315810bfaca64db1
Unified Split

Showing with 509 additions and 95 deletions.

+164 −23 comp/hdivdivfespace.cpp

+18 −0 comp/hdivfes.cpp

+7 −1 comp/python_comp.cpp

+279 −69 fem/hdivdivfe.hpp

+41 −2 ngstd/autodiffdiff.hpp
diff --git a/comp/hdivdivfespace.cpp b/comp/hdivdivfespace.cpp
@@ -14,7 +14,124 @@ namespace ngcomp
 {
  
  template<int D>
-  class DiffOpIdHDivDiv : public DiffOp<DiffOpIdHDivDiv<D> >
+  class DiffOpVecIdHDivDiv: public DiffOp<DiffOpVecIdHDivDiv<D> >
+  {
+  public:
+    enum { DIM = 1 };
+    enum { DIM_SPACE = D };
+    enum { DIM_ELEMENT = D };
+    enum { DIM_DMAT = D*(D+1)/2 };
+    enum { DIFFORDER = 0 };
+    enum { DIM_STRESS = D*(D+1)/2 };
+
+    static Array<int> GetDimensions() { return Array<int> ({D*(D+1)/2,1}); }
+
+    template <typename FEL,typename SIP>
+    static void GenerateMatrix(const FEL & bfel,const SIP & mip,
+      SliceMatrix<double,ColMajor> mat,LocalHeap & lh)
+    {
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+      fel.CalcMappedShape_Vector (mip,Trans(mat));
+    }
+
+    template <typename FEL,typename SIP,typename MAT>
+    static void GenerateMatrix(const FEL & bfel,const SIP & sip,
+      MAT & mat,LocalHeap & lh)
+    {
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+      int nd = fel.GetNDof();
+      FlatMatrix<> shape(nd,DIM_DMAT,lh);
+      fel.CalcMappedShape_Vector(sip,shape);
+      for(int i=0; i<nd; i++)
+        for(int j = 0; j <DIM_DMAT; j++)
+          mat(j,i) = shape(i,j);
+
+    }
+  };
+
+  template<int D>
+  class DiffOpIdHDivDiv: public DiffOp<DiffOpIdHDivDiv<D> >
+  {
+  public:
+    enum { DIM = 1 };
+    enum { DIM_SPACE = D };
+    enum { DIM_ELEMENT = D };
+    enum { DIM_DMAT = D*D };
+    enum { DIFFORDER = 0 };
+    enum { DIM_STRESS = D*D };
+
+    static Array<int> GetDimensions() { return Array<int> ({D,D}); }
+
+    template <typename FEL,typename SIP>
+    static void GenerateMatrix(const FEL & bfel,const SIP & mip,
+      SliceMatrix<double,ColMajor> mat,LocalHeap & lh)
+    {
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+      fel.CalcMappedShape_Matrix (mip,Trans(mat));
+    }
+
+    template <typename FEL,typename SIP,typename MAT>
+    static void GenerateMatrix(const FEL & bfel,const SIP & sip,
+      MAT & mat,LocalHeap & lh)
+    {
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+      int nd = fel.GetNDof();
+      FlatMatrix<> shape(nd,DIM_DMAT,lh);
+      fel.CalcMappedShape_Matrix(sip,shape);
+      for(int i=0; i<nd; i++)
+        for(int j = 0; j <DIM_DMAT; j++)
+          mat(j,i) = shape(i,j);
+
+    }
+  };
+
+  template<int D>
+  class DiffOpDivHDivDiv: public DiffOp<DiffOpDivHDivDiv<D> >
+  {
+  public:
+    enum { DIM = 1 };
+    enum { DIM_SPACE = D };
+    enum { DIM_ELEMENT = D };
+    enum { DIM_DMAT = D };
+    enum { DIFFORDER = 1 };
+    enum { DIM_STRESS = (D*(D+1))/2 };
+
+    template <typename FEL,typename SIP>
+    static void GenerateMatrix(const FEL & bfel,const SIP & sip,
+      SliceMatrix<double,ColMajor> mat,LocalHeap & lh)
+    {
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+
+      fel.CalcMappedDivShape (sip, Trans(mat));
+    }
+
+    template <typename FEL,typename SIP,typename MAT>
+    static void GenerateMatrix(const FEL & bfel,const SIP & sip,
+      MAT & mat,LocalHeap & lh)
+    {
+      HeapReset hr(lh);
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
+
+      int nd = fel.GetNDof();
+      FlatMatrix<> divshape(nd, D, lh);
+      fel.CalcMappedDivShape (sip, divshape);
+      for (int i=0; i<nd; i++)
+        for (int j=0; j<D; j++)
+          mat(j,i) = divshape(i,j);
+
+    }
+
+  };
+
+
+  template<int D>
+  class DiffOpIdHDivDiv_old : public DiffOp<DiffOpIdHDivDiv_old<D> >
  { 
  public:
    enum { DIM = 1 };
@@ -30,8 +147,8 @@ namespace ngcomp
    static void GenerateMatrix(const FEL & bfel, const SIP & sip,
                               MAT & mat, LocalHeap & lh)
    {
-      const HDivDivFiniteElement & fel =
-        dynamic_cast<const HDivDivFiniteElement&> (bfel);
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
      
      int nd = fel.GetNDof();
      
@@ -74,8 +191,8 @@ namespace ngcomp
  };


-    template<int D>
-  class DiffOpVecIdHDivDiv : public DiffOp<DiffOpVecIdHDivDiv<D> >
+  template<int D>
+  class DiffOpVecIdHDivDiv_old : public DiffOp<DiffOpVecIdHDivDiv_old<D> >
  { 
  public:
    enum { DIM = 1 };
@@ -91,8 +208,8 @@ namespace ngcomp
    static void GenerateMatrix(const FEL & bfel, const SIP & sip,
                               MAT & mat, LocalHeap & lh)
    {
-      const HDivDivFiniteElement & fel =
-        dynamic_cast<const HDivDivFiniteElement&> (bfel);
+      const HDivDivFiniteElement<D> & fel =
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
      
      int nd = fel.GetNDof();
      
@@ -151,7 +268,8 @@ namespace ngcomp
  };


-  template <int D> class DiffOpDivHDivDiv : public DiffOp<DiffOpDivHDivDiv<D> >
+
+  template <int D> class DiffOpDivHDivDiv_old : public DiffOp<DiffOpDivHDivDiv_old<D> >
  {
    
  public:
@@ -168,8 +286,11 @@ namespace ngcomp
    static void GenerateMatrix (const FEL & bfel, const SIP & sip,
                                MAT & mat, LocalHeap & lh)
    {
-      const HDivDivFiniteElement & fel = 
-        dynamic_cast<const HDivDivFiniteElement&> (bfel);
+      static int timer = NgProfiler::CreateTimer ("old div");
+      NgProfiler::RegionTimer reg (timer);
+
+      const HDivDivFiniteElement<D> & fel = 
+        dynamic_cast<const HDivDivFiniteElement<D>&> (bfel);
      
      int nd = fel.GetNDof();
      
@@ -230,9 +351,9 @@ namespace ngcomp
              sigma_ref(0,0) = shape(i, 0);
              sigma_ref(1,1) = shape(i, 1);
              sigma_ref(2,2) = shape(i, 2);
-              sigma_ref(0,1) = sigma_ref(1,0) = shape(i, 3);
+              sigma_ref(2,1) = sigma_ref(1,2) = shape(i, 3);
              sigma_ref(0,2) = sigma_ref(2,0) = shape(i, 4);
-              sigma_ref(2,1) = sigma_ref(1,2) = shape(i, 5);
+              sigma_ref(0,1) = sigma_ref(1,0) = shape(i, 5);
            }
          
          Vec<D> hv2;
@@ -244,22 +365,22 @@ namespace ngcomp
          
          hv2 *= iad_det.Value();
          
-          for ( int j = 0; j < D; j++ )
-            for ( int k = 0; k < D; k++ )
-              for ( int l = 0; l < D; l++ )
-                for ( int m = 0; m < D; m++ )
-                  for ( int n = 0; n < D; n++ )
-                    hv2(n) += inv_jac(m,k) *fad(n,j).Value() * sigma_ref(j,l) * fad(k,l).DValue(m);
+          // this term is zero, check why
+          //for ( int j = 0; j < D; j++ )
+          //  for ( int k = 0; k < D; k++ )
+          //    for ( int l = 0; l < D; l++ )
+          //      for ( int m = 0; m < D; m++ )
+          //        for ( int n = 0; n < D; n++ )
+          //          hv2(n) += inv_jac(m,k) *fad(n,j).Value() * sigma_ref(j,l) * fad(k,l).DValue(m);
          
          for ( int j = 0; j < D; j++)
            mat(j,i) += hv2(j);
        }
      
    }
  };
-  

-  
+
  template <int D>
  class NGS_DLL_HEADER HDivDivMassIntegrator 
    : public T_BDBIntegrator<DiffOpIdHDivDiv<D>, DiagDMat<D*D> >
@@ -407,6 +528,19 @@ namespace ngcomp
    {
      ctofdof[i] = discontinuous ? LOCAL_DOF : INTERFACE_DOF;
    }
+    if (discontinuous) return;
+
+    Array<int> innerdofs;
+    for(auto e: ma->Elements())
+    {
+      GetInnerDofNrs(e.Nr(), innerdofs);
+      for (int dof: innerdofs)
+      {
+        ctofdof[dof] = LOCAL_DOF;
+      }
+    }
+
+
  }


@@ -506,16 +640,23 @@ namespace ngcomp
    switch(ma->GetDimension())
    {
    case 2:
-      additional.Set ("vec",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv<2>>> ()); break;
+      additional.Set ("vec",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv<2>>> ());
+      additional.Set ("id_old",make_shared<T_DifferentialOperator<DiffOpIdHDivDiv_old<2>>> ());
+      additional.Set ("vec_old",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv_old<2>>> ());
+      additional.Set ("div_old",make_shared<T_DifferentialOperator<DiffOpDivHDivDiv_old<2>>> ());
      break;
    case 3:
-      additional.Set ("vec",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv<3>>> ()); break;
+      additional.Set ("vec",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv<3>>> ());
+      additional.Set ("id_old",make_shared<T_DifferentialOperator<DiffOpIdHDivDiv_old<3>>> ());
+      additional.Set ("vec_old",make_shared<T_DifferentialOperator<DiffOpVecIdHDivDiv_old<3>>> ());
+      additional.Set ("div_old",make_shared<T_DifferentialOperator<DiffOpDivHDivDiv_old<3>>> ());
+      break;
    default:
      ;
    }
    return additional;
  }
-
+  


  static RegisterFESpace<HDivDivFESpace> init ("hdivdiv");
diff --git a/comp/hdivfes.cpp b/comp/hdivfes.cpp
@@ -59,6 +59,24 @@ namespace ngcomp
      coeffs[0] = shared_ptr<CoefficientFunction> (new ConstantCoefficientFunction(1));
      integrator[VOL] = GetIntegrators().CreateBFI("masshdiv", 2, coeffs);
    }
+    if (ma->GetDimension() == 3)
+    {
+      Array<shared_ptr<CoefficientFunction>> coeffs(1);
+      coeffs[0] = shared_ptr<CoefficientFunction> (new ConstantCoefficientFunction(1));
+      integrator[VOL] = GetIntegrators().CreateBFI("masshdiv", 3, coeffs);
+    }
+    if (ma->GetDimension() == 2)
+      {
+        evaluator[VOL] = make_shared<T_DifferentialOperator<DiffOpIdHDiv<2>>>();
+        evaluator[BND] = make_shared<T_DifferentialOperator<DiffOpIdVecHDivBoundary<2>>>();
+        flux_evaluator[VOL] = make_shared<T_DifferentialOperator<DiffOpDivHDiv<2>>>();
+      }
+    else
+      {
+        evaluator[VOL] = make_shared<T_DifferentialOperator<DiffOpIdHDiv<3>>>();
+        evaluator[BND] = make_shared<T_DifferentialOperator<DiffOpIdVecHDivBoundary<3>>>();
+        flux_evaluator[VOL] = make_shared<T_DifferentialOperator<DiffOpDivHDiv<3>>>();
+      }
  }
      
    RaviartThomasFESpace :: ~RaviartThomasFESpace ()
diff --git a/comp/python_comp.cpp b/comp/python_comp.cpp
@@ -1603,7 +1603,13 @@ kwargs : For a description of the possible kwargs have a look a bit further down
         },
         py::arg("heapsize")=1000000,
         "update space after mesh-refinement")
-
+     .def("FinalizeUpdate", [](shared_ptr<FESpace> self, int heapsize)
+         { 
+           LocalHeap lh (heapsize, "FESpace::FinalizeUpdate-heap");
+           self->FinalizeUpdate(lh);
+         },
+         py::arg("heapsize")=1000000,
+         "finalize update")
    .def_property_readonly ("ndof", [](shared_ptr<FESpace> self) { return self->GetNDof(); },
                            "number of degrees of freedom")