Source code for sfepy.discrete.fem.fields_hierarchic

import numpy as nm

from sfepy.base.base import assert_
from sfepy.discrete.fem.utils import prepare_remap, prepare_translate
from sfepy.discrete.common.dof_info import expand_nodes_to_dofs
from sfepy.discrete.fem.fields_base import VolumeField, H1Mixin

[docs]class H1HierarchicVolumeField(H1Mixin, VolumeField): family_name = 'volume_H1_lobatto' def _init_econn(self): """ Initialize the extended DOF connectivity and facet orientation array. """ VolumeField._init_econn(self) for ig, ap in self.aps.iteritems(): ap.ori = nm.zeros_like(ap.econn) def _setup_facet_orientations(self): self.node_desc = self.interp.describe_nodes() def _setup_edge_dofs(self): """ Setup edge DOF connectivity. """ if self.node_desc.edge is None: return 0, None, None return self._setup_facet_dofs(1, self.node_desc.edge, self.region.get_edges, self.n_vertex_dof) def _setup_face_dofs(self): """ Setup face DOF connectivity. """ if self.node_desc.face is None: return 0, None, None return self._setup_facet_dofs(self.domain.shape.tdim - 1, self.node_desc.face, self.region.get_faces, self.n_vertex_dof + self.n_edge_dof) def _setup_facet_dofs(self, dim, facet_desc, get_facets, offset): """ Helper function to setup facet DOF connectivity, works for both edges and faces. """ facet_desc = nm.array(facet_desc) n_dof_per_facet = facet_desc.shape[1] cmesh = self.domain.cmesh facets = self.region.entities[dim] ii = nm.arange(facets.shape[0], dtype=nm.int32) all_dofs = offset + expand_nodes_to_dofs(ii, n_dof_per_facet) # Prepare global facet id remapping to field-local numbering. remap = prepare_remap(facets, cmesh.num[dim]) cconn = self.region.domain.cmesh.get_conn(self.region.tdim, dim) offs = cconn.offsets n_f = self.gel.edges.shape[0] if dim == 1 else self.gel.faces.shape[0] n_fp = 2 if dim == 1 else self.gel.surface_facet.n_vertex oris = cmesh.get_orientations(dim) for ig, ap in self.aps.iteritems(): gcells = self.region.get_cells(ig, offset=False) n_el = gcells.shape[0] indices = cconn.indices[offs[gcells[0]]:offs[gcells[-1]+1]] facets_of_cells = remap[indices] # Define global facet dof numbers. gdofs = offset + expand_nodes_to_dofs(facets_of_cells, n_dof_per_facet) # Elements of facets. iel = nm.arange(n_el, dtype=nm.int32).repeat(n_f) ies = nm.tile(nm.arange(n_f, dtype=nm.int32), n_el) # DOF columns in econn for each facet (repeating same values for # each element. iep = facet_desc[ies] ap.econn[iel[:, None], iep] = gdofs ori = oris[offs[gcells[0]]:offs[gcells[-1]+1]] if (n_fp == 2) and (ap.interp.gel.name in ['2_4', '3_8']): tp_edges = ap.interp.gel.edges ecs = ap.interp.gel.coors[tp_edges] # True = positive, False = negative edge orientation w.r.t. # reference tensor product axes. tp_edge_ori = (nm.diff(ecs, axis=1).sum(axis=2) > 0).squeeze() aux = nm.tile(tp_edge_ori, n_el) ori = nm.where(aux, ori, 1 - ori) if n_fp == 2: # Edges. # ori == 1 means the basis has to be multiplied by -1. ps = ap.interp.poly_spaces['v'] orders = ps.node_orders eori = nm.repeat(ori[:, None], n_dof_per_facet, 1) eoo = orders[iep] % 2 # Odd orders. ap.ori[iel[:, None], iep] = eori * eoo elif n_fp == 3: # Triangular faces. raise NotImplementedError else: # Quadrilateral faces. # ori encoding in 3 bits: # 0: axis swap, 1: axis 1 sign, 2: axis 2 sign # 0 = + or False, 1 = - or True # 63 -> 000 = 0 # 0 -> 001 = 1 # 30 -> 010 = 2 # 33 -> 011 = 3 # 11 -> 100 = 4 # 7 -> 101 = 5 # 52 -> 110 = 6 # 56 -> 111 = 7 # Special cases: # Both orders same and even -> 000 # Both orders same and odd -> 0?? # Bits 1, 2 are multiplied by (swapped) axial order % 2. new = nm.repeat(nm.arange(8, dtype=nm.int32), 3) translate = prepare_translate([31, 59, 63, 0, 1, 4, 22, 30, 62, 32, 33, 41, 11, 15, 43, 3, 6, 7, 20, 52, 60, 48, 56, 57], new) ori = translate[ori] eori = nm.repeat(ori[:, None], n_dof_per_facet, 1) ps = ap.interp.poly_spaces['v'] orders = ps.face_axes_nodes[iep - ps.face_indx[0]] eoo = orders % 2 eoo0, eoo1 = eoo[..., 0], eoo[..., 1] i0 = nm.where(eori < 4) i1 = nm.where(eori >= 4) eori[i0] = nm.bitwise_and(eori[i0], 2*eoo0[i0] + 5) eori[i0] = nm.bitwise_and(eori[i0], eoo1[i0] + 6) eori[i1] = nm.bitwise_and(eori[i1], eoo0[i1] + 6) eori[i1] = nm.bitwise_and(eori[i1], 2*eoo1[i1] + 5) ap.ori[iel[:, None], iep] = eori n_dof = n_dof_per_facet * facets.shape[0] assert_(n_dof == nm.prod(all_dofs.shape)) return n_dof, all_dofs, remap def _setup_bubble_dofs(self): """ Setup bubble DOF connectivity. """ if self.node_desc.bubble is None: return 0, None, None offset = self.n_vertex_dof + self.n_edge_dof + self.n_face_dof n_dof = 0 n_dof_per_cell = self.node_desc.bubble.shape[0] all_dofs = {} remaps = {} for ig, ap in self.aps.iteritems(): ii = self.region.get_cells(ig) n_cell = ii.shape[0] nd = n_dof_per_cell * n_cell group = self.domain.groups[ig] remaps[ig] = prepare_remap(ii, group.shape.n_el) aux = nm.arange(offset + n_dof, offset + n_dof + nd, dtype=nm.int32) aux.shape = (n_cell, n_dof_per_cell) iep = self.node_desc.bubble[0] ap.econn[:,iep:] = aux all_dofs[ig] = aux n_dof += nd return n_dof, all_dofs, remaps
[docs] def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None): """ Set the values of given DOFs using a function of space coordinates or value `fun`. """ if region is None: region = self.region if dpn is None: dpn = self.n_components nods = [] vals = [] for ig in self.igs: if nm.isscalar(fun): # Hack - use only vertex DOFs. gnods = self.get_dofs_in_region_group(region, ig, merge=False) n_dof = dpn * sum([nn.shape[0] for nn in gnods]) gvals = nm.zeros(n_dof, dtype=nm.dtype(type(fun))) gvals[:gnods[0].shape[0] * dpn] = fun nods.append(nm.concatenate(gnods)) vals.append(gvals) elif callable(fun): # Hack - use only vertex DOFs. gnods = self.get_dofs_in_region_group(region, ig, merge=False) n_dof = dpn * sum([nn.shape[0] for nn in gnods]) vv = fun(self.get_coor(gnods[0])) gvals = nm.zeros(n_dof, dtype=vv.dtype) gvals[:gnods[0].shape[0] * dpn] = vv nods.append(nm.concatenate(gnods)) vals.append(gvals) else: raise NotImplementedError nods, indx = nm.unique(nm.concatenate(nods), return_index=True) ii = (nm.tile(dpn * indx, dpn) + nm.tile(nm.arange(dpn, dtype=nm.int32), indx.shape[0])) vals = nm.concatenate(vals)[ii] return nods, vals
[docs] def evaluate_at(self, coors, source_vals, strategy='kdtree', close_limit=0.1, cache=None, ret_cells=False, ret_status=False, ret_ref_coors=False, verbose=True): """ Evaluate source DOF values corresponding to the field in the given coordinates using the field interpolation. Parameters ---------- coors : array The coordinates the source values should be interpolated into. source_vals : array The source DOF values corresponding to the field. strategy : str, optional The strategy for finding the elements that contain the coordinates. Only 'kdtree' is supported for the moment. close_limit : float, optional The maximum limit distance of a point from the closest element allowed for extrapolation. cache : Struct, optional To speed up a sequence of evaluations, the field mesh, the inverse connectivity of the field mesh and the KDTree instance can be cached as `cache.mesh`, `cache.offsets`, `cache.iconn` and `cache.kdtree`. Optionally, the cache can also contain the reference element coordinates as `cache.ref_coors`, `cache.cells` and `cache.status`, if the evaluation occurs in the same coordinates repeatedly. In that case the KDTree related data are ignored. ret_cells : bool, optional If True, return also the cell indices the coordinates are in. ret_status : bool, optional If True, return also the status for each point: 0 is success, 1 is extrapolation within `close_limit`, 2 is extrapolation outside `close_limit`, 3 is failure. ret_ref_coors : bool, optional If True, return also the found reference element coordinates. verbose : bool If False, reduce verbosity. Returns ------- vals : array The interpolated values. cells : array The cell indices, if `ret_cells` or `ret_status` are True. status : array The status, if `ret_status` is True. """ raise NotImplementedError