"""
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 approximation 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, trace_region=None, local=False):
"""
Get DOF connectivity of the given type in the given region.
"""
ct = conn_type[0] if isinstance(conn_type, tuple) else conn_type
nconn = self.nurbs.conn
if ct == 'cell':
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 == 'facet':
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='cell', region_name=None):
"""
Get element data dimensions.
Parameters
----------
integral : Integral instance
The integral describing used numerical quadrature.
integration : 'cell'
The term integration type. Only 'cell' 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('cell', region)
dofs = nm.unique(rconn)
else: # Facets.
rconn = self.get_econn('facet', 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 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