Merge pull request #225 from mfem/fix_eltrans_mesh

fixing eltrans (TransformBack) and mesh (unique_ptr)
mfem · Jun 30, 2024 · 9ed9f18 · 9ed9f18
2 parents 22aad8c + 0b019ce
commit 9ed9f18
Show file tree

Hide file tree

Showing 7 changed files with 102 additions and 19 deletions.
diff --git a/examples/ex23.py b/examples/ex23.py
@@ -10,6 +10,8 @@
 import numpy as np
 from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum
 
+mfem_version = mfem.MFEM_VERSION_MAJOR*10 + mfem.MFEM_VERSION_MINOR
+
 
 def run(mesh_file="",
         ref_levels=2,
@@ -30,28 +32,32 @@ def __init__(self, fespace, ess_bdr, speed):
 
             self.ess_tdof_list = mfem.intArray()
 
-            rel_tol = 1e-8
             fespace.GetEssentialTrueDofs(ess_bdr, self.ess_tdof_list)
 
             c2 = mfem.ConstantCoefficient(speed*speed)
             K = mfem.BilinearForm(fespace)
             K.AddDomainIntegrator(mfem.DiffusionIntegrator(c2))
             K.Assemble()
 
-            self.Kmat0 = mfem.SparseMatrix()
             self.Kmat = mfem.SparseMatrix()
-            dummy = mfem.intArray()
-            K.FormSystemMatrix(dummy, self.Kmat0)
-            K.FormSystemMatrix(self.ess_tdof_list, self.Kmat)
-            self.K = K
-
             self.Mmat = mfem.SparseMatrix()
+
             M = mfem.BilinearForm(fespace)
             M.AddDomainIntegrator(mfem.MassIntegrator())
             M.Assemble()
+
+            if mfem_version < 47:
+                dummy = mfem.intArray()
+                self.Kmat0 = mfem.SparseMatrix()
+                K.FormSystemMatrix(dummy, self.Kmat0)
+            K.FormSystemMatrix(self.ess_tdof_list, self.Kmat)
             M.FormSystemMatrix(self.ess_tdof_list, self.Mmat)
+
+            self.K = K
             self.M = M
 
+            rel_tol = 1e-8
+
             M_solver = mfem.CGSolver()
             M_prec = mfem.DSmoother()
             M_solver.iterative_mode = False
@@ -76,14 +82,19 @@ def __init__(self, fespace, ess_bdr, speed):
             self.T_solver = T_solver
             self.T = None
 
+            self.z = mfem.Vector(self.Height())
+
         def Mult(self, u, du_dt, d2udt2):
             # Compute:
             #    d2udt2 = M^{-1}*-K(u)
             # for d2udt2
-            z = mfem.Vector(u.Size())
-            self.Kmat.Mult(u, z)
+            z = self.z
+            self.K.FullMult(u, z)
             z.Neg()  # z = -z
+
+            z.SetSubVector(self.ess_tdof_list, 0.0)
             self.M_solver.Mult(z, d2udt2)
+            d2udt2.SetSubVector(self.ess_tdof_list, 0.0)
 
         def ImplicitSolve(self, fac0, fac1, u, dudt, d2udt2):
             # Solve the equation:
@@ -92,16 +103,21 @@ def ImplicitSolve(self, fac0, fac1, u, dudt, d2udt2):
             if self.T is None:
                 self.T = mfem.Add(1.0, self.Mmat, fac0, self.Kmat)
                 self.T_solver.SetOperator(self.T)
-            z = mfem.Vector(u.Size())
-            self.Kmat0.Mult(u, z)
+
+            z = self.z
+            self.K.FullMult(u, z)
             z.Neg()
+            z.SetSubVector(self.ess_tdof_list, 0.0)
 
             # iterate over Array<int> :D
-            for j in self.ess_tdof_list:
-                z[j] = 0.0
+            # for j in self.ess_tdof_list:
+            #    z[j] = 0.0
 
             self.T_solver.Mult(z, d2udt2)
 
+            if mfem_version >= 47:
+                d2udt2.SetSubVector(self.ess_tdof_list, 0.0)
+
         def SetParameters(self, u):
             self.T = None
 
@@ -169,7 +185,7 @@ def EvalValue(self, x):
         if (dirichlet):
             ess_bdr.Assign(1)
         else:
-            ess_bdr.Assigne(0)
+            ess_bdr.Assign(0)
 
     oper = WaveOperator(fespace, ess_bdr, speed)
 
@@ -200,7 +216,7 @@ def EvalValue(self, x):
             print("GLVis visualization disabled.")
         else:
             sout.precision(output.precision)
-            sout << "solution\n" << mesh << dudt_gf
+            sout << "solution\n" << mesh << u_gf
             sout << "pause\n"
             sout.flush()
             print(

diff --git a/mfem/_par/eltrans.i b/mfem/_par/eltrans.i
@@ -22,6 +22,8 @@ import_array();
 %import "intrules.i"
 %import "geom.i"
 
+%import "../common/mfem_config.i"
+
 %feature("shadow") mfem::ElementTransformation::Transform %{
 def Transform(self, *args):
     from .vector import Vector
@@ -35,5 +37,33 @@ def Transform(self, *args):
         return _eltrans.ElementTransformation_Transform(self, *args)
 %}
 
+%ignore mfem::ElementTransformation::TransformBack;
+%ignore mfem::IsoparametricTransformation::TransformBack;
+
 %include "fem/eltrans.hpp"
 
+//
+//  special handling for TransformBack (this is because tol_0 is protected)
+//
+namespace mfem{
+  #ifdef MFEM_USE_DOUBLE
+  %extend IsoparametricTransformation{
+     virtual int _TransformBack(const Vector &pt, IntegrationPoint &ip,
+                                const real_t phys_tol = 1e-15){
+       return self-> TransformBack(pt, ip, phys_tol);
+     }
+   };  //end of extend
+  #elif defined(MFEM_USE_SINGLE)
+  %extend IsoparametricTransformation{
+     virtual int _TransformBack(const Vector &pt, IntegrationPoint &ip,
+                                const real_t phys_tol = 1e-7){
+       return self-> TransformBack(pt, ip, phys_tol);
+     }
+   };  //end of extend
+  #endif
+ } //end of namespace
+
+%pythoncode %{
+if hasattr(IsoparametricTransformation, "_TransformBack"):
+    IsoparametricTransformation.TransformBack = IsoparametricTransformation._TransformBack
+%}
diff --git a/mfem/_par/mesh.i b/mfem/_par/mesh.i
@@ -309,6 +309,10 @@ def CartesianPartitioning(self, nxyz, return_list=False):
 %newobject mfem::Mesh::GetFaceToElementTable;
 %newobject mfem::Mesh::GetVertexToElementTable;
 
+%immutable mfem::MeshPart::mesh;
+%immutable mfem::MeshPart::nodal_fes;
+%immutable mfem::MeshPart::nodes;
+
 %include "mesh/mesh.hpp"
 %mutable;
 

diff --git a/mfem/_par/pbilinearform.i b/mfem/_par/pbilinearform.i
@@ -23,6 +23,7 @@ import_array();
 %import "handle.i"
 %import "bilinearform.i"
 %import "pfespace.i"
+%import "pgridfunc.i"
 %import "hypre.i"
 %import "../common/exception.i"
 

diff --git a/mfem/_ser/eltrans.i b/mfem/_ser/eltrans.i
@@ -22,6 +22,7 @@ import_array();
 %import "intrules.i"
 %import "geom.i"
 %import "../common/exception.i"
+%import "../common/mfem_config.i"
 
 %feature("shadow") mfem::ElementTransformation::Transform %{
 def Transform(self, *args):
@@ -39,5 +40,36 @@ def Transform(self, *args):
 %include "../common/deprecation.i"
 DEPRECATED_METHOD(mfem::IsoparametricTransformation::FinalizeTransformation())
 
+
+%ignore mfem::ElementTransformation::TransformBack;
+%ignore mfem::IsoparametricTransformation::TransformBack;
+
 %include "fem/eltrans.hpp"
 
+//
+//  special handling for TransformBack (this is because tol_0 is protected)
+//
+namespace mfem{
+  #ifdef MFEM_USE_DOUBLE
+  %extend IsoparametricTransformation{
+     virtual int _TransformBack(const Vector &pt, IntegrationPoint &ip,
+                                const real_t phys_tol = 1e-15){
+       return self-> TransformBack(pt, ip, phys_tol);
+     }
+   };  //end of extend
+  #elif defined(MFEM_USE_SINGLE)
+  %extend IsoparametricTransformation{
+     virtual int _TransformBack(const Vector &pt, IntegrationPoint &ip,
+                                const real_t phys_tol = 1e-7){
+       return self-> TransformBack(pt, ip, phys_tol);
+     }
+   };  //end of extend
+  #endif
+ } //end of namespace
+
+%pythoncode %{
+if hasattr(IsoparametricTransformation, "_TransformBack"):
+    IsoparametricTransformation.TransformBack = IsoparametricTransformation._TransformBack
+%}
+
+
diff --git a/mfem/_ser/mesh.i b/mfem/_ser/mesh.i
@@ -313,6 +313,10 @@ def CartesianPartitioning(self, nxyz, return_list=False):
 %newobject mfem::Mesh::GetFaceToElementTable;
 %newobject mfem::Mesh::GetVertexToElementTable;
 
+%immutable mfem::MeshPart::mesh;
+%immutable mfem::MeshPart::nodal_fes;
+%immutable mfem::MeshPart::nodes;
+
 %include "../common/exception.i"
 %include "mesh/mesh.hpp"
 

diff --git a/mfem/common/bilininteg_ext.i b/mfem/common/bilininteg_ext.i
@@ -255,10 +255,6 @@ namespace mfem {
     self._coeff = args
 %}
 
-%pythonappend ElasticityComponentIntegrator::ElasticityComponentIntegrator %{
-    self._coeff = parent_
-%}
-
 %pythonappend DGTraceIntegrator::DGTraceIntegrator %{
     self._coeff = args
 %}