# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
import os, bisect
import math
from typing import Tuple, Dict
import numpy as np
from scipy.optimize import fmin_l_bfgs_b
import scipy.spatial as spatial
from Muscat.Containers.Filters.FilterObjects import ElementFilter, NodeFilter
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers import ElementsContainers as EC
from Muscat.Containers import MeshModificationTools as MMT
from Muscat.Containers import MeshFieldOperations as MFO
from Muscat.Containers import MeshInspectionTools as MIT
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.FE.Fields.FEField import FEField
from Muscat.Types import MuscatFloat, ArrayLike
from Muscat.IO import XdmfWriter as XW
[docs]class VolumicMeshFromTargetSurfacesMorpher:
"""Muscat implementation of an original algorithm to morph a 3D Mesh by specifying target surfaces.
STILL LARGRLY EXPERIMENTAL
The input is a 3D mesh, with a partition of its boundary given in the form of NodalTags (intersections
of NodalTags are allowed). The target surfaces (on which the partition of the boundary is to be mapped)
are given as input. First, the algorithm extracts the input mesh boundary and maps the lines between the
partition of this boundary. This mapping is high-quality, by comparing the curvilinear abscissa of the
input and target line. The other (internal boundary points are mapped by RBF, with all the points of the
line as control points. Once these lines have been aligned, the internal points of the surfaces are
projected from the input mesh to the target surfaces using orthogonal projection. At this point, the
surfaces have been matched. The 3D mesh is then mapped by RBF using all the line points and a subset
of the (internal) boundary points as control points.
TODO: code a simpler variant where targetSurfaces are from a single mesh (and compatible elements)
"""
def __init__(self, inputMesh:Mesh, targetSurfaces:Dict[str, Mesh], tol:Dict[str, float] = None, cut_off_distances:Dict[str, float] = None,
largeTransf:bool = True, exportIntermediaryMeshes:bool = False, exportFolder:str = '.', verbose:bool = False):
"""Create an instance of VolumicMeshFromTargetSurfacesMorpher
Parameters
----------
inputMesh : Mesh
the input mesh, must a partition of its boundary given in the form of NodalTags
targetSurfaces : dict[str, Mesh]
a dict of surface target meshes, the keys must correspond to the NodalTags of inputMesh
defining the partition of its boundary
tol : dict[str, float], optional
default: None; a dict of tolerences for merging the surfaces of targetSurfaces, the keys must
correspond to the NodalTags of inputMesh defining the partition of its boundary; it is
recommended to specify a value adapted to your case
cut_off_distances : dict[str, float], optional
default: None; a dict of max distances for selected surface control points in the final 3D RBF
morphing, the keys must correspond to the NodalTags of inputMesh defining the partition of its
boundary it is recommended to specify a value adapted to your case
largeTransf : bool, optional
default: True; False if only local part of the input mesh are expected to be modified, helps
choosing a radius of the RBG morphings
exportIntermediaryMeshes : bool, optional
default: False; True for exporting intermediately produced meshes in exportFolder, in xdmf format
exportFolder : str, optional
default: '.'; folder where intermediately produced meshes are exported, if exportIntermediaryMeshes
is True
verbose : bool, optional
default: False; print additional information when True
"""
self.inputMesh = inputMesh
self.targetSurfaces = targetSurfaces
self.largeTransf = largeTransf
self.exportIntermediaryMeshes = exportIntermediaryMeshes
self.exportFolder = exportFolder
self.verbose = verbose
if tol == None:
bbox = self.inputMesh.ComputeBoundingBox()
tol_ = np.linalg.norm(bbox[1]-bbox[0])/1000.
self.tol = {}
for sn in targetSurfaces.keys():
self.tol[sn] = tol_
else:
self.tol = tol
if cut_off_distances == None:
bbox = self.inputMesh.ComputeBoundingBox()
cod_ = np.linalg.norm(bbox[1]-bbox[0])/20.
self.tol = {}
for sn in targetSurfaces.keys():
self.cut_off_distances[sn] = cod_
else:
self.cut_off_distances = cut_off_distances
self.surfaces_names = sorted(list(self.targetSurfaces.keys()))
assert set(self.surfaces_names) == set(self.tol.keys())
assert set(self.surfaces_names) == set(self.cut_off_distances.keys())
for key in self.targetSurfaces.keys():
assert len(self.inputMesh.GetNodalTag(key).GetIds()>0)
def _ComputeCurvAbsCommon(self, mesh, indices):
"""Compute a normalized curvilinear abscissa on a mesh from the indices of the nodes.
The nodes are assumed to lie on a geometric line. The path is reconstructed
to take into account cases where surface meshes are not compatible.
Parameters
----------
mesh : Mesh
considered mesh
indices : np.ndarray
indices of the nodes lying on a geometric line
Returns
-------
np.ndarray
computed normalized curvilinear abscissa
float
length of the reconstructed line
np.ndarray
indices of the nodes in mesh of the reconstructed line path
"""
candidates = mesh.nodes[indices]
dist_mat = spatial.distance_matrix(candidates, candidates)
i_, j_ = np.unravel_index(dist_mat.argmax(), dist_mat.shape)
diameter = dist_mat[i_, j_]
if np.argmax([np.linalg.norm(candidates[i_]+[-1000, -1000, -1000]), np.linalg.norm(candidates[j_]+[-1000, -1000, -1000])]) == 1:
i_, j_ = j_, i_
numbering = np.arange(candidates.shape[0])
cur_ind = i_
path = [cur_ind]
cur_pos = candidates[cur_ind]
candidates = np.delete(candidates, cur_ind, axis=0)
numbering = np.delete(numbering, cur_ind)
while len(candidates)>0:
tree = spatial.KDTree(candidates)
_, cur_ind = tree.query(cur_pos, k=1)
cur_pos = candidates[cur_ind]
path.append(numbering[cur_ind])
candidates = np.delete(candidates, cur_ind, axis=0)
numbering = np.delete(numbering, cur_ind)
path = [path[i] for i in sorted(np.unique(mesh.nodes[indices[path]], return_index=True, axis=0)[1])]
ranked_ind = indices[path]
close_loop = False
if diameter > 5*np.linalg.norm(mesh.nodes[ranked_ind[0]] - mesh.nodes[ranked_ind[-1]]):
if self.verbose:
print(">> detecting a closed loop while inspecting lines")
close_loop = True
dir = mesh.nodes[ranked_ind[1]] - mesh.nodes[ranked_ind[0]]
if np.argmax([np.linalg.norm(dir-[-1, -1, -1]), np.linalg.norm(dir-[1, 1, 1])]) == 0:
path = path[::-1]
ranked_ind = indices[path]
curv_abs = [0]
for i in range(len(path)-1):
curv_abs.append(curv_abs[i] + np.linalg.norm(mesh.nodes[ranked_ind[i+1]] - mesh.nodes[ranked_ind[i]]))
if close_loop:
curv_abs.append(curv_abs[-1] + np.linalg.norm(mesh.nodes[ranked_ind[-1]] - mesh.nodes[ranked_ind[0]]))
ranked_ind = np.hstack((ranked_ind,ranked_ind[0]))
L = curv_abs[-1]
curv_abs /= L
return curv_abs, L, ranked_ind
def _ComputeCurvAbs(self, mesh, tag):
"""Compute a normalized curvilinear abscissa on a mesh from the indices of the nodes.
The nodes are assumed to lie on a geometric line: case of compatible meshes (and unique line tag)
Parameters
----------
mesh : Mesh
considered mesh
tag : str
considered line tag
Returns
-------
np.ndarray
computed normalized curvilinear abscissa
float
length of the reconstructed line
np.ndarray
indices of the nodes in mesh of the reconstructed line path
"""
indices = mesh.GetNodalTag(tag).GetIds()
return self._ComputeCurvAbsCommon(mesh, indices)
def _ComputeCurvAbsBetween2Tags(self, mesh, tag_1, tag_2):
"""Compute a normalized curvilinear abscissa on a mesh from the indices of the nodes.
The nodes are assumed to lie on a geometric line: case of non-compatible meshes
(two line tags are considered: from surface_1 to surface_2 and surface_2 to surface_1).
Parameters
----------
mesh : Mesh
considered mesh
tag_1 : str
first considered line tag
tag_1 : str
second considered line tag
Returns
-------
np.ndarray
computed normalized curvilinear abscissa
float
length of the reconstructed line
np.ndarray
indices of the nodes in mesh of the reconstructed line path
"""
indices_1 = mesh.GetNodalTag(tag_1).GetIds()
indices_2 = mesh.GetNodalTag(tag_2).GetIds()
indices = np.unique(np.hstack((indices_1, indices_2)))
return self._ComputeCurvAbsCommon(mesh, indices)
def _MapLine(self, nodes_line_mod, curv_abscissa_base, curv_abscissa_mod, L_mod):
"""Maps the nodes of a given line of the input mesh onto the corresponding line of the
output mesh, respecting the curvilinear abscissa.
Parameters
----------
nodes_line_mod : np.ndarray
coordinates of the nodes of the considered line in the target mesh
curv_abscissa_base : np.ndarray
normalized curvilinear abscissa of the considered line in the input mesh
curv_abscissa_mod : np.ndarray
normalized curvilinear abscissa of the considered line in the target mesh
L_mod : float
length of the considered line in the target mesh
Returns
-------
np.ndarray
an array containing the position of the mapped point of the considered line of the input mesh
"""
dim_nodes = nodes_line_mod.shape[1]
nodes_lines_mapped = np.zeros((len(curv_abscissa_base), dim_nodes))
for j in range(len(curv_abscissa_base)-1):
index = min(max(bisect.bisect_left(curv_abscissa_mod, curv_abscissa_base[j])-1, 0), len(curv_abscissa_mod)-1)
assert curv_abscissa_mod[index]<=curv_abscissa_base[j]
assert curv_abscissa_mod[index+1]>curv_abscissa_base[j]
a = nodes_line_mod[index]
b = nodes_line_mod[index+1]
dl = L_mod*(curv_abscissa_base[j] - curv_abscissa_mod[index])
dir = (b-a)/np.linalg.norm(b-a)
assert dl<=np.linalg.norm(b-a)
nodes_lines_mapped[j] = a + dl * dir
nodes_lines_mapped[len(curv_abscissa_base)-1] = nodes_line_mod[-1]
return nodes_lines_mapped
def _TagTargetSurfaces(self) -> Mesh:
"""Tags the lines between the partition surfaces of the boundary of the output mesh.
Since the surface mesh may not be compatible, a decision is taken from a tolerance value
to assess where a point of the boundary of a considered surface lies on the boundary of another surface.
Each geometric line may then be sampled by two sets of points.
Returns
-------
Mesh
A mesh containing the target surfaces and the tagged lines
"""
surfaces_names = sorted(list(self.targetSurfaces.keys()))
for sn in surfaces_names:
self.targetSurfaces[sn].GetNodalTag(sn).AddToTag(np.arange(self.targetSurfaces[sn].GetNumberOfNodes()))
MMT.ComputeSkin(self.targetSurfaces[sn], md=2, inPlace=True)
nf = NodeFilter(eTag="Skin")
self.targetSurfaces[sn].GetNodalTag(f"{sn}_skin").AddToTag(nf.GetNodesIndices(self.targetSurfaces[sn]))
self.targetSurfaces[sn].DeleteElemTags("Skin")
for sn1 in surfaces_names:
space_, numberings_, _, _ = PrepareFEComputation(self.targetSurfaces[sn1])
inputFEField = FEField(name="dummy", mesh=self.targetSurfaces[sn1], space=space_, numbering=numberings_[0])
for sn2 in surfaces_names:
if sn2 != sn1:
op, _, _ = MFO.GetFieldTransferOp(inputFEField, self.targetSurfaces[sn2].nodes)
proj_nodes = op.dot(self.targetSurfaces[sn1].nodes)
res = np.linalg.norm(proj_nodes-self.targetSurfaces[sn2].nodes, axis=1)
ids = np.arange(len(res))[res<max(self.tol[sn1], self.tol[sn2])]
inds_to_add = np.intersect1d(ids, self.targetSurfaces[sn2].GetNodalTag(f"{sn2}_skin").GetIds(), assume_unique=True)
if len(inds_to_add)>0:
self.targetSurfaces[sn2].GetNodalTag(f"{sn2}_{sn1}_line").AddToTag(inds_to_add)
targetSurfaceMesh = Mesh()
for sn in surfaces_names:
self.targetSurfaces[sn].nodesTags.DeleteTags(f"{sn}_skin")
targetSurfaceMesh.Merge(self.targetSurfaces[sn])
return targetSurfaceMesh
def _TagInputMesh(self) -> Mesh:
"""Tags the lines between the partition surfaces of the boundary of the input mesh.
"""
all_surfaces_indices = np.array([], dtype=int)
for sn in self.surfaces_names:
all_surfaces_indices = np.hstack((all_surfaces_indices, self.inputMesh.GetNodalTag(sn).GetIds()))
all_surfaces_indices = np.unique(all_surfaces_indices)
self.inputMesh.nodesTags.DeleteTags("all_surfaces")
self.inputMesh.GetNodalTag("all_surfaces").AddToTag(all_surfaces_indices)
ef = ElementFilter(nTag="all_surfaces")
inputSurfaceMesh = MIT.ExtractElementsByElementFilter(self.inputMesh, ef)
inputSurfaceMesh.DeleteElemTags("Skin")
inputSurfaceMesh.nodesTags.DeleteTags('Skin')
MMT.CleanLonelyNodes(inputSurfaceMesh)
for i, sn1 in enumerate(self.surfaces_names):
ids_1 = inputSurfaceMesh.GetNodalTag(f"{sn1}").GetIds()
for j, sn2 in enumerate(self.surfaces_names):
if i<j:
ids_2 = inputSurfaceMesh.GetNodalTag(f"{sn2}").GetIds()
inds_to_add = np.intersect1d(ids_1, ids_2)
if len(inds_to_add)>0:
inputSurfaceMesh.GetNodalTag(f"{sn1}_{sn2}_line").AddToTag(inds_to_add)
for sn in self.surfaces_names:
ef = ElementFilter(nTag=sn)
local_mesh = MIT.ExtractElementsByElementFilter(inputSurfaceMesh, ef)
MMT.ComputeSkin(local_mesh, md=2, inPlace=True)
ef2 = ElementFilter(eTag="Skin", dimensionality=1)
local_mesh_boundary = MIT.ExtractElementsByElementFilter(local_mesh, ef2)
for elCont in local_mesh_boundary.elements:
bound_elements = local_mesh_boundary.GetElementsOfType(elCont.elementType)
if elCont.elementType not in inputSurfaceMesh.elements:
inputSurfaceMesh.elements.AddContainer(EC.ElementsContainer(elCont.elementType))
inputSurfaceMesh.elements[elCont.elementType].Merge(bound_elements)
inputSurfaceMesh.DeleteElemTags("Skin")
inputSurfaceMesh.nodesTags.DeleteTags(['Skin', 'all_surfaces'])
return inputSurfaceMesh
def _MorphLines(self, inputSurfaceMesh:Mesh, targetSurfaceMesh:Mesh):
"""Applies the first step: the mapping of the lines between the partition of
the boundaries from the input mesh to the target mesh (following curvilinear abscissas).
Parameters
----------
inputSurfaceMesh : Mesh
a mesh containing the surface boundaries on the input mesh, and line tags
targetSurfaceMesh : Mesh
a mesh containing the surface boundaries on the target mesh, and line tags
"""
if self.verbose:
print("Morphing using lines as control points")
inputSurfaceMesh_new_nodes = np.copy(inputSurfaceMesh.nodes)
mask_base = np.array([], dtype = int)
for i, sn1 in enumerate(self.surfaces_names):
for j, sn2 in enumerate(self.surfaces_names):
if i<j:
if f"{sn1}_{sn2}_line" in inputSurfaceMesh.nodesTags and f"{sn1}_{sn2}_line" in targetSurfaceMesh.nodesTags:
curv_abs_base, _, ind_path_base = self._ComputeCurvAbs(inputSurfaceMesh, f"{sn1}_{sn2}_line")
curv_abs_mod, L_mod, ind_path_mod = self._ComputeCurvAbsBetween2Tags(targetSurfaceMesh, f"{sn1}_{sn2}_line", f"{sn2}_{sn1}_line")
inputSurfaceMesh_new_nodes[ind_path_base] = self._MapLine(targetSurfaceMesh.nodes[ind_path_mod,:], curv_abs_base, curv_abs_mod, L_mod)
mask_base = np.hstack((mask_base, ind_path_base), dtype = int)
mask_base = np.unique(mask_base)
targetDisplacement = inputSurfaceMesh_new_nodes[mask_base,:] - inputSurfaceMesh.nodes[mask_base,:]
targetDisplacementMask = mask_base
if self.largeTransf:
bbox = inputSurfaceMesh.ComputeBoundingBox()
radius = 2.*np.linalg.norm(bbox[1]-bbox[0] + np.max(targetDisplacement, axis=0))
else:
radius = None
inputSurfaceMesh.nodes = MMT.Morphing(inputSurfaceMesh, targetDisplacement, targetDisplacementMask, radius = radius, verbose = self.verbose)
inputSurfaceMesh.nodesTags.DeleteTags("line_control_points")
inputSurfaceMesh.GetNodalTag("line_control_points").AddToTag(mask_base)
def _MorphSurfaces(self, inputSurfaceMesh:Mesh, targetSurfaceMesh:Mesh):
"""Applies the second step: the internal points of the surfaces are
projected from the input mesh to the target surfaces using orthogonal projection.
Parameters
----------
inputSurfaceMesh : Mesh
a mesh containing the surface boundaries on the input mesh, and line tags
targetSurfaceMesh : Mesh
a mesh containing the surface boundaries on the target mesh, and line tags
"""
if self.verbose:
print("Morphing using lines as control points and projecting onto target surfaces")
line_control_points = inputSurfaceMesh.GetNodalTag("line_control_points").GetIds()
for sn in self.surfaces_names:
surfaces_modified = MIT.ExtractElementsByElementFilter(targetSurfaceMesh, ElementFilter(nTag=sn, dimensionality=2))
space_, numberings_, _, _ = PrepareFEComputation(surfaces_modified)
inputFEField = FEField(name="dummy", mesh=surfaces_modified, space=space_, numbering=numberings_[0])
ids_to_treat = np.setdiff1d(inputSurfaceMesh.GetNodalTag(sn).GetIds(), line_control_points)
op, _, _ = MFO.GetFieldTransferOp(inputFEField, inputSurfaceMesh.nodes[ids_to_treat])
inputSurfaceMesh.nodes[ids_to_treat] = op.dot(surfaces_modified.nodes)
def _MorphVolume(self, inputSurfaceMesh:Mesh):
"""Applies the last step: the 3D RBF morphing using all the line points and a
subset of the (internal) boundary points as control points
Parameters
----------
inputSurfaceMesh : Mesh
a mesh containing the surface boundaries on the input mesh, and line tags
"""
if self.verbose:
print("Morphing using lines and surface points from inputSurfaceMesh as control points")
all_surfaces_indices = {}
surface_indices = np.array([], dtype=int)
for sn in self.surfaces_names:
all_surfaces_indices[sn] = self.inputMesh.GetNodalTag(sn).GetIds()
surface_indices = np.hstack((surface_indices, self.inputMesh.GetNodalTag(sn).GetIds()))
surface_indices = np.unique(surface_indices)
line_indices_in_3D_mesh = inputSurfaceMesh.GetNodalTag("line_control_points").GetIds()
all_surface_indices = np.array([], dtype=int)
for sn in self.surfaces_names:
local_surface_indices = inputSurfaceMesh.GetNodalTag(sn).GetIds()
dists = spatial.distance.cdist(inputSurfaceMesh.nodes[line_indices_in_3D_mesh], inputSurfaceMesh.nodes[local_surface_indices])
min_dist = np.min(dists, axis=0)
local_surface_indices = local_surface_indices[min_dist > self.cut_off_distances[sn]]
ii = 0
while ii < local_surface_indices.shape[0]:
tree2 = spatial.KDTree(inputSurfaceMesh.nodes[local_surface_indices])
ind = tree2.query(inputSurfaceMesh.nodes[local_surface_indices[ii]], k=2)[1][1]
d = np.linalg.norm(inputSurfaceMesh.nodes[local_surface_indices[ii]]-inputSurfaceMesh.nodes[local_surface_indices][ind])
if d < self.cut_off_distances[sn]:
local_surface_indices = np.hstack((local_surface_indices[:ind],local_surface_indices[ind+1:]))
else:
ii += 1
all_surface_indices = np.hstack((all_surface_indices, local_surface_indices))
extract_local_surface_indices = np.hstack((line_indices_in_3D_mesh, all_surface_indices))
extract_local_surface_indices = np.unique(extract_local_surface_indices)
extract_surface_indices = surface_indices[extract_local_surface_indices]
targetDisplacementMask = extract_surface_indices
targetDisplacement = inputSurfaceMesh.nodes[extract_local_surface_indices,:] - self.inputMesh.nodes[extract_surface_indices]
if self.largeTransf:
bbox = self.inputMesh.ComputeBoundingBox()
radius = 2.*np.linalg.norm(bbox[1]-bbox[0] + np.max(targetDisplacement, axis=0))
else:
radius = None
self.inputMesh.nodes = MMT.Morphing(self.inputMesh, targetDisplacement, targetDisplacementMask, radius = radius, verbose = self.verbose)
self.inputMesh.nodesTags.DeleteTags("surface_control_points")
self.inputMesh.GetNodalTag("surface_control_points").AddToTag(targetDisplacementMask)
[docs] def Morph(self) -> Mesh:
"""Applies the morphing algorithm to self.inputMesh
Returns
-------
Mesh
the morphed mesh
"""
# Tagging meshes
inputSurfaceMesh = self._TagInputMesh()
targetSurfaceMesh = self._TagTargetSurfaces()
if self.verbose == True and self.exportIntermediaryMeshes == True:
print("exporting inputSurfaceMesh.xdmf and targetSurfaceMesh.xdmf")
if self.exportIntermediaryMeshes == True:
XW.WriteMeshToXdmf(os.path.join(self.exportFolder, "inputSurfaceMesh.xdmf"), inputSurfaceMesh)
XW.WriteMeshToXdmf(os.path.join(self.exportFolder, "targetSurfaceMesh.xdmf"), targetSurfaceMesh)
# Morphing inputSurfaceMesh using lines as control points
self._MorphLines(inputSurfaceMesh, targetSurfaceMesh)
if self.verbose == True and self.exportIntermediaryMeshes == True:
print("exporting inputSurfaceMeshLineMorphed.xdmf")
if self.exportIntermediaryMeshes == True:
XW.WriteMeshToXdmf(os.path.join(self.exportFolder, "inputSurfaceMeshLineMorphed.xdmf"), inputSurfaceMesh)
# Morphing inputSurfaceMesh using lines as control points and projecting onto target surfaces
self._MorphSurfaces(inputSurfaceMesh, targetSurfaceMesh)
if self.verbose == True and self.exportIntermediaryMeshes == True:
print("exporting inputSurfaceMeshLineSurfaceMorphed.xdmf")
if self.exportIntermediaryMeshes == True:
XW.WriteMeshToXdmf(os.path.join(self.exportFolder, "inputSurfaceMeshLineSurfaceMorphed.xdmf"), inputSurfaceMesh)
# Morphing inputMesh using lines and surface points from inputSurfaceMesh as control points
self._MorphVolume(inputSurfaceMesh)
XW.WriteMeshToXdmf(os.path.join(self.exportFolder, "morphedMesh.xdmf"), self.inputMesh)
return self.inputMesh
[docs]class FoldOverFreeMaps:
"""Muscat implementation of the Fold Over Free Maps algorithm"""
def __init__(self, mesh: Mesh, boundaryEF: ElementFilter):
"""Create an instance of FoldOverFreeMaps
Parameters
----------
mesh : Mesh
The mesh to work on
boundaryEF : ElementFilter
an element filter to define the border of the mesh
"""
self.eps = 1
self.pdim = mesh.GetPointsDimensionality()
self.edim = mesh.GetElementsDimensionality()
self.nbNodes = mesh.GetNumberOfNodes()
self.mesh = mesh
self.nodeMask = np.zeros((self.nbNodes, self.pdim), dtype=bool)
for selection in boundaryEF(mesh):
nbIds = np.unique(selection.elements.connectivity[selection.indices, :])
self.nodeMask[nbIds, :] = True
from Muscat.FE.IntegrationRules import NodalEvaluationP1Quadrature
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, FESpace
self.space: FESpace = LagrangeSpaceP1
self.ip = NodalEvaluationP1Quadrature
self.spaceAtIP = {name: self.space[name].GetSpaceEvaluatedAt(ir.points, ir.weights) for name, ir in NodalEvaluationP1Quadrature.items()}
self.theta = 0.5
self.callback = None
"""callback function to be called at each interaction of the algorithm with the current mesh as argument """
[docs] def CallCallback(self):
"""Call Back wrapper"""
if self.callback is not None:
self.callback(self.mesh)
[docs] def Start(self):
"""Start the computation of the algorithm"""
n = 0
while True: # outer L-BFGS loop, Alg. 1, line 2
[Fprev, grad] = self.ComputePotential(self.mesh.nodes.flatten()) # compute $F(U^k, \varepsilon^k)$
# print(f"Fprev {Fprev}")#; exit()
[newNodes, F, _] = fmin_l_bfgs_b(self.ComputePotential, self.mesh.nodes.flatten(), factr=1e12) # inner L-BFGS loop, compute $F(U^{k+1}, \varepsilon^k)$
# newNodes = mesh.nodes.flatten() - grad/(np.linalg.norm(grad))*0.01/energy.eps
# print(newNodes)
newNodes.shape = self.mesh.nodes.shape
self.mesh.nodes = newNodes
mindet = self.MinDetJacobian()
mu = min(F / Fprev, 9e-1) * (mindet + math.sqrt(self.eps**2 + mindet**2)) / 2 # $\mu^k$, Eq. (6)
self.eps = 2 * math.sqrt(mu * (mu - mindet)) if mindet < mu else 0 # $\varepsilon^{k+1}$, Eq. (6)
self.CallCallback()
if mindet > 0 and F > 999e-3 * Fprev:
break # stopping criterion, Alg. 1, line 11
print(f"{n}, {self.eps:.5f}, {Fprev:.5f}, {mindet:.5f}")
n += 1
if n > 200: # pragma: no cover
print("Maximal number of iterations reached.")
break
[docs] def ComputePotential(self, U: ArrayLike) -> Tuple[np.number, np.ndarray]:
"""Compute the potential to be minimized
Parameters
----------
U : ArrayLike
The positions of the points
Returns
-------
Tuple[np.number, np.ndarray]
F the potential value
G the gradient of the potential
"""
# compute the energy $F$ and its gradient $\nabla F$ for the map $\vec{u}$, Eq. (5)
F = 0 # zero out the energy and the gradient
G = np.zeros((self.nbNodes, self.pdim), dtype=MuscatFloat)
nodes = U.view()
nodes.shape = self.mesh.nodes.shape
for selection in ElementFilter(dimensionality=self.edim)(self.mesh):
nbip = self.ip[selection.elementType].GetNumberOfPoints()
sAtIP = self.spaceAtIP[selection.elementType]
eps2 = self.eps**2
for id in selection.indices:
elConnectivity = selection.elements.connectivity[id, :]
for qc in range(nbip): # for every quad corner
J, det, jinv = sAtIP.GetJacobianDeterminantAndInverseAtIP(qc, nodes[elConnectivity, :])
J = J.T
chi = det / 2 + math.sqrt(eps2 + det**2) / 2 # the penalty function $\chi(\varepsilon^k, \mathrm{det}\ J)$
chip = 0.5 + det / (2 * math.sqrt(eps2 + det**2)) # its derivative $\chi'(\varepsilon^k, \mathrm{det}\ J)$
# tjtj =
tjtj = np.sum(J * J) # ♀np.trace(np.transpose(J).dot(J))
f = tjtj / (chi ** (2 / self.edim)) # quad corner shape quality
g = (det**2 + 1) / chi
F += (1 - self.theta) * f # add the term to the energy
F += self.theta * g
dfdj = (1 - self.theta) * (2 * J - det * jinv(np.eye(2)) * f * chip) / chi # $\frac{\partial f_\varepsilon}{\partial a}$: derivative w.r.t the Jacobian
if self.theta != 0:
dgdj = (self.theta) * det * jinv(((2 * det - g * chip) / chi))
dfdj += dgdj
funcDer = sAtIP.valdphidxi[qc].T
dfdu = funcDer.dot(np.transpose(dfdj)) # chain rule for the real variables (Appendix A)
G[elConnectivity, :] += dfdu[:, :]
G[self.nodeMask] = 0
return F, G.flatten()
[docs] def MinDetJacobian(self) -> np.number:
"""Compute the minimal jacobian on the mesh
Returns
-------
np.number
the min of the jacobian
"""
res = np.inf
for selection in ElementFilter(dimensionality=self.edim)(self.mesh):
ip = self.ip[selection.elementType]
sAtIP = self.spaceAtIP[selection.elementType]
lenIP = ip.GetNumberOfPoints()
for i in selection.indices:
for qc in range(lenIP): # for every quad corner
J, det, jinv = sAtIP.GetJacobianDeterminantAndInverseAtIP(qc, self.mesh.nodes[selection.elements.connectivity[i, :], :])
res = min(res, det)
return res
[docs]def CheckIntegrity_FOFM(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateSquare
n = 8
mesh = CreateSquare(dimensions=[n, n], origin=[0, 0], spacing=[0.875 / (n - 1), 1.75 / (n - 1)], ofTriangles=True)
mask = mesh.nodes[:, 1] > 0.8625
mesh.nodes[mask, 0] += 0.6
mesh.nodes[mask, 1] -= 0.6
from Muscat.IO.XdmfWriter import XdmfWriter
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
writer = XdmfWriter(fileName=TemporaryDirectory.GetTempPath() + "FoldOverFreeMaps_CheckIntegrity.xdmf")
print(writer.fileName)
writer.SetTemporal(True)
writer.SetBinary()
writer.SetHdf5(False)
writer.Open()
writer.Write(mesh)
from Muscat.Containers.Filters.FilterObjects import ElementFilter
FOFM = FoldOverFreeMaps(mesh, boundaryEF=ElementFilter(dimensionality=[0, 1, 3]))
FOFM.callback = writer.Write
FOFM.Start()
writer.Close()
if FOFM.MinDetJacobian() < 0: # pragma: no cover
raise Exception("FoldOverFreeMaps unable to unfold the mesh ")
return "ok"
[docs]def CheckIntegrity_MorphVolumicMeshFromTargetSurfaces(GUI: bool = False):
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.Containers.Filters.FilterObjects import ElementFilter, NodeFilter
from Muscat.Containers import MeshInspectionTools as MIT
from Muscat.Containers import MeshModificationTools as MMT
from Muscat.Containers import MeshCreationTools as MCT
surfaces_names_ = ['X0', 'X1', 'Y0', 'Y1', 'Z0', 'Z1']
inputMesh_ = MCT.CreateCube(dimensions=[10, 12, 14], spacing=[0.6, 0.8, 1.1], ofTetras=True)
targetMesh_ = MCT.CreateCube(dimensions=[8, 6, 9], spacing=[1., 0.85, 0.9], ofTetras=True)
for sn in surfaces_names_:
nf = NodeFilter(eTag=sn)
for mesh in [inputMesh_, targetMesh_]:
mesh.GetNodalTag(sn).AddToTag(nf.GetNodesIndices(mesh))
targetSurfaces_ = {}
tol_ = {}
cut_off_distances_ = {}
for sn in surfaces_names_:
ef = ElementFilter(eTag=sn)
inputSurfaceMesh = MIT.ExtractElementsByElementFilter(targetMesh_, ef)
MMT.CleanLonelyNodes(inputSurfaceMesh)
targetSurfaces_[sn] = inputSurfaceMesh
tol_[sn] = 1.e-4
cut_off_distances_[sn] = 0.9
morpher = VolumicMeshFromTargetSurfacesMorpher(inputMesh_, targetSurfaces_, tol_, cut_off_distances_, exportIntermediaryMeshes = False,
largeTransf = False, exportFolder = TemporaryDirectory.GetTempPath(), verbose = False)
morpher.Morph()
morpher = VolumicMeshFromTargetSurfacesMorpher(inputMesh_, targetSurfaces_, tol_, cut_off_distances_, exportIntermediaryMeshes = True,
largeTransf = True, exportFolder = TemporaryDirectory.GetTempPath(), verbose = True)
morphedMesh = morpher.Morph()
ref = np.asarray([[-1. , -1. , -1. ],
[-1. , 2.05250214, 0.70404702],
[-0.29382286, 0.43089061, 1.62046775],
[-0.29699579, 2.5237726 , 2.91439055],
[ 0.42389267, 1.86674657, 5.49586333],
[ 0.46644238, 3.25 , 4.44983021],
[ 1.23958503, 1.4405376 , 5.67279348],
[ 2.11533487, -0.19492609, -0.99932695],
[ 2.1227032 , 2.36646581, 0.35226376],
[ 2.91236193, 1.70652936, 2.40607259],
[ 2.92152694, 3.24861093, 2.32432152],
[ 3.59440583, 2.61524682, 5.88970713],
[ 4.59245849, -0.3497334 , 4.5109741 ],
[ 4.56086724, 1.95685638, 5.43289463],
[ 5.31654192, 0.50982711, -1.00011092],
[ 5.30734042, 3.19436004, 0.08784557],
[ 5.96100323, 2.4494159 , 2.59347265]])
np.testing.assert_almost_equal(morphedMesh.nodes[::100], ref)
return "ok"
[docs]def CheckIntegrity(GUI:bool=False):
functionsToTest= [
CheckIntegrity_FOFM,
CheckIntegrity_MorphVolumicMeshFromTargetSurfaces,
]
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
return RunListOfCheckIntegrities(functionsToTest)
if __name__ == '__main__':
print(CheckIntegrity(False))# pragma: no cover