Merge branch 'custom_intrules_symbolicenergy' into 'master'

custom integration rules for symbolic energy See merge request !245
NGSolve · Oct 4, 2017 · 9f00285fd957f20c290b48fd19c1ef5271438cc5 · 9f00285
2 parents b5888c4 + 1f86746
commit 9f00285fd957f20c290b48fd19c1ef5271438cc5
Unified Split

Showing with 86 additions and 71 deletions.

+5 −4 comp/python_comp.cpp

+26 −1 fem/integrator.hpp

+10 −1 fem/python_fem.cpp

+39 −53 fem/symbolicintegrator.cpp

+5 −12 fem/symbolicintegrator.hpp

+1 −0 include/ngs_stdcpp_include.hpp
diff --git a/comp/python_comp.cpp b/comp/python_comp.cpp
@@ -3113,7 +3113,7 @@ flags : dict

 	     if (ir.Size())
               {
-                 cout << IM(5) << "ir = " << ir << endl;
+                 cout << IM(1) << "WARNING: Setting the integration rule for all element types is deprecated, use LFI.SetIntegrationRule(ELEMENT_TYPE, IntegrationRule) instead!" << endl;
                 dynamic_pointer_cast<SymbolicLinearFormIntegrator>
 		   (lfi)->SetIntegrationRule(ir);                   
               }
@@ -3171,7 +3171,7 @@ flags : dict

             if (ir.Size())
               {
-                 cout << IM(5) << "ir = " << ir << endl;
+                 cout << IM(1) << "WARNING: Setting the integration rule for all element types is deprecated, use BFI.SetIntegrationRule(ELEMENT_TYPE, IntegrationRule) instead!" << endl;
                 dynamic_pointer_cast<SymbolicBilinearFormIntegrator> (bfi)
                   ->SetIntegrationRule(ir);
               }
@@ -3231,7 +3231,8 @@ flags : dict
          );
          
  m.def("SymbolicEnergy",
-          [](spCF cf, VorB vb, py::object definedon, py::object definedonelem) -> shared_ptr<BilinearFormIntegrator>
+        [](spCF cf, VorB vb, py::object definedon, py::object definedonelem)
+        -> shared_ptr<BilinearFormIntegrator>
           {
             py::extract<Region> defon_region(definedon);
             if (defon_region.check())
@@ -3249,7 +3250,7 @@ flags : dict
             return bfi;
           },
           py::arg("coefficient"), py::arg("VOL_or_BND")=VOL, py::arg("definedon")=DummyArgument(),
-           py::arg("definedonelements")=DummyArgument()
+        py::arg("definedonelements")=DummyArgument()
          );


diff --git a/fem/integrator.hpp b/fem/integrator.hpp
@@ -53,7 +53,8 @@ namespace ngfem

    /// define only on some sub-domains
    shared_ptr<BitArray> definedon_element = nullptr;
-    
+    std::array<unique_ptr<IntegrationRule>,25> userdefined_intrules;
+    std::array<unique_ptr<SIMD_IntegrationRule>,25> userdefined_simd_intrules;
  
  protected:
    void DeleteCurveIPs ( void );
@@ -103,6 +104,30 @@ namespace ngfem
      definedon_element = adefinedonelem;
    }

+    void SetIntegrationRule(ELEMENT_TYPE et, const IntegrationRule& ir)
+    {
+      userdefined_intrules[int(et)] = make_unique<IntegrationRule>(ir.Copy());
+      userdefined_simd_intrules[int(et)] = make_unique<SIMD_IntegrationRule>(*userdefined_intrules[int(et)]);
+    }
+
+    void SetIntegrationRule(const IntegrationRule& ir)
+    {
+      for(auto i : Range(25))
+        {
+          userdefined_intrules[i] = make_unique<IntegrationRule>(ir.Copy());
+          userdefined_simd_intrules[i] = make_unique<SIMD_IntegrationRule>(*userdefined_intrules[i]);
+        }
+    }
+
+    inline const IntegrationRule& GetIntegrationRule(ELEMENT_TYPE et, int order) const
+    {
+      return userdefined_intrules[et] ? *userdefined_intrules[et] : SelectIntegrationRule(et,order);
+    }
+    inline const SIMD_IntegrationRule& GetSIMDIntegrationRule(ELEMENT_TYPE et, int order) const
+    {
+      return userdefined_simd_intrules[et] ? *userdefined_simd_intrules[et] : SIMD_SelectIntegrationRule(et,order);
+    }
+
    /// defined only on some elements/facets/boundary elements
    shared_ptr<BitArray> GetDefinedOnElements () const { return definedon_element; } 
    
diff --git a/fem/python_fem.cpp b/fem/python_fem.cpp
@@ -1307,7 +1307,11 @@ void NGS_DLL_HEADER ExportNgfem(py::module &m) {

    .def("SetDefinedOnElements",  [](shared_ptr<BFI> self, shared_ptr<BitArray> ba )
                                                  { self->SetDefinedOnElements (ba); } )
-
+    .def("SetIntegrationRule", [] (shared_ptr<BFI> self, ELEMENT_TYPE et, IntegrationRule ir)
+         {
+           self -> SetIntegrationRule(et,ir);
+           return self;
+         })
    .def("CalcElementMatrix",
         [] (shared_ptr<BFI> self,
             const FiniteElement & fe, const ElementTransformation &trafo,
@@ -1404,6 +1408,11 @@ void NGS_DLL_HEADER ExportNgfem(py::module &m) {
         py::return_value_policy::reference)
    .def("SetDefinedOnElements",  [](shared_ptr<LFI> self, shared_ptr<BitArray> ba )
                                                  { self->SetDefinedOnElements (ba); } )
+    .def("SetIntegrationRule", [](shared_ptr<LFI> self, ELEMENT_TYPE et, IntegrationRule ir)
+         {
+           self->SetIntegrationRule(et,ir);
+           return self;
+         })

    .def("CalcElementVector", 
        static_cast<void(LinearFormIntegrator::*)(const FiniteElement&, const ElementTransformation&, FlatVector<double>, LocalHeap&)const>(&LinearFormIntegrator::CalcElementVector))