SfePy NTC

Source code for sfepy.discrete.fem.fields_nodal

"""
Notes
-----

Important attributes of continuous (order > 0) :class:`Field` and
:class:`SurfaceField` instances:

- `vertex_remap` : `econn[:, :n_vertex] = vertex_remap[conn]`
- `vertex_remap_i` : `conn = vertex_remap_i[econn[:, :n_vertex]]`

where `conn` is the mesh vertex connectivity, `econn` is the
region-local field connectivity.
"""
import numpy as nm

from sfepy.base.base import assert_
from sfepy.discrete.fem.utils import prepare_remap
from sfepy.discrete.common.dof_info import expand_nodes_to_dofs
from sfepy.discrete.fem.facets import get_facet_dof_permutations
from sfepy.discrete.fem.fields_base import (FEField, VolumeField, SurfaceField,
                                            H1Mixin)

[docs]class H1NodalMixin(H1Mixin): def _setup_facet_orientations(self): order = self.approx_order self.node_desc = self.poly_space.describe_nodes() edge_nodes = self.node_desc.edge_nodes if edge_nodes is not None: n_fp = self.gel.edges.shape[1] self.edge_dof_perms = get_facet_dof_permutations(n_fp, order) face_nodes = self.node_desc.face_nodes if face_nodes is not None: n_fp = self.gel.faces.shape[1] self.face_dof_perms = get_facet_dof_permutations(n_fp, order) 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.edge_dof_perms, 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.face_dof_perms, self.n_vertex_dof + self.n_edge_dof) def _setup_facet_dofs(self, dim, facet_desc, facet_perms, 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] oris = cmesh.get_orientations(dim) gcells = self.region.get_cells() n_el = gcells.shape[0] # 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) aux = offs[gcells][:, None] + ies.reshape((n_el, n_f)) indices = cconn.indices[aux] facets_of_cells = remap[indices].ravel() ori = oris[aux].ravel() perms = facet_perms[ori] # Define global facet dof numbers. gdofs = offset + expand_nodes_to_dofs(facets_of_cells, n_dof_per_facet) # DOF columns in econn for each facet. iep = facet_desc[ies] iaux = nm.arange(gdofs.shape[0], dtype=nm.int32) self.econn[iel[:, None], iep] = gdofs[iaux[:, None], perms] 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_per_cell = self.node_desc.bubble.shape[0] ii = self.region.get_cells() remap = prepare_remap(ii, self.domain.cmesh.n_el) n_cell = ii.shape[0] n_dof = n_dof_per_cell * n_cell all_dofs = nm.arange(offset, offset + n_dof, dtype=nm.int32) all_dofs.shape = (n_cell, n_dof_per_cell) iep = self.node_desc.bubble[0] self.econn[:,iep:] = all_dofs return n_dof, all_dofs, remap
[docs] def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None): """ Set the values of DOFs in a given region using a function of space coordinates or value `fun`. """ if region is None: region = self.region if dpn is None: dpn = self.n_components aux = self.get_dofs_in_region(region) nods = nm.unique(nm.hstack(aux)) if callable(fun): vals = fun(self.get_coor(nods)) elif nm.isscalar(fun): vals = nm.repeat([fun], nods.shape[0] * dpn) elif isinstance(fun, nm.ndarray): assert_(len(fun) == dpn) vals = nm.repeat(fun, nods.shape[0]) else: raise ValueError('unknown function/value type! (%s)' % type(fun)) return nods, vals
[docs] def create_basis_context(self): """ Create the context required for evaluating the field basis. """ ps = self.poly_space gps = self.gel.poly_space mesh = self.create_mesh(extra_nodes=False) ctx = ps.create_context(None, 0, 1e-15, 100, 1e-8, tdim=mesh.cmesh.tdim) geo_ctx = gps.create_context(mesh.cmesh, 0, 1e-15, 100, 1e-8) ctx.geo_ctx = geo_ctx return ctx
[docs]class H1NodalVolumeField(H1NodalMixin, VolumeField): family_name = 'volume_H1_lagrange'
[docs] def interp_v_vals_to_n_vals(self, vec): """ Interpolate a function defined by vertex DOF values using the FE geometry base (P1 or Q1) into the extra nodes, i.e. define the extra DOF values. """ if not self.node_desc.has_extra_nodes(): enod_vol_val = vec.copy() else: dim = vec.shape[1] enod_vol_val = nm.zeros((self.n_nod, dim), dtype=nm.float64) coors = self.poly_space.node_coors bf = self.gel.poly_space.eval_base(coors) bf = bf[:,0,:].copy() conn = self.econn[:, :self.gel.n_vertex] evec = nm.dot(bf, vec[conn]) enod_vol_val[self.econn] = nm.swapaxes(evec, 0, 1) return enod_vol_val
[docs]class H1DiscontinuousField(H1NodalMixin, VolumeField): family_name = 'volume_H1_lagrange_discontinuous' def _setup_global_base(self): """ Setup global DOF/base function indices and connectivity of the field. """ self._setup_facet_orientations() self._init_econn() ii = self.region.get_cells() self.bubble_remap = prepare_remap(ii, self.domain.cmesh.n_el) n_dof = nm.prod(self.econn.shape) all_dofs = nm.arange(n_dof, dtype=nm.int32) all_dofs.shape = self.econn.shape self.econn[:] = all_dofs self.n_nod = n_dof self.n_bubble_dof = n_dof self.bubble_dofs = all_dofs self.n_vertex_dof = self.n_edge_dof = self.n_face_dof = 0 self._setup_esurface()
[docs] def extend_dofs(self, dofs, fill_value=None): """ Extend DOFs to the whole domain using the `fill_value`, or the smallest value in `dofs` if `fill_value` is None. """ if self.approx_order != 0: dofs = self.average_to_vertices(dofs) new_dofs = FEField.extend_dofs(self, dofs) return new_dofs
[docs] def remove_extra_dofs(self, dofs): """ Remove DOFs defined in higher order nodes (order > 1). """ if self.approx_order != 0: dofs = self.average_to_vertices(dofs) new_dofs = FEField.remove_extra_dofs(self, dofs) return new_dofs
[docs] def average_to_vertices(self, dofs): """ Average DOFs of the discontinuous field into the field region vertices. """ data_qp, integral = self.interp_to_qp(dofs) vertex_dofs = self.average_qp_to_vertices(data_qp, integral) return vertex_dofs
[docs]class H1NodalSurfaceField(H1NodalMixin, SurfaceField): """ A field defined on a surface region. """ family_name = 'surface_H1_lagrange'
[docs] def interp_v_vals_to_n_vals(self, vec): """ Interpolate a function defined by vertex DOF values using the FE surface geometry base (P1 or Q1) into the extra nodes, i.e. define the extra DOF values. """ if not self.node_desc.has_extra_nodes(): enod_vol_val = vec.copy() else: msg = 'surface nodal fields do not support higher order nodes yet!' raise NotImplementedError(msg) return enod_vol_val