# -*- 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 Dict, List, Optional, Tuple
import numpy as np
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.FE.Fields.FEField import FEField
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Containers.Filters.FilterTools import GetFrozenFilter
from Muscat.Containers.Mesh import ElementsContainer, Mesh
from Muscat.Containers.MeshModificationTools import CleanLonelyNodes
from Muscat.Types import MuscatIndex
from collections import defaultdict
import itertools
[docs]def GetDataOverALine(startPoint: ArrayLike, stopPoint: ArrayLike, nbPoints: int, mesh: Mesh, fields: List, method: Optional[str] = None) -> Mesh:
"""Compute the values of a mesh/field over a line
startPoint and stopPoint are used to construct a line with nbPoints points.
the data in mesh.nodeFields, mesh.elemFields and fields are evaluated at every point of the line
Parameters
----------
startPoint : ArrayLike
Start point. this is passed to the CreateUniformMeshOfBars function
stopPoint : ArrayLike
Stop point. this is passed to the CreateUniformMeshOfBars function
nbPoints : int
this argument is passed to the CreateUniformMeshOfBars function
mesh : Mesh
the mesh to extract mesh.nodeFields, mesh.elemFields
fields : List
a list containing only FEFields
method : Optional[str], optional
this argument is passed to the GetFieldTransferOp function, by default None
Returns
-------
Mesh
a unstructured mesh with nodeFields populated with the data
nodeFields with contain a dictionary for every field name and the numpy vector with
the values at the points. If two or more fields have the same name only one is recovered,
the priority order staring with the hightest priority are : fields, elemFields and nodeFields
Raises
------
Exception
In the case an IPField is present in the 'fields' argument
"""
res: Dict[str, np.ndarray] = {}
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
lineMesh = CreateUniformMeshOfBars(startPoint, stopPoint, nbPoints)
from Muscat.Containers.MeshFieldOperations import GetFieldTransferOp
from Muscat.FE.Fields.FieldTools import NodeFieldToFEField, ElemFieldsToFEField
# transfer of nodal fields
nodeFieldsAsFEFields = NodeFieldToFEField(mesh)
if len(nodeFieldsAsFEFields) != 0:
firstField = next(iter(nodeFieldsAsFEFields.values()))
op, status, entities = GetFieldTransferOp(firstField, lineMesh.nodes, method=method)
for name, data in mesh.nodeFields.items():
res[name] = op.dot(data)
# transfer of cell data
cellFieldsAsFEFields = ElemFieldsToFEField(mesh)
if len(cellFieldsAsFEFields) != 0:
firstField = next(iter(cellFieldsAsFEFields.values()))
op, status, entities = GetFieldTransferOp(firstField, lineMesh.nodes, method=method)
for name, data in mesh.elemFields.items():
res[name] = op.dot(data)
# transfer of FEFields
for efField in fields:
if isinstance(efField, FEField):
op, status, status = GetFieldTransferOp(efField, lineMesh.nodes, method=method)
res[efField.name] = op.dot(efField.data)
else: # pragma: no cover
raise Exception(f"Don't know how to treat field of type ({type(efField)})")
lineMesh.nodeFields = res
return lineMesh
[docs]def GetElementsFractionInside(field: ArrayLike, mesh: Mesh, ef: Optional[ElementFilter] = None) -> np.ndarray:
"""Compute the volume fraction for each element in ids, based in the level-set field field.
Function valid for simplices (bars, triangles, tetrahedrons)
Parameters
----------
field : ArrayLike
the levelset field
mesh : Mesh
the support mesh of the point field field
ef : Optional[ElementFilter], optional
Element Filter to select elements, by default None
Returns
-------
np.ndarray
a numpy array of size number of element selected by the element filter
"""
if ef is None:
ef = ElementFilter()
frozenEF = GetFrozenFilter(ef, mesh)
res = np.empty(frozenEF.GetElementSelectionSize(), dtype=MuscatFloat)
for selection in frozenEF(mesh):
res[selection.GetSelectionSlice()] = GetElementsFractionInsideLegacy(field, mesh.nodes, selection.elementType, selection.elements, selection.indices)
return res
[docs]def GetElementsFractionInsideLegacy(field: ArrayLike, points: ArrayLike, name: ED.ElementType, elements: ElementsContainer, ids: np.ndarray) -> np.ndarray:
"""Compute the volume fraction for each element in ids, based in the level-set field field.
Function valid for simplices (bars, triangles, tetrahedrons)
Parameters
----------
field : ArrayLike
the levelset field
points : ArrayLike
the points of the mesh
name : ED.ElementType
the element name
elements : ElementsContainer
the element data
ids : np.ndarray
the ids to treat
Returns
-------
np.ndarray
the volume fraction for each element in ids
"""
from scipy.spatial import Delaunay
def triangle_lob_fraction(index, phis):
div = phis-phis[index]
div[index] = phis[index]
val = phis[index]/(div)
return np.prod(val)
def extract_signs(phis):
signs = np.sign(phis)
dominant_sign = np.sign(np.sum(signs))
if dominant_sign == -1:
signs[signs != 1] = -1
else:
signs[signs != -1] = 1
return signs
def split_signs(signs):
minuses = np.nonzero(signs == -1)[0]
pluses = np.nonzero(signs == 1)[0]
return minuses, pluses
# helper function for tets
def tetrahedron_volumic_fraction(phis):
assert phis.size == 4
assert np.count_nonzero(phis) > 0
minuses, pluses = sort_by_sign(phis)
if minuses.size == 0:
return 0.0
elif minuses.size == 1:
return tetrahedron_volumic_fraction_dominated(minuses, pluses)
elif minuses.size == 2:
return tetrahedron_volumic_fraction_non_dominated(minuses, pluses)
elif minuses.size == 3:
return 1.0 - tetrahedron_volumic_fraction_dominated(pluses, minuses)
else:
return 1.0
def tetrahedron_volumic_fraction_dominated(phi_in, phis_out):
return np.prod(phi_in / (phi_in - phis_out))
def tetrahedron_volumic_fraction_non_dominated(phis_in, phis_out):
phis_in_r = np.reshape(phis_in, (2, 1))
phis_out_r = np.reshape(phis_out, (1, 2))
ratios = phis_in_r / (phis_in_r - phis_out_r)
result = ratios[0, 0] * ratios[0, 1] * (1.0 - ratios[1, 0]) + \
ratios[1, 0] * ratios[1, 1] * (1.0 - ratios[0, 1]) + \
ratios[1, 0] * ratios[0, 1]
return result
def sort_by_sign(phis):
signs = np.sign(phis)
dominant_sign = np.sign(np.sum(signs))
if dominant_sign == -1:
minuses = phis[signs != 1]
pluses = phis[signs == 1]
else:
minuses = phis[signs == -1]
pluses = phis[signs != -1]
return minuses, pluses
# -----------------------------------------------------------------------------
res = np.zeros(len(ids)) - 1
for cpt, id in enumerate(ids):
coon = elements.connectivity[id, :]
vals = field[coon]
if np.all(vals < 0):
res[cpt] = 1
elif np.all(vals > 0):
res[cpt] = 0
else:
# elemPoints = points[coon,:]
# simplex = Delaunay(points).simplices
if ED.dimensionality[name] == 1 and len(vals) == 2:
# bar2
if vals[0] < 0:
res[cpt] = abs(vals[0])/np.sup(abs(vals))
else:
res[cpt] = 1-abs(vals[0])/np.sup(abs(vals))
elif ED.dimensionality[name] == 2 and len(vals) == 3:
# triangles
signs = extract_signs(vals)
minuses, pluses = split_signs(signs)
if minuses.size == 2:
res[cpt] = 1.0 - triangle_lob_fraction(pluses[0], vals)
elif minuses.size == 1:
res[cpt] = triangle_lob_fraction(minuses[0], vals)
elif ED.dimensionality[name] == 3 and len(vals) == 4:
# tets
res[cpt] = tetrahedron_volumic_fraction(vals)
return res
[docs]def GetVolumePerElement(inmesh: Mesh, elementFilter: Optional[ElementFilter] = None) -> np.ndarray:
"""Compute the volume (surface for 2D element and length for 1D elements) for each element selected by the elementFilter
Parameters
----------
inmesh : Mesh
the mesh to extract elements
elementFilter : Optional[ElementFilter], optional
filter to select some elements, if None the volume of all the element are computed
Returns
-------
np.ndarray
a numpy array of size number of "element selected by the elementFilter" with the volume
"""
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.SymWeakForm import GetField
from Muscat.FE.SymWeakForm import GetTestField
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Integration import IntegrateGeneral
numbering = ComputeDofNumbering(inmesh, LagrangeSpaceP0, elementFilter=elementFilter)
wform = GetField("F", 1).T*GetTestField("T", 1)
F = FEField("F", inmesh, LagrangeSpaceP0, numbering)
F.Allocate(1.)
unknownFields = [FEField("T", mesh=inmesh, space=LagrangeSpaceP0, numbering=numbering)]
_, f = IntegrateGeneral(mesh=inmesh, wform=wform, constants={}, fields=[F], unknownFields=unknownFields, elementFilter=elementFilter)
return f
[docs]def GetMeasure(inmesh: Mesh, elementFilter: Optional[ElementFilter] = None):
"""Compute the measure of the mesh (volume in 3D, surface in 2D, ...)
Only elements of the higher dimensionality are taken into account if no element filter
is supplied
Attention: if the element filter contain dimensionality mixed elements the result is the
sum of the measure (volume + surface for example)
Parameters
----------
inmesh : Mesh
the mesh to use for the computation
elementFilter : ElementFilter
element to take into account in the computation of the measure, if none
all the element of the higher dimensionality are used.
Returns
-------
MuscatFloat
the volume if the mesh contains 3D elements
the surface if the mesh contains 2D elements and no 3D elements
the length if the mesh contains 1D elements and no 3D elements nor 2D elements
"""
if elementFilter is None:
elementFilter = ElementFilter(dimensionality=inmesh.GetElementsDimensionality())
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo, ConstantSpaceGlobal
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.SymWeakForm import GetField
from Muscat.FE.SymWeakForm import GetTestField
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Integration import IntegrateGeneral
numbering = ComputeDofNumbering(inmesh, LagrangeSpaceGeo, fromConnectivity=True)
wform = GetField("F", 1).T*GetTestField("T", 1)
F = FEField("F", inmesh, LagrangeSpaceGeo, numbering)
F.Allocate(1.)
gNumbering = ComputeDofNumbering(inmesh, ConstantSpaceGlobal)
unknownFields = [FEField("T", mesh=inmesh, space=ConstantSpaceGlobal, numbering=gNumbering)]
_, f = IntegrateGeneral(mesh=inmesh, wform=wform, constants={}, fields=[F], unknownFields=unknownFields, elementFilter=elementFilter)
return f[0]
[docs]def GetVolume(inmesh: Mesh) -> MuscatFloat:
"""Compute the volume of the mesh
Only element of the bigger dimensionality are taken into account
if no elements on mesh we return zero.
for a more general function see: GetMeasure
Parameters
----------
inmesh : Mesh
the mesh to use for the computation
Returns
-------
MuscatFloat
the volume if the mesh contains 3D elements
the surface if the mesh contains 2D elements and no 3D elements
the length if the mesh contains 1D elements and no 3D elements nor 2D elements
"""
return GetMeasure(inmesh)
[docs]def ComputeNodeToElementConnectivity(inmesh: Mesh, maxNumConnections: int = 200) -> Tuple[np.ndarray, np.ndarray]:
"""Compute the node to connectivity connection for every point
Parameters
----------
inmesh : Mesh
the mesh
maxNumConnections : int, optional
the number of potential connection to a single node, normally 200 is more than enough, by default 200
Returns
-------
Tuple
numpy.ndarray, (nb point, max nb connections) item, the connectivity for each node -> element
numpy.array, the number of connections per point
"""
# generation of the dual graph
dualGraph = np.zeros((inmesh.GetNumberOfNodes(), maxNumConnections), dtype=MuscatIndex)-1
usedPoints = np.zeros(inmesh.GetNumberOfNodes(), dtype=MuscatIndex)
cpt = 0
for name, elements in inmesh.elements.items():
for i in range(elements.GetNumberOfElements()):
coon = elements.connectivity[i, :]
for j in coon:
dualGraph[j, usedPoints[j]] = cpt
usedPoints[j] += 1
cpt += 1
# we crop the output data
maxsize = np.max(np.sum(dualGraph >= 0, axis=1))
dualGraph = dualGraph[:, 0:maxsize]
return dualGraph, usedPoints
[docs]def ComputeNodeToNodeConnectivity(inmesh: Mesh, maxNumConnections=200) -> Tuple[np.ndarray, np.ndarray]:
"""Compute the node to node connection for every node in the mesh
Parameters
----------
inmesh : Mesh
the mesh
maxNumConnections : int, optional
the number of potential connection to a single node, normally 200 is more than enough, by default 200
Returns
-------
Tuple
numpy.ndarray, (nb point, max nb connections) item, the connectivity for each node -> node
numpy.array, the number of connections per point
"""
# generation of the dual graph
dualGraph = np.zeros((inmesh.GetNumberOfNodes(), maxNumConnections), dtype=MuscatIndex)-1
usedPoints = np.zeros(inmesh.GetNumberOfNodes(), dtype=MuscatIndex)
for name, elements in inmesh.elements.items():
size = elements.GetNumberOfNodesPerElement()
for i in range(elements.GetNumberOfElements()):
coon = elements.connectivity[i, :]
for j in range(size):
myIndex = coon[j]
for k in range(size):
if k == j:
continue
dualGraph[myIndex, usedPoints[myIndex]] = coon[k]
usedPoints[myIndex] += 1
# we reached the maximum number of connection
# normally we have some data duplicated
# we try to shrink the vector
if usedPoints[myIndex] == maxNumConnections: # pragma: no cover
# normally we shouldn't pas here
c = np.unique(dualGraph[myIndex, :])
dualGraph[myIndex, 0:len(c)] = c
usedPoints[myIndex] = len(c)
maxsize = 0
# we finish now we try to compact the structure
for i in range(inmesh.GetNumberOfNodes()):
c = np.unique(dualGraph[i, 0:usedPoints[i]])
dualGraph[i, 0:len(c)] = c
dualGraph[i, len(c):] = -1
usedPoints[i] = len(c)
maxsize = max(len(c), maxsize)
# we crop the output data
dualGraph = dualGraph[:, 0:maxsize]
return dualGraph, usedPoints
[docs]def ComputeElementToElementConnectivity(inmesh: Mesh, dimensionality: int = None, connectivityDimension: int = None, maxNumConnections: int = 200) -> Tuple[np.ndarray, np.ndarray]:
"""Generates the elements to element connectivity of an Mesh in the following sense.
An element is linked to another in the graph if they share:
- a face if connectivityDimension = 2
- an edge if connectivityDimension = 1
- a vertex if connectivityDimension = 0
(if connectivityDimension is initialized to None, it will be set to dimensionality - 1)
Also provides the array of number of connection per cell.
Parameters
----------
inmesh : Mesh
the input mesh
dimensionality : int, optional
dimensionality filter, by default inmesh.GetPointsDimensionality()
connectivityDimension : int, optional
type of connectivity, by default dimensionality - 1
maxNumConnections : int, optional
the number of potential connection to a single node, normally 200 is more than enough, by default 200
Returns
-------
Tuple
numpy.ndarray, (nb elements, max nb connections) item, the connections for each element -> element
numpy.array, the number of connections per elements
"""
if dimensionality == None:
dimensionality = inmesh.GetPointsDimensionality()
if connectivityDimension == None:
connectivityDimension = dimensionality - 1
assert dimensionality > connectivityDimension
elFilter = ElementFilter(dimensionality=dimensionality)
totalCpt={}
for selection in elFilter(inmesh):
ids = selection.indices
diff_dimension=dimensionality - connectivityDimension
# recovering the faces according to dimensionality discrepancy
faces=[ED.faces1,ED.faces2,ED.faces3 ][diff_dimension-1][selection.elementType]
nbElements = selection.elements.GetNumberOfElements()
for faceType,localFaceConnectivity in faces:
cpt = totalCpt.setdefault(faceType,0) + nbElements
totalCpt[faceType] = cpt
# Allocation of matrices to store all interface elements
storage = {}
for k,v in totalCpt.items():
nn = ED.numberOfNodes[k]
# storage will contain a tuple containing the connectivity of the face and the globalid of the face it relates to
storage[k] = (np.empty((v,nn),dtype=int),np.empty((v),dtype=int))
totalCpt[k] = 0
# fill storage with elements of dimensionality D
for selection in elFilter(inmesh):
elemName = selection.elements.elementType
data = selection.elements
ids = selection.indices
# recovering the faces according to dimmensionnality discrepancy
diff_dimension=dimensionality - connectivityDimension
faces=[ED.faces1,ED.faces2,ED.faces3][diff_dimension-1][selection.elementType]
nbElements = data.GetNumberOfElements()
for faceType,localFaceConnectivity in faces:
globalFaceConnectivity = data.connectivity[:,localFaceConnectivity]
cpt = totalCpt[faceType]
# Storing global connectivity in storage
storage[faceType][0][cpt:cpt+nbElements,:] = globalFaceConnectivity
# Storing the ids of the element corresponding to the connectivity above
storage[faceType][1][cpt:cpt+nbElements] = ids + selection.meshOffset
totalCpt[faceType] += nbElements
cpt_graph = 0
# dictionnary to relate faces ids (keys) and elements ids that are connected to this face (values: list of elements ids)
faceToElements=defaultdict(lambda: [])
for elemName,cpt in totalCpt.items():
if cpt == 0:
continue
store,idElems = storage[elemName]
_,index,inverse,counts = np.unique(np.sort(store,axis=1),return_index=True,return_inverse=True,return_counts=True,axis=0)
inflatedIndexes = index[inverse] # This array of the same size than store has same values on indexes that are the same elements
inflatedCounts = counts[inverse] # will be used to filter out skin elements
# Filtering out skin elements and making indexes global
for globalIdElem,globalIdFace in zip(idElems[inflatedCounts>1],inflatedIndexes[inflatedCounts>1]+cpt_graph):
# assigning the global id of the element to the current face
faceToElements[globalIdFace].append(globalIdElem)
cpt_graph+=cpt
# transforming the faceToElements dictionnary into the original elemGraph format
elemGraph = np.zeros((inmesh.GetNumberOfElements(), maxNumConnections), dtype=MuscatIndex) - 1
elemCounts= np.zeros(inmesh.GetNumberOfElements(),dtype=MuscatIndex)
# Looping on all faces (they are not skin faces since we filtered them out previously)
for connectedElements in faceToElements.values():
# for all couple of elements connected to the same face we add both of them to the element graph connectivity array
for id_elem1,id_elem2 in itertools.combinations(connectedElements,2):
elemGraph[id_elem1,elemCounts[id_elem1]]=id_elem2
elemCounts[id_elem1]+=1
elemGraph[id_elem2,elemCounts[id_elem2]]=id_elem1
elemCounts[id_elem2]+=1
maxsize = np.max(elemCounts) # crop output data
elemGraph = elemGraph[:, 0:maxsize]
return elemGraph, elemCounts
[docs]def ComputeMeshDensityAtNodes(mesh: Mesh) -> np.ndarray:
"""Function to compute the mesh size at each point
This is done using the length of the bars of each element and averaging the data at nodes.
All the element are treated, even the lower dimensionality element (not the 0D elements of course).
This means the triangles in the surface can affect the values of the nodes on the surface
No weighting is done.
Parameters
----------
mesh : Mesh
The input mesh
Returns
-------
np.ndarray
a numpy array (node field) with average of the mesh size at each point
"""
nbNodes = mesh.GetNumberOfNodes()
result = np.zeros(nbNodes, dtype=MuscatFloat)
cpt = np.zeros(nbNodes, dtype=MuscatIndex)
posOfNodes = mesh.GetPosOfNodes()
def AddBarSizesToOutput(connectivityOfBars):
dx = posOfNodes[connectivityOfBars[:, 0], 0] - posOfNodes[connectivityOfBars[:, 1], 0]
dy = posOfNodes[connectivityOfBars[:, 0], 1] - posOfNodes[connectivityOfBars[:, 1], 1]
if posOfNodes.shape[1] == 3:
dz = posOfNodes[connectivityOfBars[:, 0], 2] - posOfNodes[connectivityOfBars[:, 1], 2]
else:
dz = 0
lengths = np.sqrt(dx**2+dy**2+dz**2)
for i in range(2):
result[connectivityOfBars[:, i]] += lengths
cpt[connectivityOfBars[:, i]] += 1
for selection in ElementFilter(dimensionality=3)(mesh=mesh):
faces = ED.faces2[selection.elementType]
for fname, facesConnectivity in faces:
AddBarSizesToOutput(selection.elements.connectivity[:, facesConnectivity][:, 0:2])
for selection in ElementFilter(dimensionality=2)(mesh=mesh):
faces = ED.faces[selection.elementType]
for fname, facesConnectivity in faces:
AddBarSizesToOutput(selection.elements.connectivity[:, facesConnectivity][:, 0:2])
for selection in ElementFilter(dimensionality=1)(mesh=mesh):
AddBarSizesToOutput(selection.elements.connectivity[:, 0:2])
cpt[cpt == 0] = 1
result /= cpt
return result
[docs]def EnsureUniquenessElements(mesh: Mesh) -> None:
"""Ensure that every element if present only once on the mesh
Parameters
----------
mesh : Mesh
input mesh
Raises
------
Exception
if 2 element (event if the connectivity is permuted) a exception is raised
"""
cpt = 0
for name, data in mesh.elements.items():
dd = dict()
for el in range(data.GetNumberOfElements()):
n = tuple(np.sort(data.connectivity[el, :]))
if len(n) != len(np.unique(n)):
raise Exception("("+str(name)+") element " + str(el) + " (global["+str(el+cpt)+"])" + " use a point more than once ("+str(data.connectivity[el, :])+")")
if n in dd.keys():
raise Exception("("+str(name)+") element " + str(el) + " (global["+str(el+cpt)+"])" + " is a duplication of element " + str(dd[n]) + " (global["+str(dd[n]+cpt)+"])")
dd[n] = el
cpt = data.GetNumberOfElements()
[docs]def MeshQualityAspectRatioBeta(mesh: Mesh):
"""experimental mesh quality only available for tets
Parameters
----------
mesh : Mesh
the input mesh
Raises
------
Exception
raise if the quality is larger than 1000
"""
# https://cubit.sandia.gov/public/15.2/help_manual/WebHelp/mesh_generation/mesh_quality_assessment/tetrahedral_metrics.htm
for name, data in mesh.elements.items():
if ED.dimensionality[name] == 0:
pass
elif ED.dimensionality[name] == 1:
pass
elif ED.dimensionality[name] == 2:
pass
elif ED.dimensionality[name] == 3:
if name == ED.Tetrahedron_4:
mMax = 0
for el in range(data.GetNumberOfElements()):
n = data.connectivity[el, :]
nodes = mesh.nodes[n, :]
p0 = nodes[0, :]
p1 = nodes[1, :]
p2 = nodes[2, :]
p3 = nodes[3, :]
a = np.linalg.norm(p1-p0)
b = np.linalg.norm(p2-p0)
normal = np.cross(p1-p0, p2-p0)
# https://math.stackexchange.com/questions/128991/how-to-calculate-area-of-3d-triangle
base_area = 0.5*a*b*np.sqrt(1-(np.dot(p1-p0, p2-p0)/(a*b))**2)
normal /= np.linalg.norm(normal)
# distance to the opposite point
d = np.dot(normal, p3-p0)
volume = base_area*d
inscribed_sphere_radius = d/3
# https://math.stackexchange.com/questions/2820212/circumradius-of-a-tetrahedron
c = np.linalg.norm(p3-p0)
A = np.linalg.norm(p3-p2)
B = np.linalg.norm(p1-p3)
C = np.linalg.norm(p2-p1)
circumRadius_sphere_radius = np.sqrt((a*A+b*B+c*C) *
(a*A+b*B-c*C) *
(a*A-b*B+c*C) *
(-a*A+b*B+c*C))/(24*volume)
AspectRatioBeta = circumRadius_sphere_radius/(3.0 * inscribed_sphere_radius)
mMax = max([mMax, AspectRatioBeta])
if AspectRatioBeta > 1000:
raise Exception("Element " + str(el) + " has quality of " + str(AspectRatioBeta))
else:
raise
[docs]def ComputeMeshMinMaxLengthScale(mesh) -> Tuple[MuscatFloat, MuscatFloat]:
"""Compute a estimation of the minimal and maximal length scale of the elements
for this we compute the centroid of the element and compute the min and max distance
between the centroid and each node of the mesh.
Then return a tuple with (2*min,2*max) of all the distances on the mesh
Parameters
----------
mesh : Mesh
The input mesh
Returns
-------
Tuple[MuscatFloat,MuscatFloat]
tuple with (2*min_dist,2*max_dist) for all the elements in the mesh
"""
(resMin, resMax) = (None, None)
for name, data in mesh.elements.items():
if data.GetNumberOfNodesPerElement() < 2:
continue
if data.GetNumberOfElements() == 0:
continue
posX = mesh.nodes[data.connectivity, 0]
posY = mesh.nodes[data.connectivity, 1]
meanX = np.sum(posX, axis=1)/data.GetNumberOfNodesPerElement()
meanY = np.sum(posY, axis=1)/data.GetNumberOfNodesPerElement()
meanX.shape = (len(meanX), 1)
meanY.shape = (len(meanX), 1)
distToBarycenter2 = (posX-meanX)**2 + (posY-meanY)**2
if mesh.nodes.shape[1] == 3:
posZ = mesh.nodes[data.connectivity, 2]
meanZ = np.sum(posZ, axis=1)/data.GetNumberOfNodesPerElement()
meanZ.shape = (len(meanX), 1)
distToBarycenter2 += (posZ-meanZ)**2
distToBarycenter = np.sqrt(distToBarycenter2)
mMin = np.min(distToBarycenter)
mMax = np.max(distToBarycenter)
if resMin is None:
resMin = mMin
else:
resMin = min(resMin, mMin)
if resMax is None:
resMax = mMax
else:
resMax = min(resMax, mMax)
return (2*resMin, 2*resMax)
[docs]def ComputeNormalsAtElements(mesh: Mesh, ef:Optional[ElementFilter] = None) -> np.ndarray:
"""Compute an unique normal per element using the first 3 nodes (for 2D element in 3D)
or the first 2 nodes (for 1D elements in 2D)
Parameters
----------
mesh : Mesh
Mesh to work on
ef : ElementFilter, optional
Element Filter to select element to work on, by default None
Returns
-------
np.ndarray
A numpy array of size (mesh.GetNumberOfElements(), mesh.GetPointsDimensionality()), with the normal (in 3D or 2D) for the
selected elements
"""
if ef is None:
ef = ElementFilter(dimensionality= mesh.GetPointsDimensionality()-1)
normals = np.zeros((mesh.GetNumberOfElements(), mesh.GetPointsDimensionality()), dtype=MuscatFloat )
for selection in ef(mesh):
if ED.dimensionality[selection.elementType] in [0,3]:
continue
elif ED.dimensionality[selection.elementType] == 2:
coords = mesh.nodes[selection.elements.connectivity[selection.indices,:]]
localNormals = np.cross(
coords[..., 1, :] - coords[..., 0, :],
coords[..., 2, :] - coords[..., 0, :])
normals[selection.meshOffset:selection.meshOffset+selection.Size(),:] = localNormals
elif ED.dimensionality[selection.elementType] == 1:
coords = mesh.nodes[selection.elements.connectivity[selection.indices,:]]
axial = coords[..., 1, :] - coords[..., 0, :]
normals[selection.meshOffset:selection.meshOffset+selection.Size(),0] = axial[:, 1]
normals[selection.meshOffset:selection.meshOffset+selection.Size(),1] = -axial[:, 0]
else:
raise
norm = np.linalg.norm(normals, axis=1)
norm[norm == 0] = 1
normals /= norm[:,None]
return normals
# ------------------------- CheckIntegrity ------------------------
[docs]def CheckIntegrity_ComputeNormalsAtElements(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles, CreateMeshOf
mesh = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
normals = ComputeNormalsAtElements(mesh, ElementFilter(dimensionality = 2) )
print(normals)
np.testing.assert_equal(normals[0,:],[0.,0.,1.])
np.testing.assert_equal(normals[1,:],[1.,0.,0.])
mesh2D = CreateMeshOf([[0, 0], [1, 0], [0, 1]], [[0, 1], [1, 2]], ED.Bar_2)
normals = ComputeNormalsAtElements(mesh2D)
np.testing.assert_equal(normals[0,:],[0.,-1.])
np.testing.assert_almost_equal(normals[1,:],[np.sqrt(2)/2,np.sqrt(2)/2])
return "ok"
[docs]def CheckIntegrity_GetVolume(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf, CreateConstantRectilinearMesh, CreateSquare
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tets = [[0, 1, 2, 3],]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
vol = GetVolume(mesh)
if vol != (1./6.):
raise Exception('Error en the calculation of the volume') # pragma: no cover
vol1elem = GetVolumePerElement(mesh, ElementFilter(elementType=ED.Tetrahedron_4))
print(vol1elem)
if sum(vol1elem) != vol:
raise Exception("incompatible solution of GetVolume vs GetVolumePerElement")
myMesh = CreateConstantRectilinearMesh(dimensions=[3, 3, 3], spacing=[0.5, 0.5, 0.5], origin=[0, 0, 0])
vol = GetVolume(myMesh)
if abs(vol-1.) > 1e-8:
raise Exception('Error en the calculation of the volume') # pragma: no cover
squareMesh = CreateSquare(dimensions= [2, 2], origin = [-1.0, -1.0], spacing = [1.0, 1.0], ofTriangles = True)
from Muscat.Containers.MeshModificationTools import ComputeSkin
ComputeSkin(squareMesh, inPlace=True)
print("print volume of triangles")
vol1elem = GetVolumePerElement(squareMesh, ElementFilter(elementType=ED.Triangle_3))
print(vol1elem)
return "ok"
[docs]def CheckIntegrity_EnsureUniquenessElements(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tets = [[0, 1, 2, 3]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
# this function must work
EnsureUniquenessElements(mesh)
tets = [[0, 1, 2, 3], [3, 1, 0, 2]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
# This function must not work
try:
EnsureUniquenessElements(mesh)
return "No ok"
except:
pass
return "ok"
[docs]def CheckIntegrity_GetElementsFractionInside(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles, CreateMeshOf
meshTriangles = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
field = np.array([-1, -1, -1, 1])
for name, elements in meshTriangles.elements.items():
res = GetElementsFractionInsideLegacy(field, meshTriangles.nodes, name, elements, range(elements.GetNumberOfElements()))
print(res)
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]]
tets = [[0, 1, 2, 3],]
mesh3D = CreateMeshOf(points, tets, ED.Tetrahedron_4)
for name, elements in mesh3D.elements.items():
res = GetElementsFractionInsideLegacy(field, mesh3D.nodes, name, elements, range(elements.GetNumberOfElements()))
print(res)
res = GetElementsFractionInside(field, mesh3D, ElementFilter(dimensionality=3))
return "ok"
[docs]def CheckIntegrity_ComputeNodeToNodeConnectivity(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles
res = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
dg, unUsed = ComputeNodeToNodeConnectivity(res)
return "ok"
[docs]def CheckIntegrity_ComputeMeshMinMaxLengthScale(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tets = [[0, 1, 2, 3], [3, 1, 0, 2]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
print(ComputeMeshMinMaxLengthScale(mesh))
PrintMeshInformation(mesh)
return "ok"
[docs]def CheckIntegrity_MeshQualityAspectRatioBeta(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateCube
mesh = CreateCube(dimensions=[20, 20, 20], origin=[0.1, 0.1, 0.,], spacing=[0.9/19]*3)
MeshQualityAspectRatioBeta(mesh)
return "ok"
[docs]def CheckIntegrity_GetDataOverALine(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateCube
mesh = CreateCube(dimensions=[20, 20, 20], origin=[0.1, 0.1, 0.,], spacing=[0.9/19]*3)
mesh.nodeFields["xpos"] = mesh.nodes[:, 0]
from Muscat.FE.Fields.FieldTools import CreateFieldFromDescription
field = CreateFieldFromDescription(mesh, [(ElementFilter(zone=lambda x: x[:, 2] > 0.5), 1)])
field.name = "biMaterial"
mesh.elemFields["biMaterialAtElem"] = field.GetCellRepresentation()
print("---")
print(mesh.elemFields["biMaterialAtElem"])
print("---")
res = GetDataOverALine([0, 0, 0], [1, 1, 1], 100, mesh, [field])
print(res)
if res.GetNumberOfNodes() != 100:
raise Exception("Wrong number of point in the line")
if len(res.nodeFields) != 3:
raise Exception("Wrong number of fields")
if GUI:
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh, filename="CheckIntegrity_GetDataOverALine_bulk.xdmf", run=False)
OpenInParaView(res, filename="CheckIntegrity_GetDataOverALine_line.xdmf")
return "ok"
[docs]def CheckIntegrity(GUI: bool = False):
toTest = [
CheckIntegrity_ComputeNormalsAtElements,
CheckIntegrity_GetDataOverALine,
CheckIntegrity_MeshQualityAspectRatioBeta,
CheckIntegrity_EnsureUniquenessElements,
CheckIntegrity_ExtractElementsByMask,
CheckIntegrity_GetVolume,
CheckIntegrity_ComputeNodeToNodeConnectivity,
CheckIntegrity_ComputeMeshMinMaxLengthScale,
CheckIntegrity_GetElementsFractionInside,
CheckIntegrity_ExtractElementByTags,
]
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
return RunListOfCheckIntegrities(toTest, GUI)
if __name__ == '__main__':
print(CheckIntegrity(True)) # pragma: no cover