r""" Module for functionals defined on IndexedFaceSets.
The general form of a functional f is f(ifs, \*attributes), where
ifs is the IndexedFaceSet and attributes are the names of the attributes to
be used.
It is assumed that a functional 'knows' on what kind of cell it is defined.
Note that we also call vector-valued maps functionals.
"""
import numpy as np
from scipy.sparse import csr_matrix, diags
from . import utils as ut
# Flatness Energies
# =================
[docs]def flat_cubic_quad(ifs, attr_name):
"""Flatness energy of 4-gonal indexed face set calculated on each
face as the distance of its diagonals.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
n = np.cross(d1, d2)
norms = np.linalg.norm(n, axis=1)[:, None]
norms[np.where(norms == 0)] = 1
n = np.true_divide(n, norms)
k = np.einsum("ij, ij->i", n, v[f[:, 1]] - v[f[:, 0]])
return k
[docs]def g_flat_cubic_quad(ifs, attr_name):
"""Jacobian for flat_cubic_quad.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
scipy.sparse.csr_matrix
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
jac = np.zeros((len(f), 3 * len(v)))
n = np.cross(d1, d2)
norms = np.linalg.norm(n, axis=1)[:, None]
norms[np.where(norms == 0)] = 1
proj = v[f[:, 1]] - v[f[:, 0]]
dist = np.einsum("ij, ij->i", n, proj)[:, None] / norms
jac_v = np.cross(v[f[:, 1]] - v[f[:, 2]] - dist / norms * n, d2) / norms
jac_v1 = np.cross(v[f[:, 2]] - v[f[:, 3]] + dist / norms * n, d1) / norms
jac_v2 = np.cross(v[f[:, 0]] - v[f[:, 1]] + dist / norms * n, d2) / norms
jac_v3 = np.cross(v[f[:, 1]] - v[f[:, 0]] - dist / norms * n, d1) / norms
data = np.concatenate([jac_v, jac_v1, jac_v2, jac_v3], axis=1).ravel()
idx = ut.nflat_index(f, 3)
indptr = np.array([k * 12 for k in range(len(f) + 1)])
return csr_matrix((data, idx, indptr))
[docs]def flat_cubic_n_quad(ifs, attr_name):
"""Flatness energy of 4-gonal indexed face set calculated on each
face as the distance of its diagonals divided by the average lengths
of the diagonals.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
l = 0.5 * (np.linalg.norm(d1, axis=1) + np.linalg.norm(d2, axis=1))
n = np.cross(d1, d2)
n = np.true_divide(n, np.linalg.norm(n, axis=1)[:, None])
k = np.einsum("ij, ij->i", n, v[f[:, 1]] - v[f[:, 0]]) / l
return k
[docs]def g_flat_cubic_n_quad(ifs, attr_name):
"""Jacobian for flat_cubic_n_quad.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
scipy.sparse.csr_matrix
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
l1 = np.linalg.norm(d1, axis=-1)
l2 = np.linalg.norm(d2, axis=-1)
dist = flat_cubic_quad(ifs, attr_name)
g_dist = g_flat_cubic_quad(ifs, attr_name)
jac_v = -d1 / l1[:, None]
jac_v1 = -d2 / l2[:, None]
jac_v2 = d1 / l1[:, None]
jac_v3 = d2 / l2[:, None]
data = np.concatenate([jac_v, jac_v1, jac_v2, jac_v3], axis=1).ravel()
idx = ut.nflat_index(f, 3)
indptr = np.array([k * 12 for k in range(len(f) + 1)])
jac = csr_matrix((data, idx, indptr))
total_jac = diags(2 / (l1 + l2)) @ g_dist - diags(2 * dist / (l1 + l2) ** 2) @ jac
return total_jac
[docs]def flat_cubic_n2_quad(ifs, attr_name):
"""Flatness energy of 4-gonal indexed face set calculated on each
face as the distance of its diagonals d1, d2 mutliplied by (1/d1**2+1/d2**2).
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
l1 = np.linalg.norm(d1, axis=1)
l2 = np.linalg.norm(d2, axis=1)
fac = 1 / l1**2 + 1 / l2**2
n = np.cross(d1, d2)
n = np.true_divide(n, np.linalg.norm(n, axis=1)[:, None])
k = fac * np.einsum("ij, ij->i", n, v[f[:, 1]] - v[f[:, 0]])
return k
[docs]def g_flat_cubic_n2_quad(ifs, attr_name):
"""Jacobian for flat_cubic_n2_quad.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
scipy.sparse.csr_matrix
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
d1 = v[f[:, 2]] - v[f[:, 0]]
d2 = v[f[:, 3]] - v[f[:, 1]]
l1 = np.linalg.norm(d1, axis=-1)
l2 = np.linalg.norm(d2, axis=-1)
fac = 1 / l1**2 + 1 / l2**2
dist = flat_cubic_quad(ifs, attr_name)
g_dist = g_flat_cubic_quad(ifs, attr_name)
jac_v = -d1 / l1[:, None] ** 4
jac_v1 = -d2 / l2[:, None] ** 4
jac_v2 = d1 / l1[:, None] ** 4
jac_v3 = d2 / l2[:, None] ** 4
data = np.concatenate([jac_v, jac_v1, jac_v2, jac_v3], axis=1).ravel()
idx = ut.nflat_index(f, 3)
indptr = np.array([k * 12 for k in range(len(f) + 1)])
jac = csr_matrix((data, idx, indptr))
total_jac = diags(fac) @ g_dist - diags(2 * dist) @ jac
return total_jac
[docs]def flat_det_quad(ifs, attr_name):
"""Flatness energy of 4-gonal indexed face set calculated on each
face as the determinant of the matrix given by its horizontal, vertical and
the corresponding diagonal vector.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
diags = v[f[:, 2]] - v[f[:, 0]]
ver = v[f[:, 3]] - v[f[:, 0]]
hor = v[f[:, 1]] - v[f[:, 0]]
dets = np.linalg.det(np.stack((hor, diags, ver), axis=-1))
return dets
[docs]def g_flat_det_quad(ifs, attr_name):
"""Jacobian for flat_det_quad.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
scipy.sparse.csr_matrix
"""
a = ifs.vertex_attributes
v = a[attr_name]
f = np.int32(ifs.faces)
diags = v[f[:, 2]] - v[f[:, 0]]
ver = v[f[:, 3]] - v[f[:, 0]]
hor = v[f[:, 1]] - v[f[:, 0]]
jac_v = np.cross(hor, diags) + np.cross(diags, ver) + np.cross(ver, hor)
jac_v2 = -np.cross(hor, ver)
jac_v1 = np.cross(diags, ver)
jac_v3 = np.cross(hor, diags)
data = np.concatenate([-jac_v, jac_v1, jac_v2, jac_v3], axis=1).ravel()
idx = ut.nflat_index(f, 3)
indptr = np.array([k * 12 for k in range(len(f) + 1)])
return csr_matrix((data, idx, indptr))
# Circularity Energies
# ====================
[docs]def circ_angle_quad(ifs, attr_name):
"""Circularity energy of 4-gonal indexed face set calculated by the
angles of its faces.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
f = np.int32(ifs.faces)
attribute_dict = ifs.vertex_attributes
attr_values = attribute_dict[attr_name]
q1 = attr_values[f[:, 1]] - attr_values[f[:, 0]]
q2 = attr_values[f[:, 2]] - attr_values[f[:, 1]]
q3 = attr_values[f[:, 3]] - attr_values[f[:, 2]]
q4 = attr_values[f[:, 0]] - attr_values[f[:, 3]]
l1 = q1 / np.linalg.norm(q1, axis=1)[:, None]
l2 = q2 / np.linalg.norm(q2, axis=1)[:, None]
l3 = q3 / np.linalg.norm(q3, axis=1)[:, None]
l4 = q4 / np.linalg.norm(q4, axis=1)[:, None]
cosa = np.arccos(np.einsum("ij, ij->i", l1, -l4))
cosb = np.arccos(np.einsum("ij, ij->i", l2, -l1))
cosc = np.arccos(np.einsum("ij, ij->i", l3, -l2))
cosd = np.arccos(np.einsum("ij, ij->i", l4, -l3))
return np.concatenate([cosa + cosc - np.pi, cosb + cosd - np.pi])
[docs]def g_circ_angle_quad(ifs, attr_name):
"""Jacobian for flat_det_quad.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 3-dim. vertex attribute to use
Returns
-------
scipy.sparse.csr_matrix
"""
f = np.int32(ifs.faces)
attribute_dict = ifs.vertex_attributes
attr_values = attribute_dict[attr_name]
q1 = attr_values[f[:, 1]] - attr_values[f[:, 0]] # v1 - attr_values
q2 = attr_values[f[:, 2]] - attr_values[f[:, 1]] # v12 - v1
q3 = attr_values[f[:, 3]] - attr_values[f[:, 2]] # v2 - v12
q4 = attr_values[f[:, 0]] - attr_values[f[:, 3]] # attr_values - v2
l1 = np.linalg.norm(q1, axis=1)[:, None]
l2 = np.linalg.norm(q2, axis=1)[:, None]
l3 = np.linalg.norm(q3, axis=1)[:, None]
l4 = np.linalg.norm(q4, axis=1)[:, None]
cosa = np.einsum("ij, ij->i", q1, -q4)[:, None] / l1 / l4
cosb = np.einsum("ij, ij->i", q2, -q1)[:, None] / l2 / l1
cosc = np.einsum("ij, ij->i", q3, -q2)[:, None] / l3 / l2
cosd = np.einsum("ij, ij->i", q4, -q3)[:, None] / l4 / l3
a_v = -1 / np.sqrt(1 - cosa**2)
b_v = -1 / np.sqrt(1 - cosb**2)
c_v = -1 / np.sqrt(1 - cosc**2)
d_v = -1 / np.sqrt(1 - cosd**2)
g_ac_v = a_v * ((q4 - q1) / l1 / l4 + cosa * (q1 / l1**2 - q4 / l4**2))
g_ac_v1 = a_v * (-q4 / l1 / l4 - cosa * q1 / l1**2) + c_v * (
q3 / l2 / l3 + cosc * q2 / l2**2
)
g_ac_v2 = c_v * ((q2 - q3) / l2 / l3 + cosc * (q3 / l3**2 - q2 / l2**2))
g_ac_v3 = a_v * (q1 / l1 / l4 + cosa * q4 / l4**2) - c_v * (
q2 / l2 / l3 + cosc * q3 / l3**2
)
g_bd_v = b_v * (q2 / l1 / l2 + cosb * q1 / l1**2) - d_v * (
q3 / l3 / l4 + cosd * q4 / l4**2
)
g_bd_v1 = b_v * ((q1 - q2) / l1 / l2 - cosb * (q1 / l1**2 - q2 / l2**2))
g_bd_v2 = b_v * (-q1 / l1 / l2 - cosb * q2 / l2**2) + d_v * (
q4 / l3 / l4 + cosd * q3 / l3**2
)
g_bd_v3 = d_v * ((q3 - q4) / l3 / l4 + cosd * (q4 / l4**2 - q3 / l3**2))
data1 = np.concatenate([g_ac_v, g_ac_v1, g_ac_v2, g_ac_v3], axis=1).ravel()
data2 = np.concatenate([g_bd_v, g_bd_v1, g_bd_v2, g_bd_v3], axis=1).ravel()
data = np.concatenate([data1, data2])
idx_ac = ut.nflat_index(f, 3)
idx = np.concatenate([idx_ac, idx_ac])
indptr = np.array([k * 12 for k in range(2 * len(f) + 1)])
return csr_matrix((data, idx, indptr))
# Dirichlet Energy
# ================
[docs]def dirichlet(ifs, attr_name):
"""Dirichlet energy for indexed face sets.
Parameters
----------
ifs : ddg.datastructures.indexedfaceset.ifs.IndexedFaceSet
4-gonal indexed face set to calculate the energy on
attr_name : str
Name of the 1-dim. vertex attribute to use
Returns
-------
numpy.ndarray
"""
e = np.int32(ifs.edges)
attribute_dict = ifs.vertex_attributes
attr_values = attribute_dict[attr_name]
return attr_values[e[:, 0]] - attr_values[e[:, 1]]