Source code for sfepy.discrete.iga.fields

"""
Fields for isogeometric analysis.
"""
from __future__ import absolute_import
import numpy as nm

from sfepy.base.base import basestr, Struct
from sfepy.discrete.common.fields import parse_shape, Field
from sfepy.discrete.iga.mappings import IGMapping
from sfepy.discrete.iga.iga import get_bezier_element_entities

[docs]def parse_approx_order(approx_order): if (approx_order is None): return 0 aux = approx_order.split('+') if len(aux) == 2: return int(aux[1]) else: return 0
[docs]class IGField(Field): """ Bezier extraction based NURBS field for isogeometric analysis. Notes ----- The field has to cover the whole IGA domain. The field's NURBS basis can have higher degree than the domain NURBS basis. """ family_name = 'volume_H1_iga' def __init__(self, name, dtype, shape, region, approx_order=None, **kwargs): """ Create a Bezier element isogeometric analysis field. Parameters ---------- name : str The field name. dtype : numpy.dtype The field data type: float64 or complex128. shape : int/tuple/str The field shape: 1 or (1,) or 'scalar', space dimension (2, or (2,) or 3 or (3,)) or 'vector', or a tuple. The field shape determines the shape of the FE base functions and is related to the number of components of variables and to the DOF per node count, depending on the field kind. region : Region The region where the field is defined. approx_order : str or tuple, optional The field approximation order string or tuple with the first component in the form 'iga+<nonnegative int>'. Other components are ignored. The nonnegative int corresponds to the number of times the degree is elevated by one w.r.t. the domain NURBS description. **kwargs : dict Additional keyword arguments. """ shape = parse_shape(shape, region.domain.shape.dim) if approx_order is None: elevate_times = 0 else: if isinstance(approx_order, basestr): approx_order = (approx_order,) elevate_times = parse_approx_order(approx_order[0]) Struct.__init__(self, name=name, dtype=dtype, shape=shape, region=region, elevate_times=elevate_times) self.domain = self.region.domain self.nurbs = self.domain.nurbs.elevate(elevate_times) self._setup_kind() self.n_components = nm.prod(self.shape) self.val_shape = self.shape self.n_nod = self.nurbs.weights.shape[0] self.n_efun = nm.prod(self.nurbs.degrees + 1) self.approx_order = self.nurbs.degrees.max() self.mappings = {} self.is_surface = False def _get_facets(self, tdim): aux = get_bezier_element_entities(self.nurbs.degrees) return aux[2 - tdim]
[docs] def get_true_order(self): return nm.prod(self.nurbs.degrees)
[docs] def is_higher_order(self): """ Return True, if the field's approximation order is greater than one. """ return (self.nurbs.degrees > 1).any()
[docs] def get_econn(self, conn_type, region, is_trace=False, integration=None, local=False): """ Get DOF connectivity of the given type in the given region. """ ct = conn_type.type if isinstance(conn_type, Struct) else conn_type nconn = self.nurbs.conn if ct == 'volume': if region.name == self.region.name: conn = nconn else: cells = region.get_cells(true_cells_only=True) conn = nm.take(nconn, cells.astype(nm.int32), axis=0) elif ct == 'surface': fis = region.get_facet_indices() tdim = region.kind_tdim facets = self._get_facets(tdim) conn = [] for ii, fi in enumerate(fis): conn.append(nconn[fi[0], facets[fi[1]]]) if local: from sfepy.discrete.fem.utils import prepare_remap nods = nm.unique(conn) remap = prepare_remap(nods, nods.max() + 1) conn = [remap[ii] for ii in conn] else: raise ValueError('unsupported connectivity type! (%s)' % ct) return conn
[docs] def get_data_shape(self, integral, integration='volume', region_name=None): """ Get element data dimensions. Parameters ---------- integral : Integral instance The integral describing used numerical quadrature. integration : 'volume' The term integration type. Only 'volume' type is implemented. region_name : str The name of the region of the integral. Returns ------- data_shape : 4 ints The `(n_el, n_qp, dim, n_en)` for volume shape kind. Notes ----- - `n_el` = number of elements - `n_qp` = number of quadrature points per element/facet - `dim` = spatial dimension - `n_en` = number of element nodes """ region = self.domain.regions[region_name] _, weights = integral.get_qp(self.domain.gel.name) n_qp = weights.shape[0] data_shape = (region.shape.n_cell, n_qp, region.dim, self.n_efun) return data_shape
[docs] def get_dofs_in_region(self, region, merge=True): """ Return indices of DOFs that belong to the given region and group. Notes ----- `merge` is not used. """ idim = region.kind_tdim if idim < (self.domain.shape.tdim - 1): raise ValueError('region "%s" has no facets or cells!' % region.name) if idim == region.tdim: # Cells. rconn = self.get_econn('volume', region) dofs = nm.unique(rconn) else: # Facets. rconn = self.get_econn('surface', region) dofs = nm.unique(nm.concatenate(rconn)) if not merge: dofs = [dofs] return dofs
[docs] def get_surface_basis(self, region): from sfepy.discrete.integrals import Integral import sfepy.discrete.iga as iga from sfepy.discrete.iga.extmods.igac import eval_mapping_data_in_qp nurbs = self.nurbs facets = self._get_facets(region.kind_tdim) # Cell and face(cell) ids for each facet. fis = region.get_facet_indices() # Integral given by max. NURBS surface degree. fdegrees = iga.get_surface_degrees(nurbs.degrees) order = fdegrees.max() integral = Integral('i', order=2*order) vals, weights = integral.get_qp(self.domain.gel.surface_facet_name) # Boundary QP - use tensor product structure. bvals = iga.create_boundary_qp(vals, region.tdim) # Compute facet basis, jacobians and physical BQP. all_qp = [] all_fbfs = [] all_dets = [] for ii, (ie, ifa) in enumerate(fis): qp_coors = bvals[ifa] bfs, _, dets = eval_mapping_data_in_qp(qp_coors, nurbs.cps, nurbs.weights, nurbs.degrees, nurbs.cs, nurbs.conn, nm.array([ie])) # Facet basis. fbfs = bfs[..., facets[ifa]][0, :, 0, :] # Weight Jacobians by quadrature point weights. dets = nm.abs(dets) * weights[None, :, None, None] dets = dets[0, :, 0, :] # Physical BQP. fcps = nurbs.cps[nurbs.conn[ie, facets[ifa]]] qp = nm.dot(fbfs, fcps) all_qp.append(qp) all_fbfs.append(fbfs) all_dets.append(dets) return all_qp, all_fbfs, all_dets
[docs] def setup_extra_data(self, geometry, info, is_trace): dct = info.dc_type.type if dct != 'volume': raise ValueError('unknown dof connectivity type! (%s)' % dct)
[docs] def create_mapping(self, region, integral, integration): """ Create a new reference mapping. """ vals, weights = integral.get_qp(self.domain.gel.name) mapping = IGMapping(self.domain, region.cells, nurbs=self.nurbs) cmap = mapping.get_mapping(vals, weights) return cmap, mapping
[docs] def create_mesh(self, extra_nodes=True): """ Create a mesh corresponding to the field region. For IGA fields, this is directly the topological mesh. The `extra_nodes` argument is ignored. """ return self.domain.mesh
[docs] def create_eval_mesh(self): """ Create a mesh with the original NURBS connectivity for evaluating the field. The mesh coordinates are the NURBS control points. """ return self.domain.eval_mesh
[docs] def create_basis_context(self): """ Create the context required for evaluating the field basis. """ from sfepy.discrete.iga.extmods.igac import CNURBSContext nurbs = self.nurbs ctx = CNURBSContext(nurbs.cps, nurbs.weights, nurbs.degrees, nurbs.cs, nurbs.conn, i_max=100) return ctx
[docs] def create_output(self, dofs, var_name, dof_names=None, key=None, **kwargs): """ Convert the DOFs corresponding to the field to a dictionary of output data usable by Mesh.write(). Parameters ---------- dofs : array, shape (n_nod, n_component) The array of DOFs reshaped so that each column corresponds to one component. var_name : str The variable name corresponding to `dofs`. dof_names : tuple of str The names of DOF components. key : str, optional The key to be used in the output dictionary instead of the variable name. Returns ------- out : dict The output dictionary. """ from sfepy.discrete.iga.utils import create_mesh_and_output key = key if key is not None else var_name num = 25 if self.region.dim == 3 else 101 pars = (nm.linspace(0, 1, num),) * self.region.dim mesh, out = create_mesh_and_output(self.nurbs, pars, **{key : dofs}) out[key].var_name = var_name out[key].dofs = dof_names out['__mesh__'] = mesh return out