Merge branch 'p4est-2.8.7-upgrade' into master

Oxley p4est 2.8.7 upgrade: lnodes-based numbering, and full MPI support
(interior + boundary + Dirac; 2D/3D; scalar/system/reduced; Dirichlet/Neumann/
Robin; PCG + DIRECT) verified serial-identical on 2 and 4 ranks.

Author: Lutz Gross

.gitignore

@@ -49,3 +49,6 @@ test_*.log
 test_*.txt
 scons_py_tests*.log
 uniform_mesh_ae.silo/removed_code/
+
+# design working notes, kept out of the repo
+oxley-designs/DESIGN.md

escriptcore/py_src/linearPDEs.py

@@ -1574,11 +1574,8 @@ def addToSystem(self, op, rhs, data):
         :type data: `list`
         """
         self.getDomain().addToSystem(op, rhs, data, self.assembler)
-        if self.hasOxley():
-            self.getDomain().makeZ(self.__complex)
-            self.getDomain().makeIZ(self.__complex)
-            self.getDomain().finaliseA(op,self.__complex)
-            rhs=self.getDomain().finaliseRhs(rhs)
+        # oxley now assembles the final (element-condensed) system directly, so
+        # the old Z/IZ post-condensation hook has been removed.
 
    def addPDEToLumpedSystem(self, operator, a, b, c, hrz_lumping):
         """

oxley-designs/simple3D.py

@@ -0,0 +1,48 @@
+from .tools import *
+
+world_comm = MPI.COMM_WORLD
+domain1 = Brick(n0=20, n1=20, n2 = 10, l0=1000, l1=1000, l2=800, origin0=-500, origi1=-500, origin2=-800,  comm=world_comm, framework=None)
+# alternatively and then switch to Brick or Rectangle
+domain1 = Block(n=(20, 20, 10), l0=(1000., 1000., 1000.), origin=(-500, -500, -800.),  comm=world_comm, framework=None)
+
+refiner = Refiner(levels_max = 5)
+refner.add(Sphere(center=(0.,0, -400), radius=200, tagname="Anomaly", resolution=5))
+refiner.add(PlaneInterface(origin=(0.,0, 400), normal=(0,0,1), tagname="Deep", resolution=10))
+domain2=refiner(domain1)
+
+# now we can do things like this:
+# coordinates on all nodes:
+x_cf = ContinousFuction(domain2).getX()
+# identical to x_cf
+x_rcf = ReducedContinousFuction(domain2).getX()
+# on 8 quadrature points
+x_f = Function(domain2).getX()
+# on element center
+x_rf = ReducedFunction(domain2).getX()
+# on 4 quadrature point of faces
+x_fb = FunctionOnBoundary(domain2).getX()
+# on element center
+x_rfb = ReducedFunctionOnBoundary(domain2).getX()
+# coordinates of the DOFs = *nodes that are not hanging nodes*.
+# labeling depends on framework!!!!
+x_s = Solution(domain2).getX()
+# equivalent to x_s
+x_rs = ReducedSolution(domain2).getX()
+
+# interpolation: for instance
+interpolate(x_cf, ContinuousFunction(domain1))
+# etc
+
+# then we want to do things like:
+mypde=SingleLinearPDE(domain2)
+kappa=Scalar(1., Function(domain2))
+kappa.setTaggedValue("Anomaly",10)
+kappa.setTaggedValue("Deep",10)
+
+input=Scalar(0., FunctionOnBoundary(domain2))
+input.setTaggedValue("bottom",10)
+
+mypde.setValue(A=kappa * kronecker(mydomain),y=input, q=whereZero(x_s[2]))
+u=mypde.getSolution()
+# etc.
+

oxley-designs/tools.py

@@ -0,0 +1,132 @@
+from fontTools.misc.cython import returns
+
+
+class RefinementTask(object):
+    def __init__(self, resolution, isinterface=False, newtagname=None):
+        """
+        interface: True if only the interface is to be refined
+        tagname: the tag name to be used for elements below/inside interface
+        resolution: target resolution
+        """
+        assert resolution >0
+        self.tagname = tagname
+        self.oppositetagid = None
+        self.isinterface = isinterface
+        self.resolution = resolution
+
+    def isInterface(self):
+        return self.isinterface
+    def getTagName(self):
+        return self.tagname
+    def setTagId(self, tagid):
+        self.tagid = tagid
+    def getTagId(self):
+        return self.tagid
+    def setOppositeTageId(self, tagid):
+        self.oppositetagid = tagid
+    def getOppositeTageId(self):
+        return self.oppositetagid
+    def getResolution(self):
+        return self.resolution
+    def __call__(self, x):
+        return self.check()
+
+    def check(self, x):
+        """
+        return True if the x is inside/below the interface.
+        """
+        pass
+
+class Sphere(RefinementTask):
+    def __init__(self, center=(0.,0, 0.), radius=1., tagname="Sphere", resolution=1):
+        super().__init__(resolution, isinterface=False, tagname=tagname)
+        self.center = center
+        self.radius = radius
+    def check(self, x):
+        d = math.sqrt((x[0]-self.center[0])**2 + (x[1]-self.center[1])**2+(x[2]-self.center[2])**2)
+        return d <= self.radius
+class PlaneInterface(RefinementTask):
+    def __init__(self, origin=(0.,0, 400), normal=(0,0,1), tagname="Deep", resolution=10)
+        super().__init__(resolution, isinterface=True, tagname=tagname)
+        self.offset = inner(normal, origin)
+        self.normal = normal
+    def check(self, x):
+        return inner(normal, x) <= -self.offset
+
+
+class Refiner(object):
+    def __init__(self, levels_max = 5):
+        """
+        levels_max: maximum number of refinemnet levels
+        """
+        assert levels_max > 0
+        self.levels_max = levels_max
+        self.tasks = []
+    def add(self, tasks):
+        """
+        Add refinement tasks
+        """
+        if isinstance(tasks,list):
+            self.tasks.extend(tasks)
+        else:
+            self.tasks.append(tasks)
+    def getTasks(self):
+        return self.tasks
+
+    def refine(self, domain):
+        tags = [ t.getTagName() for t in self.tasks if t.getTagName() is not None ]
+        tagmap = domain.createTagmap(newtags=tags) -> tagmap[tagname]=tagid (chack for existing tags!)
+        [ t.setTagID(tagmap(t.getTagName()) for t in self.tasks if t.getTagName() is not None]
+
+        step = 0
+        elements0 = domain.elements
+
+        while step < self.levels_max:
+            new_elements = []
+            have_refined = False
+            for e in elements0: iterate over tree (openmp?)
+                splited_e = split(e)
+                tag_e = e.getTagID()
+                refine = False
+                for t in self.tasks:
+                    if e.size <= t.getResolution():
+                        status = [ (e2, t.check(e2.x))  for e2 in splited_e ]
+                        # subelements get tag according to side
+                        if t.getTagName() is not None:
+                            [e2.setTag(tagmap[t.getTagName()]) for e2, s in status if s ]
+                            # if all subelements are inside, change tag of parent:
+                            if all([s[1] for s in status]):
+                                tag_e = tagmap[t.getTagName()]
+                        # if we not deal with an interface then element is split if all
+                        # subelements are on the blow the interface:
+                        if  all([s[1] for s in status]):
+                            if not t.isInterface():
+                                refine = True
+                        # if there is at least one subelement is below interface we split:
+                        elif any([s[1] for s in status]):
+                            refine = True
+                if refine:
+                    have_refined = True
+                    new_elements.extend([s[1] for s in status])
+                else:
+                    new_elements.append(e.copy().setTagID(tag_e))
+            -> do we need to iterate over surface elements??? tags are not updated!!!!
+
+            if have_refined: ->MPI!!!
+                - > update p4tree from new_elements
+                elements0=new_elements
+            else:
+                break
+            step+=1
+
+    create domain and return
+
+
+
+
+
+
+
+refner.add(Sphere(center=(0.,0, -400), radius=200, tagname="Anomaly", resolution=5))
+refiner.add(PlaneInterface(origin=(0.,0, 400), normal=(0,0,1), tagname="Deep", resolution=10))
+domain2=refiner(domain1)

oxley-designs/validate_dirac3d.py

@@ -0,0 +1,16 @@
+# 3D point-source Poisson: dirac point in the domain; check the solve responds
+from esys.escript import Solution, DiracDeltaFunctions, Scalar, whereZero, kronecker, Lsup, sup, integrate, Function
+from esys.escript.linearPDEs import LinearSinglePDE, SolverOptions
+import esys.oxley as oxley, esys.ripley as ripley
+def run(name, mk):
+    d = mk([(0.5,0.5,0.5)], ["src"])
+    x = Solution(d).getX()
+    pde = LinearSinglePDE(d)
+    y_dirac = Scalar(0., DiracDeltaFunctions(d)); y_dirac.setTaggedValue("src", 1.0)
+    q = whereZero(x[0])+whereZero(x[0]-1)+whereZero(x[1])+whereZero(x[1]-1)+whereZero(x[2])+whereZero(x[2]-1)
+    pde.setValue(A=kronecker(d), y_dirac=y_dirac, q=q)
+    pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
+    u = pde.getSolution()
+    print("  %-8s sup(u)=%.6f  int(u)=%.6f"%(name, sup(u), integrate(u,Function(d))))
+run("oxley",  lambda p,t: oxley.Brick(n0=8,n1=8,n2=8, diracPoints=p, diracTags=t))
+run("ripley", lambda p,t: ripley.Brick(8,8,8, diracPoints=p, diracTags=t))

oxley-designs/validate_elasticity3d.py

@@ -0,0 +1,34 @@
+# 3D linear elasticity: -div(C:grad u) = 0, Dirichlet u; A couples components.
+from esys.escript import Solution, sup, inf, whereZero, kronecker, Lsup, Function, integrate
+from esys.escript.linearPDEs import LinearPDE, SolverOptions
+import esys.oxley as oxley, esys.ripley as ripley
+def solve(dom):
+    DIM=3
+    lam, mu = 1.0, 1.0
+    pde = LinearPDE(dom, numEquations=DIM, numSolutions=DIM)
+    A = pde.createCoefficient("A")   # A[i,j,k,l] = lam d_ij d_kl + mu(d_ik d_jl + d_il d_jk)
+    for i in range(DIM):
+        for j in range(DIM):
+            for k in range(DIM):
+                for l in range(DIM):
+                    v=0.
+                    if i==j and k==l: v+=lam
+                    if i==k and j==l: v+=mu
+                    if i==l and j==k: v+=mu
+                    A[i,j,k,l]=v
+    x = Solution(dom).getX()
+    # Dirichlet: u = (x0, 0, 0)*0.1 on all boundaries (rigid stretch), interior solves
+    onb = whereZero(x[0])+whereZero(x[0]-1)+whereZero(x[1])+whereZero(x[1]-1)+whereZero(x[2])+whereZero(x[2]-1)
+    q = pde.createCoefficient("q"); r = pde.createCoefficient("r")
+    for i in range(DIM): q[i]=onb
+    r[0]=0.1*x[0]
+    pde.setValue(A=A, q=q, r=r)
+    pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
+    return pde.getSolution()
+for n,dom in [("ripley3D",ripley.Brick(8,8,8)),("oxley3D",oxley.Block(numBlocks=(2,2,2),refine_level=2))]:
+    try:
+        u=solve(dom)
+        print("  %-9s u0[%.6f,%.6f] u1[%.6f,%.6f] u2[%.6f,%.6f] |u|int=%.6f"%(n,
+            inf(u[0]),sup(u[0]),inf(u[1]),sup(u[1]),inf(u[2]),sup(u[2]), integrate(Lsup(u)*0+u[0],Function(dom))))
+    except Exception as ex:
+        print("  %-9s FAIL %s"%(n,repr(ex)[:90]))

oxley-designs/validate_meshaccess.py

@@ -0,0 +1,114 @@
+##############################################################################
+#
+# Copyright (c) 2003-2026 by the esys.escript Group
+# https://github.com/LutzGross/esys-escript.github.io
+#
+# Primary Business: Queensland, Australia
+# Licensed under the Apache License, version 2.0
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# See CREDITS file for contributors and development history
+#
+##############################################################################
+"""
+Validation of the oxley lnodes mesh-access interface (design milestone A2).
+
+For a uniform, conforming Block domain the mesh is a regular grid, so the
+lnodes-based getMeshInfo() output can be checked analytically:
+
+  N_i          = numBlocks_i * 2**refine_level        (elements per axis)
+  numElements  = prod(N_i)
+  numNodes     = prod(N_i + 1)                         (shared nodes deduplicated)
+  nodeCoords   = the regular grid points
+  connectivity = a consistent coordinate<->index bijection
+
+The node-count / dedup checks are what distinguish correct lnodes numbering from
+the old floating-point coordinate hashing.
+"""
+import itertools
+import sys
+import numpy as np
+from esys.oxley import Block
+
+
+def _check(ok_list, name, cond):
+    print(("  PASS  " if cond else "  FAIL  ") + name)
+    ok_list.append(bool(cond))
+
+
+def _coord_set(a):
+    return set(tuple(r) for r in np.round(a, 9))
+
+
+def validate(numBlocks, length, origin, refine_level):
+    dim = len(numBlocks)
+    N = [numBlocks[i] * (2 ** refine_level) for i in range(dim)]
+    h = [length[i] / N[i] for i in range(dim)]
+    nc = 4 if dim == 2 else 8
+    exp_nelem = int(np.prod(N))
+    exp_nnode = int(np.prod([n + 1 for n in N]))
+    print("Block numBlocks=%s length=%s origin=%s refine_level=%d -> %s grid, "
+          "%d elements, %d nodes" % (numBlocks, length, origin, refine_level,
+          "x".join(map(str, N)), exp_nelem, exp_nnode))
+
+    dom = Block(numBlocks=numBlocks, length=length, origin=origin,
+                refine_level=refine_level)
+    info = dom.getMeshInfo()
+    coords = info["nodeCoords"]
+    conn = info["elementNodes"]
+
+    ok = []
+    _check(ok, "numDim", info["numDim"] == dim)
+    _check(ok, "nodesPerElement == %d" % nc, info["nodesPerElement"] == nc)
+    _check(ok, "numElements == prod(N) = %d" % exp_nelem,
+           info["numElements"] == exp_nelem)
+    _check(ok, "numNodes == prod(N+1) = %d  (dedup)" % exp_nnode,
+           info["numNodes"] == exp_nnode)
+    _check(ok, "nodeCoords shape (%d,%d)" % (exp_nnode, dim),
+           coords.shape == (exp_nnode, dim))
+    _check(ok, "elementNodes shape (%d,%d)" % (exp_nelem, nc),
+           conn.shape == (exp_nelem, nc))
+    _check(ok, "serial global ids are identity",
+           np.array_equal(info["nodeGlobalId"], np.arange(exp_nnode)))
+    _check(ok, "connectivity indices in [0,numNodes)",
+           conn.size and conn.min() >= 0 and conn.max() < exp_nnode)
+    _check(ok, "every node referenced by an element",
+           len(np.unique(conn)) == exp_nnode)
+
+    for d in range(dim):
+        lo, hi = float(coords[:, d].min()), float(coords[:, d].max())
+        _check(ok, "axis %d span == [%g,%g]" % (d, origin[d], origin[d] + length[d]),
+               np.isclose(lo, origin[d]) and np.isclose(hi, origin[d] + length[d]))
+
+    axes = [origin[d] + h[d] * np.arange(N[d] + 1) for d in range(dim)]
+    grid = np.array(list(itertools.product(*axes)))
+    _check(ok, "node coordinates == regular grid", _coord_set(coords) == _coord_set(grid))
+    _check(ok, "no duplicate node coordinates",
+           np.unique(np.round(coords, 9), axis=0).shape[0] == exp_nnode)
+
+    cells = coords[conn]                       # (nelem, nc, dim)
+    sizes = cells.max(axis=1) - cells.min(axis=1)
+    _check(ok, "every element is an axis-aligned cell of size h",
+           np.allclose(sizes, h))
+
+    _check(ok, "element tags default to 0",
+           np.array_equal(info["elementTags"], np.zeros(exp_nelem, dtype=info["elementTags"].dtype)))
+
+    print("  => %s\n" % ("OK" if all(ok) else "FAILED"))
+    return all(ok)
+
+
+CASES = [
+    ((2, 2),    (1., 1.),      (0., 0.),      0),
+    ((2, 2),    (1., 1.),      (0., 0.),      2),
+    ((3, 4),    (6., 8.),      (-1., -2.),    1),
+    ((1, 1),    (2., 3.),      (5., 5.),      3),
+    ((2, 2, 2), (1., 1., 1.),  (0., 0., 0.),  1),
+    ((3, 2, 2), (3., 2., 2.),  (-1., 0., -5.), 1),
+    ((1, 1, 1), (1., 1., 1.),  (0., 0., 0.),  2),
+]
+
+if __name__ == "__main__":
+    allok = all(validate(*c) for c in CASES)
+    print("ALL PASSED" if allok else "SOME CASES FAILED")
+    sys.exit(0 if allok else 1)