# -*- 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 List, Union, Optional
import numpy as np
from Muscat.Types import ArrayLike, MuscatIndex
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.Filters.FilterObjects import ElementFilter, NodeFilter
from Muscat.Containers.Filters.FilterOperators import FilterOperatorBase, ComplementaryFilter, FrozenFilter, FilterLike, PartialFilter
import Muscat.Containers.ElementsDescription as ED
[docs]def ElementFilterToImplicitField(elementFilter: FilterLike, mesh: Mesh, pseudoDistance: int = 2) -> np.ndarray:
"""Function to generate an iso zero level-set on the mesh to represent
the shape of the filter. This discretized iso zero on the mesh cant always
hold a 'perfect' representation of the filter, so a modified iso zero is
created. An additional parameter pseudo-distance can be increased to create
a pseudo distance. This field is created using the connectivity and not
the real distances.
Parameters
----------
elementFilter : ElementFilter
the element filter to process
mesh : Mesh
the mesh to work on
pseudoDistance : int, optional
the number of element to propagate the pseudo-distance, by default 2
Returns
-------
np.ndarray
a field over the nodes with negative values inside the domain defined by the filter
"""
def UpdateInsideOutsideNodes(mesh: Mesh, phi, insideNodes=None, outsideNodes=None, iso=0.0):
""" function to build masks (insideNodes, outsideNodes) with the information
about if a particular nodes is connected (through an element) to the
inside phi < iso or to the outside phi > iso
"""
for name, data in mesh.elements.items():
if mesh.GetPointsDimensionality() == ED.dimensionality[name]:
phis = phi[data.connectivity]
if insideNodes is not None:
elementMask = np.any(phis < iso, axis=1)
insideNodes[data.connectivity[elementMask]] = True
if outsideNodes is not None:
elementMask = np.any(phis > iso, axis=1)
outsideNodes[data.connectivity[elementMask]] = True
phi = np.zeros(mesh.GetNumberOfNodes())
insideNodes = np.zeros(mesh.GetNumberOfNodes(), dtype=bool)
insideNodeIndices = elementFilter.GetNodesIndices(mesh)
insideNodes[insideNodeIndices] = True
dim = 0
for selection in elementFilter(mesh):
dim = max(dim, ED.dimensionality[selection.elements.elementType])
outsideNodes = np.zeros(mesh.GetNumberOfNodes(), dtype=bool)
outsideNodeIndices = GetComplementaryFilter(elementFilter).GetNodesIndices(mesh)
outsideNodes[outsideNodeIndices] = True
phi[insideNodes] = -1
phi[outsideNodes] = 1
phi[np.logical_and(insideNodes, outsideNodes)] = 0
And = np.logical_and
if dim == mesh.GetPointsDimensionality():
# correction of point values
# if a point have only zeros or negative values on neighbors point then the point is set to inside
# if a point have only zeros or positive values on neighbors point then the point is set to outside
insideNodes.fill(False)
outsideNodes.fill(False)
UpdateInsideOutsideNodes(mesh, phi, insideNodes, outsideNodes, 0)
mask = phi == 0
phi[And(mask, And(np.logical_not(outsideNodes), insideNodes))] = -1/2
phi[And(mask, And(np.logical_not(insideNodes), outsideNodes))] = 1/2
insideWork = True
outsideWork = True
for i in range(1, pseudoDistance):
if insideWork:
insideNodes.fill(False)
UpdateInsideOutsideNodes(mesh, phi, insideNodes, None, float(i))
mask = phi == i
finalMask = And(mask, np.logical_not(insideNodes))
if np.any(finalMask):
phi[finalMask] = i+1
else:
insideWork = False
if outsideWork:
outsideNodes.fill(False)
UpdateInsideOutsideNodes(mesh, phi, None, outsideNodes, float(-i))
mask = phi == -i
finalMask = And(mask, np.logical_not(outsideNodes))
if np.any(finalMask):
phi[finalMask] = -(i+1)
else:
outsideWork = False
if not outsideWork and not insideWork:
break
isoZeroElements = ElementFilter(nFilter = NodeFilter(nMask = phi == 0) , nFilterTreatment="AT_LEAST_ONE", dimensionality = mesh.GetElementsDimensionality() )
from Muscat.Containers.MeshInspectionTools import ComputeMeshMinMaxLengthScale, ExtractElementsByElementFilter
isoZeroMesh = ExtractElementsByElementFilter(mesh, isoZeroElements)
if isoZeroMesh.GetNumberOfElements():
minlength, maxlength = ComputeMeshMinMaxLengthScale(isoZeroMesh)
phi *= (minlength+maxlength)/2
else: # pragma: no cover
print(isoZeroMesh)
raise RuntimeError("Empty iso zero mesh")
return phi
[docs]def GetComplementaryFilter(filter: FilterLike) -> ComplementaryFilter:
"""Create a filter with the complementary part. if filter is already a ComplementaryFilter
the original filter (filter.filters[0]) is returned
Parameters
----------
filter : FilterLike
the filter to invert
Returns
-------
ComplementaryFilter
The filter of the Complementary part of filter
"""
# the complementary of the complementary is the original filter
if isinstance(filter, ComplementaryFilter):
return filter.filters[0]
return ComplementaryFilter(filters=[filter])
[docs]def GetFrozenFilter(filter: FilterLike, mesh: Mesh, onElements: bool = True, onNodes: bool = False, withError: bool = False) -> FrozenFilter:
"""Generate a frozen filter. This is a filter with pre-evaluated ids.
This class is useful when a repeated use of a filter is needed over the same mesh.
Parameters
----------
filter : FilterLike
the filter to evaluate
mesh : Mesh
the mesh to work on
onElements : bool, optional
if True the filter will be evaluated over the elements, by default True
onNodes : bool, optional
if True the filter will be evaluated over the nodes, by default False
withError : bool, optional
if True a exception is raised if the filter has zero elements and zero nodes , by default False
Returns
-------
FrozenFilter
a FrozenFilter instance pre-evaluated using the incoming filter and the incoming mesh
if the incoming filter is already a GetFrozenFilter over the mesh, the filter is returned (no re-evaluation)
Raises
------
RuntimeError
If the user tries to freeze a frozen filter on a different mesh
"""
if isinstance(filter, FrozenFilter):
if filter.mesh is not mesh:
raise RuntimeError("Can't freeze a FrozenFilter with a new mesh")
else:
return filter
return FrozenFilter(mesh=mesh, filters=filter, onElements=onElements, onNodes=onNodes, withError=withError)
[docs]def GetListOfPartialElementFilter(elementFilter: FilterLike, nbPartitions: Union[int,MuscatIndex], meshToFrozen=None, onElements: bool = True, onNodes: bool = False) -> Union[List[PartialFilter], List[FrozenFilter]]:
"""Generate a list of PartialElementFilter for a number of partitions
Parameters
----------
elementFilter : FilterLike
_description_
nbPartitions : int
Number of partitions
meshToFrozen : Mesh, optional
if not None the filters are frozen (evaluated at the moment of creation), no more modification is allowed, by default None
Returns
-------
List[PartialElementFilter]
a list of partialElementFilter or frozen filter. the union of all the element of the list matches the original elementFilter
"""
if meshToFrozen is not None:
fFilter = GetFrozenFilter(elementFilter, mesh=meshToFrozen, onElements=onElements, onNodes=onNodes)
res = [FrozenFilter(PartialFilter(fFilter, nbPartitions, i), mesh=meshToFrozen, onElements=onElements, onNodes=onNodes) for i in range(nbPartitions)]
else:
res = [PartialFilter(elementFilter, nbPartitions, i) for i in range(nbPartitions)]
return res
[docs]def VerifyExclusiveFilters(listOfFilters: List[FilterLike], mesh: Mesh) -> bool:
"""Function to check if a list of ElementFilter is exclusive.
(each element is present at the most in one filter)
Parameters
----------
listOfFilters : List[FilterLike]
The list of filters to check
mesh : Mesh
Mesh to evaluate the filters
Returns
-------
bool
True if the filters are exclusive
"""
for elementType, data in mesh.elements.items():
mask = np.zeros(data.GetNumberOfElements(), dtype=int)-1
for fnb, f in enumerate(listOfFilters):
ids = f.GetIdsToTreat(mesh, elementType)
if np.any(mask[ids] > -1):
idd = np.where(mask[ids] > -1)[0][0]
print(f"Filter {fnb} incompatible with filter {idd} ")
return False
mask[ids] = fnb
return True
[docs]def ListOfElementFiltersFromMask(maskVector: ArrayLike) -> List[ElementFilter]:
"""Function to construct a list of filter from a element mask vector
Parameters
----------
maskVector : ArrayLike
A element size vector of ints to determine to which filter each elements belongs to.
The number of filters is calculated using np.unique
Returns
-------
List[ElementFilter]
List of ElementFilter with the associated tags as filters
"""
ids = np.unique(maskVector.flatten())
return [ElementFilter(eMask=maskVector == partition_id) for partition_id in ids]
[docs]def FilterToETag(mesh: Mesh, elementFilter: Union[ElementFilter, FilterOperatorBase], tagname: str, append:bool=False) -> None:
"""Create an Element tag with the name tagname using the elementFilter
The tag is added the mesh.
Parameters
----------
mesh : Mesh
Mesh to work on
elementFilter : Union[ElementFilter,FilterOP]
The element filter to use to select the elements
tagname : str
the name of tag to create
append: bool, default False
if False, and the tag already exist an error is raised
if True, the selection is added, a Remove Doubles is executed at the end
if the tag already exists, an exception is raised.
"""
for selection in elementFilter(mesh):
selection.elements.tags.CreateTag(tagname, not append).AddToTag(selection.indices)
if append:
selection.elements.tags[tagname].RemoveDoubles()
[docs]def ReadElementFilter(string: str) -> ElementFilter:
"""Function to read from a string all the parameter of a ElementFiler
Parameters
----------
string : str
a string in the form of `key(names) & key(options) && ...`
possible keywords are:
- eTags -> Tags(bulk1,bulk2)
- nTags -> nTags(points1,points2)
- Dim -> Dim(3,2)
- Exprs -> Exprs(x-1,y+3, x**2 +y**2 - 5)
- nTagsOn -> to select the nTags Treatment
- ExprsOn -> to select the Exprs Treatment
Example
--------
.. code-block:: python
ReadElementFilter("Tags(Inside) & Dim(2) & nTags(Toto,Tata) & Exprs(x-1, y+3, x**2 + y**2 - 5)")
Please read the documentation of ElementFilter
Returns
-------
ElementFilter
An initialized instance of an ElementFilter.
"""
from Muscat.Helpers.ParserHelper import ReadInts
tokens = [token.strip() for token in string.split("&")]
res = ElementFilter()
if len(string) == 0:
return res
# sanity check
testString = string
for k in ["nTagsOn", "nTags", "eTags", "Dim", "ExprsOn", "Exprs"]:
if testString.count(k) > 1:
raise Exception(f"Keyword '{k}' can't be used more than once")
testString = testString.replace(k, "keyword")
for t in tokens:
# sanity check
nbOpen = t.count("(")
nbClose = t.count(")")
if nbOpen != nbClose or nbOpen != 1 or nbClose != 1:
raise Exception(f"Error parsing token {t}, missing & or ( or ) ")
keyword = t.split("(")[0]
ds = t.split("(")[1].split(")")[0]
data = [d.strip() for d in ds.split(",")]
if keyword == "eTags":
res.SetETags(data)
continue
if keyword == "nTags":
res.SetNTags(data)
continue
if keyword == "Dim":
res.SetDimensionality(ReadInts(data))
continue
if keyword == "Exprs":
from Muscat.Containers.SymExpr import SymExprWithPos
res.SetZones([SymExprWithPos(t) for t in data])
continue
if keyword == "nTagsOn":
res.SetNTagsTreatment(data[0])
continue
if keyword == "ExprsOn":
res.SetZonesTreatment(data[0])
continue
raise (Exception(f"Cant read the expression '{t}', please check your syntax"))
return res
# -------------- CheckIntegrity ---------------
[docs]def CheckIntegrityGetListOfPartialElementFilter(GUI=False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Containers.ElementsDescription import Quadrangle_4
mesh = CreateSquare(dimensions=[10, 10])
ef = ElementFilter(elementType=Quadrangle_4)
listElementFilters = GetListOfPartialElementFilter(ef, 5)
frozenListElementFilters = GetListOfPartialElementFilter(ef, 5, meshToFrozen=mesh)
return "ok"
[docs]def CheckIntegrityReadElementFilter(GUI=False) -> str:
from Muscat.Helpers.CheckTools import MustFailFunction
print(ReadElementFilter("eTags(Inside) & Dim(2) & nTags(Toto,Tata) & Exprs(x-1,y+3, x**2 +y**2 - 5) "))
print(ReadElementFilter("eTags(Inside) & Dim(3)"))
print(ReadElementFilter("Dim(3) "))
print(ReadElementFilter("Exprs(x-3) & ExprsOn(ALL_NODES) & nTagsOn(ALL_NODES) "))
print(ReadElementFilter(" eTags(toto,tata,titi)"))
assert (ReadElementFilter("eTags(toto,tata,titi)").IsEquivalent(ElementFilter(eTag=["toto", "tata", "titi"])))
assert (ReadElementFilter("").IsEquivalent(ElementFilter()))
ReadElementFilter("eTags(Inside) & Dim(3,2) ")
MustFailFunction(ReadElementFilter, "eTags(Inside) & Dim(3) & Dim(2)") # double key
MustFailFunction(ReadElementFilter, "eTags(Inside) & Dim(3 ") # missing &
MustFailFunction(ReadElementFilter, "eTags(Inside) & NonExistentToken(toto)") # missing &
return "ok"
[docs]def CheckIntegrity_FilterToETag(GUI=False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Containers.ElementsDescription import Quadrangle_4
mesh = CreateSquare(dimensions=[10, 10])
ef = ElementFilter(elementType=Quadrangle_4)
FilterToETag(mesh, ef, "quads")
FilterToETag(mesh, ef, "quads", append=True)
mesh.PrepareForOutput()
if len(mesh.elements[Quadrangle_4].tags["quads"].GetIds()) != 9*9:
raise # pragma: no cover
return "ok"
[docs]def CheckIntegrity_VerifyExclusiveFilters(GUI=False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Containers.ElementsDescription import Quadrangle_4
mesh = CreateSquare(dimensions=[10, 10])
efI = ElementFilter(elementType=Quadrangle_4)
mask = np.ones(mesh.GetNumberOfElements(), dtype=bool)
efII = ElementFilter(eMask=mask)
if VerifyExclusiveFilters([efI, efII], mesh) == 1:
raise # pragma: no cover
return "ok"
[docs]def CheckIntegrity_ListOfElementFiltersFromMask(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
mesh = CreateSquare(dimensions=[10, 10])
mask = np.asarray(np.random.random(mesh.GetNumberOfElements())*5, dtype=int)
listOfFilters = ListOfElementFiltersFromMask(mask)
if VerifyExclusiveFilters(listOfFilters, mesh) == 0:
raise # pragma: no cover
tmp = TemporaryDirectory.GetTempPath()
nbElements = 0
for i, ff in enumerate(listOfFilters):
new_mesh = ExtractElementsByElementFilter(mesh, ff)
nbElements += new_mesh.GetNumberOfElements()
print(f"Number of element in the mesh {i}:" + str(new_mesh.GetNumberOfElements()))
WriteMeshToXdmf(tmp + f"CI_ListOfElementFiltersFromMask{i}.xdmf", new_mesh)
WriteMeshToXdmf(tmp + "CI_ListOfElementFiltersFromMask_base.xdmf", mesh, CellFields=[mask], CellFieldsNames=["RANK"])
print((nbElements, mesh.GetNumberOfElements()))
if nbElements != mesh.GetNumberOfElements():
return "Not OK" # pragma: no cover
return "OK"
[docs]def CheckIntegrityElementFilterToImplicitField(GUI: bool == False):
from Muscat.Containers import MeshCreationTools as UMCT
mesh = UMCT.CreateSquare(dimensions=[10, 10], origin=[0, 0], spacing=[10/9, 10/9], ofTriangles=True)
ef = ElementFilter(zone=lambda xyz: xyz[:, 0]-5.5)
field = ElementFilterToImplicitField(ef, mesh, pseudoDistance=10)
return "ok"
[docs]def CheckIntegrityGetFrozenFilter(GUI=False) -> str:
from Muscat.Helpers.CheckTools import MustFailFunction
from Muscat.Containers import MeshCreationTools as UMCT
mesh2 = UMCT.CreateSquare(dimensions=[10, 10], origin=[0, 0], spacing=[10/9, 10/9], ofTriangles=True)
mesh = UMCT.CreateSquare(dimensions=[10, 10], origin=[0, 0], spacing=[10/9, 10/9], ofTriangles=True)
ef = ElementFilter(zone=lambda xyz: xyz[:, 0]-0.5)
GetFrozenFilter(GetFrozenFilter(ef, mesh), mesh)
MustFailFunction(GetFrozenFilter, GetFrozenFilter(ef, mesh), mesh2)
return "ok"
[docs]def CheckIntegrityGetComplementaryFilter(GUI=False) -> str:
from Muscat.Containers import MeshCreationTools as UMCT
mesh = UMCT.CreateSquare(dimensions=[10, 10], origin=[0, 0], spacing=[10/9, 10/9], ofTriangles=True)
ef = ElementFilter(zone=lambda xyz: xyz[:, 0]-0.5)
assert (GetComplementaryFilter(GetComplementaryFilter(ef)) is ef)
return "ok"
[docs]def CheckIntegrity(GUI: bool = False) -> str:
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
toTest = [
CheckIntegrityElementFilterToImplicitField,
CheckIntegrityGetFrozenFilter,
CheckIntegrityGetComplementaryFilter,
CheckIntegrityReadElementFilter,
CheckIntegrity_FilterToETag,
CheckIntegrity_VerifyExclusiveFilters,
CheckIntegrity_ListOfElementFiltersFromETagList,
CheckIntegrity_ListOfElementFiltersFromMask,
CheckIntegrityGetListOfPartialElementFilter,
]
return RunListOfCheckIntegrities(toTest, GUI)
if __name__ == '__main__':
print(CheckIntegrity(GUI=True)) # pragma: no cover