"""

def __init__(self, triangulation, z, trifinder=None):

cbook._check_isinstance(Triangulation, triangulation=triangulation)

self._triangulation = triangulation

self._z = np.asarray(z)

if self._z.shape != self._triangulation.x.shape:

raise ValueError("z array must have same length as triangulation x"

" and y arrays")

cbook._check_isinstance((TriFinder, None), trifinder=trifinder)

self._trifinder = trifinder or self._triangulation.get_trifinder()

# Default scaling factors : 1.0 (= no scaling)

# Scaling may be used for interpolations for which the order of

# magnitude of x, y has an impact on the interpolant definition.

# Please refer to :meth:`_interpolate_multikeys` for details.

self._unit_x = 1.0

self._unit_y = 1.0

# Default triangle renumbering: None (= no renumbering)

# Renumbering may be used to avoid unnecessary computations

# if complex calculations are done inside the Interpolator.

# Please refer to :meth:`_interpolate_multikeys` for details.

self._tri_renum = None

# __call__ and gradient docstrings are shared by all subclasses

# (except, if needed, relevant additions).

# However these methods are only implemented in subclasses to avoid

# confusion in the documentation.

_docstring__call__ = """

_docstringgradient = r"""

def _interpolate_multikeys(self, x, y, tri_index=None,

return_keys=('z',)):

"""

# Flattening and rescaling inputs arrays x, y

# (initial shape is stored for output)

x = np.asarray(x, dtype=np.float64)

y = np.asarray(y, dtype=np.float64)

sh_ret = x.shape

if x.shape != y.shape:

raise ValueError("x and y shall have same shapes."

" Given: {0} and {1}".format(x.shape, y.shape))

x = np.ravel(x)

y = np.ravel(y)

x_scaled = x/self._unit_x

y_scaled = y/self._unit_y

size_ret = np.size(x_scaled)

# Computes & ravels the element indexes, extract the valid ones.

if tri_index is None:

tri_index = self._trifinder(x, y)

else:

if tri_index.shape != sh_ret:

raise ValueError(

"tri_index array is provided and shall"

" have same shape as x and y. Given: "

"{0} and {1}".format(tri_index.shape, sh_ret))

tri_index = np.ravel(tri_index)

mask_in = (tri_index != -1)

if self._tri_renum is None:

valid_tri_index = tri_index[mask_in]

else:

valid_tri_index = self._tri_renum[tri_index[mask_in]]

valid_x = x_scaled[mask_in]

valid_y = y_scaled[mask_in]

ret = []

for return_key in return_keys:

# Find the return index associated with the key.

try:

return_index = {'z': 0, 'dzdx': 1, 'dzdy': 2}[return_key]

except KeyError as err:

raise ValueError("return_keys items shall take values in"

" {'z', 'dzdx', 'dzdy'}") from err

# Sets the scale factor for f & df components

scale = [1., 1./self._unit_x, 1./self._unit_y][return_index]

# Computes the interpolation

ret_loc = np.empty(size_ret, dtype=np.float64)

ret_loc[~mask_in] = np.nan

ret_loc[mask_in] = self._interpolate_single_key(

return_key, valid_tri_index, valid_x, valid_y) * scale

ret += [np.ma.masked_invalid(ret_loc.reshape(sh_ret), copy=False)]

return ret

def _interpolate_single_key(self, return_key, tri_index, x, y):

"""

raise NotImplementedError("TriInterpolator subclasses" +

"""

def __init__(self, triangulation, z, trifinder=None):

TriInterpolator.__init__(self, triangulation, z, trifinder)

# Store plane coefficients for fast interpolation calculations.

self._plane_coefficients = \

self._triangulation.calculate_plane_coefficients(self._z)

def __call__(self, x, y):

return self._interpolate_multikeys(x, y, tri_index=None,

return_keys=('z',))[0]

__call__.__doc__ = TriInterpolator._docstring__call__

def gradient(self, x, y):

return self._interpolate_multikeys(x, y, tri_index=None,

return_keys=('dzdx', 'dzdy'))

gradient.__doc__ = TriInterpolator._docstringgradient

def _interpolate_single_key(self, return_key, tri_index, x, y):

if return_key == 'z':

return (self._plane_coefficients[tri_index, 0]*x +

self._plane_coefficients[tri_index, 1]*y +

self._plane_coefficients[tri_index, 2])

elif return_key == 'dzdx':

return self._plane_coefficients[tri_index, 0]

elif return_key == 'dzdy':

return self._plane_coefficients[tri_index, 1]

else:

r"""

def __init__(self, triangulation, z, kind='min_E', trifinder=None,

dz=None):

TriInterpolator.__init__(self, triangulation, z, trifinder)

# Loads the underlying c++ _triangulation.

# (During loading, reordering of triangulation._triangles may occur so

# that all final triangles are now anti-clockwise)

self._triangulation.get_cpp_triangulation()

# To build the stiffness matrix and avoid zero-energy spurious modes

# we will only store internally the valid (unmasked) triangles and

# the necessary (used) points coordinates.

# 2 renumbering tables need to be computed and stored:

# - a triangle renum table in order to translate the result from a

# TriFinder instance into the internal stored triangle number.

# - a node renum table to overwrite the self._z values into the new

# (used) node numbering.

tri_analyzer = TriAnalyzer(self._triangulation)

(compressed_triangles, compressed_x, compressed_y, tri_renum,

node_renum) = tri_analyzer._get_compressed_triangulation()

self._triangles = compressed_triangles

self._tri_renum = tri_renum

# Taking into account the node renumbering in self._z:

node_mask = (node_renum == -1)

self._z[node_renum[~node_mask]] = self._z

self._z = self._z[~node_mask]

# Computing scale factors

self._unit_x = np.ptp(compressed_x)

self._unit_y = np.ptp(compressed_y)

self._pts = np.column_stack([compressed_x / self._unit_x,

compressed_y / self._unit_y])

# Computing triangle points

self._tris_pts = self._pts[self._triangles]

# Computing eccentricities

self._eccs = self._compute_tri_eccentricities(self._tris_pts)

# Computing dof estimations for HCT triangle shape function

self._dof = self._compute_dof(kind, dz=dz)

# Loading HCT element

self._ReferenceElement = _ReducedHCT_Element()

def __call__(self, x, y):

return self._interpolate_multikeys(x, y, tri_index=None,

return_keys=('z',))[0]

__call__.__doc__ = TriInterpolator._docstring__call__

def gradient(self, x, y):

return self._interpolate_multikeys(x, y, tri_index=None,

return_keys=('dzdx', 'dzdy'))

gradient.__doc__ = TriInterpolator._docstringgradient

def _interpolate_single_key(self, return_key, tri_index, x, y):

tris_pts = self._tris_pts[tri_index]

alpha = self._get_alpha_vec(x, y, tris_pts)

ecc = self._eccs[tri_index]

dof = np.expand_dims(self._dof[tri_index], axis=1)

if return_key == 'z':

return self._ReferenceElement.get_function_values(

alpha, ecc, dof)

elif return_key in ['dzdx', 'dzdy']:

J = self._get_jacobian(tris_pts)

dzdx = self._ReferenceElement.get_function_derivatives(

alpha, J, ecc, dof)

if return_key == 'dzdx':

return dzdx[:, 0, 0]

else:

return dzdx[:, 1, 0]

else:

raise ValueError("Invalid return_key: " + return_key)

def _compute_dof(self, kind, dz=None):

"""

if kind == 'user':

if dz is None:

raise ValueError("For a CubicTriInterpolator with "

"*kind*='user', a valid *dz* "

"argument is expected.")

TE = _DOF_estimator_user(self, dz=dz)

elif kind == 'geom':

TE = _DOF_estimator_geom(self)

elif kind == 'min_E':

TE = _DOF_estimator_min_E(self)

else:

cbook._check_in_list(['user', 'geom', 'min_E'], kind=kind)

return TE.compute_dof_from_df()

@staticmethod

def _get_alpha_vec(x, y, tris_pts):

"""

ndim = tris_pts.ndim-2

a = tris_pts[:, 1, :] - tris_pts[:, 0, :]

b = tris_pts[:, 2, :] - tris_pts[:, 0, :]

abT = np.stack([a, b], axis=-1)

ab = _transpose_vectorized(abT)

OM = np.stack([x, y], axis=1) - tris_pts[:, 0, :]

metric = _prod_vectorized(ab, abT)

# Here we try to deal with the colinear cases.

# metric_inv is in this case set to the Moore-Penrose pseudo-inverse

# meaning that we will still return a set of valid barycentric

# coordinates.

metric_inv = _pseudo_inv22sym_vectorized(metric)

Covar = _prod_vectorized(ab, _transpose_vectorized(

np.expand_dims(OM, ndim)))

ksi = _prod_vectorized(metric_inv, Covar)

alpha = _to_matrix_vectorized([

[1-ksi[:, 0, 0]-ksi[:, 1, 0]], [ksi[:, 0, 0]], [ksi[:, 1, 0]]])

return alpha

@staticmethod

def _get_jacobian(tris_pts):

"""

a = np.array(tris_pts[:, 1, :] - tris_pts[:, 0, :])

b = np.array(tris_pts[:, 2, :] - tris_pts[:, 0, :])

J = _to_matrix_vectorized([[a[:, 0], a[:, 1]],

[b[:, 0], b[:, 1]]])

return J

@staticmethod

def _compute_tri_eccentricities(tris_pts):

"""

a = np.expand_dims(tris_pts[:, 2, :] - tris_pts[:, 1, :], axis=2)

b = np.expand_dims(tris_pts[:, 0, :] - tris_pts[:, 2, :], axis=2)

c = np.expand_dims(tris_pts[:, 1, :] - tris_pts[:, 0, :], axis=2)

# Do not use np.squeeze, this is dangerous if only one triangle

# in the triangulation...

dot_a = _prod_vectorized(_transpose_vectorized(a), a)[:, 0, 0]

dot_b = _prod_vectorized(_transpose_vectorized(b), b)[:, 0, 0]

dot_c = _prod_vectorized(_transpose_vectorized(c), c)[:, 0, 0]

# Note that this line will raise a warning for dot_a, dot_b or dot_c

# zeros, but we choose not to support triangles with duplicate points.

return _to_matrix_vectorized([[(dot_c-dot_b) / dot_a],

[(dot_a-dot_c) / dot_b],

"""

# 1) Loads matrices to generate shape functions as a function of

# triangle eccentricities - based on [1] p.11 '''

M = np.array([

[ 0.00, 0.00, 0.00, 4.50, 4.50, 0.00, 0.00, 0.00, 0.00, 0.00],

[-0.25, 0.00, 0.00, 0.50, 1.25, 0.00, 0.00, 0.00, 0.00, 0.00],

[-0.25, 0.00, 0.00, 1.25, 0.50, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.50, 1.00, 0.00, -1.50, 0.00, 3.00, 3.00, 0.00, 0.00, 3.00],

[ 0.00, 0.00, 0.00, -0.25, 0.25, 0.00, 1.00, 0.00, 0.00, 0.50],

[ 0.25, 0.00, 0.00, -0.50, -0.25, 1.00, 0.00, 0.00, 0.00, 1.00],

[ 0.50, 0.00, 1.00, 0.00, -1.50, 0.00, 0.00, 3.00, 3.00, 3.00],

[ 0.25, 0.00, 0.00, -0.25, -0.50, 0.00, 0.00, 0.00, 1.00, 1.00],

[ 0.00, 0.00, 0.00, 0.25, -0.25, 0.00, 0.00, 1.00, 0.00, 0.50]])

M0 = np.array([

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[-1.00, 0.00, 0.00, 1.50, 1.50, 0.00, 0.00, 0.00, 0.00, -3.00],

[-0.50, 0.00, 0.00, 0.75, 0.75, 0.00, 0.00, 0.00, 0.00, -1.50],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 1.00, 0.00, 0.00, -1.50, -1.50, 0.00, 0.00, 0.00, 0.00, 3.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.50, 0.00, 0.00, -0.75, -0.75, 0.00, 0.00, 0.00, 0.00, 1.50]])

M1 = np.array([

[-0.50, 0.00, 0.00, 1.50, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[-0.25, 0.00, 0.00, 0.75, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.50, 0.00, 0.00, -1.50, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.25, 0.00, 0.00, -0.75, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]])

M2 = np.array([

[ 0.50, 0.00, 0.00, 0.00, -1.50, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.25, 0.00, 0.00, 0.00, -0.75, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[-0.50, 0.00, 0.00, 0.00, 1.50, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[-0.25, 0.00, 0.00, 0.00, 0.75, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],

[ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]])

# 2) Loads matrices to rotate components of gradient & Hessian

# vectors in the reference basis of triangle first apex (a0)

rotate_dV = np.array([[ 1., 0.], [ 0., 1.],

[ 0., 1.], [-1., -1.],

[-1., -1.], [ 1., 0.]])

rotate_d2V = np.array([[1., 0., 0.], [0., 1., 0.], [ 0., 0., 1.],

[0., 1., 0.], [1., 1., 1.], [ 0., -2., -1.],

[1., 1., 1.], [1., 0., 0.], [-2., 0., -1.]])

# 3) Loads Gauss points & weights on the 3 sub-_triangles for P2

# exact integral - 3 points on each subtriangles.

# NOTE: as the 2nd derivative is discontinuous , we really need those 9

# points!

n_gauss = 9

gauss_pts = np.array([[13./18., 4./18., 1./18.],

[ 4./18., 13./18., 1./18.],

[ 7./18., 7./18., 4./18.],

[ 1./18., 13./18., 4./18.],

[ 1./18., 4./18., 13./18.],

[ 4./18., 7./18., 7./18.],

[ 4./18., 1./18., 13./18.],

[13./18., 1./18., 4./18.],

[ 7./18., 4./18., 7./18.]], dtype=np.float64)

gauss_w = np.ones([9], dtype=np.float64) / 9.

# 4) Stiffness matrix for curvature energy

E = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 2.]])

# 5) Loads the matrix to compute DOF_rot from tri_J at apex 0

J0_to_J1 = np.array([[-1., 1.], [-1., 0.]])

J0_to_J2 = np.array([[ 0., -1.], [ 1., -1.]])

def get_function_values(self, alpha, ecc, dofs):

"""

subtri = np.argmin(alpha, axis=1)[:, 0]

ksi = _roll_vectorized(alpha, -subtri, axis=0)

E = _roll_vectorized(ecc, -subtri, axis=0)

x = ksi[:, 0, 0]

y = ksi[:, 1, 0]

z = ksi[:, 2, 0]

x_sq = x*x

y_sq = y*y

z_sq = z*z

V = _to_matrix_vectorized([

[x_sq*x], [y_sq*y], [z_sq*z], [x_sq*z], [x_sq*y], [y_sq*x],

[y_sq*z], [z_sq*y], [z_sq*x], [x*y*z]])

prod = _prod_vectorized(self.M, V)

prod += _scalar_vectorized(E[:, 0, 0],

_prod_vectorized(self.M0, V))

prod += _scalar_vectorized(E[:, 1, 0],

_prod_vectorized(self.M1, V))

prod += _scalar_vectorized(E[:, 2, 0],

_prod_vectorized(self.M2, V))

s = _roll_vectorized(prod, 3*subtri, axis=0)

return _prod_vectorized(dofs, s)[:, 0, 0]

def get_function_derivatives(self, alpha, J, ecc, dofs):

"""

subtri = np.argmin(alpha, axis=1)[:, 0]

ksi = _roll_vectorized(alpha, -subtri, axis=0)

E = _roll_vectorized(ecc, -subtri, axis=0)

x = ksi[:, 0, 0]

y = ksi[:, 1, 0]

z = ksi[:, 2, 0]

x_sq = x*x

y_sq = y*y

z_sq = z*z

dV = _to_matrix_vectorized([

[ -3.*x_sq, -3.*x_sq],

[ 3.*y_sq, 0.],

[ 0., 3.*z_sq],

[ -2.*x*z, -2.*x*z+x_sq],

[-2.*x*y+x_sq, -2.*x*y],

[ 2.*x*y-y_sq, -y_sq],

[ 2.*y*z, y_sq],

[ z_sq, 2.*y*z],

[ -z_sq, 2.*x*z-z_sq],

[ x*z-y*z, x*y-y*z]])

# Puts back dV in first apex basis

dV = _prod_vectorized(dV, _extract_submatrices(

self.rotate_dV, subtri, block_size=2, axis=0))

prod = _prod_vectorized(self.M, dV)

prod += _scalar_vectorized(E[:, 0, 0],

_prod_vectorized(self.M0, dV))

prod += _scalar_vectorized(E[:, 1, 0],

_prod_vectorized(self.M1, dV))

prod += _scalar_vectorized(E[:, 2, 0],

_prod_vectorized(self.M2, dV))

dsdksi = _roll_vectorized(prod, 3*subtri, axis=0)

dfdksi = _prod_vectorized(dofs, dsdksi)

# In global coordinates:

# Here we try to deal with the simplest colinear cases, returning a

# null matrix.

J_inv = _safe_inv22_vectorized(J)

dfdx = _prod_vectorized(J_inv, _transpose_vectorized(dfdksi))

return dfdx

def get_function_hessians(self, alpha, J, ecc, dofs):

"""

d2sdksi2 = self.get_d2Sidksij2(alpha, ecc)

d2fdksi2 = _prod_vectorized(dofs, d2sdksi2)

H_rot = self.get_Hrot_from_J(J)

d2fdx2 = _prod_vectorized(d2fdksi2, H_rot)

return _transpose_vectorized(d2fdx2)

def get_d2Sidksij2(self, alpha, ecc):

"""

subtri = np.argmin(alpha, axis=1)[:, 0]

ksi = _roll_vectorized(alpha, -subtri, axis=0)

E = _roll_vectorized(ecc, -subtri, axis=0)

x = ksi[:, 0, 0]

y = ksi[:, 1, 0]

z = ksi[:, 2, 0]

d2V = _to_matrix_vectorized([

[ 6.*x, 6.*x, 6.*x],

[ 6.*y, 0., 0.],

[ 0., 6.*z, 0.],

[ 2.*z, 2.*z-4.*x, 2.*z-2.*x],

[2.*y-4.*x, 2.*y, 2.*y-2.*x],

[2.*x-4.*y, 0., -2.*y],

[ 2.*z, 0., 2.*y],

[ 0., 2.*y, 2.*z],

[ 0., 2.*x-4.*z, -2.*z],

[ -2.*z, -2.*y, x-y-z]])

# Puts back d2V in first apex basis

d2V = _prod_vectorized(d2V, _extract_submatrices(

self.rotate_d2V, subtri, block_size=3, axis=0))

prod = _prod_vectorized(self.M, d2V)

prod += _scalar_vectorized(E[:, 0, 0],

_prod_vectorized(self.M0, d2V))

prod += _scalar_vectorized(E[:, 1, 0],

_prod_vectorized(self.M1, d2V))

prod += _scalar_vectorized(E[:, 2, 0],

_prod_vectorized(self.M2, d2V))

d2sdksi2 = _roll_vectorized(prod, 3*subtri, axis=0)

return d2sdksi2

def get_bending_matrices(self, J, ecc):

"""

n = np.size(ecc, 0)

# 1) matrix to rotate dofs in global coordinates

J1 = _prod_vectorized(self.J0_to_J1, J)

J2 = _prod_vectorized(self.J0_to_J2, J)

DOF_rot = np.zeros([n, 9, 9], dtype=np.float64)

DOF_rot[:, 0, 0] = 1

DOF_rot[:, 3, 3] = 1

DOF_rot[:, 6, 6] = 1

DOF_rot[:, 1:3, 1:3] = J

DOF_rot[:, 4:6, 4:6] = J1

DOF_rot[:, 7:9, 7:9] = J2

# 2) matrix to rotate Hessian in global coordinates.

H_rot, area = self.get_Hrot_from_J(J, return_area=True)

# 3) Computes stiffness matrix

# Gauss quadrature.

K = np.zeros([n, 9, 9], dtype=np.float64)

weights = self.gauss_w

pts = self.gauss_pts

for igauss in range(self.n_gauss):

alpha = np.tile(pts[igauss, :], n).reshape(n, 3)

alpha = np.expand_dims(alpha, 2)

weight = weights[igauss]

d2Skdksi2 = self.get_d2Sidksij2(alpha, ecc)

d2Skdx2 = _prod_vectorized(d2Skdksi2, H_rot)

K += weight * _prod_vectorized(_prod_vectorized(d2Skdx2, self.E),

_transpose_vectorized(d2Skdx2))

# 4) With nodal (not elem) dofs

K = _prod_vectorized(_prod_vectorized(_transpose_vectorized(DOF_rot),

K), DOF_rot)

# 5) Need the area to compute total element energy

return _scalar_vectorized(area, K)

def get_Hrot_from_J(self, J, return_area=False):

"""

# Here we try to deal with the simplest colinear cases; a null

# energy and area is imposed.

J_inv = _safe_inv22_vectorized(J)

Ji00 = J_inv[:, 0, 0]

Ji11 = J_inv[:, 1, 1]

Ji10 = J_inv[:, 1, 0]

Ji01 = J_inv[:, 0, 1]

H_rot = _to_matrix_vectorized([

[Ji00*Ji00, Ji10*Ji10, Ji00*Ji10],

[Ji01*Ji01, Ji11*Ji11, Ji01*Ji11],

[2*Ji00*Ji01, 2*Ji11*Ji10, Ji00*Ji11+Ji10*Ji01]])

if not return_area:

return H_rot

else:

area = 0.5 * (J[:, 0, 0]*J[:, 1, 1] - J[:, 0, 1]*J[:, 1, 0])

return H_rot, area

def get_Kff_and_Ff(self, J, ecc, triangles, Uc):

"""

ntri = np.size(ecc, 0)

vec_range = np.arange(ntri, dtype=np.int32)

c_indices = np.full(ntri, -1, dtype=np.int32) # for unused dofs, -1

f_dof = [1, 2, 4, 5, 7, 8]

c_dof = [0, 3, 6]

# vals, rows and cols indices in global dof numbering

f_dof_indices = _to_matrix_vectorized([[

c_indices, triangles[:, 0]*2, triangles[:, 0]*2+1,

c_indices, triangles[:, 1]*2, triangles[:, 1]*2+1,

c_indices, triangles[:, 2]*2, triangles[:, 2]*2+1]])

expand_indices = np.ones([ntri, 9, 1], dtype=np.int32)

f_row_indices = _prod_vectorized(_transpose_vectorized(f_dof_indices),

_transpose_vectorized(expand_indices))

f_col_indices = _prod_vectorized(expand_indices, f_dof_indices)

K_elem = self.get_bending_matrices(J, ecc)

# Extracting sub-matrices

# Explanation & notations:

# * Subscript f denotes 'free' degrees of freedom (i.e. dz/dx, dz/dx)

# * Subscript c denotes 'condensated' (imposed) degrees of freedom

# (i.e. z at all nodes)

# * F = [Ff, Fc] is the force vector

# * U = [Uf, Uc] is the imposed dof vector

# [ Kff Kfc ]

# * K = [ ] is the laplacian stiffness matrix

# [ Kcf Kff ]

# * As F = K x U one gets straightforwardly: Ff = - Kfc x Uc

# Computing Kff stiffness matrix in sparse coo format

Kff_vals = np.ravel(K_elem[np.ix_(vec_range, f_dof, f_dof)])

Kff_rows = np.ravel(f_row_indices[np.ix_(vec_range, f_dof, f_dof)])

Kff_cols = np.ravel(f_col_indices[np.ix_(vec_range, f_dof, f_dof)])

# Computing Ff force vector in sparse coo format

Kfc_elem = K_elem[np.ix_(vec_range, f_dof, c_dof)]

Uc_elem = np.expand_dims(Uc, axis=2)

Ff_elem = - _prod_vectorized(Kfc_elem, Uc_elem)[:, :, 0]

Ff_indices = f_dof_indices[np.ix_(vec_range, [0], f_dof)][:, 0, :]

# Extracting Ff force vector in dense format

# We have to sum duplicate indices - using bincount

Ff = np.bincount(np.ravel(Ff_indices), weights=np.ravel(Ff_elem))