# -*- 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.
#
from typing import Union
from pathlib import Path
from Muscat.Types import MuscatFloat
from Muscat.IO.MeshWriter import WriteMesh
from Muscat.IO.MeshReader import ReadMesh
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.Helpers.IO.Which import Which
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.MeshTools.Remesh import RegisterRemesherClass
[docs]
class FefloRemeshing():
def __init__(self, mesh: Mesh, solution=None, levelset=None, metric=None):
self.mesh = mesh
if levelset is not None:
raise RuntimeError(' FefloRemeshing does not support the argument levelset')
self.solution = solution
self.metric = metric
self.temp_out_path = None
def __WriteMeshAndSolution(self, temp_path: Union[str, Path] = Path("temp_file"), binary: bool = True):
"""Write mesh and solution in a mesh file
Parameters
----------
temp_path : Union[str,Path]
path of the output mesh file (without extension)
binary : bool
"""
from Muscat.MeshTools.MeshTools import TagsToRefs
solutionOnOwnFile = True
self.temp_path = Path(str(temp_path))
self.temp_mesh_path = Path(str(temp_path)+".mesh"+binary*"b")
self.temp_sol_path = Path(str(temp_path)+".sol"+binary*"b")
if not self.temp_out_path:
self.temp_out_path = Path(str(temp_path)+".o.mesh"+binary*"b")
if self.solution is not None:
solution = [self.solution]
elif self.metric is not None:
solution = [self.metric]
else:
solution = None
solutionOnOwnFile = False
RefByNode, self.dictNodesRefsToTags, RefByElement, self.dictElemRefsToTags = TagsToRefs(self.mesh)
WriteMesh(str(self.temp_mesh_path), self.mesh, PointFields=solution, solutionOnOwnFile=solutionOnOwnFile, binary=binary, nodalRefNumber=RefByNode, elemRefNumber=RefByElement)
def __BuildCommandLine(self, remesher_options: dict, prepro: bool = False):
"""Build command line for feflo execution
Parameters
----------
remesher_options : dict
keys corresponding to feflo option name
and value corresponding to value if applicable
Returns
-------
CommandLine : str
the command line for feflo execution
"""
CommandLine = []
CommandLine.append(Which("feflo.a"))
CommandLine.append("-in")
CommandLine.append(str(self.temp_path))
if not prepro and self.solution is not None:
CommandLine.append("-sol")
CommandLine.append(str(self.temp_sol_path))
# Let feflo compute metric
complexity = remesher_options.get("c", None)
if complexity:
p = remesher_options.get("p", None)
if not p:
remesher_options["p"] = 10000
else:
raise ValueError("Feflo need a constraint on complexity (c) to compute the metric")
elif not prepro and self.metric is not None:
CommandLine.append("-met")
CommandLine.append(str(self.temp_sol_path))
for key, value in remesher_options.items():
CommandLine.append("-"+key)
if value:
CommandLine.append(str(value))
return " ".join(CommandLine)
def __LaunchFeflo(self, CommandLine: str):
"""Launch feflo with a subprocess using command line
Parameters
----------
CommandLine : str
the command line for feflo execution
"""
import subprocess
return subprocess.check_output(CommandLine, shell=True).decode("utf-8", errors="ignore")
def __ReadMesh(self, prefix="feflo"):
"""Read mesh file produced by feflo
Returns
-------
mesh : Mesh
feflo resulting mesh from adaption
"""
from Muscat.MeshTools.MeshTools import RefsToTags
mesh = ReadMesh(str(self.temp_out_path), ReadRefsAsField=True)
mesh = RefsToTags(mesh, self.dictNodesRefsToTags, self.dictElemRefsToTags, prefix)
return mesh
[docs]
def Prepro(self, hmin: MuscatFloat = 1e-4, geotol: MuscatFloat = 0.0001, rmax: MuscatFloat = 1):
"""Method to preprocess mesh before adaption and get a background surface mesh well oriented using INRIA Gamma3 tool feflo
Parameters
----------
hmin : MuscatFloat
Minimum mesh size allowed
geotol : MuscatFloat
Geometric tolerance to compute surface mesh
rmax : MuscatFloat
Maximum anisotropic ratio allowed
Returns
-------
backupMeshName : str
name of the file (without extension) of feflo resulting surface mesh to ensure geometry accuracy during adaptation
"""
# """
# Method to preprocess mesh before adaption and get a background surface mesh well oriented using INRIA Gamma3 tool feflo
# Args:
# hmin (float) : Name of the mesh to adapt (Default self.mesh_name)
# """
import numpy as np
hmax = (np.linalg.norm(np.max(self.mesh.nodes, axis=0) - np.min(self.mesh.nodes, axis=0)))/1000
# remesher_options = {"nordg": "", "geotol": geotol, "geoapp-allsurf-ids": "", "rmax": rmax, "hmax": hmax, "hmin": hmin, "prepro": "", "out": f"{self.temp_path}.prepro"}
remesher_options = {"nordg": "", "geotol": geotol, "geoapp-allsurf-ids": "", "prepro": "", "surf": 2, "out": f"{self.temp_path}.prepro"}
if geotol < 1e-4 or geotol > 0.1:
del remesher_options["geotol"]
CommandLine = self.__BuildCommandLine(remesher_options, prepro=True)
feflo_out = self.__LaunchFeflo(CommandLine)
with open("feflo_prepro.out", "w") as fileO:
fileO.write(feflo_out)
self.prepro_path = f"{self.temp_path}.prepro"
self.backupMeshName = f"{self.prepro_path}.back"
return self.backupMeshName
[docs]
def Runner(self, remesher_options: dict, backEndOptions: dict):
"""Run the process to write mesh in an inria mesh file,
launch feflo and read the resulting mesh
Notes
-----
The inria mesh file (called `temp_name`) will be written in a temporary directory
Parameters
----------
remesher_options : dict
keys corresponding to feflo option name
and value corresponding to value if applicable
backEndOptions: dict
Extra option for the "by file" feflo communication
"temp_name" : str
name of the output mesh file (without extension) default : "temp_file"
"binary" : bool
defaults to `True`
"prefix" : str
defaults to "feflo"
"keep_temp_files" : bool
defaults to `False`
"temp_dir" : str
directory to store temporary files
if none a Temporary Directory is used for temporary file, default None
Returns
-------
mesh : Mesh
feflo resulting mesh from adaption
"""
fefloExec = Which("feflo.a")
if fefloExec == '' or fefloExec is None:
return ("skip feflo.a is not in PATH")
if remesher_options is None:
remesher_options = {}
if backEndOptions is None:
backEndOptions = {}
temp_name = backEndOptions.get("temp_name", "temp_file")
binary = backEndOptions.get("binary", True)
prefix = backEndOptions.get("prefix", "feflo")
keep_temp_files = backEndOptions.get("keep_temp_files", False)
temp_dir = backEndOptions.get("temp_dir", None)
back = backEndOptions.get("back", False)
print_logtime = backEndOptions.get("logtime", False)
out_name = remesher_options.get("out", None)
if not temp_dir:
temp_dir = Path(TemporaryDirectory.GetTempPath())
temp_path = Path(temp_dir) / Path(temp_name)
self.__WriteMeshAndSolution(temp_path, binary)
if self.mesh.GetPointsDimensionality() == 3:
self.Prepro()
self.temp_path = self.prepro_path
if not out_name:
remesher_options["out"] = f"{temp_path}.o"
elif out_name == temp_path:
remesher_options["out"] = f"{out_name}.o"
if back:
remesher_options["back"] = self.backupMeshName
if out_name:
self.temp_out_path = Path(out_name+".mesh"+binary*"b")
if "keep-surf-ids" in remesher_options.keys():
inverse_index = {}
for key, values in self.dictElemRefsToTags.items():
for value in values:
if value not in inverse_index:
inverse_index[value] = []
inverse_index[value].append(key)
list_ref = []
tags = remesher_options["keep-surf-ids"]
for tag in tags:
list_ref.extend(inverse_index[tag])
remesher_options["keep-surf-ids"] = ','.join(map(str, list_ref))
CommandLine = self.__BuildCommandLine(remesher_options)
feflo_out = self.__LaunchFeflo(CommandLine)
if print_logtime:
print(feflo_out.split("\n")[-5].split()[1])
mesh = self.__ReadMesh(prefix)
if not keep_temp_files:
self.temp_mesh_path.unlink(missing_ok=True)
self.temp_sol_path.unlink(missing_ok=True)
self.temp_out_path.unlink(missing_ok=True)
return mesh
RegisterRemesherClass("feflo", FefloRemeshing)
RegisterRemesherClass("Feflo", FefloRemeshing)
RegisterRemesherClass("FefloByFiles", FefloRemeshing)
[docs]
def CheckIntegrity(GUI=False):
import time
import numpy as np
from Muscat.IO.GmshReader import ReadGmsh
from Muscat.TestData import GetTestDataPath
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.MeshTools.MeshModificationTools import ComputeSkin
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter, NodeFilter
tic = time.time()
fefloExec = Which("feflo.a")
if fefloExec == '' or fefloExec is None:
return ("skip")
print("Reading mesh")
filename = GetTestDataPath()+"dent3D.msh"
mesh = ReadGmsh(filename)
mesh.ConvertDataForNativeTreatment()
ComputeSkin(mesh, inPlace=True)
print(mesh)
print("Execution time:", int(time.time()-tic), "seconds")
tic = time.time()
print("Executing feflo")
obj = FefloRemeshing(mesh=mesh, solution=None)
mesh = obj.Runner({"hmin": 0.1}, backEndOptions={"binary": False})
obj = FefloRemeshing(mesh=mesh, solution=None)
mesh = obj.Runner({"hmin": 0.1}, backEndOptions={"binary": False, "temp_name": "test_feflo"})
obj = FefloRemeshing(mesh=mesh, solution=None)
mesh = obj.Runner({"hmin": 0.1}, backEndOptions={"binary": False, "temp_name": "test_feflo", "temp_dir": TemporaryDirectory.GetTempPath()})
obj = FefloRemeshing(mesh=mesh, solution=None)
mesh = obj.Runner({"hmin": 0.1}, backEndOptions={"binary": False, "temp_dir": TemporaryDirectory.GetTempPath()})
obj = FefloRemeshing(mesh=mesh, solution=np.ones(mesh.GetNumberOfNodes()))
mesh = obj.Runner({"hmin": 0.1, "c": 1000}, backEndOptions={"binary": True})
obj = FefloRemeshing(mesh=mesh, solution=np.ones(mesh.GetNumberOfNodes()))
mesh = obj.Runner({"hmin": 0.1, "c": 1000, "out": f"{TemporaryDirectory.GetTempPath()}adapted_mesh"}, backEndOptions={"binary": True})
obj = FefloRemeshing(mesh=mesh, solution=np.ones(mesh.GetNumberOfNodes()))
mesh = obj.Runner({"hmin": 0.1, "c": 1000}, backEndOptions={"binary": True, "back": True})
def AnalyticStar2D(mesh: Mesh):
siz = mesh.GetNumberOfNodes()
sol = np.zeros([siz])
for n in range(siz):
x = mesh.nodes[n][0]
y = mesh.nodes[n][1]
if x*y <= -np.pi/50:
sol[n] = 0.01*np.sin(50*x*y)
elif x*y <= 2*np.pi/50:
sol[n] = np.sin(50*x*y)
else:
sol[n] = 0.01*np.sin(50*x*y)
mesh.nodeFields["SolAtVertices0"] = sol
return sol
def Circle2D(mesh: Mesh, center=(0,0), radius=1):
siz = mesh.GetNumberOfNodes()
sol = np.zeros([siz])
for n in range(siz):
x = mesh.nodes[n][0]
y = mesh.nodes[n][1]
sol[n] = (x-center[0])**2 + (y-center[1])**2 - radius**2
mesh.nodeFields["SolAtVertices0"] = sol
return sol
mesh = ReadMesh(GetTestDataPath() + "square2D.mesh")
obj = FefloRemeshing(mesh=mesh, solution=None)
mesh = obj.Runner({"hmax": 0.1}, backEndOptions={"binary": False})
psi = AnalyticStar2D(mesh)
obj = FefloRemeshing(mesh=mesh, solution=psi)
obj.Runner({"c": 1000}, backEndOptions={"binary": False})
obj = FefloRemeshing(mesh=mesh, solution=None)
lvl = Circle2D(mesh, radius=0.5)
mesh.nodesTags.CreateTag("out", False).AddToTag(np.where(lvl > 0)[0])
ef_out = ElementFilter(elementType=ED.Triangle_3, nTag="out", nTagsTreatment="AT_LEAST_ONE")
mesh.AddElementsToTag(ef_out.GetGlobalIds(mesh), "out")
obj.Runner({"hmax": 0.01, "keep-surf-ids": ["out"]}, backEndOptions={"binary": False, "keep_temp_files": True})
print("Execution time:", int(time.time()-tic), "seconds")
return "ok"
if __name__ == '__main__':
print(CheckIntegrity())