from __future__ import absolute_import
from __future__ import print_function
import numpy as nm
from .bspline import BSpline
from sfepy.base.base import Struct
from six.moves import range
from sfepy.base.compat import lstsq
[docs]
class SplineBox(Struct):
"""
B-spline geometry parametrization. The geometry can be modified
by moving spline control points.
"""
[docs]
@staticmethod
def gen_cp_idxs(ncp):
dim = len(ncp)
if dim == 2:
idxs = nm.mgrid[0:ncp[0], 0:ncp[1]]
elif dim == 3:
idxs = nm.mgrid[0:ncp[0], 0:ncp[1], 0:ncp[2]]
else:
raise(ValueError)
cp_idxs = []
mul_idx = [1]
for ii in range(dim - 1):
mul_idx.append(mul_idx[ii] * ncp[ii])
cp_idxs = []
for ii in range(dim):
cp_idxs.append(idxs[ii].reshape(nm.prod(ncp), order='F'))
return cp_idxs, nm.array(mul_idx)
[docs]
@staticmethod
def create_spb(bbox, coors, degree=3, nsg=None):
nc, cdim = coors.shape
inside = nm.ones((nc,), dtype=bool)
nsg = nm.ones((cdim,), dtype=int) if nsg is None else nm.array(nsg)
for idim in range(cdim):
inrange = nm.logical_and(coors[:, idim] >= bbox[idim][0],
coors[:, idim] <= bbox[idim][1])
inside = nm.logical_and(inside, inrange)
ncp_tot = 1
base, uidx, ncp, cp = [], [], [], []
for idim in range(cdim):
ucoors, ucoors_idx = nm.unique(coors[inside, idim],
return_inverse=True)
ncp0 = degree + nsg[idim]
bspl = BSpline(degree, ncp=ncp0)
bspl.make_knot_vector(knot_range=(bbox[idim][0], bbox[idim][1]))
knots = bspl.get_knot_vector()
cp0 = nm.zeros((ncp0,), dtype=nm.double)
for j in range(cp0.shape[0]):
cp0[j] = nm.sum(knots[(j + 1):(j + degree + 1)]) / degree
base.append(bspl.eval_basis(t=ucoors, return_val=True))
uidx.append(ucoors_idx)
ncp.append(ncp0)
cp.append(cp0)
ncp_tot *= ncp0
cp_coors = nm.zeros((ncp_tot, cdim), dtype=nm.double)
cp_idx, mul_cp_idx = SplineBox.gen_cp_idxs(ncp)
for ii in range(cdim):
cp_coors[:, ii] = cp[ii][cp_idx[ii]]
return {'base': base,
'uidx': uidx,
'ncp': ncp,
'cp_idx': cp_idx,
'mul_cp_idx': mul_cp_idx,
'cp_coors': cp_coors,
'idxs_inside': inside}
def __init__(self, bbox, coors, nsg=None, field=None):
"""
Create a SplineBox.
Parameters
----------
bbox : array
The mesh bounding box.
coors : array
The coordinates of mesh nodes.
nsg : array
The number of segments.
"""
bbox = nm.asarray(bbox)
coors = nm.asarray(coors)
self.coors = coors.copy()
self.cdim = coors.shape[1]
if field is not None:
field = nm.asarray(field)
if len(field.shape) <= 1:
field = field[..., nm.newaxis]
self.field = field.copy()
else:
self.field = self.coors
self.__dict__.update(self.create_spb(bbox, coors, nsg=nsg))
if field is not None:
if hasattr(self, 'idxs_inside'):
b = field[self.idxs_inside, ...]
else:
b = field
a = self.get_box_matrix()
self.cp_values = lstsq(a, b)[0]
else:
self.cp_values = self.cp_coors
self.cp_values0 = self.cp_values.copy()
[docs]
def get_coors_shape(self):
"""
Get the shape of the coordinates.
"""
return self.coors.shape
[docs]
def get_control_points(self, init=False):
"""
Get the spline control points coordinates.
Returns
-------
cpt_coors : array
The coordinates of the spline control points.
init : bool
If True, return the initial state.
"""
if init:
return self.cp_values0
else:
return self.cp_values
[docs]
def set_control_points(self, cpt_coors, add=False):
"""
Set the spline control points position.
Parameters
----------
cpt_coors : array
The coordinates of the spline control points.
add : bool
If True, coors += cpt_coors
"""
if add:
self.cp_values += cpt_coors
else:
self.cp_values = cpt_coors
[docs]
def move_control_point(self, cpoint, val):
"""
Change shape of spline parametrization.
Parameters
----------
cpoint : int, list
The position (index or grid indicies) of the spline control point.
val : array
Displacement.
"""
if type(cpoint) in [list, tuple, nm.ndarray]:
idx = nm.dot(nm.array(cpoint), self.mul_cp_idx)
else:
idx = cpoint
self.cp_values[idx, :] += val
[docs]
def get_box_matrix(self):
"""
Returns:
mtx : 2D array
The matrix containing the coefficients of b-spline
basis functions.
"""
ncp, cdim = self.cp_coors.shape
mtx = nm.ones((self.uidx[0].shape[0], ncp), dtype=nm.double)
for ii in range(cdim):
mtx *= self.base[ii][self.uidx[ii], :][:, self.cp_idx[ii]]
return mtx
[docs]
def evaluate(self, cp_values=None, outside=True):
"""
Evaluate the new position of the mesh coordinates.
Parameters
----------
cp_values : array
The actual control point values. If None, use self.control_values.
outside : bool
If True, return also the coordinates outside the spline box.
Returns
-------
new_coors : array
The new position of the mesh coordinates.
"""
if cp_values is None:
cp_values = self.cp_values
mtx = self.get_box_matrix()
if outside and hasattr(self, 'idxs_inside'):
field = self.field.copy()
field[self.idxs_inside, ...] = nm.dot(mtx, cp_values)
return field
else:
return nm.dot(mtx, cp_values)
[docs]
def evaluate_derivative(self, cpoint, dirvec):
"""
Evaluate derivative of the spline
in a given control point and direction.
Parameters
----------
cpoint : int, list
The position (index or grid indicies) of the spline control point.
dirvec : array
The directional vector.
Returns
-------
diff : array
The derivative field.
"""
if type(cpoint) in [list, tuple, nm.ndarray]:
idxs = cpoint
else:
idxs = []
aux = cpoint
for ii in range(self.cdim):
idxs.append(aux // self.mul_cp_idx[-(ii + 1)])
aux = aux % self.mul_cp_idx[-(ii + 1)]
idxs = idxs[::-1]
aux = nm.ones((self.uidx[0].shape[0],), dtype=nm.double)
for ii in range(self.cdim):
aux *= self.base[ii][self.uidx[ii], idxs[ii]]
dirvec = nm.asarray(dirvec)
return nm.dot(aux[:, nm.newaxis],
nm.reshape(dirvec, (1, self.cp_values.shape[1])))
[docs]
def write_control_net(self, filename, deform_by_values=True):
"""
Write the SplineBox shape to the VTK file.
Parameters
----------
filename : str
The VTK file name.
"""
ncp = self.ncp
npt = nm.prod(ncp)
f = open(filename, 'w')
f.write("# vtk DataFile Version 2.6\nspbox file\n"
"ASCII\nDATASET UNSTRUCTURED_GRID\n\n")
f.write("POINTS %d float\n" % npt)
if self.cdim == 2:
ptformat = "%e %e 0.0\n"
elif self.cdim == 3:
ptformat = "%e %e %e\n"
cp_coors = self.cp_values if deform_by_values else self.cp_coors
for cpt in cp_coors:
f.write(ptformat % tuple(cpt))
cells = nm.array([nm.arange(0, ncp[0] - 1), nm.arange(1, ncp[0])]).T
cp = ncp[0]
nc = cp - 1
for ii in range(1, self.cdim):
cells1 = []
ncpi = ncp[ii]
for jj in range(ncpi):
cells1.append(cells + jj * cp)
nc = nc * ncpi
cells = nm.array([nm.arange(0, ncpi - 1),
nm.arange(1, ncpi)]).T * cp
for jj in range(cp):
cells1.append(cells + jj)
nc += (ncpi - 1) * cp
cells = nm.vstack(cells1)
cp *= ncp[ii]
f.write("\nCELLS %d %d\n" % (nc, 3 * nc))
for ii in cells:
f.write("2 %d %d\n" % tuple(ii))
f.write("\nCELL_TYPES %d\n" % nc)
f.write("3\n" * nc)
f.write("\nPOINT_DATA %d\n" % npt)
for ival in range(self.cp_values.shape[1]):
f.write("\nSCALARS cp_value_%d float 1\n" % (ival + 1))
f.write("LOOKUP_TABLE default\n")
f.write(str(b'\n'.join(self.cp_values[:, ival].astype('|S10'))) + '\n')
f.close()
[docs]
class SplineRegion2D(SplineBox):
"""
B-spline geometry parametrization. The boundary of the SplineRegion2D
is defined by BSpline curves.
"""
[docs]
@staticmethod
def points_in_poly(points, poly, tol=1e-6):
"""
Find which points are located inside the polygon.
"""
poly = nm.array(poly)
points = nm.array(points)
inside = nm.zeros((points.shape[0],), dtype=bool)
p1 = poly[:-1]
p2 = poly[1:]
a1 = (p2[:, 1] - p1[:, 1])
a1nz = nm.where(nm.fabs(a1) > 1e-16)[0]
a2 = (p2[a1nz, 0] - p1[a1nz, 0]) / a1[a1nz]
for jj, pt in enumerate(points):
# on edges?
if nm.any(nm.linalg.norm(p1 - pt, axis=1) +
nm.linalg.norm(p2 - pt, axis=1) -
nm.linalg.norm(p1 - p2, axis=1) < tol):
inside[jj] = True
continue
# inside?
val = nm.logical_and(
(p1[a1nz, 1] > pt[1]) != (p2[a1nz, 1] > pt[1]),
pt[0] < (a2*(pt[1] - p1[a1nz, 1]) + p1[a1nz, 0]))
if (nm.where(val)[0].shape[0] % 2) > 0:
inside[jj] = True
return nm.where(inside)[0]
[docs]
@staticmethod
def define_control_points(cp_bnd_coors, ncp):
"""
Find positions of "inner" control points depending on boundary splines.
"""
nx, ny = ncp
grid = nm.zeros(ncp, dtype=nm.int32)
grid.T.flat = nm.arange(nx * ny)
coors = nm.zeros((nx * ny, 2), dtype=nm.float64)
idxs1 = nm.arange(nx)
idxs2 = nm.arange(1, ny - 1) * nx
bcnd = nm.hstack([idxs1, idxs2 + nx - 1,
idxs1[::-1] + nx * (ny - 1), idxs2[::-1]])
coors[bcnd, :] = cp_bnd_coors
for ii in range(1, nx - 1):
for jj, t in enumerate(nm.linspace(0, 1, ny)[1:-1]):
c = (1 - t) * coors[ii, :] + t * coors[ii + (ny - 1)*nx, :]
coors[ii + nx*(jj + 1), :] = c
inside = grid[1:-1, 1:-1].flatten()
for iiter in range(5):
for ii in inside:
dx = nm.array([0., 0.])
for jj in [-1, +1, -nx, + nx]:
dx -= 0.25 * (coors[ii, :] - coors[ii + jj, :])
coors[ii] += 0.1 * dx
return coors
[docs]
@staticmethod
def create_spb(spl_bnd, coors, rho=10):
"""
Initialize SplineBox knots, control points, base functions, ...
"""
dim = 2
if coors.shape[1] != dim:
print('Only 2D SplineBoxSBnd is supported!')
raise(ValueError)
bnd_poly = []
bnd_cp = []
for s in spl_bnd:
s.set_param_n(int(rho))
bnd_poly.append(s.eval()[:-1])
bnd_cp.append(s.get_control_points()[:-1, :])
bnd_poly.append(bnd_poly[0][0, :])
ncpoints = 1
base, bspl, uidx, ncp = [], [], [], []
for si in [0, 1]:
s = spl_bnd[si]
bspl0 = BSpline(s.degree, ncp=s.ncp)
bspl0.set_knot_vector(s.knots)
bspl.append(bspl0)
base.append(None)
uidx.append(None)
ncp.append(s.ncp)
ncpoints *= s.ncp
cp_idx, mul_cp_idx = SplineBox.gen_cp_idxs(ncp)
cp_coors = SplineRegion2D.define_control_points(nm.vstack(bnd_cp), ncp)
idxs_inside = SplineRegion2D.points_in_poly(coors, nm.vstack(bnd_poly))
return {'base': base,
'bspl': bspl,
'uidx': uidx,
'ncp': ncp,
'cp_idx': cp_idx,
'mul_cp_idx': mul_cp_idx,
'cp_coors': cp_coors,
'idxs_inside': idxs_inside}
[docs]
def find_ts(self, coors):
"""
Function finds parameters (t, s) corresponding to given points (coors).
"""
from scipy.optimize import minimize
def ptdist(x, coors, spb):
for ii in range(spb.cdim):
spb.base[ii] = spb.bspl[ii].eval_basis(t=x[ii],
return_val=True)
coors_approx = spb.evaluate(outside=False)
return nm.linalg.norm(coors - coors_approx)
def gen_grid(spb, rho):
grid = nm.mgrid[0:rho, 0:rho]
t = nm.linspace(0, 1, rho)
for ii in range(spb.cdim):
spb.uidx[ii] = grid[ii, :].reshape(rho**self.cdim, order='F')
spb.base[ii] = spb.bspl[ii].eval_basis(t=t, return_val=True)
return spb.evaluate(outside=False)
rho = 100
grid = gen_grid(self, rho)
for ii in range(self.cdim):
self.uidx[ii] = nm.array([0])
ts = nm.zeros((coors.shape[0], self.cdim), dtype=nm.float64)
for ii, ic in enumerate(coors):
idx = nm.argmin(nm.linalg.norm(grid - ic, axis=1))
x0 = nm.array([idx % rho, idx // rho]) / (rho - 1.)
ts[ii] = minimize(lambda x: ptdist(x, ic, self), x0,
method='nelder-mead',
options={'xatol': 1e-5, 'disp': False}).x
return ts
def __init__(self, spl_bnd, coors, rho=1e3):
"""
Create a SplineBox which boundary is defined by B-spline curves.
Parameters
----------
spl_bnd : list
The list of BSpline objects (counterclockwise)
defining the SplineBox boundary.
coors : array
The coordinates of the mesh nodes.
rho : float
The density of points defining the boundary polygon.
"""
coors = nm.asarray(coors)
self.__dict__.update(self.create_spb(spl_bnd, coors, rho))
self.cdim = coors.shape[1]
self.coors = coors.copy()
self.field = self.coors
self.cp_values = self.cp_coors
self.ts = self.find_ts(coors[self.idxs_inside, :])
for idim in range(self.cdim):
ucoors, ucoors_idx = nm.unique(self.ts[:, idim],
return_inverse=True)
self.base[idim] = self.bspl[idim].eval_basis(t=ucoors,
return_val=True)
self.uidx[idim] = ucoors_idx