BoundaryNormalCoefficient

examples/ex11p.py

@@ -127,14 +127,15 @@
 #    serial and parallel assembly we extract the corresponding parallel
 #    matrices A and M.
 one = mfem.ConstantCoefficient(1.0)
+mone = mfem.ConstantCoefficient(-1.0)
 
 ess_bdr = mfem.intArray()
 if pmesh.bdr_attributes.Size() != 0:
     ess_bdr.SetSize(pmesh.bdr_attributes.Max())
     ess_bdr.Assign(1)
 
 a = mfem.ParBilinearForm(fespace)
-a.AddDomainIntegrator(mfem.DiffusionIntegrator(one))
+a.AddDomainIntegrator(mfem.DiffusionIntegrator(mone))
 if pmesh.bdr_attributes.Size() == 0:
     # Add a mass term if the mesh has no boundary, e.g. periodic mesh or
     # closed surface.
@@ -180,7 +181,7 @@
     precond = strumpack
 else:
     amg = mfem.HypreBoomerAMG(A)
-    amg.SetPrintLevel(0)
+    amg.SetPrintLevel(1)
     precond = amg
 
 lobpcg = mfem.HypreLOBPCG(MPI.COMM_WORLD)

mfem/_par/dist_solver.i

@@ -26,5 +26,7 @@ import_array();
 %import "pmesh.i"
 %import "solvers.i"
 
+%immutable print_level;
+
 %include "miniapps/common/dist_solver.hpp"
 

mfem/_par/setup.py

@@ -136,9 +136,11 @@ def get_extensions():
             library_dirs.append(strumpacklib)
 
     if add_cuda == '1':
-        from setup_local import cudainc
+        from setup_local import cudainc, cudatoolkit_prefix
         include_dirs.append(cudainc)
-
+        if cudatoolkit_prefix != "":
+            include_dirs.append(os.path.join(cudatoolkit_prefix, 'include'))
+
     if add_libceed == '1':
         from setup_local import libceedinc
         include_dirs.append(libceedinc)

mfem/_par/solvers.i

@@ -34,6 +34,9 @@ import_array();
 %feature("director") mfem::IterativeSolverMonitor;
 %feature("director") mfem::PyIterativeSolver;
 
+// this is to expose PrintLevel structure
+%nestedworkaround;
+
 %include "linalg/solvers.hpp"
 %include "../common/pysolvers.hpp"
 

mfem/_ser/setup.py

@@ -104,8 +104,10 @@ def get_extensions():
                "submesh", "transfermap"]
 
     if add_cuda == '1':
-        from setup_local import cudainc
+        from setup_local import cudainc, cudatoolkit_prefix
         include_dirs.append(cudainc)
+        if cudatoolkit_prefix != "":
+            include_dirs.append(os.path.join(cudatoolkit_prefix, 'include'))
     if add_libceed == '1':
         from setup_local import libceedinc
         include_dirs.append(libceedinc)

mfem/common/pycoefficient.hpp

@@ -79,6 +79,26 @@ class NumbaFunctionBase
   virtual ~NumbaFunctionBase(){}
 };
 
+// Coefficent to return the (normalized) normal vector on boundary
+class VectorBdrNormalCoefficient : public mfem::VectorCoefficient
+{
+public:
+  VectorBdrNormalCoefficient(int dim):VectorCoefficient(dim){}
+  virtual void Eval(mfem::Vector &V,
+		   mfem::ElementTransformation &T,
+		   const mfem::IntegrationPoint &ip){
+    T.SetIntPoint(&ip);
+    mfem::DenseMatrix jac = T.Jacobian();
+    V.SetSize(jac.Height());
+    CalcOrtho(jac, V);
+    double l=0.0;
+    for (int i = 0; i < vdim; i++) {
+      l = l + V[i]*V[i];
+    }
+    V /= pow(l, 0.5);
+ }
+};
+
 class NumbaCoefficientBase
 {
  private:

setup.py

@@ -109,6 +109,7 @@
 enable_cuda = False
 enable_cuda_hypre = False
 cuda_prefix = ''
+cudatoolkit_prefix=''
 cuda_arch = ''
 enable_pumi = False
 pumi_prefix = ''
@@ -921,7 +922,7 @@ def write_setup_local():
               'gslibsinc': os.path.join(gslibs_prefix, 'include'),
               'gslibpinc': os.path.join(gslibp_prefix, 'include'),
               'cxx11flag': cxx11_flag,
-              'build_mfem': '1' if build_mfem else '0'
+              'build_mfem': '1' if build_mfem else '0',
               }
 
     try:
@@ -949,6 +950,8 @@ def add_extra(xxx, inc_sub=None):
         add_extra('strumpack')
     if enable_cuda:
         add_extra('cuda')
+        params['cudatoolkit_prefix'] = cudatoolkit_prefix
+
     if enable_libceed:
         add_extra('libceed')
     if enable_suitesparse:
@@ -1220,6 +1223,10 @@ def print_config():
     print(" verbose : " + ("Yes" if verbose else "No"))
     print(" SWIG : " + swig_command)
 
+    if cuda_prefix != "":
+        print(" CUDA prefix : " + cuda_prefix)
+    if cudatoolkit_prefix != "":
+        print(" CUDAToolKit : " + cudatoolkit_prefix)
     if blas_libraries != "":
         print(" BLAS libraries : " + blas_libraries)
     if lapack_libraries != "":
@@ -1243,6 +1250,7 @@ def configure_install(self):
     global cc_command, cxx_command, mpicc_command, mpicxx_command
     global metis_64
     global enable_cuda, cuda_prefix, enable_cuda_hypre, cuda_arch
+    global cudatoolkit_prefix
     global enable_pumi, pumi_prefix
     global enable_strumpack, strumpack_prefix
     global enable_libceed, libceed_prefix, libceed_only
@@ -1381,6 +1389,7 @@ def configure_install(self):
     if enable_cuda:
         nvcc = find_command('nvcc')
         cuda_prefix = os.path.dirname(os.path.dirname(nvcc))
+        cudatoolkit_prefix = self.cudatoolkit_prefix
 
     if self.CC != '':
         cc_command = self.CC
@@ -1532,6 +1541,7 @@ class Install(_install):
         ('with-cuda', None, 'enable cuda'),
         ('with-cuda-hypre', None, 'enable cuda in hypre'),
         ('cuda-arch=', None, 'set cuda compute capability. Ex if A100, set to 80'),
+        ('cudatookilt-preifx=', None, 'set cudatoolkit prefix if cusolver.h (and others) are in differet location to cuda prefix'),
         ('with-metis64', None, 'use 64bit int in metis'),
         ('with-pumi', None, 'enable pumi (parallel only)'),
         ('pumi-prefix=', None, 'Specify locaiton of pumi'),
@@ -1574,6 +1584,7 @@ def initialize_options(self):
         self.with_cuda = False
         self.with_cuda_hypre = False
         self.cuda_arch = None
+        self.cudatoolkit_prefix = ''
         self.with_metis64 = False
 
         self.with_pumi = False

test/test_bdrnormal.py

@@ -0,0 +1,77 @@
+'''
+ test BondaryNormalCoefficient
+
+'''
+from numba import cfunc, farray, carray
+import os
+from os.path import expanduser, join
+import sys
+import io
+from mfem import path as mfem_path
+import numpy as np
+import numba
+import time
+from math import sqrt
+
+
+if len(sys.argv) > 1 and sys.argv[1] == '-p':
+    import mfem.par as mfem
+    use_parallel = True
+    from mfem.common.mpi_debug import nicePrint
+    from mpi4py import MPI
+    myid = MPI.COMM_WORLD.rank
+
+else:
+    import mfem.ser as mfem
+    use_parallel = False
+    myid = 0
+    nicePrint = print
+
+
+def check_geometry_error(filename):
+    meshfile = expanduser(join("..", 'data', filename))
+    mesh = mfem.Mesh(meshfile)
+
+    print("File: " + filename)
+    print("NE", mesh.GetNE())
+    print("NBE", mesh.GetNBE())
+    norm = mfem.VectorBdrNormalCoefficient(mesh.SpaceDimension())
+
+    @mfem.jit.scalar(dependency=(norm,))
+    def func(ptx, norm):
+        #print("point", ptx, norm)
+        #print("error", np.sum(ptx - norm)**2, sqrt(np.sum(ptx**2)))
+        #print(sqrt(np.sum(ptx**2)), np.sum(norm**2))
+        return sqrt(np.sum(ptx - norm)**2)
+
+    order = 1
+    dim = mesh.Dimension()
+    fec = mfem.H1_FECollection(order, dim)
+    fes = mfem.FiniteElementSpace(mesh, fec, 1)
+
+    gf = mfem.GridFunction(fes)
+    gf.Assign(0.0)
+
+
+    if mesh.Dimension() == 3:
+        idx = mfem.intArray(list(np.unique(mesh.GetBdrAttributeArray())))
+        gf.ProjectBdrCoefficient(func, idx)
+    else:
+        gf.ProjectCoefficient(func)
+
+    print("-- maximum geometry error", np.max(gf.GetDataArray()))
+
+    gfname = filename[:-5]+'-data.gf'
+    gf.Save(gfname)
+    print("-- error is saved as " + gfname)
+
+def run_test():
+    check_geometry_error('sphere-o2.mesh')
+    check_geometry_error('sphere-o3.mesh')
+    check_geometry_error('sphere-o4.mesh')
+    check_geometry_error('sphere-o5.mesh')
+    check_geometry_error('sphere-surface-o3.mesh')
+
+
+if __name__ == '__main__':
+    run_test()

test/test_numba.py

@@ -133,7 +133,6 @@ def check(a, b, msg):
 def run_test():
     #meshfile = expanduser(join(mfem_path, 'data', 'semi_circle.mesh'))
     mesh = mfem.Mesh(3, 3, 3, "TETRAHEDRON")
-    mesh.ReorientTetMesh()
 
     order = 1