updates

Makefile

@@ -170,7 +170,7 @@ docker-mac:
 
 
 prereq:
-	$(PYTHON) -m pip install $(PIP_FLAGS) $(PIP_INSTALL_FLAGS) Cython cython numpy scipy matplotlib pyyaml h5py pybind11 MeshPy tabulate modepy mpi4py pyamg
+	$(PYTHON) -m pip install $(PIP_FLAGS) $(PIP_INSTALL_FLAGS) wheel Cython cython numpy scipy matplotlib pyyaml h5py pybind11 MeshPy tabulate modepy mpi4py pyamg meshio
 	$(PYTHON) -m pip install $(PIP_FLAGS) $(PIP_INSTALL_FLAGS) scikit-sparse
 
 prereq-extra:

base/PyNucleus_base/LinearOperator_{SCALAR}.pxi

@@ -108,7 +108,7 @@ cdef class {SCALAR_label}LinearOperator:
                 raise NotImplementedError('Cannot multiply {} with {}:\n{}'.format(self, x, e))
 
     def __rmul__(self, x):
-        if isinstance(x, {SCALAR}):
+        if isinstance(x, (float, int, {SCALAR})):
             return {SCALAR_label}Multiply_Linear_Operator(self, x)
         else:
             raise NotImplementedError('Cannot multiply with {}'.format(x))
@@ -264,6 +264,9 @@ cdef class {SCALAR_label}LinearOperator:
 
     diagonal = property(fget=get_diagonal, fset=set_diagonal)
 
+    def getMemorySize(self):
+        return -1
+
 
 cdef class {SCALAR_label}TimeStepperLinearOperator({SCALAR_label}LinearOperator):
     def __init__(self,
@@ -356,6 +359,9 @@ cdef class {SCALAR_label}TimeStepperLinearOperator({SCALAR_label}LinearOperator)
         else:
             return super({SCALAR_label}TimeStepperLinearOperator, self).__mul__(x)
 
+    def getMemorySize(self):
+        return self.M.getMemorySize()+self.S.getMemorySize()
+
 
 cdef class {SCALAR_label}Multiply_Linear_Operator({SCALAR_label}LinearOperator):
     def __init__(self,
@@ -419,6 +425,9 @@ cdef class {SCALAR_label}Multiply_Linear_Operator({SCALAR_label}LinearOperator):
     def __repr__(self):
         return '{}*{}'.format(self.factor, self.A)
 
+    def getMemorySize(self):
+        return self.A.getMemorySize()
+
 
 cdef class {SCALAR_label}Product_Linear_Operator({SCALAR_label}LinearOperator):
     def __init__(self,
@@ -495,6 +504,9 @@ cdef class {SCALAR_label}Product_Linear_Operator({SCALAR_label}LinearOperator):
     def __repr__(self):
         return '{}*{}'.format(self.A, self.B)
 
+    def getMemorySize(self):
+        return self.A.getMemorySize()+self.B.getMemorySize()
+
 
 cdef class {SCALAR_label}VectorLinearOperator:
     def __init__(self, int num_rows, int num_columns, int vectorSize):

base/PyNucleus_base/solvers.pyx

@@ -1101,7 +1101,10 @@ cdef class complex_lu_solver(complex_solver):
         if A is not None:
             self.A = A
 
-        if isinstance(self.A, (ComplexLinearOperator, HelmholtzShiftOperator)):
+        if isinstance(self.A, ComplexDense_LinearOperator):
+            from scipy.linalg import lu_factor
+            self.lu, self.perm = lu_factor(self.A.data)
+        elif isinstance(self.A, (ComplexLinearOperator, HelmholtzShiftOperator)):
             from scipy.sparse.linalg import splu
             try:
                 if isinstance(self.A, ComplexSSS_LinearOperator):

base/PyNucleus_base/utilsFem.py

@@ -346,7 +346,7 @@ def loadDictFromHDF5(f):
                 else:
                     params[key] = LinearOperator.HDF5read(f[key])
             elif 'vertices' in f[key] and 'cells' in f[key]:
-                from PyNucleus.fem import meshNd
+                from PyNucleus.fem.mesh import meshNd
                 params[key] = meshNd.HDF5read(f[key])
             else:
                 params[key] = loadDictFromHDF5(f[key])

drivers/testDistOp.py

@@ -344,7 +344,8 @@
                     lcl_dm.inner)
     cg.maxIter = 1000
     u = lcl_dm.zeros()
-    cg(b, u)
+    with d.timer('CG solve'):
+        cg(b, u)
 
     residuals = cg.residuals
     solveGroup = d.addOutputGroup('solve', tested=True, rTol=1e-1)

fem/PyNucleus_fem/DoFMaps.pyx

@@ -798,7 +798,7 @@ cdef class DoFMap:
 
         """
         try:
-            from PyNucleus_nl.kernelsCy import RangedFractionalKernel
+            from PyNucleus_nl.kernelsCy import RangedFractionalKernel, ComplexKernel
 
             if isinstance(kernel, RangedFractionalKernel):
                 from PyNucleus_base.linear_operators import multiIntervalInterpolationOperator
@@ -831,6 +831,26 @@ cdef class DoFMap:
                     return self.scalarDM.assembleNonlocal(kernel, matrixFormat, dm2.scalarDM, returnNearField, **kwargs)
                 else:
                     return self.scalarDM.assembleNonlocal(kernel, matrixFormat, None, returnNearField, **kwargs)
+            elif isinstance(kernel, ComplexKernel):
+                from PyNucleus_nl.nonlocalLaplacian import ComplexnonlocalBuilder
+
+                builder = ComplexnonlocalBuilder(self.mesh, self, kernel, dm2=dm2, **kwargs)
+                if matrixFormat.upper() == 'DENSE':
+                    return builder.getDense()
+                elif matrixFormat.upper() == 'DIAGONAL':
+                    return builder.getDiagonal()
+                elif matrixFormat.upper() == 'SPARSIFIED':
+                    return builder.getDense(trySparsification=True)
+                elif matrixFormat.upper() == 'SPARSE':
+                    return builder.getSparse(returnNearField=returnNearField)
+                elif matrixFormat.upper() == 'H2':
+                    return builder.getH2(returnNearField=returnNearField)
+                elif matrixFormat.upper() == 'H2CORRECTED':
+                    A = builder.getH2FiniteHorizon()
+                    A.setKernel(kernel)
+                    return A
+                else:
+                    raise NotImplementedError('Unknown matrix format: {}'.format(matrixFormat))
             else:
                 from PyNucleus_nl import nonlocalBuilder
 
@@ -1153,31 +1173,57 @@ cdef class DoFMap:
         return y, dm
 
     def augmentWithBoundaryData(self,
-                                const REAL_t[::1] x,
-                                const REAL_t[::1] boundaryData):
+                                x,
+                                boundaryData):
         "Augment the finite element function with boundary data."
         cdef:
             DoFMap dm
-            fe_vector y
-            REAL_t[::1] yy
+            fe_vector yReal
+            REAL_t[::1] xReal, boundaryReal, yyReal
+            complex_fe_vector yComplex
+            COMPLEX_t[::1] xComplex, boundaryComplex, yyComplex
             INDEX_t i, k, dof, dof2, num_cells = self.mesh.num_cells
 
         if isinstance(self, Product_DoFMap):
             dm = Product_DoFMap(type(self.scalarDM)(self.mesh, tag=MAX_INT), self.numComponents)
         else:
             dm = type(self)(self.mesh, tag=MAX_INT)
-        y = dm.empty(dtype=REAL)
-        yy = y
 
-        for i in range(num_cells):
-            for k in range(self.dofs_per_element):
-                dof = self.cell2dof(i, k)
-                dof2 = dm.cell2dof(i, k)
-                if dof >= 0:
-                    yy[dof2] = x[dof]
-                else:
-                    yy[dof2] = boundaryData[-dof-1]
-        return y
+        if ((isinstance(x, fe_vector) and isinstance(boundaryData, fe_vector)) or
+            (isinstance(x, np.ndarray) and x.dtype == REAL) and (isinstance(boundaryData, np.ndarray) and boundaryData.dtype == REAL)):
+
+            xReal = x
+            boundaryReal = boundaryData
+            yReal = dm.empty(dtype=REAL)
+            yyReal = yReal
+
+            for i in range(num_cells):
+                for k in range(self.dofs_per_element):
+                    dof = self.cell2dof(i, k)
+                    dof2 = dm.cell2dof(i, k)
+                    if dof >= 0:
+                        yyReal[dof2] = xReal[dof]
+                    else:
+                        yyReal[dof2] = boundaryReal[-dof-1]
+            return yReal
+        elif ((isinstance(x, complex_fe_vector) and isinstance(boundaryData, complex_fe_vector)) or
+            (isinstance(x, np.ndarray) and x.dtype == REAL) and (isinstance(boundaryData, np.ndarray) and boundaryData.dtype == REAL)):
+            xComplex = x
+            boundaryComplex = boundaryData
+            yComplex = dm.empty(dtype=COMPLEX)
+            yyComplex = yComplex
+
+            for i in range(num_cells):
+                for k in range(self.dofs_per_element):
+                    dof = self.cell2dof(i, k)
+                    dof2 = dm.cell2dof(i, k)
+                    if dof >= 0:
+                        yyComplex[dof2] = xComplex[dof]
+                    else:
+                        yyComplex[dof2] = boundaryComplex[-dof-1]
+            return yComplex
+        else:
+            raise NotImplementedError()
 
     def getFullDoFMap(self, DoFMap complement_dm):
         cdef:

fem/PyNucleus_fem/factories.py

@@ -34,7 +34,7 @@
                          radialIndicator,
                          fractalDiffusivity, expDiffusivity,
                          componentVectorFunction)
-from . lookupFunction import lookupFunction
+from . lookupFunction import lookupFunction, vectorLookupFunction
 
 
 rhsFunSin1D = _rhsFunSin1D()
@@ -96,6 +96,7 @@ def rhsFractional2D_nonPeriodic(s):
 functionFactory.register('x1**3', monomial, params={'exponent': np.array([0., 3., 0.])})
 functionFactory.register('x2**3', monomial, params={'exponent': np.array([0., 0., 3.])})
 functionFactory.register('Lambda', Lambda)
+functionFactory.register('complexLambda', complexLambda)
 functionFactory.register('squareIndicator', squareIndicator)
 functionFactory.register('radialIndicator', radialIndicator)
 functionFactory.register('rhsBoundaryLayer2D', _rhsBoundaryLayer2D)
@@ -110,6 +111,7 @@ def rhsFractional2D_nonPeriodic(s):
 functionFactory.register('inclusionsHong', inclusionsHong)
 functionFactory.register('motorPermeability', motorPermeability)
 functionFactory.register('lookup', lookupFunction)
+functionFactory.register('vectorLookup', vectorLookupFunction)
 functionFactory.register('shiftScaleFunctor', shiftScaleFunctor)
 functionFactory.register('componentVectorFunction', componentVectorFunction, aliases=['vector'])
 
@@ -143,6 +145,7 @@ def __call__(self, mesh, *args, **kwargs):
                     circle, graded_circle, cutoutCircle, twinDisc, dumbbell, wrench,
                     Hshape, ball, rectangle, crossSquare,
                     gradedSquare, gradedBox,
+                    squareWithCircularCutout, boxWithBallCutout,
                     disconnectedInterval, disconnectedDomain,
                     double_graded_interval,
                     simpleFicheraCube, uniformSquare,
@@ -174,8 +177,11 @@ def __call__(self, mesh, *args, **kwargs):
 meshFactory.register('graded_circle', graded_circle, 2, aliases=['gradedCircle'])
 meshFactory.register('discWithInteraction', discWithInteraction, 2)
 meshFactory.register('cutoutCircle', cutoutCircle, 2, aliases=['cutoutDisc'])
+meshFactory.register('squareWithCircularCutout', squareWithCircularCutout, 2)
+meshFactory.register('boxWithBallCutout', boxWithBallCutout, 3, aliases=['boxMinusBall'])
 meshFactory.register('simpleBox', simpleBox, 3, aliases=['unitBox', 'cube', 'unitCube'])
 meshFactory.register('box', box, 3)
+meshFactory.register('ball', ball, 3)
 meshFactory.register('simpleFicheraCube', simpleFicheraCube, 3, aliases=['fichera', 'ficheraCube'])
 meshFactory.register('standardSimplex2D', standardSimplex2D, 2)
 meshFactory.register('standardSimplex3D', standardSimplex3D, 3)

fem/PyNucleus_fem/femCy.pxd

@@ -32,3 +32,34 @@ cdef class multi_function:
         public INDEX_t numInputs, numOutputs
 
     cdef void eval(self, REAL_t[::1] x, REAL_t[::1] y)
+
+
+cdef class simplexComputations:
+    cdef:
+        REAL_t[:, ::1] simplex
+    cdef void setSimplex(self, REAL_t[:, ::1] simplex)
+    cdef REAL_t evalVolume(self)
+    cdef REAL_t evalVolumeGradients(self,
+                                    REAL_t[:, ::1] gradients)
+    cdef REAL_t evalVolumeGradientsInnerProducts(self,
+                                                 REAL_t[:, ::1] gradients,
+                                                 REAL_t[::1] innerProducts)
+    cdef REAL_t evalSimplexVolumeGradientsInnerProducts(self,
+                                                        const REAL_t[:, ::1] simplex,
+                                                        REAL_t[:, ::1] gradients,
+                                                        REAL_t[::1] innerProducts)
+
+
+cdef class simplexComputations1D(simplexComputations):
+    cdef:
+        REAL_t[:, ::1] temp
+
+
+cdef class simplexComputations2D(simplexComputations):
+    cdef:
+        REAL_t[:, ::1] temp
+
+
+cdef class simplexComputations3D(simplexComputations):
+    cdef:
+        REAL_t[:, ::1] temp