# -*- 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.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.FE.IntegrationTools import IntegrateField
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]
eTags of the computational domain, if None the domain is the entire mesh
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, sdim=mesh.GetPointsDimensionality())
Grad = wf.GetField("Grad", dim, sdim=mesh.GetPointsDimensionality())
GradT = wf.GetTestField("Grad", dim, sdim=mesh.GetPointsDimensionality())
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 is None:
eFilter = ElementFilter(dimensionality=mesh.GetElementsDimensionality())
else:
eFilter = ElementFilter(dimensionality=mesh.GetElementsDimensionality(), eTag=tags)
m, f = IntegrateGeneral(mesh=mesh,
wform=Ener,
constants={},
fields=[field],
unknownFields=unknownFields,
integrationRuleName="NodalEvaluationIsoParam",
elementFilter=eFilter)
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 linear meshes
Parameters
----------
mesh : Mesh
the input mesh
psi : NDArray
a scalar nodefield
tags : List[str]
eTags of the computational domain, if None the domain is the entire mesh
Returns
-------
grad : NDarray
the gradient of psi at each node
"""
for eltype in mesh.elements.keys():
if eltype not in (ED.Point_1, ED.Bar_2, ED.Triangle_3, ED.Tetrahedron_4):
raise ValueError("ComputeGradient_L2projection works only with linear simplex meshes, use ComputeGradient instead")
siz = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
grad = np.zeros([siz, dim], dtype=MuscatFloat)
if tags is None:
eFilter = ElementFilter(dimensionality=mesh.GetElementsDimensionality())
else:
eFilter = ElementFilter(dimensionality=mesh.GetElementsDimensionality(), eTag=tags)
weightAtElement = GetVolumePerElement(mesh, elementFilter=eFilter)
selection = next(eFilter(mesh))
num = selection.Size()
conn = selection.elements.connectivity[selection.indices]
coords = mesh.nodes[conn]
data = psi[conn]
mat = np.zeros([num, dim, dim], dtype=MuscatFloat)
solmat = np.zeros([num, dim], dtype=MuscatFloat)
for i in range(dim):
mat[:, i, :] = coords[:, i+1, :] - coords[:, 0, :]
solmat[:, i] = data[:, i+1] - data[:, 0]
invmat = np.linalg.inv(mat)
res = np.squeeze(np.matmul(invmat, solmat[:, :, None]))
weight = np.zeros(siz, dtype=MuscatFloat)
for i in range(conn.shape[1]):
np.add.at(weight, conn[:, i], weightAtElement)
np.add.at(grad, conn[:, i], res * weightAtElement[:, None])
ids = np.where(weight > 0)[0]
grad[ids, :] = grad[ids, :]/weight[ids, None]
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]
eTags of the computational domain, if None the domain is the entire mesh
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, tags=tags)
for i in range(dim):
h = ComputeGradient_L2projection(mesh=mesh, psi=grad[:, i], tags=tags)
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 ComputeMetricOnSubdomain(mesh: Mesh, psi: NDArray, err: MuscatFloat = None, N: MuscatIndex = 1, p: MuscatIndex = None, gradL2: bool = False, tags: List[str] = None, isotropic: bool = False, renormalization: bool = True) -> NDArray:
"""Construct a metric field for mesh adaptation from a given scalar nodefield phi in Lp norm
From article Continuous mesh framework part II: validations and applications
by Adrien Loseille and Frédéric Alauzet, SIAM Journal of Numerical Analysis, 2011.
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 absolute interpolation error
N: MuscatIndex
desired number of elements in the adapted mesh
tags : List[str]
eTags of the computational subdomain, if None, the entire domain is considered
isotropic : bool, default False
if True, an isotropic scalar metric is returned
flag : bool, default True
if False, the global density is set to 1
Returns
-------
metric : NDarray
the metric at each node
a metric is either a vector of the sizes of the cells (if option isotropic = True)
or a vector containing the (dim+1)*dim/2 coefficients of the associated symmetric definite positive matrix
warning : order of metric component imposed by mmg3d and Inria Tools, different from order read by paraview
"""
dim = mesh.GetPointsDimensionality()
hess = ComputeHessian(mesh=mesh, psi=psi, gradL2=gradL2, tags=tags)
eigenValue, eigenVector = np.linalg.eigh(hess) # hessian is symmetric but neither positive nor definite
gamma = np.maximum(np.abs(eigenValue), 1e-5) # take absolute values of the eigenvalues
ratio = ComputeRatio(eigenValue=gamma)
density = ComputeDensity(mesh=mesh, eigenValue=gamma, ratio=ratio, err=err, N=N, p=p, renormalization=renormalization)
eigenValueOut = np.power(density[:, None]/ratio, 2./dim)
if not isotropic:
metric = BuildMet(eigenValueOut, eigenVector)
else:
metric = 1./np.sqrt(np.max(eigenValueOut, axis=1))
return metric
[docs]
def ComputeMetric(mesh: Mesh, psi: NDArray, err: MuscatFloat = None, N: MuscatIndex = 1, p: MuscatIndex = None, gradL2: bool = False, zones: List[List[str]] = None, isotropic: bool = False) -> NDArray:
"""Construct a metric field for mesh adaptation from a given scalar nodefield phi in Lp norm
From article Continuous mesh framework part II: validations and applications
by Adrien Loseille and Frédéric Alauzet, SIAM Journal of Numerical Analysis, 2011.
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 relative interpolation error (with respect to normLp(psi))
N: MuscatIndex
desired number of elements in the adapted mesh
zones : List[List[str]]
list of the eTags of the computational subdomains, if None, the entire domain is considered
isotropic : bool, default False
if True, an isotropic scalar metric is returned
Returns
-------
metric : NDarray
the metric at each node
a metric is either a vector of the sizes of the cells (if option isotropic = True)
or 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
"""
dim = mesh.GetPointsDimensionality()
if err is not None:
err = err*np.linalg.norm(psi, ord=p) # convert relative error in absolute error
if zones is None:
metric = ComputeMetricOnSubdomain(mesh=mesh, psi=psi, err=err, N=N, p=p, gradL2=gradL2, isotropic=isotropic)
else:
metrics = []
for tags in zones:
metric = ComputeMetricOnSubdomain(mesh=mesh, psi=psi, err=err, N=N, p=p, gradL2=gradL2, tags=tags, isotropic=False, renormalization=False)
metrics.append(metric)
metric = IntersectionOfMetrics(metrics)
eigenValue, eigenVector = np.linalg.eigh(fullMetric(metric))
ratio = ComputeRatio(eigenValue=eigenValue)
density = ComputeDensity(mesh=mesh, eigenValue=eigenValue, ratio=ratio, err=err, N=N, p=p)
eigenValueOut = np.power(density[:, None]/ratio, 2./dim)
if not isotropic:
metric = BuildMet(eigenValueOut, eigenVector)
else:
metric = 1./np.sqrt(np.max(eigenValueOut, axis=1))
return metric
[docs]
def ComputeRatio(eigenValue: NDArray) -> NDArray:
"""Construct at each node of the mesh the ratios r_i = h_i^3 /(h_1 h_2 h_3) where h_i is the length of the optimal metric
It is given by r_1 = gamma_1^-1 * gamma_2^1/2 * gamma_3^1/2
Parameters
----------
eigenValue : NDArray
the absolute values of the eigenValues of the hessian
Returns
-------
ratio : NDarray
"""
siz = eigenValue.shape[0]
dim = eigenValue.shape[1]
ratio = np.ones([siz, dim])
for i in range(dim):
ratio[:, i] = np.sqrt(eigenValue[:, (i+1) % dim]*eigenValue[:, (i+2) % dim])/eigenValue[:, i]
return ratio
[docs]
def ComputeDensity(mesh: Mesh, eigenValue: NDArray, ratio: NDArray, err: MuscatFloat = None, N: MuscatIndex = 1, p: MuscatIndex = None, renormalization: bool = True) -> NDArray:
"""Construct at each node of the mesh the density d = 1 /(h_1 h_2 h_3)
where h_i is the length of the optimal metric
It is given by d = N VO det / integralOfdet
where det = (sum_i ratio_i^(2/dim) * gamma_i )^(dim*p/(2*p+dim)))
Parameters
----------
mesh : Mesh
the input mesh
eigenValue : NDArray
the absolute values of the eigenValues of the hessian
err: MuscatFloat
desired level for the absolute interpolation error
N: MuscatIndex
desired number of elements in the adapted mesh
p: MuscatIndex
order of the Lp-norm (infinite by default)
Returns
-------
density : MuscatFloat
"""
dim = mesh.GetPointsDimensionality()
# compute the local part of the density
if p is None:
det = np.power(np.sum(np.power(ratio, 2./dim)*eigenValue, axis=1), dim/2.)
else:
det = np.power(np.sum(np.power(ratio, 2./dim)*eigenValue, axis=1), dim*p/(2.*p+dim))
if not renormalization: # no global renormalization at this point, it will be done after the intersection
return det
# compute the renormalization factor to have N elements in the mesh
field = FEField(name="det", mesh=mesh, data=det)
integralOfdet = IntegrateField(field)
factor = {2: 8, 3: 10}
V0 = {2: np.sqrt(3)/4, 3: np.sqrt(2)/12}
if err is not None:
if p is None:
N = np.maximum(np.power(factor[dim]*err, -dim/2.)*integralOfdet/V0[dim], 1) # on borne l'erreur afin que N soit toujours plus grand que 1, on pourrait imposer que ce soit un entier aussi
else:
N = np.maximum(np.power(factor[dim]*err, -dim/2.)*np.power(integralOfdet,(2*p+dim)/(2*p))/V0[dim], 1)
density = V0[dim]*N*det/integralOfdet
return density
[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:
metric : NDArray
a vector containing the (dim+1)*dim/2 coefficients of the associated symmetric definite positive matrix
"""
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)
metric = flattenMetric(mat_metric)
return metric
[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 dimension
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[:, None, :]
m = metrics.shape[0]
siz = metrics.shape[1]
dimF = metrics.shape[2]
dim = (dimF+1) // 2
if dim > 3:
raise ValueError("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:
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))
try:
problem.solve(solver=cp.CLARABEL)
except cp.error.SolverError:
print("change clarabel by scs, n")
problem.solve(solver=cp.SCS)
met = np.linalg.inv(W.value)
metric[n] = flattenMetric(met)*scale
if siz == 1:
return metric[0]
return metric
[docs]
def MetricFromElementsToNodes(mesh: Mesh, metricIn: NDArray) -> NDArray:
"""Compute the metric field defined at the nodes, corresponding to the geometric mean of a metric field defined at the elements
Warning : currently only works with linear simplex mesh
Parameters
----------
mesh : Mesh
the input mesh
metricsIn : NDarrays
a metric field defined at each element of the mesh
a metric is a vector containing the dimF = dim*(dim+1)/2 coefficients of the symmetric definite positive matrix
Returns
-------
metricOut : NDarrays
a metric field defined at each node of the mesh
"""
for eltype in mesh.elements.keys():
if eltype not in (ED.Point_1, ED.Bar_2, ED.Triangle_3, ED.Tetrahedron_4):
raise ValueError("MetricFromElementsToNodes works only with linear simplex meshes")
siz = mesh.GetNumberOfNodes()
dimF = metricIn.shape[1]
metricOut = np.zeros([siz, dimF, dimF], dtype=MuscatFloat)
eFilter = ElementFilter(dimensionality=mesh.GetElementsDimensionality())
weightAtElement = GetVolumePerElement(mesh, elementFilter=eFilter)
selection = next(eFilter(mesh))
conn = selection.elements.connectivity[selection.indices]
metricOut = GeometricMeanOfMatrices(metricIn, weightAtElement, conn, siz)
return metricOut
[docs]
def logm_spd(A: NDArray, eps: float = 1e-15) -> NDArray:
# A : SPD, symétrique
evals, evecs = np.linalg.eigh(A)
evals = np.clip(evals, eps, None) # évite les logs de valeurs ~0 ou négatives (erreurs num.)
return (evecs * np.log(evals)) @ evecs.T
[docs]
def expm_sym(S: NDArray) -> NDArray:
# S : symétrique (ici moyenne de logs)
evals, evecs = np.linalg.eigh(S)
return (evecs * np.exp(evals)) @ evecs.T
[docs]
def GeometricMeanOfMatrices(metricIn: NDArray, weightAtElement: NDArray = None, conn: NDArray = None, n_nodes: MuscatIndex = None) -> NDArray:
"""
This function computes the log-euclidian mean of a bunch of matrices. It is a good approximation of the Karcher mean described in
A SURVEY AND COMPARISON OF CONTEMPORARY ALGORITHMS FOR COMPUTING THE MATRIX GEOMETRIC MEAN, BEN JEURIS, RAF VANDEBRIL AND BART VANDEREYCKEN, 2012.
Parameters
----------
matrices : NDarray
a bunch of symmetric definite positive matrices of the same dimension
weights : NDarray
the corresponding weights for the mean, we renormalize the sum of the weights to 1
Returns
-------
metricOut : NDarray
the weighted log-euclidian mean
"""
from scipy import sparse
metricIn = np.asarray(metricIn, dtype=float)
if conn is None:
n_nodes = metricIn.shape[0]*(metricIn.shape[1]+1)
conn = np.arange(n_nodes).reshape(metricIn.shape[0], metricIn.shape[1]+1)
else:
conn = np.asarray(conn, dtype=np.int64)
if weightAtElement is None:
weightAtElement = np.ones(conn.shape[0], dtype=float)
else:
weightAtElement = np.asarray(weightAtElement, dtype=float).reshape(-1)
n_elem, d, _ = metricIn.shape
n_per_elem = conn.shape[1]
# 1) log(A_e) une seule fois
logA = np.empty_like(metricIn)
for e in range(n_elem):
logA[e] = logm_spd(metricIn[e])
# 2) matrice creuse d'incidence pondérée (n_nodes x n_elem)
nodes_flat = conn.reshape(-1)
elem_ids = np.repeat(np.arange(n_elem), n_per_elem)
data = weightAtElement[elem_ids]
M = sparse.csr_matrix(
(data, (nodes_flat, elem_ids)),
shape=(n_nodes, n_elem),
dtype=float
)
# 3) normalisation ligne par ligne (somme des poids par nœud)
row_sums = np.array(M.sum(axis=1)).ravel()
row_sums[row_sums == 0] = 1.0 # évite division par 0 (nœuds isolés éventuels)
Mnorm = sparse.diags(1.0 / row_sums) @ M
# 4) produit sparse : (n_nodes x n_elem) @ (n_elem x d^2) -> (n_nodes x d^2)
X = logA.reshape(n_elem, d*d)
Y = Mnorm.dot(X) # (n_nodes, d*d)
out_log = Y.reshape(n_nodes, d, d)
# 5) expm par nœud
metricOut = np.empty_like(out_log)
for i in range(n_nodes):
metricOut[i] = expm_sym(out_log[i])
return metricOut
[docs]
def ComputeMetricFieldOfExistingMesh(mesh: Mesh):
"""Compute a metric field corresponding to the input existing mesh
Parameter
----------
mesh : Mesh
the input mesh
Returns
-------
metP1 : NDarray
the metric field at each node
"""
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
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
coords = mesh.nodes[conn]
idx_i, idx_j = np.triu_indices(eltDim+1, k=1)
vij = coords[:, idx_j] - coords[:, idx_i]
met = np.einsum('...i,...j->...ij', vij, vij).sum(axis=1)
met = (eltDim+1)/2*np.linalg.inv(met)
metP1 = MetricFromElementsToNodes(mesh, met)
return metP1
[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_ComputeMetric(GUI=False):
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 = ComputeMetric(mesh, psi, err=0.01, isotropic=True)
mesh.nodeFields["sizeOfCell"] = sizeOfCell
met1 = ComputeMetric(mesh, psi, err=0.01)
met2 = ComputeMetric(mesh, psi, err=0.01, p=2)
met3 = ComputeMetric(mesh, psi, N=10000)
met4 = ComputeMetric(mesh, psi, N=1000, p=2)
print("Executing mmg")
from Muscat.MeshTools.Remesh import Remesh
mesh1 = Remesh(mesh=mesh, solution=None, metric=met1, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("2d case 1 = ", mesh1)
mesh2 = Remesh(mesh=mesh, solution=None, metric=met2, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("2d case 2 = ", mesh2)
mesh3 = Remesh(mesh=mesh, solution=None, metric=met3, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("2d case 3 = ", mesh3)
mesh4 = Remesh(mesh=mesh, solution=None, metric=met4, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("2d case 4 = ", mesh4)
# 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.01, p=2)
mesh1 = Remesh(mesh=mesh, solution=None, metric=met1, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("3d case 1 = ", mesh1)
nodes = mesh.nodes
psi2 = 1.e-4 * np.sin(50 * nodes[:, 0] * nodes[:, 1] * nodes[:, 2])
mesh.nodeFields["psi"] = psi2
met2 = ComputeMetric(mesh, psi2, N=10000, gradL2=False)
mesh2 = Remesh(mesh=mesh, solution=None, metric=met2, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print("3d case 2 = ", mesh2)
return "ok"
[docs]
def CheckIntegrity_IntersectionOfMetrics(GUI=False):
try:
import cvxpy
except:
print("Function IntersectionOfMetrics needs dependency to cvxpy, skip this test")
return "skip"
from Muscat.IO.MeshReader import ReadMesh
from Muscat.TestData import GetTestDataPath
# Article example on intersection of metrics
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] is not False:
raise ValueError("Included metric detected and should not")
# 3D Test Case
from Muscat.MeshTools.Remesh import Remesh
mesh = ReadMesh(GetTestDataPath()+"square3D.mesh")
siz = mesh.GetNumberOfNodes()
psi = np.ones([siz])
mesh.nodeFields["psi"] = psi
met1 = ComputeMetric(mesh, psi, err=0.01, p=2)
nodes = mesh.nodes
psi2 = 1.e-4 * np.sin(50 * nodes[:, 0] * nodes[:, 1] * nodes[:, 2])
mesh.nodeFields["psi"] = psi2
met2 = ComputeMetric(mesh, psi2, N=1000, gradL2=True)
met = IntersectionOfMetrics(metrics=[met1, met2], preprocessing=False)
for i in range(siz):
fig = PlotMetric([met1[i], met2[i], met[i]])
if GUI:
import matplotlib.pyplot as plt
plt.show()
nmesh = Remesh(mesh=mesh, solution=None, metric=met, remesher_options={}, backEnd="MmgInMemory", backEndOptions=None)
print(nmesh)
return "ok"
[docs]
def CheckIntegrity_GeoMeanOfMatrices(GUI=False):
from Muscat.IO.MeshReader import ReadMesh
from Muscat.TestData import GetTestDataPath
# Article example on mean of metrics
met1 = np.array([25., 4., 1.])
met2 = np.array([20., 1., 1.])
met3 = np.array([1., 1., 20.])
mat = GeometricMeanOfMatrices(np.array((fullMetric(met1), fullMetric(met2), fullMetric(met3))))
print(mat)
# 2d example on mean of metrics
met1 = np.array([1., 0., 1.])
met2 = np.array([1., 0., 1.])
met3 = np.array([1., 0., 1.])
mat = GeometricMeanOfMatrices(np.array((fullMetric(met1), fullMetric(met2), fullMetric(met3))))
print(mat)
check_equality = np.round(flattenMetric(mat), decimals=5) == np.array([1., 0., 1.])
if not np.all(check_equality):
raise ValueError("Wrong values for geometric mean of metrics")
# 3d example on mean of metrics
met1 = np.array([1., 0., 1., 0., 0., 1.])
met2 = np.array([1., 0., 1., 0., 0., 1.])
met3 = np.array([1., 0., 1., 0., 0., 1.])
print(fullMetric(met1))
mat = GeometricMeanOfMatrices(np.array((fullMetric(met1), fullMetric(met2), fullMetric(met3))))
print(mat)
check_equality = np.round(flattenMetric(mat), decimals=5) == np.array([1., 0., 1., 0., 0., 1.])
if not np.all(check_equality):
raise ValueError("Wrong values for geometric mean of metrics")
# 3D Test Case
mesh = ReadMesh(GetTestDataPath()+"square3D.mesh")
siz = mesh.GetNumberOfElements(dim=mesh.GetElementsDimensionality())
metricIn = np.zeros([siz, 3, 3], dtype=MuscatFloat)
metricIn[:] = fullMetric(np.array([1., 0., 1., 0., 0., 1.]))
_ = MetricFromElementsToNodes(mesh, metricIn)
ComputeMetricFieldOfExistingMesh(mesh)
return "ok"
[docs]
def CheckIntegrity(GUI: bool = False) -> str:
print(CheckIntegrity_ComputeMetric())
print(CheckIntegrity_IntersectionOfMetrics())
print(CheckIntegrity_GeoMeanOfMatrices())
print("done")
return "ok"
if __name__ == '__main__':
print(CheckIntegrity())