Source code for sfepy.discrete.fem.fe_surface

import numpy as nm

from sfepy.base.base import get_default, Struct
from sfepy.discrete.fem.facets import build_orientation_map
from sfepy.discrete.fem.utils import prepare_remap

[docs] class FESurface(Struct): """Description of a surface of a finite element domain.""" def __init__(self, name, region, efaces, volume_econn, volume_region=None): """nodes[leconn] == econn""" """nodes are sorted by node number -> same order as region.vertices""" self.name = get_default(name, 'surface_data_%s' % region.name) face_indices = region.get_facet_indices() faces = efaces[face_indices[:, 1]] if faces.size == 0 and not region.is_empty: raise ValueError('region with no faces! (%s)' % region.name) if volume_region is None: ii = face_indices[:, 0] elif hasattr(volume_region, 'get_cell_indices'): ii = volume_region.get_cell_indices(face_indices[:, 0]) else: ii = volume_region try: ee = volume_econn[ii] except: raise ValueError('missing region face indices! (%s)' % region.name) econn = nm.empty(faces.shape, dtype=nm.int32) for ir, face in enumerate(faces): econn[ir] = ee[ir, face] nodes = nm.unique(econn) if len(nodes): remap = prepare_remap(nodes, nodes.max() + 1) leconn = remap[econn].copy() else: leconn = econn.copy() n_fa, n_fp = econn.shape face_type = 's%d' % n_fp # Store bkey in SurfaceData, so that base function can be # queried later. bkey = 'b%s' % face_type[1:] self.econn = econn self.fis = nm.ascontiguousarray(face_indices.astype(nm.int32)) self.n_fa, self.n_fp = n_fa, n_fp self.nodes = nodes self.leconn = leconn self.face_type = face_type self.bkey = bkey self.meconn = {} self.mleconn = {} self.set_orientation_map()
[docs] @staticmethod def from_region(name, region, ret_gel=False): face_indices = region.get_facet_indices() cells = face_indices[:, 0] econn, gel = region.domain.get_conn(ret_gel=True, tdim=region.tdim, cells=cells) surface = FESurface(name, region, gel.get_surface_entities(), econn, slice(0, econn.shape[0])) if ret_gel: return surface, gel.surface_facet else: return surface
[docs] def set_orientation_map(self): n_fp = self.n_fp if n_fp <= 4: oo, _, _ = build_orientation_map(n_fp) ori_map = nm.zeros((nm.max(list(oo.keys())) + 1, n_fp), dtype=nm.int32) ori_map[list(oo.keys())] = nm.array([ii[1] for ii in oo.values()]) self.ori_map = ori_map else: self.ori_map = None
[docs] def setup_mirror_connectivity(self, region, mirror_name): """ Setup mirror surface connectivity required to integrate over a mirror region. 1. Get orientation of the faces: a) for own elements -> ooris b) for mirror elements -> moris 2. orientation -> permutation. """ def get_omap(reg, ori_map): if reg.tdim == (reg.dim - 1): cells = reg.get_cells() conn = reg.domain.get_conn(tdim=reg.tdim)[cells, :] return nm.argsort(conn) else: conn = reg.cmesh.get_conn_as_graph(reg.dim, reg.dim - 1) oris = reg.cmesh.facet_oris fis = reg.get_facet_indices() return ori_map[oris[conn.indptr[fis[:, 0]] + fis[:, 1]]] if mirror_name in self.meconn: return mregion = region.get_mirror_region(mirror_name) mmap = get_omap(mregion, self.ori_map) omap = get_omap(region, self.ori_map) econn = self.econn if mirror_name in region.mirror_maps\ and region.mirror_maps[mirror_name] is not None: mirror_map = region.mirror_maps[mirror_name] omap = omap[mirror_map] econn = econn[mirror_map] n_el, n_ep = econn.shape ii = nm.repeat(nm.arange(n_el)[:, None], n_ep, 1) meconn = nm.empty_like(econn) meconn[ii, mmap] = econn[ii, omap] nodes = nm.unique(meconn) remap = prepare_remap(nodes, nodes.max() + 1) self.meconn[mirror_name] = meconn self.mleconn[mirror_name] = remap[meconn].copy()
[docs] def get_connectivity(self, local=False, trace_region=None): """ Return the surface element connectivity. Parameters ---------- local : bool If True, return local connectivity w.r.t. surface nodes, otherwise return global connectivity w.r.t. all mesh nodes. trace_trace : None or str If not None, return mirror connectivity according to `local`. """ if trace_region is None: if local: return self.leconn else: return self.econn else: if local: return self.mleconn[trace_region] else: return self.meconn[trace_region]
[docs] class FEPhantomSurface(FESurface): """A phantom surface of the region with tdim=2.""" def __init__(self, name, region, volume_econn): self.name = get_default(name, 'surface_data_%s' % region.name) ii = region.get_cells() econn = volume_econn[ii] nodes = nm.unique(econn) if len(nodes): remap = prepare_remap(nodes, nodes.max() + 1) leconn = remap[econn].copy() else: leconn = econn.copy() n_fa, n_fp = econn.shape face_type = 's%d' % n_fp # Store bkey in SurfaceData, so that base function can be # queried later. bkey = 'b%s' % face_type[1:] self.econn = econn self.fis = -nm.ones((n_fa, 2), dtype=nm.int32) self.fis[:, 0] = ii self.n_fa, self.n_fp = n_fa, n_fp self.nodes = nodes self.leconn = leconn self.face_type = face_type self.bkey = bkey self.meconn = {} self.mleconn = {} self.set_orientation_map()