Source code for nutils.topology

# 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 topology module defines the topology objects, notably the
:class:`StructuredTopology` and :class:`UnstructuredTopology`. Maintaining
strict separation of topological and geometrical information, the topology
represents a set of elements and their interconnectivity, boundaries,
refinements, subtopologies etc, but not their positioning in physical space. The
dimension of the topology represents the dimension of its elements, not that of
the the space they are embedded in.

The primary role of topologies is to form a domain for :mod:`nutils.function`
objects, like the geometry function and function bases for analysis, as well as
provide tools for their construction. It also offers methods for integration and
sampling, thus providing a high level interface to operations otherwise written
out in element loops. For lower level operations topologies can be used as
:mod:`nutils.element` iterators.

from . import element, function, util, numpy, parallel, matrix, log, config, numeric, cache, transform, warnings, _
import functools,, itertools, functools, operator

_identity = lambda x: x

[docs]class Topology: 'topology base class' # subclass needs to implement: .elements def __init__(self, ndims): 'constructor' assert numeric.isint(ndims) and ndims >= 0 self.ndims = ndims def __str__(self): 'string representation' return '{}(#{})'.format(self.__class__.__name__, len(self)) def __len__(self): return len(self.elements) def __iter__(self): return iter(self.elements) def getitem(self, item): return EmptyTopology(self.ndims) def __getitem__(self, item): if not isinstance(item, tuple): item = item, if all(it in (...,slice(None)) for it in item): return self topo = self.getitem(item) if len(item) != 1 or not isinstance(item[0],str) \ else functools.reduce(operator.or_, map(self.getitem, item[0].split(',')), EmptyTopology(self.ndims)) if not topo: raise KeyError(item) return topo def __invert__(self): return OppositeTopology(self) def __or__(self, other): assert isinstance(other, Topology) and other.ndims == self.ndims return other if not self \ else self if not other \ else NotImplemented if isinstance(other, UnionTopology) \ else UnionTopology((self,other)) __ror__ = lambda self, other: self.__or__(other) def __and__(self, other): # Strategy: loop over combined elements sorted by .transform while keeping # track of the origin (mine=True for self, mine=False for other), and # select an element if it is equal to or a refinement of the previous # (hold) element and it originates from the other topology (mine == need). # Hold is not updated in case of a match because it might match multiple # children. elems = [] need = None for elem, mine in sorted([(elem, True) for elem in self] + [(elem, False) for elem in other], key=lambda v: v[0].transform): if mine == need and elem.transform[:len(hold.transform)] == hold.transform: assert elem.opposite[:len(hold.opposite)] == hold.opposite elems.append(elem) else: hold = elem need = not mine return UnstructuredTopology(self.ndims, elems) __rand__ = lambda self, other: self.__and__(other) def __add__(self, other): return self | other def __contains__(self, element): ielem = self.edict.get(element.transform) return ielem is not None and self.elements[ielem] == element def __sub__(self, other): assert isinstance(other, Topology) and other.ndims == self.ndims return other.__rsub__(self) def __rsub__(self, other): assert isinstance(other, Topology) and other.ndims == self.ndims return other - other.subset(self, newboundary=getattr(self,'boundary',None)) def __mul__(self, other): return ProductTopology(self, other) def edict(self): '''transform -> ielement mapping''' return {elem.transform: ielem for ielem, elem in enumerate(self)} def border_transforms(self): border_transforms = set() for belem in self.boundary: try: ielem, tail = transform.lookup_item(belem.transform, self.edict) except KeyError: pass else: border_transforms.add(self.elements[ielem].transform) return border_transforms @property def refine_iter(self): topo = self for irefine in log.count('refinement level'): yield topo topo = topo.refined stdfunc = lambda self, *args, **kwargs: self.basis('std', *args, **kwargs) linearfunc = lambda self, *args, **kwargs: self.basis('std', degree=1, *args, **kwargs) splinefunc = lambda self, *args, **kwargs: self.basis('spline', *args, **kwargs) bubblefunc = lambda self, *args, **kwargs: self.basis('bubble', *args, **kwargs) discontfunc = lambda self, *args, **kwargs: self.basis('discont', *args, **kwargs) def basis(self, name, *args, **kwargs): if self.ndims == 0: return function.asarray([1]) f = getattr(self, 'basis_' + name) return f(*args, **kwargs)
[docs] @log.title @util.single_or_multiple def elem_eval(self, funcs, ischeme, separate=False, geometry=None, asfunction=False, edit=_identity, *, arguments=None): 'element-wise evaluation' fcache = cache.WrapperCache() assert not separate or not asfunction, '"separate" and "asfunction" are mutually exclusive' if geometry: iwscale = function.J(geometry, self.ndims) npoints = len(self) slices = range(npoints) else: iwscale = 1 slices = [] npoints = 0 for elem in log.iter('elem', self): ipoints, iweights = ischeme[elem] if isinstance(ischeme, else fcache[elem.reference.getischeme](ischeme) np = len(ipoints) slices.append(slice(npoints,npoints+np)) npoints += np nprocs = min(config.nprocs, len(self)) zeros = parallel.shzeros if nprocs > 1 else numpy.zeros retvals = [] idata = [] for ifunc, func in enumerate(funcs): func = function.asarray(edit(func * iwscale)) func = function.zero_argument_derivatives(func) retval = zeros((npoints,)+func.shape, dtype=func.dtype) idata.extend(function.Tuple([ifunc, function.Tuple(ind), f.simplified]) for ind, f in function.blocks(func)) retvals.append(retval) idata = function.Tuple(idata) if idata.graphviz() if arguments is None: arguments = {} for ielem, elem in parallel.pariter(log.enumerate('elem', self), nprocs=nprocs): ipoints, iweights = ischeme[elem] if isinstance(ischeme, else fcache[elem.reference.getischeme](ischeme) s = slices[ielem], try: for ifunc, index, data in idata.eval(_transforms=(elem.transform, elem.opposite), _points=ipoints, _cache=fcache, **arguments):[ifunc], s+numpy.ix_(*[ind for (ind,) in index]),,data) if geometry else data) except function.EvaluationError: warnings.warn('not all functions evaluated successfully') for ifunc, indexfunc, datafunc in idata: index = indexfunc.eval(_transforms=(elem.transform, elem.opposite), _points=ipoints, _cache=fcache, **arguments) try: data = datafunc.eval(_transforms=(elem.transform, elem.opposite), _points=ipoints, _cache=fcache, **arguments) except function.EvaluationError: retvals[ifunc][s+numpy.ix_(*[ind for (ind,) in index])] = numpy.nan else:[ifunc], s+numpy.ix_(*[ind for (ind,) in index]),,data) if geometry else data) log.debug('cache', fcache.stats)'created', ', '.join('{}({})'.format(retval.__class__.__name__, ','.join(str(n) for n in retval.shape)) for retval in retvals)) if asfunction: if geometry: retvals = [function.elemwise({elem.transform: value for elem, value in zip(self, retval)}, shape=retval.shape[1:]) for retval in retvals] else: tsp = [(elem.transform, s, fcache[elem.reference.getischeme](ischeme)[0]) for elem, s in zip(self, slices)] retvals = [function.sampled({trans: (numeric.const(retval[s], copy=False), points) for trans, s, points in tsp}, self.ndims) for retval in retvals] elif separate: retvals = [[retval[s] for s in slices] for retval in retvals] return retvals
[docs] @log.title @util.single_or_multiple def elem_mean(self, funcs, geometry, ischeme, *, arguments=None): 'element-wise average' retvals = self.elem_eval((1,)+funcs, geometry=geometry, ischeme=ischeme, arguments=arguments) return [v / retvals[0][(slice(None),)+(_,)*(v.ndim-1)] for v in retvals[1:]]
def _integrate(self, funcs, ischeme, fcache=None, arguments=None): if arguments is None: arguments = {} # Functions may consist of several blocks, such as originating from # chaining. Here we make a list of all blocks consisting of triplets of # argument id, evaluable index, and evaluable values. blocks = [(ifunc, function.Tuple(ind), f.simplified) for ifunc, func in enumerate(funcs) for ind, f in function.blocks(function.zero_argument_derivatives(func))] block2func, indices, values = zip(*blocks) if blocks else ([],[],[]) log.debug('integrating {} distinct blocks'.format('+'.join( str(block2func.count(ifunc)) for ifunc in range(len(funcs))))) if function.Tuple(values).graphviz() if fcache is None: fcache = cache.WrapperCache() # To allocate (shared) memory for all block data we evaluate indexfunc to # build an nblocks x nelems+1 offset array, and nblocks index lists of # length nelems. offsets = numpy.zeros((len(blocks), len(self)+1), dtype=int) if blocks: sizefunc = function.stack([f.size for ifunc, ind, f in blocks]).simplified for ielem, elem in enumerate(self): n, = sizefunc.eval(_transforms=(elem.transform, elem.opposite), _cache=fcache, **arguments) offsets[:,ielem+1] = offsets[:,ielem] + n # Since several blocks may belong to the same function, we post process the # offsets to form consecutive intervals in longer arrays. The length of # these arrays is captured in the nfuncs-array nvals. nvals = numpy.zeros(len(funcs), dtype=int) for iblock, ifunc in enumerate(block2func): offsets[iblock] += nvals[ifunc] nvals[ifunc] = offsets[iblock,-1] # The data_index list contains shared memory index and value arrays for # each function argument. nprocs = min(config.nprocs, len(self)) empty = parallel.shempty if nprocs > 1 else numpy.empty data_index = [ (empty(n, dtype=float), empty((funcs[ifunc].ndim,n), dtype=int)) for ifunc, n in enumerate(nvals) ] # In a second, parallel element loop, valuefunc is evaluated to fill the # data part of data_index using the offsets array for location. Each # element has its own location so no locks are required. The index part of # data_index is filled in the same loop. It does not use valuefunc data but # benefits from parallel speedup. valueindexfunc = function.Tuple(function.Tuple([value]+list(index)) for value, index in zip(values, indices)) for ielem, elem in parallel.pariter(log.enumerate('elem', self), nprocs=nprocs): ipoints, iweights = ischeme[elem] if isinstance(ischeme, else fcache[elem.reference.getischeme](ischeme) assert iweights is not None, 'no integration weights found' for iblock, (intdata, *indices) in enumerate(valueindexfunc.eval(_transforms=(elem.transform, elem.opposite), _points=ipoints, _cache=fcache, **arguments)): s = slice(*offsets[iblock,ielem:ielem+2]) data, index = data_index[block2func[iblock]] w_intdata =, intdata) data[s] = w_intdata.ravel() si = (slice(None),) + (_,) * (w_intdata.ndim-1) for idim, (ii,) in enumerate(indices): index[idim,s].reshape(w_intdata.shape)[...] = ii[si] si = si[:-1] log.debug('cache', fcache.stats) return data_index
[docs] @log.title @util.single_or_multiple def integrate(self, funcs, ischeme='gauss', degree=None, geometry=None, force_dense=False, fcache=None, edit=_identity, *, arguments=None): 'integrate' if degree is not None: ischeme += str(degree) iwscale = function.J(geometry, self.ndims) if geometry else 1 integrands = [function.asarray(edit(func * iwscale)) for func in funcs] data_index = self._integrate(integrands, ischeme, fcache, arguments) return [matrix.assemble(data, index, integrand.shape, force_dense) for integrand, (data,index) in zip(integrands, data_index)]
[docs] @log.title def integral(self, func, ischeme='gauss', degree=None, geometry=None, edit=_identity): 'integral' if degree is not None: ischeme += str(degree) iwscale = function.J(geometry, self.ndims) if geometry else 1 integrand = edit(func * iwscale) from . import solver return solver.Integral([((self, ischeme), integrand)])
[docs] def projection(self, fun, onto, geometry, **kwargs): 'project and return as function' weights = self.project(fun, onto, geometry, **kwargs) return
[docs] @log.title def project(self, fun, onto, geometry, tol=0, ischeme='gauss', degree=None, droptol=1e-12, exact_boundaries=False, constrain=None, verify=None, ptype='lsqr', precon='diag', edit=_identity, *, arguments=None, **solverargs): 'L2 projection of function onto function space' log.debug('projection type:', ptype) if degree is not None: ischeme += str(degree) if constrain is None: constrain = util.NanVec(onto.shape[0]) else: constrain = constrain.copy() if exact_boundaries: constrain |= self.boundary.project(fun, onto, geometry, constrain=constrain, title='boundaries', ischeme=ischeme, tol=tol, droptol=droptol, ptype=ptype, edit=edit, arguments=arguments) assert isinstance(constrain, util.NanVec) assert constrain.shape == onto.shape[:1] avg_error = None # setting this depends on projection type if ptype == 'lsqr': assert ischeme is not None, 'please specify an integration scheme for lsqr-projection' fun2 = function.asarray(fun)**2 if len(onto.shape) == 1: Afun = function.outer(onto) bfun = onto * fun elif len(onto.shape) == 2: Afun = function.outer(onto).sum(2) bfun = function.sum(onto * fun, -1) if fun2.ndim: fun2 = fun2.sum(-1) else: raise Exception assert fun2.ndim == 0 A, b, f2, area = self.integrate([Afun,bfun,fun2,1], geometry=geometry, ischeme=ischeme, edit=edit, arguments=arguments, title='building system') N = A.rowsupp(droptol) if numpy.equal(b, 0).all(): constrain[~constrain.where&N] = 0 avg_error = 0. else: solvecons = constrain.copy() solvecons[~(constrain.where|N)] = 0 u = A.solve(b, solvecons, tol=tol, symmetric=True, precon=precon, **solverargs) constrain[N] = u[N] err2 = f2 -*b-A.matvec(u), u) # can be negative ~zero due to rounding errors avg_error = numpy.sqrt(err2) / area if err2 > 0 else 0 elif ptype == 'convolute': assert ischeme is not None, 'please specify an integration scheme for convolute-projection' if len(onto.shape) == 1: ufun = onto * fun afun = onto elif len(onto.shape) == 2: ufun = function.sum(onto * fun, axis=-1) afun = function.norm2(onto) else: raise Exception u, scale = self.integrate([ufun, afun], geometry=geometry, ischeme=ischeme, edit=edit, arguments=arguments) N = ~constrain.where & (scale > droptol) constrain[N] = u[N] / scale[N] elif ptype == 'nodal': ## data = function.Tuple([fun, onto]) ## F = W = 0 ## for elem in self: ## f, w = data(elem, 'bezier2') ## W += w.sum(axis=-1).sum(axis=0) ## F += numeric.contract(f[:,_,:], w, axis=[0,2]) ## I = (W!=0) F = numpy.zeros(onto.shape[0]) W = numpy.zeros(onto.shape[0]) I = numpy.zeros(onto.shape[0], dtype=bool) fun = function.zero_argument_derivatives(function.asarray(fun)) data = function.Tuple(function.Tuple([fun, onto_f.simplified, function.Tuple(onto_ind)]) for onto_ind, onto_f in function.blocks(function.zero_argument_derivatives(onto))) for elem in self: ipoints, iweights = elem.getischeme('bezier2') for fun_, onto_f_, onto_ind_ in data.eval(_transforms=(elem.transform, elem.opposite), _points=ipoints, **arguments or {}): onto_f_ = onto_f_.swapaxes(0,1) # -> dof axis, point axis, ... indfun_ = fun_[(slice(None),)+numpy.ix_(*onto_ind_[1:])] assert onto_f_.shape[0] == len(onto_ind_[0]) assert onto_f_.shape[1:] == indfun_.shape W[onto_ind_[0]] += onto_f_.reshape(onto_f_.shape[0],-1).sum(1) F[onto_ind_[0]] += (onto_f_ * indfun_).reshape(onto_f_.shape[0],-1).sum(1) I[onto_ind_[0]] = True I[constrain.where] = False constrain[I] = F[I] / W[I] else: raise Exception('invalid projection {!r}'.format(ptype)) numcons = constrain.where.sum() info = 'constrained {}/{} dofs'.format(numcons, constrain.size) if avg_error is not None: info += ', error {:.2e}/area'.format(avg_error) if verify is not None: assert numcons == verify, 'number of constraints does not meet expectation: {} != {}'.format(numcons, verify) return constrain def simplex(self): simplices = [simplex for elem in self for simplex in elem.simplices] return UnstructuredTopology(self.ndims, simplices)
[docs] def refined_by(self, refine): 'create refined space by refining dofs in existing one' refine = set(item.transform if isinstance(item,element.Element) else item for item in refine) refined = [] for elem in self: if elem.transform in refine: refined.extend(elem.children) else: refined.append(elem) return self.hierarchical(refined, precise=True)
def hierarchical(self, refined, precise=False): return HierarchicalTopology(self, refined, precise) @property def refined(self): return RefinedTopology(self)
[docs] def refine(self, n): 'refine entire topology n times' if numpy.iterable(n): assert len(n) == self.ndims assert all(ni == n[0] for ni in n) n = n[0] return self if n <= 0 else self.refined.refine(n-1)
[docs] @log.title def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None, *, arguments=None): 'trim element along levelset' if arguments is None: arguments = {} fcache = cache.WrapperCache() levelset = function.zero_argument_derivatives(levelset).simplified if leveltopo is None: ischeme = 'vertex{}'.format(maxrefine) refs = [elem.reference.trim(levelset.eval(_transforms=(elem.transform, elem.opposite), _points=fcache[elem.reference.getischeme](ischeme)[0], _cache=fcache, **arguments), maxrefine=maxrefine, ndivisions=ndivisions) for elem in log.iter('elem', self)] else:'collecting leveltopo elements') bins = [[] for ielem in range(len(self))] for elem in leveltopo: ielem, tail = transform.lookup_item(elem.transform, self.edict) bins[ielem].append(tail) refs = [] for elem, ctransforms in'elem', self, bins): levels = numpy.empty(elem.reference.nvertices_by_level(maxrefine)) cover = list(fcache[elem.reference.vertex_cover](tuple(sorted(ctransforms)), maxrefine)) # confirm cover and greedily optimize order mask = numpy.ones(len(levels), dtype=bool) while mask.any(): imax = numpy.argmax([mask[indices].sum() for trans, points, indices in cover]) trans, points, indices = cover.pop(imax) levels[indices] = levelset.eval(_transforms=(elem.transform + trans,), _points=points, _cache=fcache, **arguments) mask[indices] = False refs.append(elem.reference.trim(levels, maxrefine=maxrefine, ndivisions=ndivisions)) log.debug('cache', fcache.stats) return SubsetTopology(self, refs, newboundary=name)
[docs] def subset(self, elements, newboundary=None, strict=False): 'intersection' refs = [element.EmptyReference(self.ndims)] * len(self) for elem in elements: try: ielem = self.edict[elem.transform] except KeyError: assert not strict, 'elements do not form a strict subset' else: ref = self.elements[ielem].reference & elem.reference if strict: assert ref == elem.reference, 'elements do not form a strict subset' refs[ielem] = ref if not any(refs): return EmptyTopology(self.ndims) return SubsetTopology(self, refs, newboundary)
def withgroups(self, vgroups={}, bgroups={}, igroups={}, pgroups={}): return WithGroupsTopology(self, vgroups, bgroups, igroups, pgroups) if vgroups or bgroups or igroups or pgroups else self withsubdomain = lambda self, **kwargs: self.withgroups(vgroups=kwargs) withboundary = lambda self, **kwargs: self.withgroups(bgroups=kwargs) withinterfaces = lambda self, **kwargs: self.withgroups(igroups=kwargs) withpoints = lambda self, **kwargs: self.withgroups(pgroups=kwargs) @log.title @util.single_or_multiple def elem_project(self, funcs, degree, ischeme=None, check_exact=False, *, arguments=None): if arguments is None: arguments = {} if ischeme is None: ischeme = 'gauss{}'.format(degree*2) blocks = function.Tuple([function.Tuple([function.Tuple((function.Tuple(ind), f.simplified)) for ind, f in function.blocks(function.zero_argument_derivatives(func))]) for func in funcs]) bases = {} extractions = [[] for ifunc in range(len(funcs))] for elem in log.iter('elem', self): try: points, projector, basis = bases[elem.reference] except KeyError: points, weights = elem.reference.getischeme(ischeme) coeffs = elem.reference.get_poly_coeffs('bernstein', degree=degree) basis = numeric.poly_eval(coeffs[_], points) npoints, nfuncs = basis.shape A =, basis[:,:,_] * basis[:,_,:]) projector = numpy.linalg.solve(A, basis.T * weights) bases[elem.reference] = points, projector, basis for ifunc, ind_val in enumerate(blocks.eval(_transforms=(elem.transform, elem.opposite), _points=points, **arguments)): if len(ind_val) == 1: (allind, sumval), = ind_val else: allind, where = zip(*[numpy.unique([i for ind, val in ind_val for i in ind[iax]], return_inverse=True) for iax in range(funcs[ifunc].ndim)]) sumval = numpy.zeros([len(n) for n in (points,) + allind]) for ind, val in ind_val: I, where = zip(*[(w[:len(n)], w[len(n):]) for w, n in zip(where, ind)]), numpy.ix_(range(len(points)), *I), val) assert not any(where) ex =, sumval) if check_exact: numpy.testing.assert_almost_equal(sumval,, ex), decimal=15) extractions[ifunc].append((allind, ex)) return extractions @log.title def volume(self, geometry, ischeme='gauss', degree=1, *, arguments=None): return self.integrate(1, geometry=geometry, ischeme=ischeme, degree=degree, arguments=arguments) @log.title def check_boundary(self, geometry, elemwise=False, ischeme='gauss', degree=1, tol=1e-15, print=print, *, arguments=None): if elemwise: for elem in self: elem.reference.check_edges(tol=tol, print=print) volume = self.volume(geometry, ischeme=ischeme, degree=degree, arguments=arguments) zeros, volumes = self.boundary.integrate([geometry.normal(), geometry * geometry.normal()], geometry=geometry, ischeme=ischeme, degree=degree, arguments=arguments) if numpy.greater(abs(zeros), tol).any(): print('divergence check failed: {} != 0'.format(zeros)) if numpy.greater(abs(volumes - volume), tol).any(): print('divergence check failed: {} != {}'.format(volumes, volume)) def volume_check(self, geometry, ischeme='gauss', degree=1, decimal=15, *, arguments=None): warnings.deprecation('volume_check will be removed in future, us check_boundary instead') self.check_boundary(geometry=geometry, ischeme=ischeme, degree=degree, tol=10**-decimal, arguments=arguments) def indicator(self, subtopo): if isinstance(subtopo, str): subtopo = self[subtopo] transforms = tuple(sorted(elem.transform for elem in self)) values = numeric.const([int(trans in subtopo.edict) for trans in transforms]) assert len(subtopo) == values.sum(0), '{} is not a proper subtopology of {}'.format(subtopo, self) return function.Get(values, axis=0, item=function.FindTransform(transforms, function.Promote(self.ndims, trans=function.TRANS))) def select(self, indicator, ischeme='bezier2', *, arguments=None): values = self.elem_eval(indicator, ischeme, separate=True, arguments=arguments) selected = [elem for elem, value in zip(self, values) if numpy.greater(value, 0).any()] return UnstructuredTopology(self.ndims, selected) def prune_basis(self, basis): used = numpy.zeros(len(basis), dtype=bool) for axes, func in function.blocks(basis): dofmap = axes[0] for elem in self: dofs = dofmap.eval(_transforms=(elem.transform, elem.opposite)) used[dofs] = True return function.mask(basis, used) def locate(self, geom, points, ischeme='vertex', scale=1, tol=1e-12, eps=0, maxiter=100, *, arguments=None): nprocs = min(config.nprocs, len(self)) if arguments is None: arguments = {} if geom.ndim == 0: geom = geom[_] points = points[...,_] assert geom.shape == (self.ndims,) points = numpy.asarray(points, dtype=float) assert points.ndim == 2 and points.shape[1] == self.ndims vertices = self.elem_eval(geom, ischeme=ischeme, separate=True, arguments=arguments) bboxes = numpy.array([numpy.mean(v,axis=0) * (1-scale) + numpy.array([numpy.min(v,axis=0), numpy.max(v,axis=0)]) * scale for v in vertices]) # nelems x {min,max} x ndims vref = element.getsimplex(0) ielems = parallel.shempty(len(points), dtype=int) xis = parallel.shempty((len(points),len(geom)), dtype=float) for ipoint, point in parallel.pariter(log.enumerate('point', points), nprocs=nprocs): ielemcandidates, = numpy.logical_and(numpy.greater_equal(point, bboxes[:,0,:]), numpy.less_equal(point, bboxes[:,1,:])).all(axis=-1).nonzero() for ielem in sorted(ielemcandidates, key=lambda i: numpy.linalg.norm(bboxes[i].mean(0)-point)): converged = False elem = self.elements[ielem] xi, w = elem.reference.getischeme('gauss1') xi = (,xi) / w.sum())[_] if len(xi) > 1 else xi.copy() J = function.localgradient(geom, self.ndims) geom_J = function.Tuple((function.zero_argument_derivatives(geom), function.zero_argument_derivatives(J))).simplified for iiter in range(maxiter): point_xi, J_xi = geom_J.eval(_transforms=(elem.transform, elem.opposite), _points=xi, **arguments) err = numpy.linalg.norm(point - point_xi) if err < tol: converged = True break if iiter and err > prev_err: break prev_err = err xi += numpy.linalg.solve(J_xi, point - point_xi) if converged and elem.reference.inside(xi[0], eps=eps): ielems[ipoint] = ielem xis[ipoint], = xi break else: raise LocateError('failed to locate point: {}'.format(point)) pelems = [] for ielem, xi in zip(ielems, xis): elem = self.elements[ielem] trans = transform.Matrix(numpy.empty((len(xi), 0)), xi) pelems.append(element.Element(vref, elem.transform + (trans,), elem.opposite and elem.opposite + (trans,), oriented=True)) return UnstructuredTopology(0, pelems) def supp(self, basis, mask=None): if mask is None: mask = numpy.ones(len(basis), dtype=bool) elif isinstance(mask, list) or numeric.isarray(mask) and mask.dtype == int: tmp = numpy.zeros(len(basis), dtype=bool) tmp[mask] = True mask = tmp else: assert numeric.isarray(mask) and mask.dtype == bool and mask.shape == basis.shape[:1] indfunc = function.Tuple([ind[0] for ind, f in function.blocks(basis)]) subset = [] for elem in self: try: ind, = numpy.concatenate(indfunc.eval(_transforms=(elem.transform, elem.opposite)), axis=1) except function.EvaluationError: pass else: if mask[ind].any(): subset.append(elem) if not subset: return EmptyTopology(self.ndims) return self.subset(subset, newboundary='supp', strict=True) def revolved(self, geom): assert geom.ndim == 1 revdomain = self * RevolutionTopology() angle, = function.rootcoords(1) geom, angle = function.bifurcate(geom, angle) revgeom = function.concatenate([geom[0] * function.trignormal(angle), geom[1:]]) simplify = cache.replace(lambda op: function.zeros(()) if op is angle else None) return revdomain, revgeom, simplify def extruded(self, geom, nelems, periodic=False, bnames=('front','back')): assert geom.ndim == 1 root = transform.RootTrans('extrude', shape=[nelems if periodic else 0]) extopo = self * StructuredLine(root, i=0, j=nelems, periodic=periodic, bnames=bnames) exgeom = function.concatenate(function.bifurcate(geom, function.rootcoords(1))) return extopo, exgeom
class LocateError(Exception): pass
[docs]class WithGroupsTopology(Topology): 'item topology' def __init__(self, basetopo, vgroups={}, bgroups={}, igroups={}, pgroups={}): assert vgroups or bgroups or igroups or pgroups self.basetopo = basetopo self.vgroups = vgroups.copy() self.bgroups = bgroups.copy() self.igroups = igroups.copy() self.pgroups = pgroups.copy() super().__init__(basetopo.ndims) assert all(topo is Ellipsis or isinstance(topo, str) or isinstance(topo, Topology) and topo.ndims == basetopo.ndims and set(self.basetopo.edict).issuperset(topo.edict) for topo in self.vgroups.values()) def withgroups(self, vgroups={}, bgroups={}, igroups={}, pgroups={}): args = [] for groups, newgroups in (self.vgroups,vgroups), (self.bgroups,bgroups), (self.igroups,igroups), (self.pgroups,pgroups): groups = groups.copy() groups.update(newgroups) args.append(groups) return WithGroupsTopology(self.basetopo, *args) def __iter__(self): return iter(self.basetopo) def __len__(self): return len(self.basetopo) def getitem(self, item): if not isinstance(item, str): return self.basetopo.getitem(item) try: itemtopo = self.vgroups[item] except KeyError: return self.basetopo.getitem(item) else: return itemtopo if isinstance(itemtopo, Topology) else self.basetopo[itemtopo] @property def edict(self): return self.basetopo.edict @property def border_transforms(self): return self.basetopo.border_transforms @property def connectivity(self): return self.basetopo.connectivity @property def structure(self): return self.basetopo.structure @property def elements(self): return self.basetopo.elements @property def boundary(self): return self.basetopo.boundary.withgroups(self.bgroups) @property def interfaces(self): baseitopo = self.basetopo.interfaces # last minute orientation fix igroups = {name: UnstructuredTopology(self.ndims-1, [elem if elem.transform in baseitopo.edict else elem.flipped for elem in elems]) for name, elems in self.igroups.items()} return baseitopo.withgroups(igroups) @property def points(self): return self.basetopo.points.withgroups(self.pgroups) def basis(self, name, *args, **kwargs): return self.basetopo.basis(name, *args, **kwargs) def refined(self): groups = [{name: topo.refined if isinstance(topo,Topology) else topo for name, topo in groups.items()} for groups in (self.vgroups,self.bgroups,self.igroups,self.pgroups)] return self.basetopo.refined.withgroups(*groups)
[docs]class OppositeTopology(Topology): 'opposite topology' def __init__(self, basetopo): self.basetopo = basetopo super().__init__(basetopo.ndims) def getitem(self, item): return ~(self.basetopo.getitem(item)) def __iter__(self): return (elem.flipped for elem in self.basetopo) def __len__(self): return len(self.basetopo) def elements(self): return tuple(self) def __invert__(self): return self.basetopo
[docs]class EmptyTopology(Topology): 'empty topology' def __iter__(self): return iter([]) def __len__(self): return 0 def __or__(self, other): assert self.ndims == other.ndims return other def __rsub__(self, other): return other @property def elements(self): return ()
[docs]class Point(Topology): 'point' def __init__(self, trans, opposite=None): assert trans[-1].fromdims == 0 self.elem = element.Element(element.getsimplex(0), trans, opposite, oriented=True) super().__init__(ndims=0) def __iter__(self): yield self.elem @property def elements(self): return self.elem,
[docs]class StructuredLine(Topology): 'structured topology' def __init__(self, root, i, j, periodic=False, bnames=None): 'constructor' assert isinstance(i,int) and isinstance(j,int) and j > i assert not bnames or len(bnames) == 2 and all(isinstance(bname,str) for bname in bnames) assert isinstance(root, transform.TransformItem) self.root = root self.i = i self.j = j self.periodic = periodic self.bnames = bnames or () super().__init__(ndims=1) def _transforms(self): # one extra left and right for opposites, even if periodic=True return tuple((self.root, transform.Shift([float(offset)])) for offset in range(self.i-1, self.j+1)) def __iter__(self): reference = element.getsimplex(1) return (element.Element(reference, trans) for trans in self._transforms[1:-1]) def __len__(self): return self.j - self.i def elements(self): return tuple(self) def boundary(self): if self.periodic: return EmptyTopology(ndims=0) transforms = self._transforms right, left = element.LineReference().edge_transforms bnd = Point(transforms[1] + (left,), transforms[0] + (right,)), Point(transforms[-2] + (right,), transforms[-1] + (left,)) return UnionTopology(bnd, self.bnames) def interfaces(self): transforms = self._transforms right, left = element.LineReference().edge_transforms points = [Point(trans + (left,), opp + (right,)) for trans, opp in zip(transforms[2:-1], transforms[1:-2])] if self.periodic: points.append(Point(transforms[1] + (left,), transforms[-2] + (right,))) return UnionTopology(points) @classmethod def _bernstein_poly(cls, degree): 'bernstein polynomial coefficients' @classmethod def _spline_coeffs(cls, p, n): 'spline polynomial coefficients' assert p >= 0, 'invalid polynomial degree {}'.format(p) if p == 0: assert n == -1 return numpy.array([[[1.]]]) assert 1 <= n < 2*p extractions = numpy.empty((n, p+1, p+1)) extractions[0] = numpy.eye(p+1) for i in range(1, n): extractions[i] = numpy.eye(p+1) for j in range(2, p+1): for k in reversed(range(j, p+1)): alpha = 1. / min(2+k-j, n-i+1) extractions[i-1,:,k] = alpha * extractions[i-1,:,k] + (1-alpha) * extractions[i-1,:,k-1] extractions[i,-j-1:-1,-j-1] = extractions[i-1,-j:,-1] # magic bernstein triangle poly = numpy.zeros([p+1,p+1], dtype=int) for k in range(p//2+1): poly[k,k] = root = (-1)**p if k == 0 else (poly[k-1,k] * (k*2-1-p)) / k for i in range(k+1,p+1-k): poly[i,k] = poly[k,i] = root = (root * (k+i-p-1)) / i poly = poly[::-1].astype(float) return numeric.const(numeric.contract(extractions[:,_,:,:], poly[_,:,_,:], axis=-1).transpose(0,2,1), copy=False)
[docs] def basis_spline(self, degree, periodic=None, removedofs=None): 'spline from vertices' if periodic is None: periodic = self.periodic if numpy.iterable(degree): degree, = degree if numpy.iterable(removedofs): removedofs, = removedofs strides = 1, 1 shape = len(self), degree+1 ndofs = sum(s*(n-1) for s, n in zip(strides, shape))+1 dofs = numpy.arange(ndofs) if periodic and degree > 0: assert ndofs >= 2 * degree dofs[-degree:] = dofs[:degree] ndofs -= degree dofs = numpy.lib.stride_tricks.as_strided(dofs, shape=shape, strides=tuple(s*dofs.strides[0] for s in strides)) dofs = numeric.const(dofs, copy=False) p = degree n = 2*p-1 nelems = len(self) if periodic: if nelems == 1: # periodicity on one element can only mean a constant coeffs = [self._spline_coeffs(0, n)] dofs = numeric.const([[0]], copy=False) else: coeffs = list(self._spline_coeffs(p, n)[p-1:p]) * nelems else: coeffs = list(self._spline_coeffs(p, min(nelems,n))) if len(coeffs) < nelems: coeffs = coeffs[:p-1] + coeffs[p-1:p] * (nelems-2*(p-1)) + coeffs[p:] coeffs = numeric.const(coeffs, copy=False) func = function.polyfunc(coeffs, dofs, ndofs, self._transforms[1:-1], issorted=False) if not removedofs: return func mask = numpy.ones(ndofs, dtype=bool) mask[list(removedofs)] = False return function.mask(func, mask)
[docs] def basis_discont(self, degree): 'discontinuous shape functions' ref = element.LineReference() coeffs = [ref.get_poly_coeffs('bernstein', degree=degree)]*len(self) ndofs = ref.get_ndofs(degree) dofs = numeric.const(numpy.arange(ndofs*len(self), dtype=int).reshape(len(self), ndofs), copy=False) return function.polyfunc(coeffs, dofs, ndofs*len(self), self._transforms[1:-1], issorted=False)
[docs] def basis_std(self, degree, periodic=None, removedofs=None): 'spline from vertices' if periodic is None: periodic = self.periodic strides = max(1, degree), 1 shape = len(self), degree+1 ndofs = sum(s*(n-1) for s, n in zip(strides, shape))+1 dofs = numpy.arange(ndofs) if periodic and degree > 0: dofs[-1] = dofs[0] ndofs -= 1 dofs = numpy.lib.stride_tricks.as_strided(dofs, shape=shape, strides=tuple(s*dofs.strides[0] for s in strides)) dofs = numeric.const(dofs, copy=False) coeffs = [element.LineReference().get_poly_coeffs('bernstein', degree=degree)]*len(self) func = function.polyfunc(coeffs, dofs, ndofs, self._transforms[1:-1], issorted=False) if not removedofs: return func mask = numpy.ones(ndofs, dtype=bool) mask[list(removedofs)] = False return function.mask(func, mask)
def __str__(self): 'string representation' return '{}({}:{})'.format(self.__class__.__name__, self.i, self.j)
[docs]class StructuredTopology(Topology): 'structured topology' def __init__(self, root, axes, nrefine=0, bnames=None): 'constructor' self.root = root self.axes = tuple(axes) self.nrefine = nrefine self.shape = tuple(axis.j - axis.i for axis in self.axes if axis.isdim) if bnames is None: assert len(self.axes) <= 3 bnames = ('left', 'right'), ('bottom', 'top'), ('front', 'back') bnames = itertools.chain.from_iterable(n for axis, n in zip(self.axes, bnames) if axis.isdim and not axis.isperiodic) self._bnames = tuple(bnames) assert len(self._bnames) == sum(2 for axis in self.axes if axis.isdim and not axis.isperiodic) assert all(isinstance(bname,str) for bname in self._bnames) super().__init__(len(self.shape)) def __iter__(self): reference = util.product(element.getsimplex(1 if axis.isdim else 0) for axis in self.axes) return (element.Element(reference, trans, opp, oriented=True) for trans, opp in zip(self._transform.flat, self._opposite.flat)) def __len__(self): return, dtype=int) def getitem(self, item): if not isinstance(item, tuple): return EmptyTopology(self.ndims) assert all(isinstance(it,slice) for it in item) and len(item) <= self.ndims if all(it == slice(None) for it in item): # shortcut return self axes = [] idim = 0 for axis in self.axes: if axis.isdim and idim < len(item): s = item[idim] start, stop, stride = s.indices(axis.j - axis.i) assert stride == 1 assert stop > start if start > 0 or stop < axis.j - axis.i: axis = DimAxis(axis.i+start, axis.i+stop, isperiodic=False) idim += 1 axes.append(axis) return StructuredTopology(self.root, axes, self.nrefine, bnames=self._bnames) def elements(self): return tuple(self) @property def periodic(self): dimaxes = (axis for axis in self.axes if axis.isdim) return tuple(idim for idim, axis in enumerate(dimaxes) if axis.isdim and axis.isperiodic) @staticmethod def mktransforms(axes, root, nrefine): assert nrefine >= 0 updim = [] rmdims = numpy.zeros(len(axes), dtype=bool) for order, side, idim in sorted((axis.ibound, axis.side, idim) for idim, axis in enumerate(axes) if not axis.isdim): ref = util.product(element.getsimplex(0 if rmdim else 1) for rmdim in rmdims) iedge = (idim - rmdims[:idim].sum()) * 2 + 1 - side updim.append(ref.edge_transforms[iedge]) rmdims[idim] = True grid = [numpy.arange(axis.i>>nrefine, ((axis.j-1)>>nrefine)+1) if axis.isdim else numpy.array([(axis.i-1 if axis.side else axis.j)>>nrefine]) for axis in axes] indices = numeric.broadcast(*numeric.ix(grid)) transforms = numeric.asobjvector([transform.Shift(numpy.array(index, dtype=float))] for index in log.iter('elem', indices, indices.size)).reshape(indices.shape) if nrefine: scales = numeric.asobjvector([trans] for trans in (element.LineReference()**len(axes)).child_transforms).reshape((2,)*len(axes)) for irefine in log.range('level', nrefine-1, -1, -1): offsets = numpy.array([r[0] for r in grid]) grid = [numpy.arange(axis.i>>irefine,((axis.j-1)>>irefine)+1) if axis.isdim else numpy.array([(axis.i-1 if axis.side else axis.j)>>irefine]) for axis in axes] A = transforms[numpy.broadcast_arrays(*numeric.ix(r//2-o for r, o in zip(grid, offsets)))] B = scales[numpy.broadcast_arrays(*numeric.ix(r%2 for r in grid))] transforms = A + B shape = tuple(axis.j - axis.i for axis in axes if axis.isdim) return numeric.asobjvector(transform.canonical([root] + trans + updim) for trans in log.iter('canonical', transforms.flat)).reshape(shape) @log.title def _transform(self): return self.mktransforms(self.axes, self.root, self.nrefine) @log.title def _opposite(self): nbounds = len(self.axes) - self.ndims if nbounds == 0: return self._transform axes = [BndAxis(axis.i, axis.j, axis.ibound, not axis.side) if not axis.isdim and axis.ibound==nbounds-1 else axis for axis in self.axes] return self.mktransforms(axes, self.root, self.nrefine) @property def structure(self): warnings.deprecation('topology.structure will be removed in future') reference = util.product(element.getsimplex(1 if axis.isdim else 0) for axis in self.axes) return numeric.asobjvector(element.Element(reference, trans, opp, oriented=True) for trans, opp in numpy.broadcast(self._transform, self._opposite)).reshape(self.shape) def connectivity(self): connectivity = numpy.empty(self.shape+(self.ndims,2), dtype=int) connectivity[...] = -1 ielems = numpy.arange(len(self)).reshape(self.shape) for idim in range(self.ndims): s = (slice(None),)*idim s1 = s + (slice(1,None),) s2 = s + (slice(0,-1),) connectivity[s2+(...,idim,0)] = ielems[s1] connectivity[s1+(...,idim,1)] = ielems[s2] if idim in self.periodic: connectivity[s+(-1,...,idim,0)] = ielems[s+(0,)] connectivity[s+(0,...,idim,1)] = ielems[s+(-1,)] return connectivity.reshape(len(self), self.ndims*2) def boundary(self): 'boundary' nbounds = len(self.axes) - self.ndims btopo = EmptyTopology(self.ndims-1) jdim = 0 for idim, axis in enumerate(self.axes): if not axis.isdim or axis.isperiodic: continue btopos = [ StructuredTopology( root=self.root, axes=self.axes[:idim] + (BndAxis(n,n if not axis.isperiodic else 0,nbounds,side),) + self.axes[idim+1:], nrefine=self.nrefine, bnames=self._bnames[:jdim*2]+self._bnames[jdim*2+2:]) for side, n in enumerate((axis.i,axis.j)) ] btopo |= UnionTopology(btopos, self._bnames[jdim*2:jdim*2+2]) jdim += 1 return btopo def interfaces(self): 'interfaces' assert self.ndims > 0, 'zero-D topology has no interfaces' itopos = [] nbounds = len(self.axes) - self.ndims for idim, axis in enumerate(self.axes): if not axis.isdim: continue bndprops = [BndAxis(i, i, ibound=nbounds, side=True) for i in range(axis.i+1, axis.j)] if axis.isperiodic: assert axis.i == 0 bndprops.append(BndAxis(axis.j, 0, ibound=nbounds, side=True)) itopos.append(EmptyTopology(self.ndims-1) if not bndprops else UnionTopology(StructuredTopology(self.root, self.axes[:idim] + (axis,) + self.axes[idim+1:], self.nrefine) for axis in bndprops)) assert len(itopos) == self.ndims return UnionTopology(itopos, names=['dir{}'.format(idim) for idim in range(self.ndims)]) def _basis_spline(self, degree, knotvalues=None, knotmultiplicities=None, periodic=None): 'spline with structure information' if periodic is None: periodic = self.periodic if numeric.isint(degree): degree = [degree]*self.ndims assert len(degree)==self.ndims if knotvalues is None: knotvalues = [None]*self.ndims if knotmultiplicities is None: knotmultiplicities = [None]*self.ndims vertex_structure = numpy.array(0) stdelems = [] dofshape = [] slices = [] cache = {} for idim in range(self.ndims): p = degree[idim] n = self.shape[idim] isperiodic = idim in periodic k = knotvalues[idim] if k is None: #Defaults to uniform spacing k = numpy.arange(n+1) else: k = numpy.asarray(k) m = knotmultiplicities[idim] if m is None: #Defaults to open spline without internal repetitions m = numpy.array([p+1]+[1]*(n-1)+[p+1]) if not isperiodic else numpy.ones(n+1,dtype=int) else: m = numpy.array(m) #Make copy to prevent overwriting of data assert min(m)>0 and max(m)<=p+1, 'Incorrect multiplicity encountered' assert len(k)==len(m), 'Length mismatch between knots vector and knot multiplicities vector' assert len(k)==n+1, 'Knot vector size does not match the topology size' if not isperiodic: nd = sum(m)-p-1 npre = p+1-m[0] #Number of knots to be appended to front npost = p+1-m[-1] #Number of knots to be appended to rear m[0] = m[-1] = p+1 else: assert m[0]==m[-1], 'Periodic spline multiplicity expected' assert m[0]<p+1, 'Endpoint multiplicity for periodic spline should be p or smaller' nd = sum(m[:-1]) npre = npost = 0 k = numpy.concatenate([k[-p-1:-1]+k[0]-k[-1], k, k[1:1+p]-k[0]+k[-1]]) m = numpy.concatenate([m[-p-1:-1], m, m[1:1+p]]) km = numpy.array([ki for ki,mi in zip(k,m) for cnt in range(mi)],dtype=float) assert len(km)==sum(m) assert nd>0, 'No basis functions defined. Knot vector too short.' stdelems_i = [] slices_i = [] offsets = numpy.cumsum(m[:-1])-p if isperiodic: offsets = offsets[p:-p] offset0 = offsets[0]+npre for offset in offsets: start = max(offset0-offset,0) #Zero unless prepending influence stop = p+1-max(offset-offsets[-1]+npost,0) #Zero unless appending influence slices_i.append(slice(offset-offset0+start,offset-offset0+stop)) lknots = km[offset:offset+2*p] - km[offset] #Copy operation required if p: #Normalize for optimized caching lknots /= lknots[-1] key = (tuple(numeric.round(lknots*numpy.iinfo(numpy.int32).max)), p) try: coeffs = cache[key] except KeyError: coeffs = numeric.const(self._localsplinebasis(lknots, p).T, copy=False) cache[key] = coeffs stdelems_i.append(coeffs[start:stop]) stdelems.append(stdelems_i) numbers = numpy.arange(nd) if isperiodic: numbers = numpy.concatenate([numbers,numbers[:p]]) vertex_structure = vertex_structure[...,_]*nd+numbers dofshape.append(nd) slices.append(slices_i) #Cache effectivity log.debug('Local knot vector cache effectivity: {}'.format(100*(1.-len(cache)/float(sum(self.shape))))) # deduplicate stdelems and compute tensorial products `unique` with indices `index` # such that unique[index[i,j]] == poly_outer_product(stdelems[0][i], stdelems[1][j]) index = numpy.array(0) for stdelems_i in stdelems: unique_i = tuple(set(stdelems_i)) unique = unique_i if not index.ndim \ else [numeric.poly_outer_product(a, b) for a in unique for b in unique_i] index = index[...,_] * len(unique_i) + tuple(map(unique_i.index, stdelems_i)) coeffs = [unique[i] for i in index.flat] dofmap = [numeric.const(vertex_structure[S].ravel(), copy=False) for S in itertools.product(*slices)] return coeffs, dofmap, dofshape
[docs] def basis_spline(self, degree, knotvalues=None, knotmultiplicities=None, periodic=None, removedofs=None): 'spline basis' if removedofs == None: removedofs = [None] * self.ndims else: assert len(removedofs) == self.ndims coeffs, dofmap, dofshape = self._basis_spline(degree=degree, knotvalues=knotvalues, knotmultiplicities=knotmultiplicities, periodic=periodic) func = function.polyfunc(coeffs, dofmap, util.product(dofshape), (elem.transform for elem in self), issorted=False) if not any(removedofs): return func mask = numpy.ones((), dtype=bool) for idofs, ndofs in zip(removedofs, dofshape): mask = mask[...,_].repeat(ndofs, axis=-1) if idofs: mask[..., [numeric.normdim(ndofs,idof) for idof in idofs]] = False assert mask.shape == tuple(dofshape) return function.mask(func, mask.ravel())
def basis_bspline(self, *args, **kwargs): warnings.deprecation('basis "bspline" has been merged with "spline"') return self.basis_spline(*args, **kwargs) @staticmethod def _localsplinebasis (lknots, p): assert numeric.isarray(lknots), 'Local knot vector should be numpy array' assert len(lknots)==2*p, 'Expected 2*p local knots' #Based on Algorithm A2.2 Piegl and Tiller N = [None]*(p+1) N[0] = numpy.poly1d([1.]) if p > 0: assert numpy.less(lknots[:-1]-lknots[1:], numpy.spacing(1)).all(), 'Local knot vector should be non-decreasing' assert lknots[p]-lknots[p-1]>numpy.spacing(1), 'Element size should be positive' lknots = lknots.astype(float) xi = numpy.poly1d([lknots[p]-lknots[p-1],lknots[p-1]]) left = [None]*p right = [None]*p for i in range(p): left[i] = xi - lknots[p-i-1] right[i] = -xi + lknots[p+i] saved = 0. for r in range(i+1): temp = N[r]/(lknots[p+r]-lknots[p+r-i-1]) N[r] = saved+right[r]*temp saved = left[i-r]*temp N[i+1] = saved assert all(Ni.order==p for Ni in N) return numeric.const([Ni.coeffs for Ni in N]).T[::-1]
[docs] def basis_discont(self, degree): 'discontinuous shape functions' ref = util.product([element.LineReference()]*self.ndims) coeffs = [ref.get_poly_coeffs('bernstein', degree=degree)]*len(self) ndofs = ref.get_ndofs(degree) dofs = numeric.const(numpy.arange(ndofs*len(self), dtype=int).reshape(len(self), ndofs), copy=False) return function.polyfunc(coeffs, dofs, ndofs*len(self), (elem.transform for elem in self), issorted=False)
[docs] def basis_std(self, degree, removedofs=None, periodic=None): 'spline from vertices' if periodic is None: periodic = self.periodic if numeric.isint(degree): degree = (degree,) * self.ndims if removedofs == None: removedofs = [None] * self.ndims else: assert len(removedofs) == self.ndims dofshape = [] slices = [] vertex_structure = numpy.array(0) for idim in range(self.ndims): periodic_i = idim in periodic n = self.shape[idim] p = degree[idim] nd = n * p + 1 numbers = numpy.arange(nd) if periodic_i and p > 0: numbers[-1] = numbers[0] nd -= 1 vertex_structure = vertex_structure[...,_] * nd + numbers dofshape.append(nd) slices.append([slice(p*i,p*i+p+1) for i in range(n)]) lineref = element.LineReference() coeffs = [functools.reduce(numeric.poly_outer_product, (lineref.get_poly_coeffs('bernstein', degree=p) for p in degree))]*len(self) dofs = [numeric.const(vertex_structure[S].ravel(), copy=False) for S in numpy.broadcast(*numpy.ix_(*slices))] func = function.polyfunc(coeffs, dofs, numpy.product(dofshape), self._transform.ravel(), issorted=False) if not any(removedofs): return func mask = numpy.ones((), dtype=bool) for idofs, ndofs in zip(removedofs, dofshape): mask = mask[...,_].repeat(ndofs, axis=-1) if idofs: mask[..., [numeric.normdim(ndofs,idof) for idof in idofs]] = False assert mask.shape == tuple(dofshape) return function.mask(func, mask.ravel())
@property def refined(self): 'refine non-uniformly' axes = [DimAxis(i=axis.i*2,j=axis.j*2,isperiodic=axis.isperiodic) if axis.isdim else BndAxis(i=axis.i*2,j=axis.j*2,ibound=axis.ibound,side=axis.side) for axis in self.axes] return StructuredTopology(self.root, axes, self.nrefine+1, bnames=self._bnames) def __str__(self): 'string representation' return '{}({})'.format(self.__class__.__name__, 'x'.join(str(n) for n in self.shape))
[docs]class UnstructuredTopology(Topology): 'unstructured topology' def __init__(self, ndims, elements): self.elements = tuple(elements) assert all(elem.ndims == ndims for elem in self.elements) super().__init__(ndims) @log.title def connectivity(self): edges = {} connectivity = [] for ielem, elem in log.enumerate('elem', self): connectivity.append(-numpy.ones(elem.nedges, dtype=int)) for iedge, belem in enumerate(elem.edges): belemcoords = belem.vertices edgekey = tuple(sorted(belemcoords)) try: jelem, jedge = edges.pop(edgekey) except KeyError: edges[edgekey] = ielem, iedge else: # TODO assert transformation equivalence connectivity[ielem][iedge] = jelem connectivity[jelem][jedge] = ielem return tuple(connectivity) def boundary(self): elements = [elem.edge(iedge) for elem, ioppelems in zip(self, self.connectivity) for iedge in numpy.where(ioppelems == -1)[0]] return UnstructuredTopology(self.ndims-1, elements) def interfaces(self): seen = set() elements = [] for ielem, ioppelems in enumerate(self.connectivity): elem = self.elements[ielem] for iedge, ioppelem in enumerate(ioppelems): if ioppelem == -1: continue try: seen.remove((ielem, iedge)) except KeyError: ioppedge = tuple(self.connectivity[ioppelem]).index(ielem) oppelem = self.elements[ioppelem] elements.append(elem.edge(iedge).withopposite(oppelem.edge(ioppedge), oriented=False)) seen.add((ioppelem, ioppedge)) assert not seen return UnstructuredTopology(self.ndims-1, elements)
[docs] def basis_bubble(self): 'bubble from vertices' assert self.ndims == 2 coeffs = numeric.const([[[1,-1, 0, 0], [-1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 0, 0, 0], [ 1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 1, 0, 0], [ 0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 0, 0, 0], [ 0, 27,-27, 0], [0,-27, 0, 0], [0, 0, 0, 0]]]) nmap = [] dofmap = {} for ielem, elem in enumerate(self): assert isinstance(elem.reference, element.TriangleReference) assert elem.nverts == 3 dofs = numpy.empty(4, dtype=int) for i, v in enumerate(elem.vertices): dof = dofmap.get(v) if dof is None: dof = len(self) + len(dofmap) dofmap[v] = dof dofs[i] = dof dofs[3] = ielem nmap.append(numeric.const(dofs, copy=False)) ndofs = len(self)+len(dofmap) return function.polyfunc([coeffs]*len(self), nmap, ndofs, (elem.transform for elem in self), issorted=False)
def basis_spline(self, degree): assert degree == 1 return self.basis('std', degree)
[docs] def basis_discont(self, degree): 'discontinuous shape functions' assert numeric.isint(degree) and degree >= 0 coeffs = [] nmap = [] ndofs = 0 for elem in self: elemcoeffs = elem.reference.get_poly_coeffs('bernstein', degree=degree) coeffs.append(elemcoeffs) nmap.append(numeric.const(ndofs + numpy.arange(len(elemcoeffs)), copy=False)) ndofs += len(elemcoeffs) degrees = set(n-1 for c in coeffs for n in c.shape[1:]) return function.polyfunc(coeffs, nmap, ndofs, (elem.transform for elem in self), issorted=False)
def _basis_c0_structured(self, name, degree): 'C^0-continuous shape functions with lagrange stucture' assert numeric.isint(degree) and degree >= 0 if degree == 0: raise ValueError('Cannot build a C^0-continuous basis of degree 0. Use basis \'discont\' instead.') nlocaldofs = 0 elem_slices = [] coeffs = [] for elem in self: elem_coeffs = elem.reference.get_poly_coeffs(name, degree=degree) coeffs.append(elem_coeffs) elem_slices.append(slice(nlocaldofs, nlocaldofs+len(elem_coeffs))) nlocaldofs += len(elem_coeffs) dofmap = -numpy.ones([nlocaldofs], dtype=int) for ielem, elem in enumerate(self): # Loop over all neighbors of elem and merge dofs. for iedge, jelem in enumerate(self.connectivity[ielem]): if jelem < 0: continue jedge = self.connectivity[jelem].tolist().index(ielem) #if ielem < jelem or ielem == jelem and iedge < jedge: # continue idofs = elem.reference.get_edge_dofs(degree, iedge) verts = tuple(elem.edge(iedge).vertices) neighbor = self.elements[jelem] neighbor_verts = tuple(neighbor.edge(jedge).vertices) jdofs = neighbor.reference.get_edge_dofs(degree, jedge) #jdofs = jdofs[neighbor.reference.edges[jedge][1].get_dof_transpose_map(degree, tuple(map(neighbor_verts.index, verts)))] foo = neighbor.reference.edges[jedge][1].get_dof_transpose_map(degree, tuple(map(neighbor_verts.index, verts))) for idof, j in zip(idofs, foo): ridof = idof = elem_slices[ielem].start+idof # Resolve idof. while dofmap[ridof] >= 0: ridof = dofmap[ridof] # Resolve jdof. rjdof = jdof = elem_slices[jelem].start+jdofs[j] while dofmap[rjdof] >= 0: rjdof = dofmap[rjdof] if ridof != rjdof: dofmap[max(ridof,rjdof)] = min(ridof,rjdof) if ridof != idof: dofmap[idof] = min(ridof,rjdof) # Assign dof numbers. ndofs = 0 for i in range(len(dofmap)): if dofmap[i] < 0: dofmap[i] = ndofs ndofs += 1 else: dofmap[i] = dofmap[dofmap[i]] dofs = tuple(numeric.const(dofmap[s]) for s in elem_slices) return function.polyfunc(coeffs, dofs, ndofs, (elem.transform for elem in self), issorted=False)
[docs] def basis_lagrange(self, degree): 'lagrange shape functions' return self._basis_c0_structured('lagrange', degree)
[docs] def basis_bernstein(self, degree): 'bernstein shape functions' return self._basis_c0_structured('bernstein', degree)
basis_std = basis_bernstein
[docs]class UnionTopology(Topology): 'grouped topology' def __init__(self, topos, names=()): self._topos = tuple(topos) assert all(isinstance(topo, Topology) for topo in self._topos) self._names = tuple(names)[:len(self._topos)] assert all(isinstance(name,str) for name in self._names) assert len(set(self._names)) == len(self._names), 'duplicate name' ndims = self._topos[0].ndims assert all(topo.ndims == ndims for topo in self._topos) super().__init__(ndims) def getitem(self, item): topos = [topo if name == item else topo.getitem(item) for topo, name in itertools.zip_longest(self._topos, self._names)] return functools.reduce(operator.or_, topos, EmptyTopology(self.ndims)) def __or__(self, other): if not isinstance(other, UnionTopology): return UnionTopology(self._topos + (other,), self._names) return UnionTopology(self._topos[:len(self._names)] + other._topos + self._topos[len(self._names):], self._names + other._names) def elements(self): elements = [] for trans, elems in util.gather((elem.transform, elem) for topo in self._topos for elem in topo): if len(elems) == 1: elements.append(elems[0]) else: refs = [elem.reference for elem in elems] while len(refs) > 1: # sweep all possible unions until a single reference is left nrefs = len(refs) iref = 0 while iref < len(refs)-1: for jref in range(iref+1, len(refs)): try: unionref = refs[iref] | refs[jref] except TypeError: pass else: refs[iref] = unionref del refs[jref] break iref += 1 assert len(refs) < nrefs, 'incompatible elements in union' unionref, = refs opposite = elems[0].opposite assert all(elem.opposite == opposite for elem in elems[1:]) elements.append(element.Element(unionref, trans, opposite, oriented=True)) return elements @property def refined(self): return UnionTopology([topo.refined for topo in self._topos], self._names)
[docs]class SubsetTopology(Topology): 'trimmed' def __init__(self, basetopo, refs, newboundary=None): if newboundary is not None: assert isinstance(newboundary, str) or isinstance(newboundary, Topology) and newboundary.ndims == basetopo.ndims-1 assert len(refs) == len(basetopo) self.refs = tuple(refs) self.basetopo = basetopo self.newboundary = newboundary super().__init__(basetopo.ndims) def getitem(self, item): return self.basetopo.getitem(item).subset(self.elements, strict=False) def __rsub__(self, other): if self.basetopo == other: refs = [elem.reference - ref for elem, ref in zip(self.basetopo, self.refs)] return SubsetTopology(self.basetopo, refs, ~self.newboundary if isinstance(self.newboundary,Topology) else self.newboundary) return super().__rsub__(other) def __or__(self, other): if not isinstance(other, SubsetTopology) or self.basetopo != other.basetopo: return super().__or__(other) refs = [ref1 | ref2 for ref1, ref2 in zip(self.refs, other.refs)] if all(elem.reference == ref for elem, ref in zip(self.basetopo, refs)): return self.basetopo return SubsetTopology(self.basetopo, refs) # TODO boundary def connectivity(self): mask = numpy.array([bool(ref) for ref in self.refs] + [False]) # trailing false serves to map -1 to -1 renumber = numpy.cumsum(mask)-1 renumber[~mask] = -1 connectivity = tuple(tuple(renumber[ioppelems[iedge]] if edgeref and iedge < len(ioppelems) else -1 for iedge, edgeref in enumerate(ref.edge_refs)) for ref, ioppelems in zip(self.refs, self.basetopo.connectivity) if ref) return connectivity def elements(self): return tuple(element.Element(ref, elem.transform, elem.opposite) for elem, ref in zip(self.basetopo, self.refs) if ref) @property def refined(self): elems = [child for elem in self for child in elem.children if child] return self.basetopo.refined.subset(elems, self.newboundary.refined if isinstance(self.newboundary,Topology) else self.newboundary, strict=True) @log.title def boundary(self): brefs = [element.EmptyReference(self.ndims-1)] * len(self.basetopo.boundary) # subset of original boundary newbelems = [] # newly formed boundary elements connectivity = self.basetopo.connectivity elements = self.basetopo.elements baseboundary = self.basetopo.boundary for ielem, ioppelems in enumerate(connectivity): ref = self.refs[ielem] if not ref: continue elem = elements[ielem] newbelems.extend(element.Element(edge, elem.transform + (etrans,), elem.transform + (etrans.flipped,)) for etrans, edge in ref.edges[elem.reference.nedges:]) for iedge, ioppelem in enumerate(ioppelems): bref = ref.edge_refs[iedge] if not bref: continue if ioppelem == -1: index = baseboundary.edict[transform.canonical(elem.transform + (ref.edge_transforms[iedge],))] brefs[index] = bref # by construction, bref must be equal or subset of original else: ioppedge = tuple(connectivity[ioppelem]).index(ielem) oppref = self.refs[ioppelem] if oppref: bref -= oppref.edge_refs[ioppedge] if bref: newbelems.append(element.Element(bref, elem.transform + (ref.edge_transforms[iedge],), elements[ioppelem].edge(ioppedge).transform)) origboundary = SubsetTopology(baseboundary, brefs) trimboundary = OrientedGroupsTopology(self.newboundary if isinstance(self.newboundary,Topology) else self.basetopo.interfaces, newbelems) return UnionTopology([trimboundary, origboundary], names=[self.newboundary] if isinstance(self.newboundary,str) else []) @log.title def interfaces(self): irefs = [element.EmptyReference(self.ndims-1)] * len(self.basetopo.interfaces) # subset of original interfaces connectivity = self.basetopo.connectivity elements = self.basetopo.elements baseinterfaces = self.basetopo.interfaces for ielem, ioppelems in enumerate(connectivity): ref = self.refs[ielem] if not ref: continue # edge is empty elem = elements[ielem] for iedge, ioppelem in enumerate(ioppelems): if ioppelem == -1: continue # edge lies on the boundary oppref = self.refs[ioppelem] if not oppref: continue # edge is empty index = baseinterfaces.edict.get(transform.canonical(elem.transform + (ref.edge_transforms[iedge],))) if index is None: continue # edge is not oriented with an interface; rely on connectivity to also yield flipped element ioppedge = tuple(connectivity[ioppelem]).index(ielem) irefs[index] = ref.edge_refs[iedge] & oppref.edge_refs[ioppedge] return SubsetTopology(baseinterfaces, irefs) @log.title def basis(self, name, *args, **kwargs): if isinstance(self.basetopo, HierarchicalTopology): warnings.warn('basis may be linearly dependent; a linearly indepent basis is obtained by trimming first, then creating hierarchical refinements') basis = self.basetopo.basis(name, *args, **kwargs) return self.prune_basis(basis)
[docs]class OrientedGroupsTopology(UnstructuredTopology): 'unstructured topology with undirected semi-overlapping basetopology' def __init__(self, basetopo, elems): self.basetopo = basetopo super().__init__(basetopo.ndims, elems) def getitem(self, item): elements = [] for elem in self.basetopo.getitem(item): try: ielem, tail = transform.lookup_item(elem.transform, self.edict) except KeyError: elem = elem.flipped try: ielem, tail = transform.lookup_item(elem.transform, self.edict) except KeyError: continue if tail: raise NotImplementedError ref = self.elements[ielem].reference & elem.reference elements.append(element.Element(ref, elem.transform, elem.opposite, oriented=True)) return UnstructuredTopology(self.ndims, elements)
[docs]class RefinedTopology(Topology): 'refinement' def __init__(self, basetopo): self.basetopo = basetopo super().__init__(basetopo.ndims) def getitem(self, item): return self.basetopo.getitem(item).refined def elements(self): return tuple([child for elem in self.basetopo for child in elem.children]) def boundary(self): return self.basetopo.boundary.refined
[docs]class TrimmedTopologyItem(Topology): 'trimmed topology item' def __init__(self, basetopo, refdict): self.basetopo = basetopo self.refdict = refdict super().__init__(basetopo.ndims) def elements(self): elements = [] for elem in self.basetopo: ref = elem.reference & self.refdict[elem.transform] if ref: elements.append(element.Element(ref, elem.transform, elem.opposite, oriented=True)) return elements
[docs]class TrimmedTopologyBoundaryItem(Topology): 'trimmed topology boundary item' def __init__(self, btopo, trimmed, othertopo): self.btopo = btopo self.trimmed = trimmed self.othertopo = othertopo super().__init__(btopo.ndims) def elements(self): belems = [elem for elem in self.trimmed if elem.opposite in self.btopo.edict] if self.othertopo: belems.extend(self.othertopo) return belems
[docs]class HierarchicalTopology(Topology): 'collection of nested topology elments' def __init__(self, basetopo, allelements, precise): 'constructor' assert not isinstance(basetopo, HierarchicalTopology) self.basetopo = basetopo self.allelements = tuple(allelements) if precise: self.elements = self.allelements super().__init__(basetopo.ndims) def getitem(self, item): return self.basetopo.getitem(item).hierarchical(self.allelements, precise=False) def hierarchical(self, elements, precise=False): return self.basetopo.hierarchical(elements, precise) def elements(self): itemelems = [] for elem in self.allelements: try: ielem, tail = transform.lookup_item(elem.transform, self.basetopo.edict) except KeyError: continue itemelem = self.basetopo.elements[ielem] ref = itemelem.reference for trans in tail: index = ref.child_transforms.index(trans) ref = ref.child_refs[index] if not ref: break else: ref &= elem.reference if ref: itemelems.append(element.Element(ref, elem.transform, elem.opposite, oriented=True)) return itemelems @log.title def levels(self): levels = [self.basetopo] for elem in self: try: ielem, tail = transform.lookup_item(elem.transform, self.basetopo.edict) except KeyError: raise Exception('element is not a refinement of basetopo') else: nrefine = len(tail) while nrefine >= len(levels): levels.append(levels[-1].refined) assert elem.transform in levels[nrefine].edict, 'element is not a refinement of basetopo' return tuple(levels) def refined(self): elements = [child for elem in self for child in elem.children] return self.basetopo.hierarchical(elements, precise=True) @log.title def boundary(self): 'boundary elements' basebtopo = self.basetopo.boundary edgepool = [edge for elem in self if transform.lookup(elem.transform, self.basetopo.border_transforms) for edge in elem.edges if edge is not None] belems = [] for edge in edgepool: # superset of boundary elements try: iedge, tail = transform.lookup_item(edge.transform, basebtopo.edict) except KeyError: pass else: opptrans = basebtopo.elements[iedge].opposite + tail belems.append(element.Element(edge.reference, edge.transform, opptrans, oriented=True)) return basebtopo.hierarchical(belems, precise=True) @log.title def interfaces(self): 'interfaces' # Build a lookup table for level and element indices given elements in this # topology. elem_index_level = { elem: (ielem, ilevel) for ilevel, level in enumerate(self.levels) for ielem, elem in enumerate(level) } oriented = isinstance(self.basetopo, StructuredTopology) edict = self.edict interfaces = [] for elem in log.iter('elem', self): # Get `level`, element number at `level` of `elem`. ielem, ilevel = elem_index_level[elem] level = self.levels[ilevel] # Loop over neighbours of `elem`. for ielemedge, ineighbor in enumerate(level.connectivity[ielem]): if ineighbor < 0: # Not an interface. continue neighbor = level.elements[ineighbor] # Lookup `neighbor` (from the same `level` as `elem`) in this topology. head, tail = transform.lookup(neighbor.transform, edict) or (None, None) if not head: # `neighbor` not found, hence refinements of `neighbor` are present. # The interface of this edge will be added when we encounter the # refined elements. continue # Find the edge of `neighbor` between `neighbor` and `elem`. ineighboredge = tuple(level.connectivity[ineighbor]).index(ielem) if not tail and (ielem, ielemedge) > (ineighbor, ineighboredge): # `neighbor` itself, not a parent of, exists in this topology (`tail` # is empty). To make sure we add this interface only once we # continue here if the current element has a higher index (in # `level`) than the neighbor (or a higher edge number if the elements # are equal, which might occur when there is only one element in a # periodic dimension). continue # Create and add the interface between `elem` and `neighbor`. elemedge = elem.edges[ielemedge] neighboredge = neighbor.edges[ineighboredge] interfaces.append(element.Element(elemedge.reference, elemedge.transform, neighboredge.transform, oriented=oriented)) return UnstructuredTopology(self.ndims-1, interfaces)
[docs] @log.title def basis(self, name, *args, **kwargs): 'build hierarchical function space' # The law: a basis function is retained if all elements of self can # evaluate it through cascade, and at least one element of self can # evaluate it directly. # Procedure: per refinement level, track which basis functions have at # least one supporting element coinsiding with self ('touched') and no # supporting element finer than self ('supported'). dofs_coeffs = [] renumber = [] supports = [] length = 0 for topo in log.iter('level', self.levels): basis = topo.basis(name, *args, **kwargs) # shape functions for current level supported = numpy.ones(len(basis), dtype=bool) # True if dof is fully contained in self or parents touchtopo = numpy.zeros(len(basis), dtype=bool) # True if dof touches at least one elem in self (axes,func), = function.blocks(basis) dofmap, = axes if isinstance(func, function.Polyval): coeffs = func.coeffs assert coeffs.ndim == 1+self.ndims elif func.isconstant: assert func.ndim == 1 coeffs = func[(slice(None),*(_,)*self.ndims)] else: raise ValueError for elem in topo: trans = elem.transform idofs, = dofmap.eval(_transforms=(elem.transform, elem.opposite)) if trans in self.edict: touchtopo[idofs] = True elif transform.lookup(trans, self.edict): supported[idofs] = False support = supported & touchtopo supports.append(support) cumsum_support = numpy.cumsum(support) renumber.append(cumsum_support+(length-1)) length += cumsum_support[-1] dofs_coeffs.append(function.Tuple((dofmap, coeffs))) dofs = [] coeffs = [] transforms = tuple(sorted(elem.transform for elem in self)) for trans in transforms: hcoeffs = [] hdofs = [] ibase, tail = transform.lookup_item(trans, self.basetopo.edict) for ilevel in range(len(tail)+1): (idofs,), (icoeffs,) = dofs_coeffs[ilevel].eval(_transforms=(trans,)) isupport = supports[ilevel][idofs] if not isupport.any(): continue hdofs.extend(map(renumber[ilevel].__getitem__, idofs[isupport])) hcoeffs.extend(transform.transform_poly(tail[ilevel:], icoeffs[isupport])) dofs.append(numeric.const(hdofs)) coeffs.append(numeric.poly_stack(hcoeffs)) return function.polyfunc(coeffs, dofs, length, transforms, issorted=True)
[docs]class ProductTopology(Topology): 'product topology' def __init__(self, topo1, topo2): self.topo1 = topo1 self.topo2 = topo2 super().__init__(topo1.ndims+topo2.ndims) def __len__(self): return len(self.topo1) * len(self.topo2) def structure(self): return self.topo1.structure[(...,)+(_,)*self.topo2.ndims] * self.topo2.structure def elements(self): return (numpy.array(self.topo1.elements, dtype=object)[:,_] * numpy.array(self.topo2.elements, dtype=object)[_,:]).ravel() def __iter__(self): return self.elements.flat @property def refined(self): return self.topo1.refined * self.topo2.refined def refine(self, n): if numpy.iterable(n): assert len(n) == self.ndims else: n = (n,)*self.ndims return self.topo1.refine(n[:self.topo1.ndims]) * self.topo2.refine(n[self.topo1.ndims:]) def getitem(self, item): return self.topo1.getitem(item) * self.topo2 | self.topo1 * self.topo2.getitem(item) if isinstance(item, str) \ else self.topo1[item[:self.topo1.ndims]] * self.topo2[item[self.topo1.ndims:]] def basis(self, name, *args, **kwargs): def _split(arg): if isinstance(arg, (list,tuple)): assert len(arg) == self.ndims return arg[:self.topo1.ndims], arg[self.topo1.ndims:] else: return arg, arg splitargs = [_split(arg) for arg in args] splitkwargs = [(name,)+_split(arg) for name, arg in kwargs.items()] basis1, basis2 = function.bifurcate( self.topo1.basis(name, *[arg1 for arg1, arg2 in splitargs], **{name: arg1 for name, arg1, arg2 in splitkwargs}), self.topo2.basis(name, *[arg2 for arg1, arg2 in splitargs], **{name: arg2 for name, arg1, arg2 in splitkwargs})) return function.ravel(function.outer(basis1,basis2), axis=0) def boundary(self): return self.topo1 * self.topo2.boundary + self.topo1.boundary * self.topo2 def interfaces(self): return self.topo1 * self.topo2.interfaces + self.topo1.interfaces * self.topo2
[docs]class RevolutionTopology(Topology): 'topology consisting of a single revolution element' def __init__(self): self.elements = element.Element(element.RevolutionReference(), [transform.RootTrans('angle',(1,))]), self.boundary = EmptyTopology(ndims=0) super().__init__(ndims=1) def basis(self, name, *args, **kwargs): return function.asarray([1])
[docs]class MultipatchTopology(Topology): 'multipatch topology' Patch = collections.namedtuple('Patch', ['topo', 'verts', 'boundaries']) Patch.__qualname__ = 'MultipatchTopology.Patch'
[docs] @staticmethod def build_boundarydata(connectivity): 'build boundary data based on connectivity' boundarydata = [] for patch in connectivity: ndims = len(patch.shape) patchboundarydata = [] for dim, side in itertools.product(range(ndims), [-1, 0]): # ignore vertices at opposite face verts = numpy.array(patch) opposite = tuple({0:-1, -1:0}[side] if i == dim else slice(None) for i in range(ndims)) verts[opposite] = verts.max()+1 if len(set(verts.flat)) != 2**(ndims-1)+1: raise NotImplementedError('Cannot compute canonical boundary if vertices are used more than once.') # reverse axes such that lowest vertex index is at first position reverse = tuple(slice(None, None, -1) if i else slice(None) for i in numpy.unravel_index(verts.argmin(), verts.shape)) verts = verts[reverse] # transpose such that second lowest vertex connects to lowest vertex in first dimension, third in second dimension, et cetera k = [verts[tuple(1 if i == j else 0 for j in range(ndims))] for i in range(ndims)] transpose = tuple(sorted(range(ndims), key=k.__getitem__)) verts = verts.transpose(transpose) # boundarid boundaryid = tuple(verts[...,0].flat) patchboundarydata.append((boundaryid,dim,side,reverse,transpose)) boundarydata.append(tuple(patchboundarydata)) # TODO: boundary sanity checks return boundarydata
def __init__(self, patches): 'constructor' self.patches = tuple(self.Patch(*patch) for patch in patches) self._patchinterfaces = {} for patch in self.patches: for boundaryid, dim, side, reverse, transpose in patch.boundaries: self._patchinterfaces.setdefault(boundaryid, []).append((patch.topo, dim, side, reverse, transpose)) self._patchinterfaces = { boundaryid: tuple(data) for boundaryid, data in self._patchinterfaces.items() if len(data) > 1 } super().__init__(self.patches[0].topo.ndims) def elements(self): return tuple(itertools.chain.from_iterable(patch.topo for patch in self.patches)) def getitem(self, key): for i in range(len(self.patches)): if key == 'patch{}'.format(i): return self.patches[i].topo else: return UnionTopology(patch.topo.getitem(key) for patch in self.patches)
[docs] def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultiplicities=None): '''spline from vertices Create a spline basis with degree ``degree`` per patch. If ``patchcontinuous``` is true the basis is $C^0$-continuous at patch interfaces. ''' if knotvalues is None: knotvalues = {None: None} else: knotvalues, _knotvalues = {}, knotvalues for edge, k in _knotvalues.items(): if k is None: rk = None else: k = tuple(k) rk = k[::-1] if edge is None: knotvalues[edge] = k else: l, r = edge assert (l,r) not in knotvalues assert (r,l) not in knotvalues knotvalues[(l,r)] = k knotvalues[(r,l)] = rk if knotmultiplicities is None: knotmultiplicities = {None: None} else: knotmultiplicities, _knotmultiplicities = {}, knotmultiplicities for edge, k in _knotmultiplicities.items(): if k is None: rk = None else: k = tuple(k) rk = k[::-1] if edge is None: knotmultiplicities[edge] = k else: l, r = edge assert (l,r) not in knotmultiplicities assert (r,l) not in knotmultiplicities knotmultiplicities[(l,r)] = k knotmultiplicities[(r,l)] = rk missing = object() coeffs = [] dofmap = [] transforms = [] dofcount = 0 commonboundarydofs = {} for ipatch, patch in enumerate(self.patches): transforms.extend(elem.transform for elem in patch.topo) # build structured spline basis on patch `patch.topo` patchknotvalues = [] patchknotmultiplicities = [] for idim in range(self.ndims): left = tuple(0 if j == idim else slice(None) for j in range(self.ndims)) right = tuple(1 if j == idim else slice(None) for j in range(self.ndims)) dimknotvalues = set() dimknotmultiplicities = set() for edge in zip(patch.verts[left].flat, patch.verts[right].flat): v = knotvalues.get(edge, knotvalues.get(None, missing)) m = knotmultiplicities.get(edge, knotmultiplicities.get(None, missing)) if v is missing: raise 'missing edge' dimknotvalues.add(v) if m is missing: raise 'missing edge' dimknotmultiplicities.add(m) if len(dimknotvalues) != 1: raise 'ambiguous knot values for patch {}, dimension {}'.format(ipatch, idim) if len(dimknotmultiplicities) != 1: raise 'ambiguous knot multiplicities for patch {}, dimension {}'.format(ipatch, idim) patchknotvalues.append(next(iter(dimknotvalues))) patchknotmultiplicities.append(next(iter(dimknotmultiplicities))) patchcoeffs, patchdofmap, patchdofcount = patch.topo._basis_spline(degree, knotvalues=patchknotvalues, knotmultiplicities=patchknotmultiplicities) coeffs.extend(patchcoeffs) dofmap.extend(numeric.const(dofs+dofcount, copy=False) for dofs in patchdofmap) if patchcontinuous: # reconstruct multidimensional dof structure dofs = dofcount + numpy.arange(, dtype=int).reshape(patchdofcount) for boundaryid, dim, side, reverse, transpose in patch.boundaries: # get patch boundary dofs and reorder to canonical form boundarydofs = dofs[reverse].transpose(transpose)[...,0].ravel() # append boundary dofs to list (in increasing order, automatic by outer loop and dof increment) commonboundarydofs.setdefault(boundaryid, []).append(boundarydofs) dofcount += if patchcontinuous: # build merge mapping: merge common boundary dofs (from low to high) pairs = itertools.chain(*(zip(*dofs) for dofs in commonboundarydofs.values() if len(dofs) > 1)) merge = {} for dofs in sorted(pairs): dst = merge.get(dofs[0], dofs[0]) for src in dofs[1:]: merge[src] = dst # build renumber mapping: renumber remaining dofs consecutively, starting at 0 remainder = set(merge.get(dof, dof) for dof in range(dofcount)) renumber = dict(zip(sorted(remainder), range(len(remainder)))) # apply mappings dofmap = tuple(numeric.const(tuple(renumber[merge.get(dof, dof)] for dof in v.flat), dtype=int).reshape(v.shape) for v in dofmap) dofcount = len(remainder) return function.polyfunc(coeffs, dofmap, dofcount, transforms, issorted=False)
[docs] def basis_discont(self, degree): 'discontinuous shape functions' bases = [patch.topo.basis('discont', degree=degree) for patch in self.patches] coeffs = [] dofs = [] ndofs = 0 for patch in self.patches: basis = patch.topo.basis('discont', degree=degree) (axes,func), = function.blocks(basis) patch_dofmap, = axes if isinstance(func, function.Polyval): patch_coeffs = func.coeffs assert patch_coeffs.ndim == 1+self.ndims elif func.isconstant: assert func.ndim == 1 patch_coeffs = func[(slice(None),*(_,)*self.ndims)] else: raise ValueError patch_coeffs_dofs = function.Tuple((patch_coeffs, patch_dofmap)) for elem in patch.topo: (elem_coeffs,), (elem_dofs,) = patch_coeffs_dofs.eval(_transforms=(elem.transform,)) coeffs.append(elem_coeffs) dofs.append(numeric.const(ndofs+elem_dofs, copy=False)) ndofs += len(basis) return function.polyfunc(coeffs, dofs, ndofs, (elem.transform for patch in self.patches for elem in patch.topo), issorted=False)
[docs] def basis_patch(self): 'degree zero patchwise discontinuous basis' npatches = len(self.patches) coeffs = [numeric.const(1, dtype=int).reshape(1, *(1,)*self.ndims)]*npatches dofs = numeric.const(range(npatches), dtype=int)[:,_] return function.polyfunc(coeffs, dofs, npatches, ((patch.topo.root,) for patch in self.patches), issorted=False) def boundary(self): 'boundary' subtopos = [] subnames = [] for i, patch in enumerate(self.patches): names = dict(zip(itertools.product(range(self.ndims), [0,-1]), patch.topo._bnames)) for boundaryid, dim, side, reverse, transpose in patch.boundaries: if boundaryid in self._patchinterfaces: continue subtopos.append(patch.topo.boundary[names[dim,side]]) subnames.append('patch{}-{}'.format(i, names[dim,side])) if len(subtopos) == 0: return EmptyTopology(self.ndims-1) else: return UnionTopology(subtopos, subnames) def interfaces(self): '''interfaces Return a topology with all element interfaces. The patch interfaces are accessible via the group ``'interpatch'`` and the interfaces *inside* a patch via ``'intrapatch'``. ''' intrapatchtopo = EmptyTopology(self.ndims-1) if not self.patches else \ UnionTopology(patch.topo.interfaces for patch in self.patches) btopos = [] bconnectivity = [] for boundaryid, patchdata in self._patchinterfaces.items(): if len(patchdata) > 2: raise ValueError('Cannot create interfaces of multipatch topologies with more than two interface connections.') pairs = [] for topo, dim, side, reverse, transpose in patchdata: names = dict(zip(itertools.product(range(self.ndims), [0,-1]), topo._bnames)) # get structured set of boundary elements elems = topo.boundary[names[dim, side]].structure # add singleton axis elems = elems[tuple(_ if i == dim else slice(None) for i in range(self.ndims))] # apply canonical transformation elems = elems[reverse].transpose(transpose)[..., 0] shape = elems.shape pairs.append(elems.flat) # join element pairs elems = [ element.Element(elem_a.reference, elem_a.transform, elem_b.transform, oriented=True) for elem_a, elem_b in zip(*pairs) ] # create structured topology of joined element pairs bpatch = numpy.array(boundaryid).reshape((2,)*(self.ndims-1)) #btopos.append(StructuredTopology(numpy.array(elems).reshape(shape))) btopos.append(UnstructuredTopology(self.ndims-1, elems)) bconnectivity.append(bpatch) # create multipatch topology of interpatch boundaries interpatchtopo = MultipatchTopology(zip(btopos, bconnectivity, self.build_boundarydata(bconnectivity))) return UnionTopology((intrapatchtopo, interpatchtopo), ('intrapatch', 'interpatch')) def refined(self): 'refine' return MultipatchTopology((patch.topo.refined, patch.verts, patch.boundaries) for patch in self.patches)
# UTILITY FUNCTIONS DimAxis = collections.namedtuple('DimAxis', ['i','j','isperiodic']) DimAxis.isdim = True BndAxis = collections.namedtuple('BndAxis', ['i','j','ibound','side']) BndAxis.isdim = False def common_refine(topo1, topo2): warnings.deprecation('common_refine(a, b) will be removed in future; use a & b instead') return topo1 & topo2 # vim:shiftwidth=2:softtabstop=2:expandtab:foldmethod=indent:foldnestmax=2