updates

base/PyNucleus_base/LinearOperator_{SCALAR}.pxi

@@ -220,6 +220,13 @@ cdef class {SCALAR_label}LinearOperator:
             return Dense_LinearOperator.HDF5read(node)
         elif node.attrs['type'] == 'diagonal':
             return diagonalOperator.HDF5read(node)
+        elif node.attrs['type'] == 'interpolationOperator':
+            return interpolationOperator.HDF5read(node)
+        elif node.attrs['type'] == 'multiIntervalInterpolationOperator':
+            return multiIntervalInterpolationOperator.HDF5read(node)
+        elif node.attrs['type'] == 'h2':
+            from PyNucleus_nl.clusterMethodCy import H2Matrix
+            return H2Matrix.HDF5read(node)
         else:
             raise NotImplementedError(node.attrs['type'])
 

base/PyNucleus_base/linear_operators.pyx

@@ -1351,12 +1351,37 @@ cdef class interpolationOperator(sumMultiplyOperator):
     def __repr__(self):
         return '<%dx%d %s with %d interpolation nodes>' % (self.num_rows, self.num_columns, self.__class__.__name__, self.numInterpolationNodes)
 
+    def HDF5write(self, node):
+        node.attrs['type'] = 'interpolationOperator'
+        node.attrs['left'] = self.left
+        node.attrs['right'] = self.right
+        node.create_dataset('nodes', data=np.array(self.nodes, copy=False))
+        for i in range(len(self.ops)):
+            grp = node.create_group(str(i))
+            self.ops[i].HDF5write(grp)
+
+    @staticmethod
+    def HDF5read(node):
+        left = node.attrs['left']
+        right = node.attrs['right']
+        nodes = np.array(node['nodes'], dtype=REAL)
+        ops = []
+        for i in range(nodes.shape[0]):
+            ops.append(LinearOperator.HDF5read(node[str(i)]))
+        return interpolationOperator(ops, nodes, left, right)
+
+    cpdef void assure_constructed(self):
+        for i in range(len(self.ops)):
+            if isinstance(self.ops[i], delayedConstructionOperator):
+                self.ops[i].assure_constructed()
+
 
 cdef class multiIntervalInterpolationOperator(LinearOperator):
     cdef:
         public list ops
         INDEX_t selected
-        REAL_t left, right
+        readonly REAL_t left
+        readonly REAL_t right
 
     def __init__(self, list intervals, list nodes, list ops):
         shape = ops[0][0].shape
@@ -1371,6 +1396,12 @@ cdef class multiIntervalInterpolationOperator(LinearOperator):
             self.ops.append(interpolationOperator(ops[k], nodes[k], left, right))
         self.selected = -1
 
+    def get(self):
+        if self.selected != -1:
+            return self.ops[self.selected].val
+        else:
+            return np.nan
+
     def set(self, REAL_t val, BOOL_t derivative=False):
         cdef:
             interpolationOperator op
@@ -1434,6 +1465,29 @@ cdef class multiIntervalInterpolationOperator(LinearOperator):
     def isSparse(self):
         return self.getSelectedOp().isSparse()
 
+    def HDF5write(self, node):
+        node.attrs['type'] = 'multiIntervalInterpolationOperator'
+        for i in range(len(self.ops)):
+            grp = node.create_group(str(i))
+            self.ops[i].HDF5write(grp)
+
+    @staticmethod
+    def HDF5read(node):
+        numOps = len(node)
+        ops = []
+        nodes = []
+        intervals = []
+        for i in range(numOps):
+            op = LinearOperator.HDF5read(node[str(i)])
+            ops.append(op.ops)
+            nodes.append(op.nodes)
+            intervals.append((op.left, op.right))
+        return multiIntervalInterpolationOperator(intervals, nodes, ops)
+
+    cpdef void assure_constructed(self):
+        for i in range(len(self.ops)):
+            self.ops[i].assure_constructed()
+
 
 cdef class delayedConstructionOperator(LinearOperator):
     def __init__(self, INDEX_t numRows, INDEX_t numCols):
@@ -1490,3 +1544,7 @@ cdef class delayedConstructionOperator(LinearOperator):
     def isSparse(self):
         self.assure_constructed()
         return self.A.isSparse()
+
+    def HDF5write(self, node):
+        self.assure_constructed()
+        self.A.HDF5write(node)

base/PyNucleus_base/performanceLogger.pyx

@@ -30,7 +30,7 @@ cpdef void endMemRegion(str key):
 
 
 cdef class FakeTimer:
-    def __init__(self):
+    def __init__(self, str key=''):
         pass
 
     cdef void start(self):
@@ -176,7 +176,13 @@ cdef class PLogger(FakePLogger):
         s = ''
         for key in sorted(self.values.keys()):
             if totalsOnly:
-                s += '{}: {} ({} calls)\n'.format(str(key), sum(self.values[key]), len(self.values[key]))
+                try:
+                    s += '{}: {} ({} calls)\n'.format(str(key), sum(self.values[key]), len(self.values[key]))
+                except TypeError:
+                    if len(self.values[key]):
+                        s += str(key) +': ' + self.values[key][0].__repr__() + '\n'
+                    else:
+                        s += str(key) +': ' + self.values[key].__repr__() + '\n'
             else:
                 s += str(key) +': ' + self.values[key].__repr__() + '\n'
         return s

base/PyNucleus_base/tupleDict_{VALUE}.pxi

@@ -264,3 +264,13 @@ cdef class tupleDict{VALUE}:
                     print(i, j, self.indexL[i][j], self.indexL[i][j+1])
                     return False
         return True
+
+    def toDict(self):
+        cdef:
+            INDEX_t e[2]
+            {VALUE_t} val
+            dict d = {}
+        self.startIter()
+        while self.next(e, &val):
+            d[(e[0], e[1])] = val
+        return d

base/PyNucleus_base/utilsFem.py

@@ -102,7 +102,10 @@ def mergeOrdered(a_list, b_list):
             val = data[key]
             # (number of calls, min over calls, mean over calls, med over calls, max over calls)
             # on rank
-            data2[key] = (len(val), np.min(val), np.mean(val), np.median(val), np.max(val))
+            try:
+                data2[key] = (len(val), np.min(val), np.mean(val), np.median(val), np.max(val))
+            except:
+                pass
         data = data2
         # gather data for all ranks
         if self.comm is not None:
@@ -1475,15 +1478,21 @@ def __call__(self):
                 newValue = getattr(self.baseObj, prop)
                 oldValue = self.cached_args.get(prop, None)
                 args.append(newValue)
-                cached_args[prop] = newValue
                 # TODO: keep hash?
                 try:
                     if isinstance(newValue, np.ndarray):
+                        cached_args[prop] = newValue.copy()
                         if (newValue != oldValue).any():
+                            dependencyLogger.log(self.logLevel, 'Values for {} differ: \'{}\' != \'{}\', calling \'{}\''.format(prop, oldValue, newValue, self.fun.__name__))
                             needToBuild = True
+                        else:
+                            dependencyLogger.log(self.logLevel, 'Values for {} are identical: \'{}\' == \'{}\''.format(prop, oldValue, newValue))
                     elif newValue != oldValue:
+                        cached_args[prop] = newValue
                         dependencyLogger.log(self.logLevel, 'Values for {} differ: \'{}\' != \'{}\', calling \'{}\''.format(prop, oldValue, newValue, self.fun.__name__))
                         needToBuild = True
+                    else:
+                        dependencyLogger.log(self.logLevel, 'Values for {} are identical: \'{}\' == \'{}\''.format(prop, oldValue, newValue))
                 except Exception as e:
                     dependencyLogger.log(logging.WARN, 'Cannot compare values {}, {} for property \'{}\', exception {}, force call \'{}\''.format(oldValue, newValue, prop, e, self.fun.__name__))
                     needToBuild = True

docs/example2.rst

@@ -23,7 +23,7 @@ The singularity :math:`\beta` of the kernel depends on the family of kernels:
 
 - fractional type: :math:`\beta(x,y)=d+2s(x,y)`, where :math:`d` is the spatial dimension and :math:`s(x,y)` is the fractional order.
 - constant type :math:`\beta(x,y)=0`
-- peridynamic type :math:`\beta(x,y)=-1`
+- peridynamic type :math:`\beta(x,y)=1`
 
 At present, the only implemented interaction regions are balls in the 2-norm:
 

drivers/runFractional.py

@@ -16,6 +16,7 @@
 
 d = driver(MPI.COMM_WORLD)
 d.add('saveOperators', False)
+d.add('vtkOutput', "")
 p = fractionalLaplacianProblem(d, False)
 discrProblem = discretizedNonlocalProblem(d, p)
 
@@ -60,9 +61,12 @@
     if p.element != 'P0':
         plotDefaults['shading'] = 'gouraud'
 
-if d.startPlot('solution'):
+if p.dim < 3 and d.startPlot('solution'):
     mS.plotSolution()
-if mS.error is not None and d.startPlot('error'):
+if p.dim < 3 and mS.error is not None and d.startPlot('error'):
     mS.error.plot(**plotDefaults)
 
+if d.vtkOutput != "":
+    mS.exportVTK(d.vtkOutput)
+
 d.finish()

drivers/runFractionalHeat.py

@@ -53,9 +53,9 @@
     if p.element != 'P0':
         plotDefaults['shading'] = 'gouraud'
 
-if d.startPlot('solution'):
+if p.dim < 3 and d.startPlot('solution'):
     mS.plotSolution()
-if mS.error is not None and d.startPlot('error'):
+if p.dim < 3 and mS.error is not None and d.startPlot('error'):
     mS.error.plot(**plotDefaults)
 
 d.finish()

drivers/testDistOp.py

@@ -207,8 +207,7 @@
 if d.buildDistributedH2Bcast:
     with d.timer('distributed, bcast build'):
         A_distributedH2Bcast = dm.assembleNonlocal(nPP.kernel, matrixFormat='H2', comm=d.comm,
-                                                   params={'assembleOnRoot': False,
-                                                           'forceUnsymmetric': True})
+                                                   params={'assembleOnRoot': False})
     with d.timer('distributed, bcast matvec'):
         print('Distributed:     ', A_distributedH2Bcast)
         y_distributedH2Bcast = A_distributedH2Bcast*x
@@ -220,7 +219,6 @@
     with d.timer('distributed, halo build'):
         A_distributedH2 = dm.assembleNonlocal(nPP.kernel, matrixFormat='H2', comm=d.comm,
                                               params={'assembleOnRoot': False,
-                                                      'forceUnsymmetric': True,
                                                       'localFarFieldIndexing': True},
                                               PLogger=tm.PLogger)
     t = d.addOutputGroup('TimersH2', timerOutputGroup(driver=d))
@@ -345,11 +343,12 @@
     cg.maxIter = 1000
     u = lcl_dm.zeros()
     with d.timer('CG solve'):
-        cg(b, u)
+        iterCG = cg(b, u)
 
     residuals = cg.residuals
     solveGroup = d.addOutputGroup('solve', tested=True, rTol=1e-1)
     solveGroup.add('residual norm', residuals[-1])
+    solveGroup.add('CG iterations', iterCG)
 
     # pure Neumann condition -> add nullspace components to match analytic solution
     if nPP.boundaryCondition in (NEUMANN, HOMOGENEOUS_NEUMANN) and nPP.analyticSolution is not None:

fem/PyNucleus_fem/DoFMaps.pxd

@@ -19,7 +19,9 @@ cdef class DoFMap:
         public meshBase mesh  #: The underlying mesh
         readonly INDEX_t dim  #: The spatial dimension of the underlying mesh
         BOOL_t reordered
-        public list localShapeFunctions  #: List of local shape functions
+        public list _localShapeFunctions  #: List of local shape functions
+        void** _localShapeFunctionsPtr
+        BOOL_t vectorValued
         public REAL_t[:, ::1] nodes  #: The barycentric coordinates of the DoFs
         public REAL_t[:, ::1] dof_dual
         public INDEX_t num_dofs  #: The number of DoFs of the finite element space
@@ -46,6 +48,8 @@ cdef class DoFMap:
     cpdef void getVertexDoFs(self, INDEX_t[:, ::1] v2d)
     cpdef void resetUsingIndicator(self, function indicator)
     cpdef void resetUsingFEVector(self, REAL_t[::1] ind)
+    cdef shapeFunction getLocalShapeFunction(self, INDEX_t dofNo)
+    cdef vectorShapeFunction getLocalVectorShapeFunction(self, INDEX_t dofNo)
 
 
 cdef class P1_DoFMap(DoFMap):

fem/PyNucleus_fem/DoFMaps.pyx

@@ -10,6 +10,7 @@ cimport numpy as np
 from libc.math cimport isnan
 from PyNucleus_base.myTypes import INDEX, REAL, COMPLEX, BOOL
 from cpython cimport Py_buffer
+from libc.stdlib cimport malloc, free
 from PyNucleus_base.blas cimport assign, assign3, assignScaled, matmat
 from PyNucleus_base.ip_norm cimport vector_t, ip_serial, norm_serial, wrapRealInnerToComplex, wrapRealNormToComplex
 from PyNucleus_base import uninitialized
@@ -330,6 +331,42 @@ cdef class DoFMap:
         ind = self.interpolate(indicator)
         self.resetUsingFEVector(ind)
 
+    @property
+    def localShapeFunctions(self):
+        return self._localShapeFunctions
+
+    @localShapeFunctions.setter
+    def localShapeFunctions(self, list localShapeFunctions):
+        cdef:
+            INDEX_t i
+            shapeFunction sf
+            vectorShapeFunction vsf
+        if isinstance(localShapeFunctions[0], shapeFunction):
+            self.vectorValued = False
+            self._localShapeFunctions = localShapeFunctions
+            self._localShapeFunctionsPtr = <void**>malloc(len(localShapeFunctions)*sizeof(void*))
+            for i in range(len(localShapeFunctions)):
+                sf = self._localShapeFunctions[i]
+                self._localShapeFunctionsPtr[i] = <void*>sf
+        elif isinstance(localShapeFunctions[0], vectorShapeFunction):
+            self.vectorValued = True
+            self._localShapeFunctions = localShapeFunctions
+            self._localShapeFunctionsPtr = <void**>malloc(len(localShapeFunctions)*sizeof(void*))
+            for i in range(len(localShapeFunctions)):
+                vsf = self._localShapeFunctions[i]
+                self._localShapeFunctionsPtr[i] = <void*>vsf
+        else:
+            raise NotImplementedError()
+
+    cdef shapeFunction getLocalShapeFunction(self, INDEX_t dofNo):
+        return <shapeFunction>self._localShapeFunctionsPtr[dofNo]
+
+    cdef vectorShapeFunction getLocalVectorShapeFunction(self, INDEX_t dofNo):
+        return <vectorShapeFunction>self._localShapeFunctionsPtr[dofNo]
+
+    def __del__(self):
+        free(self._localShapeFunctionsPtr)
+
     cpdef void resetUsingFEVector(self, REAL_t[::1] ind):
         cdef:
             INDEX_t[:, ::1] new_dofs = uninitialized((self.mesh.num_cells,
@@ -2506,18 +2543,19 @@ cdef class Product_DoFMap(DoFMap):
         self.num_dofs = self.numComponents*self.scalarDM.num_dofs
         self.num_boundary_dofs = self.numComponents*self.scalarDM.num_boundary_dofs
 
-        self.localShapeFunctions = []
+        localShapeFunctions = []
         self.nodes = uninitialized((self.dofs_per_element, self.mesh.dim+1), dtype=REAL)
         self.dof_dual = np.zeros((self.dofs_per_element, numComponents), dtype=REAL)
         i = 0
         for dofNo in range(self.scalarDM.dofs_per_element):
             for component in range(numComponents):
                 phi = self.scalarDM.localShapeFunctions[dofNo]
-                self.localShapeFunctions.append(productSpaceShapeFunction(phi, numComponents, component, self.mesh.dim))
+                localShapeFunctions.append(productSpaceShapeFunction(phi, numComponents, component, self.mesh.dim))
                 for j in range(dim+1):
                     self.nodes[i, j] = self.scalarDM.nodes[dofNo, j]
                 self.dof_dual[i, component] = 1.
                 i += 1
+        self.localShapeFunctions = localShapeFunctions
 
     def __repr__(self):
         return '({})^{} with {} DoFs and {} boundary DoFs.'.format(type(self.scalarDM).__name__,

fem/PyNucleus_fem/lookupFunction.pyx

@@ -33,7 +33,7 @@ cdef class lookupFunction(function):
         for k in range(self.dm.dofs_per_element):
             dof = self.dm.cell2dof(cellNo, k)
             if dof >= 0:
-                shapeFun = self.dm.localShapeFunctions[k]
+                shapeFun = self.dm.getLocalShapeFunction(k)
                 val += shapeFun.eval(self.cellFinder.bary)*self.u[dof]
         return val
 
@@ -75,7 +75,7 @@ cdef class vectorLookupFunction(vectorFunction):
         for k in range(self.dm.dofs_per_element):
             dof = self.dm.cell2dof(cellNo, k)
             if dof >= 0:
-                shapeFun = self.dm.localShapeFunctions[k]
+                shapeFun = self.dm.getLocalVectorShapeFunction(k)
                 shapeFun.setCell(self.mesh.cells[cellNo, :])
                 shapeFun.eval(self.cellFinder.bary, self.gradients, self.temp)
                 for componentNo in range(self.mesh.dim):

fem/PyNucleus_fem/mesh.py

@@ -2024,14 +2024,14 @@ class mesh0d(meshNd):
 
 
 class mesh1d(meshNd):
-    def plot(self, boundary=None, info=False):
+    def plot(self, vertices=True, boundary=None, info=False):
         import matplotlib.pyplot as plt
         X = np.array([v[0] for v in self.vertices])
         if self.vertices.shape[1] == 1:
             Y = np.zeros_like(X)
             lenX = X.max()-X.min()
             plt.xlim([X.min()-lenX*0.1, X.max()+lenX*0.1])
-            plt.plot(X, Y, 'o-', zorder=1)
+            plt.plot(X, Y, 'o-' if vertices else '-', zorder=1)
         else:
             v = self.vertices_as_array
             c = self.cells_as_array
@@ -2040,7 +2040,8 @@ def plot(self, boundary=None, info=False):
                      [v[c[:, 0], 1],
                       v[c[:, 1], 1]],
                      c='k')
-            plt.scatter(self.vertices_as_array[:, 0], self.vertices_as_array[:, 1])
+            if vertices:
+                plt.scatter(self.vertices_as_array[:, 0], self.vertices_as_array[:, 1])
             lenX = v[:, 0].max()-v[:, 0].min()
             plt.xlim([v[:, 0].min()-lenX*0.1, v[:, 0].max()+lenX*0.1])
             lenY = v[:, 1].max()-v[:, 1].min()