Source code for nutils.element

# Copyright (c) 2014 Evalf
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

"""
The element module defines reference elements such as the :class:`QuadElement`
and :class:`TriangularElement`, but also more exotic objects like the
:class:`TrimmedElement`. A set of (interconnected) elements together form a
:mod:`nutils.topology`. Elements have edges and children (for refinement), which
are in turn elements and map onto self by an affine transformation. They also
have a well defined reference coordinate system, and provide pointsets for
purposes of integration and sampling.
"""

from . import log, util, numpy, config, numeric, function, cache, transform, warnings, _
import re, math, itertools, operator, functools


## ELEMENT

[docs]class Element: 'element class' __slots__ = 'reference', 'transform', 'opposite' def __init__(self, reference, trans, opptrans=None, oriented=False): assert isinstance(reference, Reference) trans = transform.canonical(trans) assert trans[-1].fromdims == reference.ndims and trans[0].todims == None if opptrans is not None: opptrans = transform.canonical(opptrans) assert opptrans[-1].fromdims == reference.ndims and opptrans[0].todims == None if not oriented: vtx1 = transform.apply(trans, reference.vertices) if vtx1 != transform.apply(opptrans, reference.vertices): for ptrans in reference.permutation_transforms: if vtx1 == transform.apply(opptrans + (ptrans,), reference.vertices): opptrans += ptrans, break else: raise Exception('Did not find a conforming permutation for the opposing transformation') self.reference = reference self.transform = trans self.opposite = opptrans or trans def withopposite(self, opp, oriented=False): if isinstance(opp, tuple): return Element(self.reference, self.transform, opp, oriented) assert isinstance(opp, Element) and opp.reference == self.reference return Element(self.reference, self.transform, opp.transform, oriented or opp.opposite==self.transform) def __mul__(self, other): self_is_iface = self.opposite != self.transform other_is_iface = other.opposite != other.transform trans = transform.Bifurcate(self.transform, other.transform), if self_is_iface != other_is_iface: opptrans = transform.Bifurcate(self.opposite, other.opposite), else: opptrans = None return Element(self.reference * other.reference, trans, opptrans, oriented=True) def __getnewargs__(self): return self.reference, self.transform, self.opposite, True def __hash__(self): return hash((self.reference, self.transform, self.opposite)) def __eq__(self, other): return self is other or isinstance(other,Element) \ and self.reference == other.reference \ and self.transform == other.transform \ and self.opposite == other.opposite @property def vertices(self): return transform.apply(self.transform, self.reference.vertices) @property def ndims(self): return self.reference.ndims @property def nverts(self): return self.reference.nverts @property def nedges(self): return self.reference.nedges @property def edges(self): return [self.edge(i) for i in range(self.nedges)] def edge(self, iedge): trans, edge = self.reference.edges[iedge] return Element(edge, self.transform + (trans,), self.opposite and self.opposite + (trans,), oriented=True) if edge else None @property def children(self): return [Element(child, self.transform + (trans,), self.opposite and self.opposite + (trans,), oriented=True) for trans, child in self.reference.children if child] @property def flipped(self): assert self.opposite, 'element does not define an opposite' return Element(self.reference, self.opposite, self.transform, oriented=True) @property def simplices(self): return [Element(reference, self.transform + (trans,), self.opposite and self.opposite + (trans,), oriented=True) for trans, reference in self.reference.simplices] def __str__(self): return 'Element({})'.format(self.vertices)
## REFERENCE ELEMENTS
[docs]class Reference(cache.Immutable): 'reference element' def __init__(self, ndims:int): self.ndims = ndims @property def nverts(self): return len(self.vertices) __and__ = lambda self, other: self if self == other else NotImplemented __or__ = lambda self, other: self if self == other else NotImplemented __rand__ = lambda self, other: self.__and__(other) __ror__ = lambda self, other: self.__or__(other) __sub__ = __rsub__ = lambda self, other: self.empty if self == other else NotImplemented __bool__ = __nonzero__ = lambda self: bool(self.volume) @property def empty(self): return EmptyReference(self.ndims) def __mul__(self, other): assert isinstance(other, Reference) return TensorReference(self, other) def __pow__(self, n): assert numeric.isint(n) and n >= 0 return getsimplex(0) if n == 0 \ else self if n == 1 \ else self * self**(n-1) @property def nedges(self): return len(self.edge_transforms) @property def nchildren(self): return len(self.child_transforms) @property def edges(self): return list(zip(self.edge_transforms, self.edge_refs)) @property def children(self): return list(zip(self.child_transforms, self.child_refs)) @cache.property def connectivity(self): # Nested tuple with connectivity information about edges of children: # connectivity[ichild][iedge] = ioppchild (interface) or -1 (boundary). childmap = [[-1] * child.nedges for child in self.child_refs] vmap = {} for ichild, (ctrans, cref) in enumerate(self.children): for iedge, etrans in enumerate(cref.edge_transforms): v = ctrans * etrans try: jchild, jedge = vmap.pop(v.flipped) except KeyError: vmap[v] = ichild, iedge else: childmap[jchild][jedge] = ichild childmap[ichild][iedge] = jchild for etrans, eref in self.edges: for ctrans in eref.child_transforms: vmap.pop(etrans * ctrans, None) assert not any(self.child_refs[ichild].edge_refs[iedge] for ichild, iedge in vmap.values()), 'not all boundary elements recovered' return tuple(map(tuple, childmap)) @cache.property def ribbons(self): # tuples of (iedge1,jedge1), (iedge2,jedge2) pairs assert self.ndims >= 2 transforms = {} ribbons = [] for iedge1, (etrans1,edge1) in enumerate(self.edges): if edge1: for iedge2, (etrans2,edge2) in enumerate(edge1.edges): if edge2: key = tuple(sorted(tuple(p) for p in (etrans1 * etrans2).apply(edge2.vertices))) try: jedge1, jedge2 = transforms.pop(key) except KeyError: transforms[key] = iedge1, iedge2 else: assert self.edge_refs[jedge1].edge_refs[jedge2] == edge2 ribbons.append(((iedge1,iedge2), (jedge1,jedge2))) assert not transforms return tuple(ribbons) permutation_transforms = () def getischeme(self, ischeme): match = re.match('([a-zA-Z]+)(.*)', ischeme) assert match, 'cannot parse integration scheme {!r}'.format(ischeme) ptype, args = match.groups() get = getattr(self, 'getischeme_'+ptype) ipoints, iweights = get(eval(args)) if args else get() return numeric.const(ipoints, copy=False), numeric.const(iweights, copy=False) if iweights is not None else None @classmethod def register(cls, ptype, func): setattr(cls, 'getischeme_{}'.format(ptype), func) def with_children(self, child_refs): child_refs = tuple(child_refs) if not any(child_refs): return self.empty if child_refs == self.child_refs: return self return WithChildrenReference(self, child_refs) @cache.property def volume(self): ipoints, iweights = self.getischeme('gauss{}'.format(1)) return iweights.sum() @cache.property def centroid(self): ipoints, iweights = self.getischeme('gauss{}'.format(1)) return ipoints.T.dot(iweights)/iweights.sum()
[docs] def trim(self, levels, maxrefine, ndivisions): 'trim element along levelset' assert len(levels) == self.nvertices_by_level(maxrefine) return self if not self or numpy.greater_equal(levels, 0).all() \ else self.empty if numpy.less_equal(levels, 0).all() \ else self.with_children(cref.trim(clevels, maxrefine-1, ndivisions) for cref, clevels in zip(self.child_refs, self.child_divide(levels,maxrefine))) if maxrefine > 0 \ else self.slice(lambda vertices: numeric.dot(numeric.poly_eval(self._linear_bernstein[_], vertices), levels), ndivisions)
@cache.property def _linear_bernstein(self): return self.get_poly_coeffs('bernstein', degree=1) def slice(self, levelfunc, ndivisions): # slice along levelset by recursing over dimensions levels = levelfunc(self.vertices) assert numeric.isint(ndivisions) assert len(levels) == self.nverts if numpy.greater_equal(levels, 0).all(): return self if numpy.less_equal(levels, 0).all(): return self.empty nbins = 2**ndivisions if self.ndims == 1: l0, l1 = levels xi = numpy.round(l0/(l0-l1) * nbins) if xi in (0,nbins): return self.empty if xi == 0 and l1 < 0 or xi == nbins and l0 < 0 else self v0, v1 = self.vertices midpoint = v0 + (xi/nbins) * (v1-v0) refs = [edgeref if levelfunc(edgetrans.apply(numpy.zeros((1,0)))) > 0 else edgeref.empty for edgetrans, edgeref in self.edges] else: refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges] if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) <= 1: return self if sum(bool(ref) for ref in refs) <= 1: return self.empty clevel = levelfunc(self.centroid[_])[0] select = clevel*levels<=0 if clevel!=0 else levels!=0 levels = levels[select] vertices = self.vertices[select] xi = numpy.round(levels/(levels-clevel) * nbins) midpoint = numpy.mean(vertices + (self.centroid-vertices)*(xi/nbins)[:,_], axis=0) if tuple(refs) == tuple(self.edge_refs): return self if not any(refs): return self.empty mosaic = MosaicReference(self, refs, midpoint) return self.empty if mosaic.volume == 0 else mosaic if mosaic.volume < self.volume else self def cone(self, trans, tip): return Cone(self, trans, tip) def check_edges(self, tol=1e-15, print=print): volume = 0 zero = 0 for trans, edge in self.edges: if edge: xe, we = edge.getischeme('gauss1') w_normal = we[:,_] * trans.ext zero += w_normal.sum(0) volume += numeric.contract(trans.apply(xe), w_normal, axis=0) if numpy.greater(abs(zero), tol).any(): print('divergence check failed: {} != 0'.format(zero)) if numpy.greater(abs(volume - self.volume), tol).any(): print('divergence check failed: {} != {}'.format(volume, self.volume)) def vertex_cover(self, ctransforms, maxrefine): if maxrefine < 0: raise Exception('maxrefine is too low') npoints = self.nvertices_by_level(maxrefine) allindices = numpy.arange(npoints) if len(ctransforms) == 1: assert not ctransforms[0] return ((), self.getischeme('vertex{}'.format(maxrefine))[0], allindices), if maxrefine == 0: raise Exception('maxrefine is too low') cbins = [[] for ichild in range(self.nchildren)] for ctrans in ctransforms: ichild = self.child_transforms.index(ctrans[0]) cbins[ichild].append(ctrans[1:]) if not all(cbins): raise Exception('transformations to not form an element cover') fcache = cache.WrapperCache() return tuple(((ctrans,) + trans, points, cindices[indices]) for ctrans, cref, cbin, cindices in zip(self.child_transforms, self.child_refs, cbins, self.child_divide(allindices,maxrefine)) for trans, points, indices in fcache[cref.vertex_cover](tuple(sorted(cbin)), maxrefine-1)) def __str__(self): return self.__class__.__name__ __repr__ = __str__ def get_ndofs(self, degree): raise NotImplementedError def get_poly_coeffs(self, basis, **kwargs): raise NotImplementedError def get_edge_dofs(self, degree, iedge): raise NotImplementedError def get_dof_transpose_map(self, degree, vertex_transpose_map): raise NotImplementedError
[docs]class EmptyReference(Reference): 'inverse reference element' volume = 0 edge_transforms = () edge_refs = () child_transforms = () child_refs = () @property def vertices(self): return numeric.const(numpy.zeros((0, self.ndims)), copy=False) __and__ = __sub__ = lambda self, other: self if other.ndims == self.ndims else NotImplemented __or__ = lambda self, other: other if other.ndims == self.ndims else NotImplemented __rsub__ = lambda self, other: other if other.ndims == self.ndims else NotImplemented
[docs] def trim(self, levels, maxrefine, ndivisions): return self
def inside(self, point, eps=0): return False
[docs]class RevolutionReference(Reference): 'modify gauss integration to always return a single point' def __init__(self): super().__init__(ndims=1) def vertices(self): return numeric.const([[0.]]) # NOTE unclear if this is the desired outcome @property def edge_transforms(self): # only used in check_edges return transform.Updim(numpy.zeros((1,0)), [-numpy.pi], isflipped=True), transform.Updim(numpy.zeros((1,0)), [+numpy.pi], isflipped=False) @property def edge_refs(self): # idem edge_transforms return PointReference(), PointReference() @property def simplices(self): return (transform.Identity(self.ndims), self), def getischeme(self, ischeme): return numeric.const([[0.]]), numeric.const([2 * numpy.pi]) def inside(self, point, eps=0): return True
[docs]class SimplexReference(Reference): 'simplex reference' @property def vertices(self): return numeric.const(numpy.concatenate([numpy.zeros(self.ndims)[_,:], numpy.eye(self.ndims)], axis=0), copy=False) @cache.property def edge_refs(self): assert self.ndims > 0 return (getsimplex(self.ndims-1),) * (self.ndims+1) @cache.property def edge_transforms(self): assert self.ndims > 0 return tuple(transform.SimplexEdge(self.ndims, i) for i in range(self.ndims+1)) @property def child_refs(self): return tuple([self] * (2**self.ndims)) @property def child_transforms(self): return tuple(transform.SimplexChild(self.ndims, ichild) for ichild in range(2**self.ndims)) @cache.property def permutation_transforms(self): transforms = [] for verts in itertools.permutations(tuple(v for v in self.vertices)): offset = verts[0] linear = verts[1:]-verts[0] transforms.append(transform.Square(linear.T, offset)) return tuple(transforms) @cache.property def ribbons(self): return tuple(((iedge1,iedge2),(iedge2+1,iedge1)) for iedge1 in range(self.ndims+1) for iedge2 in range(iedge1,self.ndims)) def getischeme_vtk(self): return self.vertices, None def getischeme_vertex(self, n=0): if n == 0: return self.vertices, None return self.getischeme_bezier(2**n+1) @property def simplices(self): return (transform.Identity(self.ndims), self), def get_ndofs(self, degree): prod = lambda start, stop: functools.reduce(operator.mul, range(start, stop), 1) return prod(degree+1, degree+1+self.ndims) // prod(1, self.ndims+1) def get_poly_coeffs(self, basis, **kwargs): f = getattr(self, '_get_poly_coeffs_{}'.format(basis), None) if f: return f(**kwargs) else: raise ValueError('basis {!r} undefined on {}'.format(basis, type(self).__qualname__)) def _integer_barycentric_coordinates(self, degree): return ( (degree-sum(i),*i[::-1]) for i in itertools.product(*[range(degree+1)]*self.ndims) if sum(i) <= degree) def _get_poly_coeffs_bernstein(self, degree): ndofs = self.get_ndofs(degree) if self.ndims == 0: return numpy.ones((ndofs,), dtype=int) coeffs = numpy.zeros((ndofs,)+(degree+1,)*self.ndims, dtype=int) for i, p in enumerate(self._integer_barycentric_coordinates(degree)): p = p[1:] for q in itertools.product(*[range(degree+1)]*self.ndims): if sum(p+q) <= degree: coeffs[(i,)+tuple(map(operator.add, p, q))] = (-1)**sum(q)*math.factorial(degree)//(math.factorial(degree-sum(p+q))*util.product(map(math.factorial, p+q))) assert i == ndofs - 1 return numeric.const(coeffs, copy=False) def _get_poly_coeffs_lagrange(self, degree): if self.ndims == 0: coeffs = numpy.ones((1,)) elif degree == 0: coeffs = numpy.ones((1,*[1]*self.ndims)) else: P = numpy.array(tuple(self._integer_barycentric_coordinates(degree)), dtype=int)[:,1:] coeffs_ = numpy.linalg.inv(((P[:,_,:]/degree)**P[_,:,:]).prod(-1)) coeffs = numpy.zeros((len(P),*[degree+1]*self.ndims), dtype=float) for i, p in enumerate(P): coeffs[(slice(None),*p)] = coeffs_[i] return numeric.const(coeffs, copy=False) def get_edge_dofs(self, degree, iedge): return numeric.const(tuple(i for i, j in enumerate(self._integer_barycentric_coordinates(degree)) if j[iedge] == 0), dtype=int) def get_dof_transpose_map(self, degree, vertex_transpose_map): vertex_transpose_map = tuple(vertex_transpose_map) if len(vertex_transpose_map) != self.nverts or set(vertex_transpose_map) != set(range(self.nverts)): raise ValueError('invalid vertex indices: {!r}'.format(vertex_transpose_map)) return numeric.const(tuple(i for i, j in sorted(enumerate(self._integer_barycentric_coordinates(degree)), key=lambda ij: tuple(map(ij[1].__getitem__, vertex_transpose_map[::-1])))), dtype=int)
[docs]class PointReference(SimplexReference): '0D simplex' def __init__(self): super().__init__(ndims=0) def getischeme(self, ischeme): return numeric.const(numpy.empty([1,0])), numeric.const([1.]) def inside(self, point, eps=0): return True def nvertices_by_level(self, n): return 1 def child_divide(self, vals, n): return vals,
[docs]class LineReference(SimplexReference): '1D simplex' def __init__(self): self._bernsteincache = [] # TEMPORARY super().__init__(ndims=1) def getischeme_gauss(self, degree): assert isinstance(degree, int) and degree >= 0 x, w = gauss(degree) return x[:,_], w def getischeme_uniform(self, n): return (numpy.arange(.5,n) / n)[:,_], numeric.const.full([n], 1/n) def getischeme_bezier(self, np): return numpy.linspace(0, 1, np)[:,_], None def nvertices_by_level(self, n): return 2**n + 1 def child_divide(self, vals, n): assert n > 0 assert len(vals) == self.nvertices_by_level(n) m = (len(vals)+1) // 2 return vals[:m], vals[m-1:] def inside(self, point, eps=0): x, = point return -eps <= x <= 1+eps
[docs]class TriangleReference(SimplexReference): '2D simplex' def __init__(self): super().__init__(ndims=2)
[docs] def getischeme_gauss(self, degree): '''get integration scheme http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf''' if isinstance(degree, tuple): assert len(degree) == self.ndims degree = sum(degree) assert isinstance(degree, int) and degree >= 0 I = [0,0], J = [1,1],[0,1],[1,0] K = [1,2],[2,0],[0,1],[2,1],[1,0],[0,2] icw = [ (I, [1/3], 1) ] if degree <= 1 else [ (J, [2/3,1/6], 1/3) ] if degree == 2 else [ (I, [1/3], -9/16), (J, [3/5,1/5], 25/48), ] if degree == 3 else [ (J, [0.816847572980458,0.091576213509771], 0.109951743655322), (J, [0.108103018168070,0.445948490915965], 0.223381589678011), ] if degree == 4 else [ (I, [1/3], 0.225), (J, [0.797426985353088,0.101286507323456], 0.125939180544827), (J, [0.059715871789770,0.470142064105115], 0.132394152788506), ] if degree == 5 else [ (J, [0.873821971016996,0.063089014491502], 0.050844906370207), (J, [0.501426509658180,0.249286745170910], 0.116786275726379), (K, [0.636502499121399,0.310352451033785,0.053145049844816], 0.082851075618374), ] if degree == 6 else [ (I, [1/3.], -0.149570044467671), (J, [0.479308067841924,0.260345966079038], 0.175615257433204), (J, [0.869739794195568,0.065130102902216], 0.053347235608839), (K, [0.638444188569809,0.312865496004875,0.048690315425316], 0.077113760890257), ] if degree > 7: warnings.warn('inexact integration for polynomial of degree {}'.format(degree)) return numpy.concatenate([numpy.take(c,i) for i, c, w in icw], axis=0), \ numpy.concatenate([[w*.5] * len(i) for i, c, w in icw])
def getischeme_uniform(self, n): points = numpy.arange(1./3, n) / n nn = n**2 C = numpy.empty([2,n,n]) C[0] = points[:,_] C[1] = points[_,:] coords = C.reshape(2, nn) flip = coords.sum(0) > 1 coords[:,flip] = 1 - coords[::-1,flip] weights = numeric.const.full([nn], .5/nn) return coords.T, weights def getischeme_bezier(self, np): points = numpy.linspace(0, 1, np) return numpy.array([[x,y] for i, y in enumerate(points) for x in points[:np-i]]), None def nvertices_by_level(self, n): m = 2**n + 1 return ((m+1)*m) // 2 def child_divide(self, vals, n): assert len(vals) == self.nvertices_by_level(n) np = 1 + 2**n # points along parent edge mp = 1 + 2**(n-1) # points along child edge cvals = [] for i in range(mp): j = numpy.arange(mp-i) cvals.append([vals[b+a*np-(a*(a-1))//2] for a, b in [(i,j),(i,mp-1+j),(mp-1+i,j),(i+j,mp-1-j)]]) return numpy.concatenate(cvals, axis=1) def inside(self, point, eps=0): x, y = point return x >= -eps and y >= -eps and 1-x-y >= -eps
[docs]class TetrahedronReference(SimplexReference): '3D simplex' # TETRAHEDRON: # c\d # a-b # # EDGES: # d\ d\ d\ c\ # b-c a-c a-b a-b # SUBDIVIDED TETRAHEDRON: # f\ i\j # d-e\g-h # a-b-c # # SUBDIVIDED EDGES: # j\ j\ j\ f\ # h-i\ g-i\ g-h\ d-e\ # c-e-f a-d-f a-b-c a-b-c # # CHILDREN: # d\g e\h f\i i\j e\g g\h g\i h\i # a-b b-c d-e g-h b-d b-e d-e e-g _children_vertices = [0,1,3,6], [1,2,4,7], [3,4,5,8], [6,7,8,9], [1,3,4,6], [1,4,6,7], [3,4,6,8], [4,6,7,8] def __init__(self): super().__init__(ndims=3) def getindices_vertex(self, n): m = 2**n+1 indis = numpy.arange(m) return numpy.array([[i,j,k] for k in indis for j in indis[:m-k] for i in indis[:m-j-k]]) def getischeme_vertex(self, n): return self.getindices_vertex(n)/(2**n), None def nvertices_by_level(self, n): m = 2**n+1 return ((m+2)*(m+1)*m)//6 def child_divide(self, vals, n): assert len(vals) == self.nvertices_by_level(n) child_indices = self.getindices_vertex(1) offset = numpy.array([1,0,0,0]) linear = numpy.array([[-1,-1,-1],[1,0,0],[0,1,0],[0,0,1]]) m = 2**n+1 cvals = [] for child_ref, child_vertices in zip(self.child_refs,self._children_vertices): V = child_indices[child_vertices] child_offset = (2**(n-1))*V.T.dot(offset) child_linear = V.T.dot(linear) original = child_ref.getindices_vertex(n-1) transformed = original.dot(child_linear.T) + child_offset i, j, k = transformed.T cvals.append(vals[((k-1)*k*(2*k-1)//6 - (1+2*m)*(k-1)*k//2 + m*(m+1)*k)//2 + ((2*(m-k)+1)*j-j**2)//2 + i]) return numpy.array(cvals)
[docs] def getischeme_gauss(self, degree): '''get integration scheme http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf''' if isinstance(degree, tuple): assert len(degree) == 3 degree = sum(degree) assert isinstance(degree, int) and degree >= 0 I = [0,0,0], J = [1,1,1],[0,1,1],[1,1,0],[1,0,1] K = [0,1,1],[1,0,1],[1,1,0],[1,0,0],[0,1,0],[0,0,1] L = [0,1,1],[1,0,1],[1,1,0],[2,1,1],[1,2,1],[1,1,2],[1,0,2],[0,2,1],[2,1,0],[1,2,0],[0,1,2],[2,0,1] icw = [ (I, [1/4], 1), ] if degree == 1 else [ (J, [0.5854101966249685,0.1381966011250105], 1/4), ] if degree == 2 else [ (I, [.25], -.8), (J, [.5,1/6], .45), ] if degree == 3 else [ (I, [.25], -.2368/3), (J, [0.7857142857142857,0.0714285714285714], .1372/3), (K, [0.1005964238332008,0.3994035761667992], .448/3), ] if degree == 4 else [ (I, [.25], 0.1817020685825351), (J, [0,1/3.], 0.0361607142857143), (J, [8/11.,1/11.], 0.0698714945161738), (K, [0.4334498464263357,0.0665501535736643], 0.0656948493683187), ] if degree == 5 else [ (J, [0.3561913862225449,0.2146028712591517], 0.0399227502581679), (J, [0.8779781243961660,0.0406739585346113], 0.0100772110553207), (J, [0.0329863295731731,0.3223378901422757], 0.0553571815436544), (L, [0.2696723314583159,0.0636610018750175,0.6030056647916491], 0.0482142857142857), ] if degree == 6 else [ (I, [.25], 0.1095853407966528), (J, [0.7653604230090441,0.0782131923303186], 0.0635996491464850), (J, [0.6344703500082868,0.1218432166639044], -0.3751064406859797), (J, [0.0023825066607383,0.3325391644464206], 0.0293485515784412), (K, [0,.5], 0.0058201058201058), (L, [.2,.1,.6], 0.1653439153439105) ] if degree == 7 else [ (I, [.25], -0.2359620398477557), (J, [0.6175871903000830,0.1274709365666390], 0.0244878963560562), (J, [0.9037635088221031,0.0320788303926323], 0.0039485206398261), (K, [0.4502229043567190,0.0497770956432810], 0.0263055529507371), (K, [0.3162695526014501,0.1837304473985499], 0.0829803830550589), (L, [0.0229177878448171,0.2319010893971509,0.5132800333608811], 0.0254426245481023), (L, [0.7303134278075384,0.0379700484718286,0.1937464752488044], 0.0134324384376852), ] if degree > 8: warnings.warn('inexact integration for polynomial of degree {}'.format(degree)) return numpy.concatenate([numpy.take(c,i) for i, c, w in icw], axis=0), \ numpy.concatenate([[w/6] * len(i) for i, c, w in icw])
[docs]class TensorReference(Reference): 'tensor reference' _re_ischeme = re.compile('([a-zA-Z]+)(.*)') def __init__(self, ref1, ref2): assert not isinstance(ref1, TensorReference) self.ref1 = ref1 self.ref2 = ref2 super().__init__(ref1.ndims + ref2.ndims) def __mul__(self, other): assert isinstance(other, Reference) return TensorReference(self.ref1, self.ref2 * other) @cache.property def vertices(self): vertices = numpy.empty((self.ref1.nverts, self.ref2.nverts, self.ndims), dtype=float) vertices[:,:,:self.ref1.ndims] = self.ref1.vertices[:,_] vertices[:,:,self.ref1.ndims:] = self.ref2.vertices[_,:] return numeric.const(vertices.reshape(self.ref1.nverts*self.ref2.nverts, self.ndims), copy=False) @property def centroid(self): return numpy.concatenate([self.ref1.centroid, self.ref2.centroid]) def nvertices_by_level(self, n): return self.ref1.nvertices_by_level(n) * self.ref2.nvertices_by_level(n) def child_divide(self, vals, n): np1 = self.ref1.nvertices_by_level(n) np2 = self.ref2.nvertices_by_level(n) return [v2.swapaxes(0,1).reshape((-1,)+vals.shape[1:]) for v1 in self.ref1.child_divide(vals.reshape((np1,np2)+vals.shape[1:]), n) for v2 in self.ref2.child_divide(v1.swapaxes(0,1), n)] def __str__(self): return '{}*{}'.format(self.ref1, self.ref2) def getischeme_vtk(self): if self.ref1.ndims == 0: assert self.ref1.nverts == 1 points, weights = self.ref2.getischeme_vtk() elif self.ref1.ndims == self.ref2.ndims == 1: points = numpy.empty([2, 2, 2]) points[...,:1] = self.ref1.vertices[:,_] points[0,:,1:] = self.ref2.vertices points[1,:,1:] = self.ref2.vertices[::-1] elif self.ref1.ndims <= 1 and self.ref2.ndims >= 1: points = numpy.empty([self.ref1.nverts, self.ref2.nverts, self.ndims]) points[...,:self.ref1.ndims] = self.ref1.vertices[:,_] points[...,self.ref1.ndims:] = self.ref2.vertices[_,:] elif self.ref1.ndims >= 1 and self.ref2.ndims <= 1: points = numpy.empty([self.ref2.nverts, self.ref1.nverts, self.ndims]) points[...,:self.ref1.ndims] = self.ref1.vertices[_,:] points[...,self.ref1.ndims:] = self.ref2.vertices[:,_] else: raise NotImplementedError return points.reshape(self.nverts, self.ndims), None def getischeme(self, ischeme): if '*' in ischeme: ischeme1, ischeme2 = ischeme.split('*', 1) else: match = self._re_ischeme.match(ischeme) assert match, 'cannot parse integration scheme {!r}'.format(ischeme) ptype, args = match.groups() get = getattr(self, 'getischeme_'+ptype, None) if get: return get(eval(args)) if args else get() if args and ',' in args: args = eval(args) assert len(args) == self.ndims ischeme1 = ptype+','.join(str(n) for n in args[:self.ref1.ndims]) ischeme2 = ptype+','.join(str(n) for n in args[self.ref1.ndims:]) else: ischeme1 = ischeme2 = ischeme ipoints1, iweights1 = self.ref1.getischeme(ischeme1) ipoints2, iweights2 = self.ref2.getischeme(ischeme2) ipoints = numpy.empty((len(ipoints1), len(ipoints2), self.ndims)) ipoints[:,:,0:self.ref1.ndims] = ipoints1[:,_,:self.ref1.ndims] ipoints[:,:,self.ref1.ndims:self.ndims] = ipoints2[_,:,:self.ref2.ndims] iweights = numeric.const((iweights1[:,_] * iweights2[_,:]).ravel(), copy=False) if iweights1 is not None and iweights2 is not None else None return numeric.const(ipoints.reshape(len(ipoints1) * len(ipoints2), self.ndims), copy=False), iweights @cache.property def edge_transforms(self): edge_transforms = [] if self.ref1.ndims: edge_transforms.extend(transform.TensorEdge1(trans1, self.ref2.ndims) for trans1 in self.ref1.edge_transforms) if self.ref2.ndims: edge_transforms.extend(transform.TensorEdge2(self.ref1.ndims, trans2) for trans2 in self.ref2.edge_transforms) return tuple(edge_transforms) @property def edge_refs(self): edge_refs = [] if self.ref1.ndims: edge_refs.extend(edge1 * self.ref2 for edge1 in self.ref1.edge_refs) if self.ref2.ndims: edge_refs.extend(self.ref1 * edge2 for edge2 in self.ref2.edge_refs) return tuple(edge_refs) @cache.property def ribbons(self): if self.ref1.ndims == 0: return self.ref2.ribbons if self.ref2.ndims == 0: return self.ref1.ribbons ribbons = [] for iedge1 in range(self.ref1.nedges): #iedge = self.ref1.edge_refs[iedge] * self.ref2 for iedge2 in range(self.ref2.nedges): #jedge = self.ref1 * self.ref2.edge_refs[jedge] jedge1 = self.ref1.nedges + iedge2 jedge2 = iedge1 if self.ref1.ndims > 1: iedge2 += self.ref1.edge_refs[iedge1].nedges ribbons.append(((iedge1,iedge2),(jedge1,jedge2))) if self.ref1.ndims >= 2: ribbons.extend(self.ref1.ribbons) if self.ref2.ndims >= 2: ribbons.extend(((iedge1+self.ref1.nedges,iedge2+self.ref1.nedges), (jedge1+self.ref1.nedges,jedge2+self.ref1.nedges)) for (iedge1,iedge2), (jedge1,jedge2) in self.ref2.ribbons) return tuple(ribbons) @cache.property def child_transforms(self): return tuple(transform.TensorChild(trans1, trans2) for trans1 in self.ref1.child_transforms for trans2 in self.ref2.child_transforms) @property def child_refs(self): return tuple(child1 * child2 for child1 in self.ref1.child_refs for child2 in self.ref2.child_refs) def inside(self, point, eps=0): return self.ref1.inside(point[:self.ref1.ndims],eps) and self.ref2.inside(point[self.ref1.ndims:],eps) @property def simplices(self): return tuple((transform.TensorChild(trans1, trans2), TensorReference(simplex1, simplex2)) for trans1, simplex1 in self.ref1.simplices for trans2, simplex2 in self.ref2.simplices) def get_ndofs(self, degree): return self.ref1.get_ndofs(degree)*self.ref2.get_ndofs(degree) def get_poly_coeffs(self, basis, **kwargs): return numeric.poly_outer_product(self.ref1.get_poly_coeffs(basis, **kwargs), self.ref2.get_poly_coeffs(basis, **kwargs)) def get_edge_dofs(self, degree, iedge): if not numeric.isint(iedge) or iedge < 0 or iedge >= self.nedges: raise IndexError('edge index out of range') nd2 = self.ref2.get_ndofs(degree) if iedge < self.ref1.nedges: dofs1 = self.ref1.get_edge_dofs(degree, iedge) dofs2 = range(self.ref2.get_ndofs(degree)) else: dofs1 = range(self.ref1.get_ndofs(degree)) dofs2 = self.ref2.get_edge_dofs(degree, iedge-self.ref1.nedges) return numeric.const(tuple(d1*nd2+d2 for d1, d2 in itertools.product(dofs1, dofs2)), dtype=int) @property def _flat_refs(self): for ref in self.ref1, self.ref2: if isinstance(ref, TensorReference): yield from ref._flat_refs else: yield ref def get_dof_transpose_map(self, degree, vertex_transpose_map): vertex_transpose_map = tuple(vertex_transpose_map) if len(vertex_transpose_map) != self.nverts: raise ValueError('invalid vertex indices: {!r}'.format(vertex_transpose_map)) refs = tuple(ref for ref in self._flat_refs if ref.nverts > 1) # Let `ref_verts[i]` be a permutation of `range(refs[i].nverts)`. The # `vertex_transpose_map` should be the tensor product of the # `ref_verts[i]*vertex_strides[i]` for all `i`, permuted by `perm` and # flattened. The `ref_strides` recovers the original structure from the # permuted and flattened `vertex_transpose_map`. We reverse engineer the # per ref vertices, `ref_verts`, and permutation of the references, `perm`, # and apply the same permutation and flattening to the tensor product of # the dofs. stride = 1 vertex_strides = [] for ref in refs[::-1]: vertex_strides.insert(0, stride) stride *= ref.nverts verts = numpy.array(0, dtype=int) dofs = numpy.array(0, dtype=int) ref_strides = [] for ref, stride in zip(refs, vertex_strides): ref_idx = [vertex_transpose_map.index(i*stride) for i in range(ref.nverts)] ref_verts = numpy.argsort(ref_idx) verts = verts[...,None]*len(ref_verts)+ref_verts ref_dofs = ref.get_dof_transpose_map(degree, ref_verts) dofs = dofs[...,None]*len(ref_dofs)+ref_dofs ref_strides.append(ref_idx[ref_verts[1]]-ref_idx[ref_verts[0]]) perm = numpy.argsort(ref_strides)[::-1] # Verify that `vertex_transpose_map` is in fact a tensor product of the # `ref_verts`. if not numpy.all(numpy.equal(numpy.transpose(verts, perm).ravel(), vertex_transpose_map)): raise ValueError('invalid transformation: {!r}'.format(vertex_transpose_map)) return numeric.const(numpy.transpose(dofs, perm).ravel())
[docs]class Cone(Reference): 'cone' def __init__(self, edgeref, etrans, tip:numeric.const): assert etrans.fromdims == edgeref.ndims assert etrans.todims == len(tip) super().__init__(len(tip)) self.edgeref = edgeref self.etrans = etrans self.tip = tip ext = etrans.ext self.extnorm = numpy.linalg.norm(ext) self.height = numpy.dot(etrans.offset - tip, ext) / self.extnorm assert self.height >= 0, 'tip is positioned at the negative side of edge' @cache.property def vertices(self): return numeric.const(numpy.vstack([[self.tip], self.etrans.apply(self.edgeref.vertices)]), copy=False) @cache.property def edge_transforms(self): edge_transforms = [self.etrans] if self.edgeref.ndims > 0: for trans, edge in self.edgeref.edges: if edge: b = self.etrans.apply(trans.offset) A = numpy.hstack([numpy.dot(self.etrans.linear, trans.linear), (self.tip-b)[:,_]]) newtrans = transform.Updim(A, b, isflipped=self.etrans.isflipped^trans.isflipped^(self.ndims%2==1)) # isflipped logic tested up to 3D edge_transforms.append(newtrans) else: edge_transforms.append(transform.Updim(numpy.zeros((1,0)), self.tip, isflipped=not self.etrans.isflipped)) return edge_transforms @cache.property def edge_refs(self): edge_refs = [self.edgeref] if self.edgeref.ndims > 0: extrudetrans = transform.Updim(numpy.eye(self.ndims-1)[:,:-1], numpy.zeros(self.ndims-1), isflipped=self.ndims%2==0) tip = numpy.array([0]*(self.ndims-2)+[1], dtype=float) edge_refs.extend(edge.cone(extrudetrans, tip) for edge in self.edgeref.edge_refs if edge) else: edge_refs.append(getsimplex(0)) return edge_refs def getischeme_gauss(self, degree): if self.nverts == self.ndims+1: # simplex spoints, sweights = getsimplex(self.ndims).getischeme_gauss(degree) offset = self.vertices[0,:] linear = self.vertices[1:,:] - offset points = numpy.dot(spoints, linear) + offset weights = sweights * abs(numpy.linalg.det(linear)) else: epoints, eweights = self.edgeref.getischeme('gauss{}'.format(degree)) tpoints, tweights = getsimplex(1).getischeme_gauss(degree + self.ndims - 1) tx, = tpoints.T points = (tx[:,_,_] * (self.etrans.apply(epoints)-self.tip)[_,:,:] + self.tip).reshape(-1, self.ndims) wx = tx**(self.ndims-1) * tweights * self.extnorm * self.height weights = (eweights[_,:] * wx[:,_]).ravel() return points, weights def getischeme_bezier(self, degree): assert self.nverts == self.ndims+1 spoints, none = getsimplex(self.ndims).getischeme_bezier(degree) offset = self.vertices[0,:] linear = self.vertices[1:,:] - offset return numpy.dot(spoints, linear) + offset, None def getischeme_vtk(self): if self.nverts == 4 and self.ndims==3: # tetrahedron I = slice(None) elif self.nverts == 5 and self.ndims==3: # pyramid I = numpy.array([1,2,4,3,0]) elif self.nverts == 3 and self.ndims==2: # triangle I = slice(None) else: raise Exception('invalid number of points: {}'.format(self.nverts)) return self.vertices[I], None @property def simplices(self): if self.nverts == self.ndims+1 or self.edgeref.ndims == 2 and self.edgeref.nverts == 4: # simplices and square-based pyramids are ok return [(transform.Identity(self.ndims), self)] return tuple((transform.Identity(self.ndims), ref.cone(self.etrans*trans,self.tip)) for trans, ref in self.edgeref.simplices) def inside(self, point, eps=0): # point = etrans.apply(epoint) * xi + tip * (1-xi) => etrans.apply(epoint) = tip + (point-tip) / xi xi = numpy.dot(self.etrans.ext, point-self.tip) / numpy.dot(self.etrans.ext, self.etrans.offset-self.tip) return 0 < xi <= 1+eps and self.edgeref.inside(numpy.linalg.solve( numpy.dot(self.etrans.linear.T, self.etrans.linear), numpy.dot(self.etrans.linear.T, self.tip + (point-self.tip)/xi - self.etrans.offset)), eps=eps)
[docs]class NeighborhoodTensorReference(TensorReference): 'product reference element' def __init__(self, ref1, ref2, neighborhood, transf): '''Neighborhood of elem1 and elem2 and transformations to get mutual overlap in right location. Returns 3-element tuple: * neighborhood, as given by Element.neighbor(), * transf1, required rotation of elem1 map: {0:0, 1:pi/2, 2:pi, 3:3*pi/2}, * transf2, required rotation of elem2 map (is indep of transf1 in Topology.''' super().__init__(ref1, ref2) self.neighborhood = neighborhood self.transf = transf def singular_ischeme_quad(self, points): transfpoints = numpy.empty(points.shape) def transform(points, transf): x, y = points[:,0], points[:,1] tx = x if transf in (0,1,6,7) else 1-x ty = y if transf in (0,3,4,7) else 1-y return function.stack((ty, tx) if transf%2 else (tx, ty), axis=1) transfpoints[:,:2] = transform(points[:,:2], self.transf[0]) transfpoints[:,2:] = transform(points[:,2:], self.transf[1]) return transfpoints
[docs] def get_tri_bem_ischeme(self, ischeme): 'Some cached quantities for the singularity quadrature scheme.' points, weights = (getsimplex(1)**4).getischeme(ischeme) eta1, eta2, eta3, xi = points.T if self.neighborhood == 0: temp = xi*eta1*eta2*eta3 pts0 = xi*eta1*(1 - eta2) pts1 = xi - pts0 pts2 = xi - temp pts3 = xi*(1 - eta1) pts4 = pts0 + temp pts5 = xi*(1 - eta1*eta2) pts6 = xi*eta1 - temp points = numpy.array( [[1-xi, 1-pts2, 1-xi, 1-pts5, 1-pts2, 1-xi ], [pts1, pts3, pts4, pts0, pts6, pts0], [1-pts2, 1-xi, 1-pts5, 1-xi, 1-xi, 1-pts2], [pts3, pts1, pts0, pts4, pts0, pts6]]).reshape(4, -1).T points = points * [-1,1,-1,1] + [1,0,1,0] # flipping in x -GJ weights = numpy.concatenate(6*[xi**3*eta1**2*eta2*weights]) elif self.neighborhood == 1: A = xi*eta1 B = A*eta2 C = A*eta3 D = B*eta3 E = xi - B F = A - B G = xi - D H = B - D I = A - D points = numpy.array( [[1-xi, 1-xi, 1-E, 1-G, 1-G ], [C, G, F, H, I ], [1-E, 1-G, 1-xi, 1-xi, 1-xi], [F, H, D, A, B ]]).reshape(4, -1).T temp = xi*A weights = numpy.concatenate([A*temp*weights] + 4*[B*temp*weights]) elif self.neighborhood == 2: A = xi*eta2 B = A*eta3 C = xi*eta1 points = numpy.array( [[1-xi, 1-A ], [C, B ], [1-A, 1-xi], [B, C ]]).reshape(4, -1).T weights = numpy.concatenate(2*[xi**2*A*weights]) else: assert self.neighborhood == -1, 'invalid neighborhood {!r}'.format(self.neighborhood) points = numpy.array([eta1*eta2, 1-eta2, eta3*xi, 1-xi]).T weights = eta2*xi*weights return points, weights
[docs] def get_quad_bem_ischeme(self, ischeme): 'Some cached quantities for the singularity quadrature scheme.' quad = getsimplex(1)**4 points, weights = quad.getischeme(ischeme) eta1, eta2, eta3, xi = points.T if self.neighborhood == 0: xe = xi*eta1 A = (1 - xi)*eta3 B = (1 - xe)*eta2 C = xi + A D = xe + B points = numpy.array( [[A, B, A, D, B, C, C, D], [B, A, D, A, C, B, D, C], [C, D, C, B, D, A, A, B], [D, C, B, C, A, D, B, A]]).reshape(4, -1).T weights = numpy.concatenate(8*[xi*(1-xi)*(1-xe)*weights]) elif self.neighborhood == 1: ox = 1 - xi A = xi*eta1 B = xi*eta2 C = ox*eta3 D = C + xi E = 1 - A F = E*eta3 G = A + F points = numpy.array( [[D, C, G, G, F, F ], [B, B, B, xi, B, xi], [C, D, F, F, G, G ], [A, A, xi, B, xi, B ]]).reshape(4, -1).T weights = numpy.concatenate(2*[xi**2*ox*weights] + 4*[xi**2*E*weights]) elif self.neighborhood == 2: A = xi*eta1 B = xi*eta2 C = xi*eta3 points = numpy.array( [[xi, A, A, A ], [A, xi, B, B ], [B, B, xi, C ], [C, C, C, xi]]).reshape(4, -1).T weights = numpy.concatenate(4*[xi**3*weights]) else: assert self.neighborhood == -1, 'invalid neighborhood {!r}'.format(self.neighborhood) return points, weights
[docs] def getischeme_singular(self, n): 'get integration scheme' gauss = 'gauss{}'.format(n*2-2) assert self.ref1 == self.ref2 == getsimplex(1)**2 points, weights = self.get_quad_bem_ischeme(gauss) return self.singular_ischeme_quad(points), weights
[docs]class OwnChildReference(Reference): 'forward self as child' def __init__(self, baseref): self.baseref = baseref self.child_refs = baseref, self.child_transforms = transform.Identity(baseref.ndims), super().__init__(baseref.ndims) @property def vertices(self): return self.baseref.vertices @property def edge_transforms(self): return self.baseref.edge_transforms @property def edge_refs(self): return [OwnChildReference(edge) for edge in self.baseref.edge_refs] def getischeme(self, ischeme): return self.baseref.getischeme(ischeme) @property def simplices(self): return self.baseref.simplices def get_ndofs(self, degree): return self.baseref.get_ndofs(degree) def get_poly_coeffs(self, basis, **kwargs): return self.baseref.get_poly_coeffs(basis, **kwargs) def get_edge_dofs(self, degree, iedge): return self.baseref.get_edge_dofs(degree, iedge) def get_dof_transpose_map(self, degree, vertex_transpose_map): return self.baseref.get_dof_transpose_map(degree, vertex_transpose_map)
[docs]class WithChildrenReference(Reference): 'base reference with explicit children' def __init__(self, baseref, child_refs:tuple): assert len(child_refs) == baseref.nchildren and any(child_refs) and child_refs != baseref.child_refs assert all(isinstance(child_ref,Reference) for child_ref in child_refs) assert all(child_ref.ndims == baseref.ndims for child_ref in child_refs) self.baseref = baseref self.child_transforms = baseref.child_transforms self.child_refs = child_refs super().__init__(baseref.ndims) def check_edges(self, tol=1e-15, print=print): super().check_edges(tol=tol, print=print) for cref in self.child_refs: cref.check_edges(tol=tol, print=print) @property def vertices(self): return self.baseref.vertices @property def permutation_transforms(self): return self.baseref.permutation_transforms def nvertices_by_level(self, n): return self.baseref.nvertices_by_level(n) def child_divide(self, vals, n): return self.baseref.child_divide(vals, n) __sub__ = lambda self, other: self.empty if other in (self,self.baseref) else self.baseref.with_children(self_child-other_child for self_child, other_child in zip(self.child_refs, other.child_refs)) if isinstance(other, WithChildrenReference) and other.baseref in (self,self.baseref) else NotImplemented __rsub__ = lambda self, other: self.baseref.with_children(other_child - self_child for self_child, other_child in zip(self.child_refs, other.child_refs)) if other == self.baseref else NotImplemented __and__ = lambda self, other: self if other == self.baseref else other if isinstance(other,WithChildrenReference) and self == other.baseref else self.baseref.with_children(self_child & other_child for self_child, other_child in zip(self.child_refs, other.child_refs)) if isinstance(other, WithChildrenReference) and other.baseref == self.baseref else NotImplemented __or__ = lambda self, other: other if other == self.baseref else self.baseref.with_children(self_child | other_child for self_child, other_child in zip(self.child_refs, other.child_refs)) if isinstance(other, WithChildrenReference) and other.baseref == self.baseref else NotImplemented @cache.property def __extra_edges(self): extra_edges = [(ichild, iedge, cref.edge_refs[iedge]) for ichild, cref in enumerate(self.child_refs) if cref for iedge in range(self.baseref.child_refs[ichild].nedges, cref.nedges)] for ichild, edges in enumerate(self.baseref.connectivity): cref = self.child_refs[ichild] if not cref: continue # new child is empty for iedge, jchild in enumerate(edges): if jchild == -1: continue # iedge already is an external boundary coppref = self.child_refs[jchild] if coppref == self.baseref.child_refs[jchild]: continue # opposite is complete, so iedge cannot form a new external boundary eref = cref.edge_refs[iedge] if coppref: # opposite new child is not empty eref -= coppref.edge_refs[self.baseref.connectivity[jchild].index(ichild)] if eref: extra_edges.append((ichild, iedge, eref)) return extra_edges def subvertex(self, ichild, i): assert 0<=ichild<self.nchildren npoints = 0 for childindex, child in enumerate(self.child_refs): if child: points, weights = child.getischeme('vertex{}'.format(i-1)) assert weights is None if childindex == ichild: rng = numpy.arange(npoints, npoints+len(points)) npoints += len(points) elif ichild==childindex: rng = numpy.array([],dtype=int) return npoints, rng
[docs] def getischeme(self, ischeme): 'get integration scheme' if ischeme.startswith('vertex'): return self.baseref.getischeme(ischeme) allcoords = [] allweights = [] for trans, ref in self.children: if ref: points, weights = ref.getischeme(ischeme) allcoords.append(trans.apply(points)) if weights is not None: allweights.append(weights * abs(float(trans.det))) coords = numpy.concatenate(allcoords, axis=0) weights = numpy.concatenate(allweights, axis=0) \ if len(allweights) == len(allcoords) else None return coords, weights
@property def simplices(self): return [(trans2*trans1, simplex) for trans2, child in self.children for trans1, simplex in (child.simplices if child else [])] @cache.property def edge_transforms(self): return tuple(self.baseref.edge_transforms) \ + tuple(transform.ScaledUpdim(self.child_transforms[ichild], self.child_refs[ichild].edge_transforms[iedge]) for ichild, iedge, ref in self.__extra_edges) @cache.property def edge_refs(self): refs = [] for etrans, eref in self.baseref.edges: children = [] if eref: for ctrans, cref in eref.children: ctrans_, etrans_ = etrans.swapup(ctrans) ichild = self.baseref.child_transforms.index(ctrans_) cref = self.child_refs[ichild] children.append(cref.edge_refs[cref.edge_transforms.index(etrans_)] if cref else EmptyReference(self.ndims-1)) refs.append(eref.with_children(children)) for ichild, iedge, ref in self.__extra_edges: refs.append(OwnChildReference(ref)) return tuple(refs) @cache.property def connectivity(self): # same as base implementation but cheaper return tuple(tuple(edges[iedge] if iedge < len(edges) and edges[iedge] != -1 and self.child_refs[edges[iedge]] else -1 for iedge in range(self.child_refs[ichild].nedges)) for ichild, edges in enumerate(self.baseref.connectivity)) def inside(self, point, eps=0): return any(cref.inside(ctrans.invapply(point), eps=eps) for ctrans, cref in self.children) def get_ndofs(self, degree): return self.baseref.get_ndofs(degree) def get_poly_coeffs(self, basis, **kwargs): return self.baseref.get_poly_coeffs(basis, **kwargs) def get_edge_dofs(self, degree, iedge): return self.baseref.get_edge_dofs(degree, iedge) def get_dof_transpose_map(self, degree, vertex_transpose_map): return self.baseref.get_dof_transpose_map(degree, vertex_transpose_map)
[docs]class MosaicReference(Reference): 'triangulation' def __init__(self, baseref, edge_refs:tuple, midpoint:numeric.const): assert len(edge_refs) == baseref.nedges assert edge_refs != tuple(baseref.edge_refs) self.baseref = baseref self._edge_refs = edge_refs self._midpoint = midpoint self.edge_refs = list(edge_refs) self.edge_transforms = list(baseref.edge_transforms) if baseref.ndims == 1: assert any(edge_refs) and not all(edge_refs), 'invalid 1D mosaic: exactly one edge should be non-empty' iedge, = [i for i, edge in enumerate(edge_refs) if edge] self.edge_refs.append(getsimplex(0)) self.edge_transforms.append(transform.Updim(linear=numpy.zeros((1,0)), offset=midpoint, isflipped=not baseref.edge_transforms[iedge].isflipped)) else: newedges = [(etrans1, etrans2, edge) for (etrans1,orig), new in zip(baseref.edges, edge_refs) for etrans2, edge in new.edges[orig.nedges:]] for (iedge1,iedge2), (jedge1,jedge2) in baseref.ribbons: Ei = edge_refs[iedge1] ei = Ei.edge_refs[iedge2] if Ei else EmptyReference(Ei.ndims-1) Ej = edge_refs[jedge1] ej = Ej.edge_refs[jedge2] if Ej else EmptyReference(Ej.ndims-1) ejsubi = ej - ei if ejsubi: newedges.append((self.edge_transforms[jedge1], Ej.edge_transforms[jedge2], ejsubi)) eisubj = ei - ej if eisubj: newedges.append((self.edge_transforms[iedge1], Ei.edge_transforms[iedge2], eisubj)) extrudetrans = transform.Updim(numpy.eye(baseref.ndims-1)[:,:-1], numpy.zeros(baseref.ndims-1), isflipped=baseref.ndims%2==0) tip = numpy.array([0]*(baseref.ndims-2)+[1], dtype=float) for etrans, trans, edge in newedges: b = etrans.apply(trans.offset) A = numpy.hstack([numpy.dot(etrans.linear, trans.linear), (midpoint-b)[:,_]]) newtrans = transform.Updim(A, b, isflipped=etrans.isflipped^trans.isflipped^(baseref.ndims%2==1)) # isflipped logic tested up to 3D self.edge_transforms.append(newtrans) self.edge_refs.append(edge.cone(extrudetrans, tip)) super().__init__(baseref.ndims) @cache.property def vertices(self): vertices = [] for etrans, eref in self.edges: indices = [] for vertex in etrans.apply(eref.vertices).tolist(): try: index = vertices.index(vertex) except ValueError: index = len(vertices) vertices.append(vertex) indices.append(index) return numeric.const(vertices) def __and__(self, other): if other in (self,self.baseref): return self if isinstance(other, MosaicReference) and other.baseref == self: return other if isinstance(other, MosaicReference) and self.baseref == other.baseref and numpy.equal(other._midpoint, self._midpoint).all(): isect_edge_refs = [selfedge & otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs)] if not any(isect_edge_refs): return self.empty return MosaicReference(self.baseref, isect_edge_refs, self._midpoint) return NotImplemented def __or__(self, other): if other in (self,self.baseref): return other if isinstance(other, MosaicReference) and self.baseref == other.baseref and numpy.equal(other._midpoint, self._midpoint).all(): union_edge_refs = [selfedge | otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs)] if tuple(union_edge_refs) == tuple(self.baseref.edge_refs): return self.baseref return MosaicReference(self.baseref, union_edge_refs, self._midpoint) return NotImplemented def __sub__(self, other): if other in (self,self.baseref): return self.empty if isinstance(other, MosaicReference) and other.baseref == self: inv_edge_refs = [baseedge - edge for baseedge, edge in zip(self.edge_refs, other._edge_refs)] return MosaicReference(self, inv_edge_refs, other._midpoint) return NotImplemented def __rsub__(self, other): if other == self.baseref: inv_edge_refs = [baseedge - edge for baseedge, edge in zip(other.edge_refs, self._edge_refs)] return MosaicReference(other, inv_edge_refs, self._midpoint) return NotImplemented def nvertices_by_level(self, n): return self.baseref.nvertices_by_level(n) @cache.property def subrefs(self): return [ref.cone(trans,self._midpoint) for trans, ref in zip(self.baseref.edge_transforms, self._edge_refs) if ref] @property def simplices(self): return [simplex for subvol in self.subrefs for simplex in subvol.simplices]
[docs] def getischeme(self, ischeme): 'get integration scheme' if ischeme.startswith('vertex'): return self.baseref.getischeme(ischeme) allpoints, allweights = zip(*[subvol.getischeme(ischeme) for subvol in self.subrefs]) points = numpy.concatenate(allpoints, axis=0) weights = None if any(w is None for w in allweights) else numpy.concatenate(allweights, axis=0) return points, weights
def inside(self, point, eps=0): return any(subref.inside(point, eps=eps) for subref in self.subrefs) def get_ndofs(self, degree): return self.baseref.get_ndofs(degree) def get_poly_coeffs(self, basis, **kwargs): return self.baseref.get_poly_coeffs(basis, **kwargs) def get_edge_dofs(self, degree, iedge): return self.baseref.get_edge_dofs(degree, iedge) def get_dof_transpose_map(self, degree, vertex_transpose_map): return self.baseref.get_dof_transpose_map(degree, vertex_transpose_map)
# UTILITY FUNCTIONS _gauss = [] def gauss(degree): n = degree // 2 while len(_gauss) <= n: _gauss.append(None) gaussn = _gauss[n] if gaussn is None: k = numpy.arange(n) + 1 d = k / numpy.sqrt(4*k**2-1) x, w = numpy.linalg.eigh(numpy.diagflat(d,-1)) # eigh operates (by default) on lower triangle _gauss[n] = gaussn = numeric.const((x+1) * .5, copy=False), numeric.const(w[0]**2, copy=False) return gaussn def getsimplex(ndims): Simplex_by_dim = PointReference, LineReference, TriangleReference, TetrahedronReference return Simplex_by_dim[ndims]() def index_or_append(items, item): try: index = items.index(item) except ValueError: index = len(items) items.append(item) return index def arglexsort(triangulation): return numpy.argsort(numeric.asobjvector(tuple(tri) for tri in triangulation)) # vim:shiftwidth=2:softtabstop=2:expandtab:foldmethod=indent:foldnestmax=2