# -*- 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 Dict, List, Tuple
from Muscat.Types import MuscatFloat, MuscatIndex
from numpy.typing import NDArray
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Containers.Mesh import Mesh
import Muscat.Containers.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.Containers.MeshInspectionTools import GetVolumePerElement
from Muscat.Containers.InriaMeshTools import flattenMetric, fullMetric
import time
[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, tags: List[str] = None) -> 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 el in range(num):
grad[conn[el, :], :] += res[el,:]*weightAtElement[el]
weight[ conn[el, :] ] += weightAtElement[el]
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
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, psi, tags = tags)
for i in range(dim):
h = ComputeGradient_L2projection(mesh, grad[:,i], tags = tags)
hess[:, i, :] += h
hess[:, :, i] += h
else:
grad = ComputeGradient(mesh, psi, tags = tags)
for i in range(dim):
h = ComputeGradient(mesh, 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, hmin: MuscatFloat = 0.001, hmax: MuscatFloat = 1000, 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
hmin, hmax : MuscatFloat
bounds for the size of the 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
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
dimF = (dim+1)*dim // 2
metric = np.zeros([siz, dimF], dtype = MuscatFloat)
det = np.zeros(siz, dtype = MuscatFloat)
if zones == None:
hess = ComputeHessian(mesh, psi, gradL2)
else:
hess = 0
for tags in zones:
h = ComputeHessian(mesh, psi, gradL2, tags = tags)
hess = hess*(h == 0) + (hess == 0)*h
for n in range(siz):
if np.isclose(hess[n], np.zeros([dim,dim]), rtol=1e-14, atol=1e-14).all():
mat_metric = np.zeros([dim, dim], dtype = MuscatFloat)
# order of metric component imposed by mmg3d, different from order read by paraview
np.fill_diagonal(mat_metric, 1./hmax**2)
metric[n] = flattenMetric(mat_metric)
det[n] = 1./hmax**(2*dim)
else:
eigenValue, eigenVector = np.linalg.eigh(hess[n]) # computing eigenvalues EigenValue and matrix of eigenvectors EigenVector of the hessian matrix
eigenValue = abs(eigenValue) # eigenValue = 1/h^2 where h is the local mesh size
eigenValue = np.minimum(np.maximum(eigenValue, 1./hmax**2), 1./hmin**2)
det[n] = np.prod(eigenValue)
mat_metric = eigenVector @ np.diag(eigenValue) @ np.linalg.inv(eigenVector) # absolute value of the eigenvalues
metric[n] = flattenMetric(mat_metric) # order of metric component imposed by mmg3d, different from order read by paraview
# on pourra enlever la dependance en le maillage tetrahedrique lineaire ici
# et calculer l'integrale du determinant grace a l integrateur de muscat
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
weightAtElement = GetVolumePerElement(mesh, ElementFilter(elementType = eltype))/(dim+1)
integralOfDet = 0
if p is None:
for el in range(cell.GetNumberOfElements()):
integralOfDet += weightAtElement[el]*np.sum(np.power(det[conn[el,:]], 1./2))
else:
for el in range(cell.GetNumberOfElements()):
integralOfDet += weightAtElement[el]*np.sum(np.power(det[conn[el,:]], p/(2*p+dim)))
if N is None:
metric *= dim/err
if p is not None:
metric *= np.power(integralOfDet, 1./p)
else:
metric *= np.power(N,2./dim)*np.power(integralOfDet, -2./dim)
if p is not None:
for n in range(siz):
metric[n] *= np.power(det[n], -1./(2*p+dim))
return metric
[docs]def ObviousSolutionOrScale(node_metrics: NDArray, active_scaling=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 : List of NDarrays
a list of metric at same node of the same shape
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 node 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.eigh(fmet)[0]
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=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_scaling : bool
Activate or not the scaling of node metrics to improve clarabel convergence
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] # number of input metric fields in the List
siz = metrics.shape[1] # number of nodes in the fields
dimF = metrics.shape[2] # size of the flatten metric vector
if dimF == 1:
dim = 1
inverse1 = lambda M: [1/M[0]]
inverse2 = lambda W: [1/W[0, 0]]
elif dimF == 3:
dim = 2
def inverse1(M: NDArray):
det = M[0] * M[2] - M[1] * M[1]
inv = [[M[2], -M[1]], [-M[1], M[0]]] / det
return inv
def inverse2(W: NDArray):
det = W[0, 0] * W[1, 1] - W[0, 1] * W[0, 1]
inv = [W[1, 1], -W[0, 1], W[0, 0]] / det
return inv
elif dimF == 6:
dim = 3
def inverse1(M: NDArray): # invert vector M as a square symmetric matrix
det = M[0] * (M[2] * M[5] - M[4] * M[4]) + M[1] * (M[4] * M[3] - M[5] * M[1]) + M[3] * (M[1] * M[4] - M[3] * M[2])
inv = [[M[2] * M[5] - M[4] * M[4], M[4] * M[3] - M[5] * M[1], M[1] * M[4] - M[3] * M[2]],
[M[4] * M[3] - M[5] * M[1], M[5] * M[0] - M[3] * M[3], M[3] * M[1] - M[4] * M[0]],
[M[1] * M[4] - M[3] * M[2], M[3] * M[1] - M[4] * M[0], M[0] * M[2] - M[1] * M[1]]] / det
return inv
def inverse2(W: NDArray): # invert and flatten matrix W (order of metric components imposed by mmg3d)
det = W[0, 0] * (W[1, 1] * W[2, 2] - W[1, 2] * W[2, 1]) + W[0, 1] * (W[1, 2] * W[2, 0] - W[2, 2] * W[1, 0]) + W[0, 2] * (W[1, 0] * W[2, 1] - W[2, 0] * W[1, 1])
inv = [W[1, 1] * W[2, 2] - W[2, 1] * W[1, 2], W[1, 2] * W[2, 0] - W[1, 0] * W[2, 2], W[2, 2] * W[0, 0] - W[0, 2]* W[2, 0],
W[1, 0] * W[2, 1] - W[2, 0] * W[1, 1], W[2, 0] * W[0, 1] - W[0, 0] * W[2, 1], W[0, 0] * W[1, 1] - W[1, 0] * W[0, 1]] / det
return inv
else:
raise ValueError(f"dimension can only be 1, 2 or 3")
metric = np.empty_like(metrics[0], dtype=MuscatFloat)
W = cp.Variable((dim, dim), PSD = True, name = "W")
l = cp.Variable(m, nonneg = True, name = "l")
for n in range(siz):
constraints = [l <= 1.]
obvious, idx, scale = ObviousSolutionOrScale(metrics[:,n], active_scaling)
if obvious:
metric[n] = metrics[idx, n]
else:
for k in range(m):
inv = inverse1(metrics[k, n] / scale)
constraints.append(W << l[k] * inv)
objective = cp.Maximize(cp.atoms.log_det(W))
problem = cp.Problem(objective, constraints)
problem.solve(solver='CLARABEL')
metric[n] = inverse2(W.value) * scale
if siz == 1:
return metric[0]
return metric
[docs]def ComputeSizeOfCell(mesh: Mesh, psi: NDArray, err:MuscatFloat = 0.01, hmin:MuscatFloat = 0.001, hmax:MuscatFloat = 100, 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
hmin, hmax : MuscatFloat
bounds for the size of the mesh
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
"""
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
sizeOfCell = np.zeros(siz, dtype = MuscatFloat)
if zones == None:
hess = ComputeHessian(mesh, psi)
else:
hess = 0
for tags in zones:
h = ComputeHessian(mesh, psi, tags = tags)
hess = hess*(h == 0) + (hess == 0)*h
for n in range(siz):
if np.isclose(hess[n], np.zeros([dim,dim]), rtol=1e-14, atol=1e-14).all():
sizeOfCell[n] = 1./hmax**2
else:
eigenValue, eigenVector = np.linalg.eigh(hess[n]) # computing eigenvalues EigenValue and matrix of eigenvectors EigenVector of the hessian matrix
eigenValue = abs(eigenValue) # eigenValue = 1/h^2 where h is the local mesh size
eigenValue = np.minimum(np.maximum(eigenValue, 1./hmax**2), 1./hmin**2)
sizeOfCell[n] = 1./np.sqrt(np.max(eigenValue))
sizeOfCell *= 1./np.sqrt(dim/err)
return sizeOfCell
[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):
#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)
check_equality = np.round(met, decimals=8) == np.array([3.21058607, -1.03509769, 6.4211717 ])
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.0001)
mesh.nodeFields["sizeOfCell"] = sizeOfCell
met1 = ComputeMetric(mesh, psi, err = 0.0001)
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.1, p=2)
met = IntersectionOfMetrics(metrics = [met1, met2])
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.Containers.MmgRemeshing import MmgRemeshing
obj = MmgRemeshing(mesh = mesh, solution = None, metric = met1)
nmesh = obj.Runner(options = {}, temp_name = "star_problem", binary = True)
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, err = 0.0001, 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.0001, gradL2=True)
met = IntersectionOfMetrics(metrics = [met1, met2])
return "ok"
if __name__ == '__main__':
print(CheckIntegrity())