fixed ex26, ex26p. adjusted to use a constrained prolongation operator

mfem · Mar 23, 2024 · 03a0b1d · 03a0b1d
1 parent 5d8ee99
commit 03a0b1d
Show file tree

Hide file tree

Showing 9 changed files with 95 additions and 81 deletions.
diff --git a/examples/ex26.py b/examples/ex26.py
@@ -18,34 +18,30 @@ def run(order_refinements=2,
 
     class DiffusionMultigrid(mfem.PyGeometricMultigrid):
         def __init__(self, fespaces, ess_bdr):
-            mfem.PyGeometricMultigrid.__init__(self, fespaces)
+            mfem.PyGeometricMultigrid.__init__(self, fespaces, ess_bdr)
 
             self.smoothers = []
-            self.ConstructCoarseOperatorAndSolver(fespaces.GetFESpaceAtLevel(0),
-                                                  ess_bdr)
+            self.ConstructCoarseOperatorAndSolver(
+                fespaces.GetFESpaceAtLevel(0),)
 
             for level in range(1, fespaces.GetNumLevels()):
                 self.ConstructOperatorAndSmoother(fespaces.GetFESpaceAtLevel(level),
-                                                  ess_bdr)
+                                                  level)
 
-        def ConstructBilinearForm(self, fespace, ess_bdr):
+        def ConstructBilinearForm(self, fespace):
             form = mfem.BilinearForm(fespace)
             form.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL)
             form.AddDomainIntegrator(mfem.DiffusionIntegrator(one))
             form.Assemble()
             self.AppendBilinearForm(form)
 
-            ess = mfem.intArray()
-            self.AppendEssentialTDofs(ess)
-            fespace.GetEssentialTrueDofs(ess_bdr, ess)
+        def ConstructCoarseOperatorAndSolver(self, coarse_fespace):
 
-        def ConstructCoarseOperatorAndSolver(self, coarse_fespace, ess_bdr):
-
-            self.ConstructBilinearForm(coarse_fespace, ess_bdr)
+            self.ConstructBilinearForm(coarse_fespace)
 
             opr = mfem.OperatorPtr()
             opr.SetType(mfem.Operator.ANY_TYPE)
-            self.bfs[-1].FormSystemMatrix(self.essentialTrueDofs[-1], opr)
+            self.bfs[0].FormSystemMatrix(self.essentialTrueDofs[0], opr)
             opr.SetOperatorOwner(False)
 
             pcg = mfem.CGSolver()
@@ -58,20 +54,20 @@ def ConstructCoarseOperatorAndSolver(self, coarse_fespace, ess_bdr):
             self.smoothers.append((pcg, opr))
             self.AddLevel(opr.Ptr(), pcg, False, False)
 
-        def ConstructOperatorAndSmoother(self, fespace, ess_bdr):
-            self.ConstructBilinearForm(fespace, ess_bdr)
+        def ConstructOperatorAndSmoother(self, fespace, level):
+            self.ConstructBilinearForm(fespace)
 
             opr = mfem.OperatorPtr()
             opr.SetType(mfem.Operator.ANY_TYPE)
-            self.bfs[-1].FormSystemMatrix(self.essentialTrueDofs[-1], opr)
+
+            self.bfs[level].FormSystemMatrix(
+                self.essentialTrueDofs[level], opr)
             opr.SetOperatorOwner(False)
 
             diag = mfem.Vector(fespace.GetTrueVSize())
-            self.bfs[-1].AssembleDiagonal(diag)
-
+            self.bfs[level].AssembleDiagonal(diag)
             smoother = mfem.OperatorChebyshevSmoother(opr.Ptr(), diag,
-                                                      self.essentialTrueDofs[-1], 2)
-
+                                                      self.essentialTrueDofs[level], 2)
             self.smoothers.append((smoother, opr))
             self.AddLevel(opr.Ptr(), smoother, False, False)
 

diff --git a/examples/ex26p.py b/examples/ex26p.py
@@ -23,36 +23,31 @@ def run(order_refinements=2,
 
     class DiffusionMultigrid(mfem.PyGeometricMultigrid):
         def __init__(self, fespaces, ess_bdr):
-            mfem.PyGeometricMultigrid.__init__(self, fespaces)
+            mfem.PyGeometricMultigrid.__init__(self, fespaces, ess_bdr)
 
             self.smoothers = []
-            self.ConstructCoarseOperatorAndSolver(fespaces.GetFESpaceAtLevel(0),
-                                                  ess_bdr)
+            self.ConstructCoarseOperatorAndSolver(
+                fespaces.GetFESpaceAtLevel(0),)
 
             for level in range(1, fespaces.GetNumLevels()):
                 self.ConstructOperatorAndSmoother(fespaces.GetFESpaceAtLevel(level),
-                                                  ess_bdr)
+                                                  level)
 
-        def ConstructBilinearForm(self, fespace, ess_bdr, partial_assembly):
+        def ConstructBilinearForm(self, fespace, partial_assembly):
             form = mfem.ParBilinearForm(fespace)
             if partial_assembly:
                 form.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL)
             form.AddDomainIntegrator(mfem.DiffusionIntegrator(one))
             form.Assemble()
             self.AppendBilinearForm(form)
 
-            ess = mfem.intArray()
-            self.AppendEssentialTDofs(ess)
-            fespace.GetEssentialTrueDofs(ess_bdr, ess)
+        def ConstructCoarseOperatorAndSolver(self, coarse_fespace,):
 
-        def ConstructCoarseOperatorAndSolver(self, coarse_fespace,
-                                             ess_bdr):
-
-            self.ConstructBilinearForm(coarse_fespace, ess_bdr, False)
+            self.ConstructBilinearForm(coarse_fespace, False)
 
             hypreCoarseMat = mfem.HypreParMatrix()
-            self.bfs[-1].FormSystemMatrix(self.essentialTrueDofs[-1],
-                                          hypreCoarseMat)
+            self.bfs[0].FormSystemMatrix(self.essentialTrueDofs[0],
+                                         hypreCoarseMat)
 
             amg = mfem.HypreBoomerAMG(hypreCoarseMat)
             amg.SetPrintLevel(-1)
@@ -69,20 +64,21 @@ def ConstructCoarseOperatorAndSolver(self, coarse_fespace,
             self.AddLevel(hypreCoarseMat, pcg, False, False)
             self.amg = amg
 
-        def ConstructOperatorAndSmoother(self, fespace, ess_bdr):
-            self.ConstructBilinearForm(fespace, ess_bdr, True)
+        def ConstructOperatorAndSmoother(self, fespace, level):
+            self.ConstructBilinearForm(fespace, True)
 
             opr = mfem.OperatorPtr()
             opr.SetType(mfem.Operator.ANY_TYPE)
 
-            self.bfs[-1].FormSystemMatrix(self.essentialTrueDofs[-1], opr)
+            self.bfs[level].FormSystemMatrix(
+                self.essentialTrueDofs[level], opr)
             opr.SetOperatorOwner(False)
 
             diag = mfem.Vector(fespace.GetTrueVSize())
-            self.bfs[-1].AssembleDiagonal(diag)
+            self.bfs[level].AssembleDiagonal(diag)
 
             smoother = mfem.OperatorChebyshevSmoother(opr.Ptr(), diag,
-                                                      self.essentialTrueDofs[-1],
+                                                      self.essentialTrueDofs[level],
                                                       2,
                                                       fespace.GetParMesh().GetComm())
 

diff --git a/mfem/_par/array.i b/mfem/_par/array.i
@@ -174,3 +174,6 @@ INSTANTIATE_ARRAY_NUMPYARRAY(uint, unsigned int, NPY_UINT)       // 32bit
 %template(intArray2D) mfem::Array2D<int>;
 %template(doubleArray2D) mfem::Array2D<double>;
 
+/* Array< Array<int> *> */
+IGNORE_ARRAY_METHODS(mfem::Array<int> *)
+INSTANTIATE_ARRAY2(Array<int> *, Array<int>, intArray, 1)
diff --git a/mfem/_par/bilinearform.i b/mfem/_par/bilinearform.i
@@ -119,6 +119,10 @@ namespace mfem {
 
 } //end of namespace
 
+%import "../common/array_instantiation_macro.i"
+IGNORE_ARRAY_METHODS(mfem::BilinearForm *)
+INSTANTIATE_ARRAY0(BilinearForm *, BilinearForm, 1)
+
 %include "fem/bilinearform.hpp"
 
 // instatitate template methods

diff --git a/mfem/_par/multigrid.i b/mfem/_par/multigrid.i
@@ -30,27 +30,24 @@ LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array<mfem::Operator*>&,
 LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array<mfem::Solver*>&,
 			   mfem::Solver*)
 LIST_TO_MFEMOBJ_BOOLARRAY_IN(const mfem::Array<bool>& )
-//%feature("director") mfem::PyGeometricMultigrid;
 
 %pythonprepend mfem::PyGeometricMultigrid::AppendBilinearForm %{
-   if not hasattr(self, "_forms"): self._forms = []
-   self._forms.append(form)
-   form.thisown = 0
-%}
-%pythonprepend mfem::PyGeometricMultigrid::AppendEssentialTDofs %{
-   if not hasattr(self, "_esss"): self._esss = []
-   self._esss.append(ess)
-   ess.thisown = 0
+  if not hasattr(self, "forms"):
+      self.forms=[]
+  self.forms.append(form)
+  form.thisown=0
 %}
+
 %feature("shadow") mfem::PyGeometricMultigrid::_pybfs %{
   @property
   def bfs(self):
-     return self._forms
- %}
+     return self._GetBilinearFormArray()
+%}
+
 %feature("shadow") mfem::PyGeometricMultigrid::_pyess %{
   @property
   def essentialTrueDofs(self):
-     return self._esss
+     return self._GetEssentialTrueDofs()
 %}
 
 %include "fem/multigrid.hpp"
@@ -60,15 +57,25 @@ LIST_TO_MFEMOBJ_BOOLARRAY_IN(const mfem::Array<bool>& )
 class PyGeometricMultigrid : public GeometricMultigrid
 {
 public:
- PyGeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_)
+  PyGeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_)
    : GeometricMultigrid(fespaces_){}
 
+  PyGeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_,
+                      const Array<int> &ess_bdr_)
+   : GeometricMultigrid(fespaces_, ess_bdr_){}
+
   void AppendBilinearForm(BilinearForm *form){
     bfs.Append(form);
   }
-  void AppendEssentialTDofs(Array<int> *ess){
-      essentialTrueDofs.Append(ess);
+
+  Array<BilinearForm*> *_GetBilinearFormArray(){
+    return &bfs;
+  }
+
+  Array<Array<int> *> *_GetEssentialTrueDofs(){
+    return &essentialTrueDofs;
   }
+
   void _pybfs(void){}
   void _pyess(void){}
 };

diff --git a/mfem/_ser/array.i b/mfem/_ser/array.i
@@ -186,3 +186,7 @@ INSTANTIATE_ARRAY_NUMPYARRAY(uint, unsigned int, NPY_UINT)       // 32bit
 %template(intArray2D) mfem::Array2D<int>;
 %template(doubleArray2D) mfem::Array2D<double>;
 
+/* Array< Array<int> *> */
+IGNORE_ARRAY_METHODS(mfem::Array<int> *)
+INSTANTIATE_ARRAY2(Array<int> *, Array<int>, intArray, 1)
+
diff --git a/mfem/_ser/multigrid.i b/mfem/_ser/multigrid.i
@@ -31,33 +31,26 @@ LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array<mfem::Operator*>&,
 LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array<mfem::Solver*>&,
 			   mfem::Solver*)
 LIST_TO_MFEMOBJ_BOOLARRAY_IN(const mfem::Array<bool>& )
- /*
-ObjectArrayInput(mfem::Solver *);
-ObjectArrayInput(mfem::Operator *);
-BoolArrayInput(bool);
- */
-//%feature("director") mfem::PyGeometricMultigrid;
 
 %pythonprepend mfem::PyGeometricMultigrid::AppendBilinearForm %{
-   if not hasattr(self, "_forms"): self._forms = []
-   self._forms.append(form)
-   form.thisown = 0
-%}
-%pythonprepend mfem::PyGeometricMultigrid::AppendEssentialTDofs %{
-   if not hasattr(self, "_esss"): self._esss = []
-   self._esss.append(ess)
-   ess.thisown = 0
+  if not hasattr(self, "forms"):
+      self.forms=[]
+  self.forms.append(form)
+  form.thisown=0
 %}
+
 %feature("shadow") mfem::PyGeometricMultigrid::_pybfs %{
   @property
   def bfs(self):
-     return self._forms
- %}
+     return self._GetBilinearFormArray()
+%}
+
 %feature("shadow") mfem::PyGeometricMultigrid::_pyess %{
   @property
   def essentialTrueDofs(self):
-     return self._esss
+     return self._GetEssentialTrueDofs()
 %}
+
 %include "fem/multigrid.hpp"
 
 %inline %{
@@ -67,18 +60,24 @@ class PyGeometricMultigrid : public GeometricMultigrid
 public:
  PyGeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_)
    : GeometricMultigrid(fespaces_){}
+ PyGeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_,
+                      const Array<int> &ess_bdr_)
+   : GeometricMultigrid(fespaces_, ess_bdr_){}
 
   void AppendBilinearForm(BilinearForm *form){
     bfs.Append(form);
   }
-  void AppendEssentialTDofs(Array<int> *ess){
-      essentialTrueDofs.Append(ess);
+
+  Array<BilinearForm*> *_GetBilinearFormArray(){
+    return &bfs;
+  }
+
+  Array<Array<int> *> *_GetEssentialTrueDofs(){
+    return &essentialTrueDofs;
   }
+
   void _pybfs(void){}
   void _pyess(void){}
 };
   } /* end of namespace */
 %}
-
-
-
diff --git a/mfem/common/array_instantiation_macro.i b/mfem/common/array_instantiation_macro.i
@@ -1,8 +1,8 @@
-%define INSTANTIATE_ARRAY0(XXX, YYY, USEPTR)
+%define INSTANTIATE_ARRAY2(XXX, YYY, ZZZ, USEPTR)
 #if USEPTR == 1
-%template(##YYY##Ptr##Array) mfem::Array<mfem::XXX>;
+%template(##ZZZ##Ptr##Array) mfem::Array<mfem::XXX>;
 #else
-%template(##YYY##Array) mfem::Array<mfem::XXX>;
+%template(##ZZZ##Array) mfem::Array<mfem::XXX>;
 #endif
 %extend mfem::Array<mfem::XXX> {
 
@@ -63,6 +63,11 @@ PyObject * __getitem__(PyObject* param) {
 INSTANTIATE_ARRAY0(XXX, XXX, 0)
 %enddef
 
+%define INSTANTIATE_ARRAY0(XXX, YYY, USEPTR)
+INSTANTIATE_ARRAY2(XXX, YYY, YYY, USEPTR)
+%enddef
+
+
 %define INSTANTIATE_ARRAY_INT
 %template(intArray) mfem::Array<int>;
 %extend mfem::Array<int> {

diff --git a/setup.py b/setup.py
@@ -61,10 +61,10 @@
 repos_sha = {
     # "mfem": "00b2a0705f647e17a1d4ffcb289adca503f28d42", # version 4.5.2
     # "mfem": "962774d5ffa84ceed3bc670e52388250ee028da1",  # version 4.5.2 + distsolve
-    # "mfem": "69fbae732d5279c8d0f42c5430c4fd5656731d00", # version 4.6
+    # "mfem": "69fbae732d5279c8d0f42c5430c4fd5656731d00",  # version 4.6
     # "mfem": "8bb929c2ff86cdf2ee9bb058cc75e59acb07bb94",  # doftrans simplification (Nov. 15. 2023)
-    # after socket connection fix (Nov. 29 2023)
-    "mfem": "4a45c70d1269d293266b77a3a025a9756d10ed8f",
+    # "mfem": "4a45c70d1269d293266b77a3a025a9756d10ed8f",  # socket connection fix (Nov. 29 2023)
+    "mfem": "68873fa4d403c7c94a653c7bc781815ff5b2734d",    # use constrained prolongation operator in ex26, ex26p (Nov. 30 2023)
     "gklib": "a7f8172703cf6e999dd0710eb279bba513da4fec",
     "metis": "94c03a6e2d1860128c2d0675cbbb86ad4f261256", }