# -*- coding: utf-8 -*-UnstructuredFeaSym
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
import numpy as np
from typing import List, Tuple
from Muscat.Types import MuscatFloat, MuscatIndex
from numpy.typing import NDArray
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Filters.FilterBase import Intersect1D
from Muscat.MeshContainers.Mesh import Mesh
import Muscat.MeshContainers.ElementsDescription as ED
import Muscat.FE.SymWeakForm as wf
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Fields.FEField import FEField
from Muscat.MeshTools.MeshInspectionTools import GetVolumePerElement
from Muscat.MeshTools.InriaMeshTools import flattenMetric, fullMetric
[docs]
def ComputeGradient(mesh: Mesh, psi: NDArray, tags: List[str] = None) -> NDArray:
"""Compute the gradient at each node from a scalar nodefield using a variational formulation:
Find Grad such that int Grad . GradT dx = int grad(Psi) . GradT dx for all GradT
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
tags : List[str]
tags of the computational domain
Returns
-------
grad : NDarray
the gradient of psi at each node
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
grad = np.zeros([siz, dim], dtype = MuscatFloat)
Psi = wf.GetField("Psi", 1)
Grad = wf.GetField("Grad", dim)
GradT = wf.GetTestField("Grad", dim)
Ener = Grad.T*GradT + wf.Gradient(Psi, sdim = dim).T*GradT
numbering = ComputeDofNumbering(mesh, LagrangeSpaceP1, fromConnectivity = True)
psi = np.ascontiguousarray(psi)
field = FEField(name = "Psi",
mesh = mesh,
space = LagrangeSpaceP1,
numbering = numbering,
data = psi
)
dofs = ["Grad_0", "Grad_1", "Grad_2"]
unknownFields = [FEField(name = name,
mesh = mesh,
space = LagrangeSpaceP1,
numbering = numbering
) for name in dofs[:dim] ]
if tags == None:
ff = ElementFilter(dimensionality = mesh.GetElementsDimensionality())
else:
ff = ElementFilter(dimensionality = mesh.GetElementsDimensionality(), eTag = tags)
m, f = IntegrateGeneral(mesh = mesh,
wform = Ener,
constants = {},
fields = [field],
unknownFields = unknownFields,
integrationRuleName = "NodalEvaluationIsoParam",
elementFilter = ff)
d = m.diagonal()
g = f/(d + (d == 0))
for i in range(dim):
grad[:, i] = g[i*siz:(i+1)*siz]
return grad
[docs]
def ComputeGradient_L2projection(mesh: Mesh, psi: NDArray) -> NDArray:
"""Compute the gradient at each node from a scalar nodefield psi using formula:
grad_i(psi) = sum_{K in S_i} \|K\| grad_K(psi) / sum_{K in S_i} \|K\|
Warning : works only with simplex meshes
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
tags : List[str]
tags of the computational domain
Returns
-------
grad : NDarray
the gradient of psi at each node
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
grad = np.zeros([siz, dim], dtype = MuscatFloat)
eltDim = mesh.GetElementsDimensionality()
if eltDim == 3:
eltype = ED.Tetrahedron_4
elif eltDim == 2:
eltype = ED.Triangle_3
else:
eltype = ED.Bar
cell = mesh.GetElementsOfType(eltype)
conn = cell.connectivity
coords = mesh.nodes[conn]
num = cell.GetNumberOfElements()
mat = np.zeros([num, dim, dim], dtype = MuscatFloat)
for i in range(dim):
mat[:, i, :] = coords[:, i+1, :]-coords[:, 0, :]
invmat = np.linalg.inv(mat)
data = psi[conn]
solmat = np.zeros([num, dim], dtype = MuscatFloat)
for i in range(dim):
solmat[:, i] = data[:, i+1] - data[:, 0]
res = np.squeeze(np.matmul(invmat, solmat[:, :, np.newaxis]))
weight = np.zeros(siz, dtype = MuscatFloat)
weightAtElement = GetVolumePerElement(mesh, ElementFilter(elementType = eltype))
for i in range(conn.shape[1]):
weight[conn[:, i]] += weightAtElement
grad[conn[:, i]] += res * weightAtElement[:, np.newaxis]
ids = np.where(weight > 0)[0]
grad[ids, :] = grad[ids, :]/weight[ids, np.newaxis]
return grad
[docs]
def ComputeHessian(mesh: Mesh, psi: NDArray, gradL2: bool = False, tags: List[str] = None) -> NDArray:
"""Compute the hessian at each node from a scalar nodefield phi using twice the ComputeGradient function
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
tags : List[str]
tags of the computational domain
gradL2 : bool
if True, use a L2 projection to compute the gradients
Returns
-------
hess : NDarray
the symmetric matrix hessian of psi at each node
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
hess = np.zeros([siz, dim, dim], dtype = MuscatFloat)
if gradL2:
grad = ComputeGradient_L2projection(mesh = mesh, psi = psi)
for i in range(dim):
h = ComputeGradient_L2projection(mesh = mesh, psi = grad[:,i])
hess[:, i, :] += h
hess[:, :, i] += h
else:
grad = ComputeGradient(mesh = mesh, psi = psi, tags = tags)
for i in range(dim):
h = ComputeGradient(mesh = mesh, psi = grad[:,i], tags = tags)
hess[:, i, :] += h
hess[:, :, i] += h
return hess/2
[docs]
def ComputeMetric(mesh: Mesh, psi: NDArray, p: MuscatIndex = None, err: MuscatFloat = 0.01, N: MuscatIndex = None, gradL2: bool = False, zones: List[List[str]] = None) -> NDArray:
"""Construct a metric field for mesh adaptation from a given scalar nodefield phi in Lp norm
Warning : currently works only with simplex meshes
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
p: MuscatIndex
order of the Lp-norm (infinite by default)
err: MuscatFloat
desired level for the interpolation error
N: MuscatIndex
desired number of nodes in the adapted mesh
zones : List[List[str]]
list of the various computational domains (to allow a remeshing by part in the presence of interfaces)
Returns
-------
metric : NDarray
the metric at each node
a metric is a vector containing the (dim+1)*dim/2 coefficients of the associated symmetric definite positive matrix
warning : order of metric component imposed by mmg3d, different from order read by paraview
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
hmax = np.linalg.norm(np.max(mesh.nodes, axis=0) - np.min(mesh.nodes, axis=0))
if zones is None:
hess = ComputeHessian(mesh=mesh, psi=psi, gradL2=gradL2)
else:
hess = 0
nodes_idx = {}
for tags in zones:
h = ComputeHessian(mesh=mesh, psi=psi, gradL2=gradL2, tags=tags)
node_idx = np.array([i if val.any() != 0 else -1 for i, val in enumerate(h)]) # A ameliorer avec filtres
idx = np.delete(node_idx, np.where(node_idx == -1)) # A ameliorer avec filtres
n_intersect = []
for value in nodes_idx.values():
n_intersect.extend(Intersect1D(idx, value))
nodes_idx[tags[0]] = idx
if n_intersect:
eigenValue, eigenVector = np.linalg.eigh(hess[n_intersect])
eigenValue = abs(eigenValue)
met1 = BuildMet(eigenValue, eigenVector)
eigenValue, eigenVector = np.linalg.eigh(h[n_intersect])
eigenValue = abs(eigenValue)
met2 = BuildMet(eigenValue, eigenVector)
int_met = IntersectionOfMetrics(metrics=[flattenMetric(met1), flattenMetric(met2)], preprocessing=True) # Autoriser l'entrée au format plein ou flat
int_met = fullMetric(int_met) # Autoriser la sortie au format plein ou flat
hess = hess * (h == 0) + (hess == 0) * h # Utiliser node_idx
if n_intersect:
hess[n_intersect] = int_met
eigenValue, eigenVector = np.linalg.eigh(hess)
eigenValue = abs(eigenValue)
if p is None:
norm = np.linalg.norm(psi, ord=None)
else:
norm = np.linalg.norm(psi, ord=p)
if N is None:
coef = dim/(err * norm)
else:
coef = 1.
coef = ComputeCoef(mesh, eigenValue, hmax, p, N, coef)
eigenValue *= coef
eigenValue = np.maximum(eigenValue, 1./hmax**2)
# print(eigenValue)
diag = np.zeros([siz, dim, dim])
for i in range(dim):
diag[:, i, i] = eigenValue[:, i]
mat_metric = eigenVector @ diag @ np.linalg.inv(eigenVector)
metric = flattenMetric(mat_metric)
return metric
[docs]
def BuildMet(eigenValue: NDArray, eigenVector: NDArray) -> NDArray:
"""Build metric from eigenValue and eigenVector
Args:
eigenValue (NDArray): computed eigenValues for metric on nodes
eigenVector (NDArray): eigenVectors for metric on nodes
Returns:
NDArray: metric n*n
"""
siz = eigenValue.shape[0]
dim = eigenValue.shape[1]
diag = np.zeros([siz, dim, dim])
for i in range(dim):
diag[:, i, i] = eigenValue[:, i]
mat_metric = eigenVector @ diag @ np.linalg.inv(eigenVector)
return mat_metric
[docs]
def ComputeCoef(mesh: Mesh, eigenValue: NDArray, hmax: MuscatFloat, p: MuscatIndex, N: MuscatIndex, coef: MuscatFloat = 1.) -> MuscatFloat:
"""Coefficient computation to adapt eigenValues on constraints
Returns:
MuscatFloat: Coefficient multiplicator to apply to eigenValues
"""
dim = mesh.GetPointsDimensionality()
eltDim = mesh.GetElementsDimensionality()
if eltDim == 3:
eltype = ED.Tetrahedron_4
elif eltDim == 2:
eltype = ED.Triangle_3
else:
eltype = ED.Bar_2
cell = mesh.GetElementsOfType(eltype)
conn = cell.connectivity
weightAtElement = GetVolumePerElement(mesh, ElementFilter(elementType=eltype))/(dim+1)
det = np.prod(eigenValue, axis=1)[:, np.newaxis]
if max(det) < 1e-15: # cas pathologique tous les det sont nuls
eigenValue = np.clip(eigenValue, 1./(hmax**2), None)
det = np.prod(eigenValue, axis=1)[:, np.newaxis]
det = np.clip(det, 1e-15, None)
if p is None:
if N is not None:
integralOfDet = np.sum(weightAtElement*np.ravel(np.sum(np.power(det[conn], 1./2), axis=1)))
else:
integralOfDet = np.sum(weightAtElement*np.ravel(np.sum(np.power(det[conn], p/(2*p+dim)), axis=1)))
if N is None:
if p is not None:
coef *= np.power(integralOfDet, 1./p)
else:
coef = np.power(N, 2./dim)*np.power(integralOfDet, -2./dim)
if p is not None:
coef *= np.power(det, -1./(2*p+dim))
return coef
[docs]
def ComputeSizeOfCell(mesh: Mesh, psi: NDArray, err: MuscatFloat = 0.01, zones: List[List[str]] = None) -> NDArray:
"""Construct a map of size of mesh cells from a given scalar nodefield phi in L^inf norm
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
err: MuscatFloat
desired level for the interpolation error
zones : List[List[str]]
list of the various computational domains (to allow a remeshing by part in the presence of interfaces)
Returns
-------
sizeOfCell : NDarray
the resulting size map at each node
"""
dim = mesh.GetPointsDimensionality()
hmax = np.linalg.norm(np.max(mesh.nodes, axis = 0) - np.min(mesh.nodes, axis = 0))
norm = np.linalg.norm(psi, ord = None)
if zones == None:
hess = ComputeHessian(mesh = mesh, psi = psi)
else:
hess = 0
for tags in zones:
h = ComputeHessian(mesh = mesh, psi = psi, tags = tags)
hess = hess*(h == 0) + (hess == 0)*h
eigenValue = np.linalg.eigvalsh(hess)
eigenValue = abs(eigenValue)
eigenValue = np.maximum(eigenValue, 1./hmax**2)
sizeOfCell = 1./np.sqrt(np.max(eigenValue, axis = 1))
sizeOfCell *= 1./np.sqrt(dim/(err * norm))
return sizeOfCell
[docs]
def ObviousSolutionOrScale(node_metrics: NDArray, active_scaling: bool = True) -> Tuple[bool, MuscatIndex, MuscatFloat]:
"""Look if one of the metric to intersect at node is obviously solution,
otherwise compute scaling to apply to help solver CLARABEL to compute intersection problem
Parameter
----------
node_metrics : NDarrays
an array of metrics of the same shape at a given node
a metric is a vector containing the dimF = dim*(dim+1)/2 coefficients of the symmetric definite positive matrix
active_scaling : bool
Activate or not the scaling of the metrics to improve clarabel convergence
Returns
-------
obvious : bool
Is there a metric obviously solution
index : int
Index of the metric solution
scale : float
Scaling factor to apply to the metrics to help solver CLARABEL converge
"""
node_mets = np.abs(node_metrics)
node_mets = np.ma.masked_values(node_mets, 0, rtol = 1e-16, atol = 1e-16)
min_by_met = np.min(node_mets, axis = 1)
max_by_met = np.max(node_mets, axis = 1)
max_min = np.max(min_by_met)
index = np.argmax(min_by_met)
max_other = np.max(np.delete(max_by_met, index))
if max_other < max_min:
fmet = fullMetric(node_metrics[index])
valp = np.linalg.eigvalsh(fmet)
if np.min(valp) > max_other:
obvious = True
return obvious, index, 1.
obvious = False
scale_max = np.max(max_by_met)
scale_min = np.min(min_by_met)
scale = 10**((np.log10(scale_max) + np.log10(scale_min))/2)
if 0.01 < scale < 10 or not active_scaling:
scale = 1.
return obvious, None, scale
[docs]
def IntersectionOfMetrics(metrics: List, active_scaling: bool = True, preprocessing: bool = True) -> NDArray:
"""Construct the Löwner-John metric from a set of metrics
From article LMI approximations for the radius of the intersection of ellipsoids: a survey, D. Henrion, S. Tarbouriech, D. Arzelier, 1998
Parameters
----------
metrics : List of NDarrays
a list of metric nodefields of the same shape
a metric is a vector containing the dimF = dim*(dim+1)/2 coefficients of the symmetric definite positive matrix
active_preprocessing : bool
Activate or not the pre processing of the metrics to improve clarabel convergence
active_scaling : bool
Activate or not the scaling of the metrics
Returns
-------
metric : NDarray
at each node, the Löwner-John metric corresponding to the largest ellipsoid included in the intersection of the input metrics
"""
import cvxpy as cp
metrics = np.copy(metrics)
if len(metrics.shape) == 2:
metrics = metrics[:, np.newaxis, :]
m = metrics.shape[0]
siz = metrics.shape[1]
dimF = metrics.shape[2]
dim = (dimF + 1) // 2
if dim > 3:
raise ValueError(f"dimension can only be 1, 2 or 3")
metric = np.empty_like(metrics[0], dtype = MuscatFloat)
W = cp.Variable(shape = (dim, dim), PSD = True)
l = cp.Variable(shape = m, nonneg = True)
objective = cp.Maximize(cp.atoms.log_det(W))
constraints = [l <= 1]
P = []
for k in range(m):
Pk = cp.Parameter(shape = (dim, dim), PSD = True)
P.append(Pk)
constraints += [W << l[k] * Pk]
problem = cp.Problem(objective, constraints)
for n in range(siz):
scale = 1.
if preprocessing == True:
obvious, idx, scale = ObviousSolutionOrScale(metrics[:, n], active_scaling)
if obvious:
metric[n] = metrics[idx, n]
continue
for k in range(m):
P[k].value = np.linalg.inv(fullMetric(metrics[k, n] / scale))
problem.solve(solver = cp.CLARABEL)
met = np.linalg.inv(W.value)
metric[n] = flattenMetric(met) * scale
if siz == 1:
return metric[0]
return metric
[docs]
def plotMetric(list_met: List):
"""Function to generate a matplolib figure of ellipses/ellipsoids representing metrics
Args:
list_met (List): List of metrics to plot as ellipses,
metrics written under flat mmg format
Returns:
fig (Matplotlib figure): figure that can be plot or saved
"""
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
fig = plt.figure()
dimF = list_met[0].shape[0]
if dimF == 3:
ax = fig.add_subplot()
elif dimF == 6:
ax = fig.add_subplot(111, projection = '3d')
cmap = plt.get_cmap('hsv', len(list_met)+1)
max_radius = 0.
for i, met in enumerate(list_met):
label = "metric {}".format(i)
if i == len(list_met)-1:
label = "Intersection"
fmet = fullMetric(met)
eigenValue, eigenVector = np.linalg.eigh(fmet)
eigenVector = np.transpose(eigenVector)
if dimF == 3:
theta = np.rad2deg(np.arccos(np.dot(eigenVector[0], [1, 0])))*np.sign(eigenVector[0][1])
shape = Ellipse((0,0), 2./np.sqrt(eigenValue[0]), 2./np.sqrt(eigenValue[1]), angle = theta, fill = False, color = cmap(i), label = label)
ax.add_patch(shape)
elif dimF == 6:
coefs = (eigenValue[0], eigenValue[1], eigenValue[2]) # a0/c x**2 + a1/c y**2 + a2/c z**2 = 1
# Radii corresponding to the coefficients:
rx, ry, rz = 1/np.sqrt(coefs)
max_radius = max(rx, ry, rz, max_radius)
# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v))
y = ry * np.outer(np.sin(u), np.sin(v))
z = rz * np.outer(np.ones_like(u), np.cos(v))
# Align ellipsoid axes with eigenvectors
for ij, _ in np.ndenumerate(x):
[x[ij], y[ij], z[ij]] = np.dot([x[ij], y[ij], z[ij]], eigenVector)
# Plot Ellipsoid:
ax.plot_surface(x, y, z, rstride = 4, cstride = 4, color = cmap(i), label = label, alpha = 0.25)
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.axis('on')
plt.legend()
return fig
[docs]
def CheckIntegrity(GUI=False):
try:
import cvxpy
except:
print("Function IntersectionOfMetrics needs dependency to cvxpy, skip this test")
return "skip"
#Article example
met1 = flattenMetric(np.array([[2, 1], [1, 3]]))
met2 = flattenMetric(np.array([[3, -2], [-2, 2]]))
met3 = flattenMetric(np.array([[1, -2], [-2, 6]]))
list_met2D = [met1, met2, met3]
met = IntersectionOfMetrics(list_met2D, preprocessing = True)
check_equality = np.round(met, decimals=5) == np.array([3.21059, -1.0351, 6.42117 ])
if not np.all(check_equality):
raise ValueError ("Wrong values for intersected metric solution")
# Metric included in the second one
met1 = np.array([89880.67130345, -11007.81972081, 1767.87697063])/500000.
met2 = np.array([0.0002, 0., 0.0002])
list_met2D = [met1, met2]
toto = ObviousSolutionOrScale(np.array(list_met2D))
if toto != (True, 0, 1.):
raise ValueError ("Obvious included metric solution not detected")
# Intersection of a metric with a very high anisotropic ratio and a more isotropic one
met1 = np.array([8756.98123322, -295.69639205, 8769.95285002])
met2 = np.array([928636.13929854, -837467.06190469, 755250.03119965])
list_met2D = [met1, met2]
toto = ObviousSolutionOrScale(np.array(list_met2D))
if toto[0] != False:
raise ValueError ("Included metric detected and should not")
print("Star problem")
from Muscat.IO.MeshReader import ReadMesh
from Muscat.TestData import GetTestDataPath
# Compute analytic star
mesh = ReadMesh(GetTestDataPath()+"square2D.mesh")
mesh.ConvertDataForNativeTreatment()
print(mesh)
siz = mesh.GetNumberOfNodes()
psi = np.zeros([siz])
for n in range(siz):
x = mesh.nodes[n][0]
y = mesh.nodes[n][1]
if x*y <= -np.pi/50:
psi[n] = 0.01*np.sin(50*x*y)
elif x*y <= 2*np.pi/50:
psi[n] = np.sin(50*x*y)
else:
psi[n] = 0.01*np.sin(50*x*y)
mesh.nodeFields["psi"] = psi
sizeOfCell = ComputeSizeOfCell(mesh, psi, err = 0.01)
mesh.nodeFields["sizeOfCell"] = sizeOfCell
met1 = ComputeMetric(mesh, psi, err = 0.01)
try:
import cvxpy
except:
print("Function IntersectionOfMetrics needs dependency to cvxpy, skip this test")
else:
# y = x function
psi = np.zeros([siz])
for n in range(siz):
x = mesh.nodes[n][0]
y = mesh.nodes[n][1]
psi[n] = np.exp(-1*abs(x-y))
mesh.nodeFields["psi"] = psi
met2 = ComputeMetric(mesh, psi, err = 0.01, p = 2)
met = IntersectionOfMetrics(metrics = [met1, met2], preprocessing = True)
mesh.nodeFields["met"] = met
for i in range(siz):
fig = plotMetric([met1[i], met2[i], met[i]])
if GUI:
import matplotlib.pyplot as plt
plt.show()
print("Executing mmg")
from Muscat.MeshTools.Remesh import Remesh
nmesh = Remesh(mesh=mesh, solution=None, metric=met1, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print(nmesh)
# 3D Test Case
mesh = ReadMesh(GetTestDataPath()+"square3D.mesh")
siz = mesh.GetNumberOfNodes()
psi = np.ones([siz])
mesh.nodeFields["psi"] = psi
met1 = ComputeMetric(mesh, psi, p = 2, N = 1000)
nodes = mesh.nodes
psi2 = 1.e-4 * np.sin(50 * nodes[:, 0] * nodes[:, 1] * nodes[:, 2])
mesh.nodeFields["psi"] = psi2
met2 = ComputeMetric(mesh, psi2, err = 0.001, gradL2 = True)
met = IntersectionOfMetrics(metrics = [met1, met2], preprocessing = False)
return "ok"
if __name__ == '__main__':
print(CheckIntegrity())