Fix bug in sparse assembly, add lookupTwoPoint

Author: cgcgcg

fem/PyNucleus_fem/DoFMaps.pyx

@@ -2264,6 +2264,21 @@ def str2DoFMap(element):
         raise NotImplementedError('Unknown DoFMap: {}'.format(element))
 
 
+def DoFMap2str(dm):
+    if isinstance(dm, P0_DoFMap):
+        return 'P0'
+    elif isinstance(dm, P1_DoFMap):
+        return 'P1'
+    elif isinstance(dm, P2_DoFMap):
+        return 'P2'
+    elif isinstance(dm, P3_DoFMap):
+        return 'P3'
+    elif isinstance(dm, N1e_DoFMap):
+        return 'N1e'
+    else:
+        raise NotImplementedError('Unknown DoFMap: {}'.format(type(dm)))
+
+
 def getAvailableDoFMaps():
     return ['P0', 'P1', 'P2', 'P3', 'N1e']
 

nl/PyNucleus_nl/nonlocalAssembly.pyx

@@ -18,6 +18,7 @@ from PyNucleus_base.ip_norm cimport (ip_distributed_nonoverlapping,
                                      norm_distributed_nonoverlapping)
 from PyNucleus_fem.functions cimport function, constant
 from PyNucleus_fem.DoFMaps cimport P0_DoFMap, P1_DoFMap, Product_DoFMap
+from PyNucleus_fem.DoFMaps import DoFMap2str
 from PyNucleus_fem.meshCy cimport sortEdge, encode_edge, decode_edge
 from PyNucleus_fem.femCy cimport local_matrix_t
 from PyNucleus_fem.femCy import assembleMatrix, mass_1d_sym_scalar_anisotropic, mass_2d_sym_scalar_anisotropic

nl/PyNucleus_nl/nonlocalAssembly_{SCALAR}.pxi

@@ -1105,7 +1105,7 @@ cdef class {SCALAR_label}nonlocalBuilder:
         # We want to capture all element x element interactions.
         # We set up a temporary dofmap and construct a near field wrt that.
         if self.dm2 is None:
-            treeDM = dofmapFactory('P1', self.dm.mesh, -1)
+            treeDM = dofmapFactory(DoFMap2str(self.dm), self.dm.mesh, -1)
         else:
             treeDM = self.dm.combine(self.dm2)
 

nl/PyNucleus_nl/nonlocalProblems.py

@@ -37,7 +37,8 @@
                                  leftRightTwoPoint,
                                  interfaceTwoPoint,
                                  smoothedLeftRightTwoPoint,
-                                 lambdaTwoPoint)
+                                 lambdaTwoPoint,
+                                 lookupTwoPoint)
 from . interactionDomains import (fullSpace,
                                   ball1_retriangulation,
                                   ball1_barycenter,
@@ -103,6 +104,7 @@ def build(self, name, *args, **kwargs):
 twoPointFunctionFactory.register('leftRight', leftRightTwoPoint, aliases=['leftRightTwoPoint'])
 twoPointFunctionFactory.register('interface', interfaceTwoPoint, aliases=['interfaceTwoPoint'])
 twoPointFunctionFactory.register('lambda', lambdaTwoPoint)
+twoPointFunctionFactory.register('lookup', lookupTwoPoint)
 
 interactionFactory = factory()
 interactionFactory.register('fullSpace', fullSpace, aliases=['full'])

nl/PyNucleus_nl/twoPointFunctions.pxd

@@ -8,6 +8,8 @@
 cimport numpy as np
 from PyNucleus_base.myTypes cimport INDEX_t, REAL_t, COMPLEX_t, BOOL_t
 from PyNucleus_fem.functions cimport function, complexFunction
+from PyNucleus_fem.DoFMaps cimport DoFMap, shapeFunction
+from PyNucleus_fem.meshCy cimport cellFinder2
 
 
 include "twoPointFunctions_decl_REAL.pxi"

nl/PyNucleus_nl/twoPointFunctions_decl_{SCALAR}.pxi

@@ -23,6 +23,14 @@ cdef class {SCALAR_label}constantTwoPoint({SCALAR_label}twoPointFunction):
         public {SCALAR}_t value
 
 
+cdef class {SCALAR_label}lookupTwoPoint({SCALAR_label}twoPointFunction):
+    cdef:
+        public DoFMap dm
+        cellFinder2 cellFinder
+        {SCALAR}_t[::1] vals1, vals2
+        {SCALAR}_t[:, ::1] A
+
+
 cdef class {SCALAR_label}parametrizedTwoPointFunction({SCALAR_label}twoPointFunction):
     cdef:
         void *params

nl/PyNucleus_nl/twoPointFunctions_{SCALAR}.pxi

@@ -151,6 +151,101 @@ cdef class {SCALAR_label}constantTwoPoint({SCALAR_label}twoPointFunction):
         return {SCALAR_label}constantTwoPoint, (self.value, )
 
 
+cdef class {SCALAR_label}lookupTwoPoint({SCALAR_label}twoPointFunction):
+    def __init__(self, DoFMap dm, {SCALAR}_t[:, ::1] A, BOOL_t symmetric, cellFinder2 cF=None):
+        super({SCALAR_label}lookupTwoPoint, self).__init__(symmetric, 1)
+        self.dm = dm
+        self.A = A
+        if cF is None:
+            self.cellFinder = cellFinder2(self.dm.mesh)
+        else:
+            self.cellFinder = cF
+        self.vals1 = uninitialized((self.dm.dofs_per_element), dtype={SCALAR})
+        self.vals2 = uninitialized((self.dm.dofs_per_element), dtype={SCALAR})
+
+    cdef void eval(self, REAL_t[::1] x, REAL_t[::1] y, {SCALAR}_t[::1] value):
+        cdef:
+            shapeFunction shapeFun
+            REAL_t val
+            INDEX_t cellNo1, cellNo2, dof1, dof2, k1, k2
+        cellNo1 = self.cellFinder.findCell(x)
+        if cellNo1 == -1:
+            value[0] = 0.
+            return
+        for k1 in range(self.dm.dofs_per_element):
+            dof2 = self.dm.cell2dof(cellNo1, k1)
+            if dof2 >= 0:
+                shapeFun = self.dm.getLocalShapeFunction(k1)
+                shapeFun.evalPtr(&self.cellFinder.bary[0], NULL, &val)
+                self.vals1[k1] = val
+            else:
+                self.vals1[k1] = 0.
+
+        cellNo2 = self.cellFinder.findCell(y)
+        if cellNo2 == -1:
+            value[0] = 0.
+            return
+        for k2 in range(self.dm.dofs_per_element):
+            dof1 = self.dm.cell2dof(cellNo2, k2)
+            if dof1 >= 0:
+                shapeFun = self.dm.getLocalShapeFunction(k2)
+                shapeFun.evalPtr(&self.cellFinder.bary[0], NULL, &val)
+                self.vals2[k2] = val
+            else:
+                self.vals2[k2] = 0.
+
+        value[0] = 0.
+        for k1 in range(self.dm.dofs_per_element):
+            dof1 = self.dm.cell2dof(cellNo1, k1)
+            for k2 in range(self.dm.dofs_per_element):
+                dof2 = self.dm.cell2dof(cellNo2, k2)
+                value[0] += self.vals1[k1]*self.A[dof1, dof2]*self.vals2[k2]
+
+    cdef void evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, {SCALAR}_t* value):
+        cdef:
+            shapeFunction shapeFun
+            REAL_t val
+            INDEX_t cellNo1, cellNo2, dof, dof1, dof2, k1, k2
+        cellNo1 = self.cellFinder.findCellPtr(x)
+        if cellNo1 == -1:
+            value[0] = 0.
+            return
+        for k1 in range(self.dm.dofs_per_element):
+            dof = self.dm.cell2dof(cellNo1, k1)
+            if dof >= 0:
+                shapeFun = self.dm.getLocalShapeFunction(k1)
+                shapeFun.evalPtr(&self.cellFinder.bary[0], NULL, &val)
+                self.vals1[k1] = val
+            else:
+                self.vals1[k1] = 0.
+
+        cellNo2 = self.cellFinder.findCellPtr(y)
+        if cellNo2 == -1:
+            value[0] = 0.
+            return
+        for k2 in range(self.dm.dofs_per_element):
+            dof = self.dm.cell2dof(cellNo2, k2)
+            if dof >= 0:
+                shapeFun = self.dm.getLocalShapeFunction(k2)
+                shapeFun.evalPtr(&self.cellFinder.bary[0], NULL, &val)
+                self.vals2[k2] = val
+            else:
+                self.vals2[k2] = 0.
+
+        value[0] = 0.
+        for k1 in range(self.dm.dofs_per_element):
+            dof1 = self.dm.cell2dof(cellNo1, k1)
+            for k2 in range(self.dm.dofs_per_element):
+                dof2 = self.dm.cell2dof(cellNo2, k2)
+                value[0] += self.vals1[k1]*self.A[dof1, dof2]*self.vals2[k2]
+
+    def __repr__(self):
+        return 'lookup {}'.format(self.dm)
+
+    def __reduce__(self):
+        return {SCALAR_label}lookupTwoPoint, (self.dm, np.array(self.A), self.symmetric)
+
+
 cdef class {SCALAR_label}parametrizedTwoPointFunction({SCALAR_label}twoPointFunction):
     def __init__(self, BOOL_t symmetric, INDEX_t valueSize):
         super({SCALAR_label}parametrizedTwoPointFunction, self).__init__(symmetric, valueSize)