"""
Computational domain for isogeometric analysis.
"""
from __future__ import absolute_import
import os.path as op
import numpy as nm
from sfepy.base.base import assert_, Struct
from sfepy.discrete.common.domain import Domain
from sfepy.discrete.iga import iga
from sfepy.discrete.iga import io
from sfepy.discrete.iga.extmods.igac import eval_in_tp_coors
import six
from six.moves import range
[docs]
class NurbsPatch(Struct):
"""
Single NURBS patch data.
"""
def __init__(self, knots, degrees, cps,
weights, cs, conn):
degrees = nm.asarray(degrees, dtype=nm.int32)
cs = [nm.asarray(cc, dtype=nm.float64) for cc in cs]
if cs[0].ndim == 3:
cs = [nm.ascontiguousarray(cc[:, None, ...]) for cc in cs]
Struct.__init__(self, name='nurbs', knots=knots, degrees=degrees,
cps=cps, weights=weights, cs=cs, conn=conn)
self.n_els = [len(ii) for ii in cs]
self.dim = len(self.n_els)
def _get_ref_coors_1d(self, pars, axis):
uk = nm.unique(self.knots[axis])
indices = nm.searchsorted(uk[1:], pars)
ref_coors = nm.empty_like(pars)
for ii in range(len(uk) - 1):
ispan = nm.where(indices == ii)[0]
pp = pars[ispan]
ref_coors[ispan] = (pp - uk[ii]) / (uk[ii+1] - uk[ii])
return uk, indices, ref_coors
def __call__(self, u=None, v=None, w=None, field=None):
"""
Igakit-like interface for NURBS evaluation.
"""
pars = [u]
if v is not None: pars += [v]
if w is not None: pars += [w]
indices = []
rcs = []
for ia, par in enumerate(pars):
uk, indx, rc = self._get_ref_coors_1d(par, ia)
indices.append(indx.astype(nm.uint32))
rcs.append(rc)
out = eval_in_tp_coors(field, indices,
rcs, self.cps, self.weights,
self.degrees,
self.cs, self.conn)
return out
[docs]
def evaluate(self, field, u=None, v=None, w=None):
"""
Igakit-like interface for NURBS evaluation.
"""
return self(u, v, w, field)
def _to_igakit(self):
import igakit.cad as cad
n_efuns = self.degrees + 1
nks = nm.array([len(ii) for ii in self.knots])
shape = tuple(nks - n_efuns)
cps = self.cps.reshape(shape + (-1,))
weights = self.weights.reshape(shape)
return cad.NURBS(self.knots, cps, weights=weights)
def _from_igakit(self, inurbs):
cs = iga.compute_bezier_extraction(inurbs.knots, inurbs.degree)
n_els = [len(ii) for ii in cs]
conn, bconn = iga.create_connectivity(n_els, inurbs.knots,
inurbs.degree)
cps = inurbs.points[..., :self.dim].copy()
cps = cps.reshape((-1, self.dim))
return NurbsPatch(inurbs.knots, inurbs.degree, cps,
inurbs.weights.ravel(), cs, conn)
[docs]
def elevate(self, times=0):
"""
Elevate the patch degrees several `times` by one.
Returns
-------
nurbs : NurbsPatch instance
Either `self` if `times` is zero, or a new instance.
"""
if times == 0: return self
aux = self._to_igakit()
for ia in range(self.dim):
aux.elevate(ia, times)
assert_(nm.isfinite(aux.points).all(),
'igakit degree elevation failed for axis %d!' % ia)
return self._from_igakit(aux)
[docs]
class IGDomain(Domain):
"""
Bezier extraction based NURBS domain for isogeometric analysis.
"""
[docs]
@staticmethod
def from_file(filename):
"""
filename : str
The name of the IGA domain file.
"""
data = io.read_iga_data(filename)
name = op.splitext(filename)[0]
return IGDomain.from_data(*(data + (name,)))
[docs]
@staticmethod
def read_domain_from_hdf5(fd, group):
"""
Create a domain from the given hdf5 data group.
fd: tables.File
HDF5 file handle to read the mesh from.
group: tables.group.Group
HDF5 data group (of file fd) to read the mesh from.
"""
data = io.read_iga_data(fd, group)
return IGDomain.from_data(*data)
[docs]
def write_domain_to_hdf5(self, fd, group):
"""
Save the domain to a hdf5 file.
fd: tables.File
HDF5 file handle to write the mesh to.
group: tables.group.Group
HDF5 data group (of file fd) to write the mesh to.
"""
io.write_iga_data(fd, group, *(self._get_io_data() + (self.name,)))
def _get_io_data(self):
"""
Return the data describing the domain for storing the domain in a hdf5
file.
TODO - data for regions recreating
"""
knots = self.nurbs.knots
degrees = self.nurbs.degrees
cps = self.nurbs.cps
weights = self.nurbs.weights
cs = self.nurbs.cs
conn = self.nurbs.conn
bcps = self.bmesh.cps
bweights = self.bmesh.weights
bconn = self.bmesh.conn
return (knots, degrees, cps, weights, cs, conn,
bcps, bweights, bconn, self.vertex_set_bcs)
[docs]
@staticmethod
def from_data(knots, degrees, cps, weights, cs, conn,
bcps, bweights, bconn, regions, name='iga_domain_from_data'):
"""
Create the IGA domain from the given data.
"""
nurbs = NurbsPatch(knots, degrees, cps, weights, cs, conn)
bmesh = Struct(name='bmesh', cps=bcps, weights=bweights, conn=bconn)
domain = IGDomain(name, nurbs=nurbs, bmesh=bmesh, regions=regions)
return domain
def __init__(self, name, nurbs, bmesh, regions=None, **kwargs):
"""
Create an IGA domain.
Parameters
----------
name : str
The domain name.
"""
Domain.__init__(self, name, nurbs=nurbs, bmesh=bmesh, regions=regions,
**kwargs)
from sfepy.discrete.fem.geometry_element import create_geometry_elements
from sfepy.discrete.fem import Mesh
from sfepy.discrete.fem.utils import prepare_remap
tconn = iga.get_bezier_topology(bmesh.conn, nurbs.degrees)
itc = nm.unique(tconn)
remap = prepare_remap(itc, bmesh.conn.max() + 1)
ltcoors = bmesh.cps[itc]
ltconn = remap[tconn]
n_nod, dim = ltcoors.shape
n_el = ltconn.shape[0]
self.shape = Struct(n_nod=n_nod, dim=dim, tdim=0, n_el=n_el)
desc = '%d_%d' % (dim, bmesh.conn.shape[1])
mat_id = nm.zeros(bmesh.conn.shape[0], dtype=nm.int32)
eval_mesh = Mesh.from_data(self.name + '_eval', nurbs.cps, None,
[nurbs.conn], [mat_id], [desc])
self.eval_mesh = eval_mesh
desc = '%d_%d' % (dim, 2**dim)
mat_id = nm.zeros(ltconn.shape[0], dtype=nm.int32)
self.mesh = Mesh.from_data(self.name + '_topo', ltcoors, None, [ltconn],
[mat_id], [desc])
self.cmesh = self.mesh.cmesh
self.cmesh_tdim = self.mesh.cmesh_tdim
gels = create_geometry_elements()
self.cmesh.set_local_entities(gels)
self.cmesh.setup_entities()
self.shape.tdim = self.cmesh.tdim
self.gel = gels[desc]
if regions is not None:
self.vertex_set_bcs = {}
for key, val in six.iteritems(self.regions):
self.vertex_set_bcs[key] = remap[val]
self.reset_regions()