# -*- 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 Tuple, Optional, Union, List, Dict
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
from Muscat.Helpers.ProgressBar import PrintProgressBar
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.MeshCreationTools import QuadToLin
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.LinAlg.Transform import Transform
from Muscat.FE.Fields.FEField import FEField
from Muscat.Containers.FieldTransferOpPython import GetFieldTransferOpPython
from Muscat.Helpers.Logger import Error
[docs]def ApplyRotationMatrixTensorField(fields: Dict, fieldsToTreat: List[List], baseNames: List = ["v1", "v2"], inPlace: bool = False, prefix: str = "new_", inverse: bool = False) -> Dict:
"""Apply a rotation operation on the fields
Parameters
----------
fields : dict[ArrayLike]
dictionary of field. Keys are names, values are data
fieldsToTreat : ArrayLike
is a 3x3 list of list with the names of the fields to create the local tensor.
For example [["S00","S01","S02"]["S01","S11","S12"]["S02","S12","S22"]]
baseNames : List, optional
the names of the fields to extract the direction (each field must be on size (nb entries, 3) ), by default ["v1","v2"]
inPlace : bool, optional
True to put the output back to the field dictionary, by default False
prefix : str, optional
prefix to prepend to the new field, by default `new_`
inverse : bool, optional
True to apply the inverse transform, by default False
Returns
-------
dict
a dictionary containing the field after the rotation operation
"""
nbEntries = fields[fieldsToTreat[0][0]].shape[0]
bs = Transform()
bs.keepNormalized = True
v1 = fields[baseNames[0]]
v2 = fields[baseNames[1]]
tempData = np.zeros((nbEntries, 3, 3))
for i in range(3):
for j in range(3):
tempData[:, i, j] = fields[fieldsToTreat[i][j]][:]
for i in range(nbEntries):
bs.SetFirst(v1[i, :])
bs.SetSecond(v2[i, :])
if inverse:
tempData[i, :, :] = bs.ApplyInvTransformTensor(tempData[i, :, :])
else:
tempData[i, :, :] = bs.ApplyTransformTensor(tempData[i, :, :])
output = {}
for i in range(3):
for j in range(3):
output[fieldsToTreat[i][j]] = tempData[:, i, j]
output = {(prefix + x): v for x, v in output.items()}
if inPlace:
fields.update(output)
return output
[docs]def CopyFieldsFromOriginalMeshToTargetMesh(inMesh: Mesh, outMesh: Mesh) -> None:
"""Function to copy fields (nodeFields and elemFields) for the original mesh to the
derived mesh ( f(inMesh) -> outMesh )
Parameters
----------
inMesh : Mesh
the source mesh, we extract the fields from inMesh.nodeFields and inMesh.elemFields.
outMesh : Mesh
The target mesh, we push the new fields into outMesh.nodeFields and outMesh.elemFields.
"""
def Work(inDict, outDict, oid):
for k, d in inDict.items():
outDict[k] = d[oid]
Work(inMesh.nodeFields, outMesh.nodeFields, outMesh.originalIDNodes)
# Compute the transfer array for the elemFields
cpt1 = 0
cpt2 = 0
oid = np.empty(outMesh.GetNumberOfElements(), dtype=MuscatIndex)
for name, data in inMesh.elements.items():
if name in outMesh.elements:
elements = outMesh.elements[name]
oid[cpt2 : cpt2 + elements.GetNumberOfElements()] = elements.originalIds + cpt1
cpt2 += elements.GetNumberOfElements()
cpt1 += data.GetNumberOfElements()
Work(inMesh.elemFields, outMesh.elemFields, oid)
[docs]def GetFieldTransferOp(
inputField: FEField, targetPoints: ArrayLike, method: Union[str, None] = None, verbose: bool = False, elementFilter: Optional[ElementFilter] = None
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
try:
return GetFieldTransferOpCpp(inputField, targetPoints, method, verbose=verbose, elementFilter=elementFilter)
except ImportError: # pragma: no cover
Error("Error in the cpp version GetFieldTransferOp. Using Python (slow) version ")
return GetFieldTransferOpPython(inputField, targetPoints, method, verbose=verbose, elementFilter=elementFilter)
except: # pragma: no cover
raise
[docs]def GetFieldTransferOpCpp(
inputField: FEField, targetPoints: ArrayLike, method: Union[str, None] = None, verbose: bool = False, elementFilter: Optional[ElementFilter] = None, options: Optional[Dict] = None, nbThreads=None
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Compute the transfer operator from the inputField to the target points so, Using the cpp implementation:
valueAtTargetPoints = op.dot(FEField.data)
The methods available are:
- "Interp/Nearest"
- "Nearest/Nearest"
- "Interp/Clamp"
- "Interp/Extrap"
- "Interp/ZeroFill"
Possible values (int) for the returned status are:
- 0: "Nearest"
- 1: "Interp"
- 2: "Extrap"
- 3: "Clamp"
- 4: "ZeroFill"
- 5: "MultipleExtrapPoint"
Parameters
----------
inputField : FEField
the FEField to be transferred
targetPoints : ArrayLike
Numpy array of the target points. Position to extract the values
method : Union[str,None], optional
A couple for the algorithm used when the point is inside/outside the mesh
- "Interp" -> to use the interpolation of the FEField to extract the values
- "Nearest" -> to use the closest node to extract the values
- "Clamp" -> to use the closest point to the surface if the point is outside
- "Extrap" -> to use shape function of the closest element to compute a extrapolation
- "ZeroFill" -> fill with zero in the case the point is outside
If None is provided then "Interp/Clamp" is used
verbose : bool, optional
Print a progress bar, by default False
elementFilter : Optional[ElementFilter], optional
ElementFilter to extract the information from only a subdomain, by default None
options : Optional[Dict], optional
Dictionary containing the one or more key:values pairs:
- "usePointSearch": bool , to activate the research of potential element by point, default (True)
- "useElementSearch": bool , to activate the research of potential element by a center of the closest element and it neighbors, default (False)
- "useElementSearchFast": bool , to activate the research of potential element by a center of the closest element only (option useElementSearch must be true) default (False)
- "useEdgeSearch": bool , to activate the research of potential element using the closest edge. default (True)
- "DifferentialOperator": int, to activate the computation of the derivative of the field (not the value), default -1.
possible values are; -1 for value transfer; 0 to compute the d/dx; 1 to compute d/dy; 2 to compute d/dz
nbThreads : Optional[int], optional
Set the maximum number of thread to use
Returns
-------
Tuple [np.ndarray,np.ndarray]
a tuple with 2 object containing:
- op, sparse matrix with the operator to make the transfer
- status: vector of ints with the status transfer for each target point:
- entities: vector of the ids of the entities (node or element) used to compute the transfer operator.
If status = 0, then entity is a node id. if status :math:`\in [1,2,3,5]` the status is an element id.
return op, status, entities
"""
from Muscat.Containers.NativeTransfer import NativeTransfer
nt = NativeTransfer()
if nbThreads is not None:
nt.SetMaxConcurrency(nbThreads)
nt.SetVerbose(verbose)
nt.SetTargetPoints(targetPoints)
if method is None:
method = "Interp/Clamp"
nt.SetTransferMethod(method)
defaultOptions = {
"usePointSearch": True,
"useElementSearch": False,
"useElementSearchFast": False,
"useEdgeSearch": True,
}
if options is None:
options = {}
defaultOptions.update(options)
dispatch = {
"usePointSearch": nt.SetUsePointSearch,
"useElementSearch": nt.SetUseElementSearch,
"useElementSearchFast": nt.SetUseElementSearchFast,
"useEdgeSearch": nt.SetUseEdgeSearch,
"DifferentialOperator": nt.SetDifferentialOperator,
}
for k, v in defaultOptions.items():
if k in dispatch.keys():
dispatch[k](v)
else: # pragma: no cover
raise RuntimeError(f"Option {k} not valid")
nt.SetSourceFEField(inputField, elementFilter)
nt.Compute()
op = nt.GetOperator()
status = nt.GetStatus()
entities = nt.GetEntitiesIds()
return op, status, entities
[docs]def TransportPosToPoints(mesh: Mesh, points: ArrayLike, method: str = "Interp/Clamp", verbose: bool = False) -> np.ndarray:
"""For each point in points compute the position of the source data by transporting the position field
Parameters
----------
mesh : Mesh
the input mesh from where we
points : ArrayLike
the target points
method : str, optional
Transfer method, by default "Interp/Clamp"
verbose : bool, optional
True to print more output during the computation of the transfer operator, by default None
Returns
-------
np.ndarray
for each point in points the position of the source data
"""
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, offset, NGauss = PrepareFEComputation(mesh, numberOfComponents=1)
field = FEField("", mesh=mesh, space=space, numbering=numberings[0])
op, status, entities = GetFieldTransferOp(field, points, method=method, verbose=verbose, elementFilter=ElementFilter())
return op.dot(mesh.nodes)
[docs]def TransportPos(imesh: Mesh, tmesh: Mesh, tspace, tnumbering, method: str = "Interp/Clamp", verbose: bool = False) -> np.ndarray:
"""Function to transport the position from the input mesh (imesh)
to a target FEField target mesh, target space, tnumbering
Parameters
----------
imesh : Mesh
input mesh
tmesh : Mesh
target mesh
tspace : _type_
target space
tnumbering : _type_
target numbering
method : str, optional
method for the interpolation/extrapolation, by default "Interp/Clamp"
verbose : bool, optional
True to print more output during the computation of the transfer operator, by default None
Returns
-------
np.ndarray
a numpy.array with 3 FEField for each component of the position
"""
from Muscat.FE.Fields.FEField import FEField
pos = TransportPosToPoints(imesh, tmesh.nodes, method=method, verbose=verbose)
names = ["x", "y", "z"][0 : imesh.GetPointsDimensionality()]
posFields = np.array([FEField("ipos_" + names[x], tmesh, space=tspace, numbering=tnumbering, data=pos[:, x]) for x in range(len(names))])
return posFields
[docs]def PointToCellData(mesh: Mesh, pointfield: np.ndarray, dim: Optional[MuscatIndex] = None) -> np.ndarray:
"""Convert a point field to cell field. the cell values are compute the average
of the values on the nodes for each element
Parameters
----------
mesh : Mesh
mesh support of the field
pointfield : np.ndarray
a field defined on the points
dim : int, optional
dimensionality filter. if dim != the result array contain only values for the selected element
( len(result) is smaller than mesh.GetNumberOfElements() ) , by default None
Returns
-------
np.ndarray
a numpy array with the field on the elements
"""
from Muscat.Containers.Filters.FilterTools import GetFrozenFilter
if dim is not None:
filt = GetFrozenFilter(ElementFilter(dimensionality=dim), mesh)
else:
filt = GetFrozenFilter(ElementFilter(), mesh)
nbelemtns = filt.GetElementSelectionSize()
if len(pointfield.shape) == 2:
ncols = pointfield.shape[1]
res = np.zeros((nbelemtns, ncols), dtype=float)
else:
ncols = 1
res = np.zeros((nbelemtns), dtype=float)
# cpt = 0
for selection in filt(mesh):
# name,data,ids
if len(pointfield.shape) == 1:
valAtCenter = (np.sum(pointfield[selection.elements.connectivity], axis=1) / selection.elements.GetNumberOfNodesPerElement()).flatten()
# res[cpt:cpt+selection.elements.GetNumberOfElements()] = valAtCenter
res[selection.GetSelectionSlice()] = valAtCenter
else:
for i in range(ncols):
valAtCenter = (np.sum(pointfield[selection.elements.connectivity, i], axis=1) / selection.elements.GetNumberOfNodesPerElement()).flatten()
# res[cpt:cpt+data.GetNumberOfElements(),i] = valAtCenter
res[selection.GetSelectionSlice(), i] = valAtCenter
# cpt += len(ids)
return res
[docs]def QuadFieldToLinField(quadMesh: Mesh, quadField: np.ndarray, linMesh: Optional[Mesh] = None) -> np.ndarray:
"""Extract the linear part of the field quadField.
Parameters
----------
quadMesh : Mesh
the quadratic mesh supporting the quad field
quadField : np.ndarray
the field to be converted to linear
linMesh : Mesh, optional
the support of the linear field, if None a linear conversion of the quadMesh is done (with QuandTolin), by default None
Returns
-------
np.ndarray
_description_
"""
if linMesh == None:
linMesh = QuadToLin(quadMesh)
extractIndices = np.arange(quadMesh.GetNumberOfNodes())[linMesh.originalIDNodes]
return quadField[extractIndices]
# ------------------------- CheckIntegrity ------------------------
[docs]def CheckIntegrityApplyRotationMatrixTensorField(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
inputmesh = CreateUniformMeshOfBars(0, 1, 10)
nn = inputmesh.GetNumberOfNodes()
inputmesh.nodeFields["Sxx"] = np.ones(nn)
inputmesh.nodeFields["Syy"] = np.full(nn, 2)
inputmesh.nodeFields["Szz"] = np.full(nn, 3)
inputmesh.nodeFields["Sxy"] = np.full(nn, 0)
inputmesh.nodeFields["Sxz"] = np.full(nn, 0)
inputmesh.nodeFields["Syz"] = np.full(nn, 0)
inputmesh.nodeFields["v1"] = np.zeros((nn, 3))
inputmesh.nodeFields["v1"][:, 0] = 1
inputmesh.nodeFields["v2"] = np.zeros((nn, 3))
inputmesh.nodeFields["v2"][:, 2] = 1
res = ApplyRotationMatrixTensorField(inputmesh.nodeFields, [["Sxx", "Sxy", "Sxz"], ["Sxy", "Syy", "Syz"], ["Sxz", "Syz", "Szz"]], baseNames=["v1", "v2"], inPlace=False, prefix="", inverse=False)
print(list(res.keys()))
res["v1"] = inputmesh.nodeFields["v1"]
res["v2"] = inputmesh.nodeFields["v2"]
res2 = ApplyRotationMatrixTensorField(res, [["Sxx", "Sxy", "Sxz"], ["Sxy", "Syy", "Syz"], ["Sxz", "Syz", "Szz"]], baseNames=["v1", "v2"], inPlace=True, prefix="new_", inverse=True)
# for k, v in res.items():
# print(k, v)
return "ok"
[docs]def RunTransfer(inputFEField, data, outmesh, methods = ["Interp/Nearest", "Nearest/Nearest", "Interp/Clamp", "Interp/Extrap"]) -> None:
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tempdir = TemporaryDirectory.GetTempPath()
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
WriteMeshToXdmf(tempdir + "GetFieldTransferOp_Original" + inputFEField.name + ".xdmf", inputFEField.mesh, PointFields=[data], PointFieldsNames=["OriginalData"])
PointFields = []
PointFieldsNames = []
from Muscat.Helpers.Logger import logger
import time
from Muscat.Helpers.TextFormatHelper import TFormat
def PrintRow(datarow):
res = "|"
for d in datarow:
if type(d) is str:
res += TFormat.Center(str(d), fill=" ", width=20) + "|"
else:
res += TFormat.Center("%.4f" % d, fill=" ", width=20) + "|"
print(res)
PrintRow(["method", "cpp [ms]", "python [ms]", "speedup "])
for method in methods:
cppStartTime = time.time()
opCpp, statusCpp, entitiesIdsCpp = GetFieldTransferOpCpp(inputFEField, outmesh.nodes, method=method, verbose=logger.getEffectiveLevel() <= 20, nbThreads=1)
cppStopTime = time.time()
newdataCpp = opCpp.dot(data)
PointFieldsNames.append(method + "Cpp")
PointFields.append(newdataCpp)
PointFieldsNames.append(method + "Status" + "Cpp")
PointFields.append(statusCpp)
newPosCpp = opCpp.dot(inputFEField.mesh.nodes)
PointFieldsNames.append(method + "Pos" + "Cpp")
PointFields.append(newPosCpp)
PointFieldsNames.append(method + "Entities" + "Cpp")
PointFields.append(entitiesIdsCpp)
pythonStartTime = time.time()
opPython, statusPython, entitiesIdsPython = GetFieldTransferOpPython(inputFEField, outmesh.nodes, method=method, verbose=logger.getEffectiveLevel() <= 20)
pythonStopTime = time.time()
newdataPython = opPython.dot(data)
PointFieldsNames.append(method + "Python")
PointFields.append(newdataPython)
PointFieldsNames.append(method + "Status" + "Python")
PointFields.append(statusPython)
newPosPython = opPython.dot(inputFEField.mesh.nodes)
PointFieldsNames.append(method + "Pos" + "Python")
PointFields.append(newPosPython)
PointFieldsNames.append(method + "Entities" + "Python")
PointFields.append(entitiesIdsPython)
PrintRow([method, cppStopTime - cppStartTime, pythonStopTime - pythonStartTime, (pythonStopTime - pythonStartTime) / ((cppStopTime - cppStartTime) if cppStopTime - cppStartTime > 0 else 1)])
PointFields.append(newdataCpp-newdataPython)
PointFieldsNames.append(method+"_Diff_cpp_python")
try:
np.testing.assert_allclose(newdataCpp, newdataPython)
except: # pragma: no cover
print(method)
print(newdataCpp.shape)
print(newdataPython.shape)
print(newdataCpp)
print(newdataPython)
WriteMeshToXdmf(tempdir + "GetFieldTransferOp_TargetMesh" + inputFEField.name + ".xdmf", outmesh, PointFields=PointFields, PointFieldsNames=PointFieldsNames, Binary = True, HDF5= False)
print(f"Last file with the data : \n {tempdir}GetFieldTransferOp_TargetMesh{inputFEField.name}.xdmf")
print(f"Cpp and python version of the field transfer gives differents solutions for method : {method}")
raise
WriteMeshToXdmf(tempdir + "GetFieldTransferOp_TargetMesh" + inputFEField.name + ".xdmf", outmesh, PointFields=PointFields, PointFieldsNames=PointFieldsNames, Binary = True, HDF5= False)
[docs]def CheckIntegrity1D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
inputmesh = CreateUniformMeshOfBars(0, 1, 10)
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
b = CreateUniformMeshOfBars(-0.1, 1.1, 15)
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="", mesh=inputmesh, space=space, numbering=numberings[0])
# for el,data in inputmesh.elements.items():
# print(data.connectivity)
if GUI: # pragma: no cover
import matplotlib.pyplot as plt
fig, axs = plt.subplots(nrows=2, ncols=3, constrained_layout=True)
axis = axs.flat
else:
axis = [None] * 5
data = (inputmesh.nodes[:, 0] - 0.5) ** 2
from Muscat.Helpers.Logger import logger
for method, ax in zip(["Interp/Nearest", "Nearest/Nearest", "Interp/Clamp", "Interp/Extrap", "Interp/ZeroFill"], axis):
op = GetFieldTransferOp(inputFEField, b.nodes, method=method, verbose=logger.getEffectiveLevel() <= 20)[0]
result = op.dot(data)
if GUI: # pragma: no cover
ax.plot(inputmesh.nodes[:, 0], data, "X-", label="Original Data")
ax.plot(b.nodes[:, 0], result, "o:", label=method)
legend = ax.legend(loc="upper center", shadow=True, fontsize="x-large")
ax.set(xlabel="x", ylabel="data", title=method)
if GUI: # pragma: no cover
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tempdir = TemporaryDirectory.GetTempPath()
fig.savefig(tempdir + "test.png")
plt.show()
return "ok"
[docs]def CheckIntegrity2D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
inputmesh = CreateSquare(dimensions=[5, 5], origin=[0, 0], spacing=[1.0 / 4, 1.0 / 4.0])
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
N = 10
b = CreateSquare(dimensions=[N, N], origin=[-0.1, -0.1], spacing=[1.2 / (N - 1), 1.2 / (N - 1)])
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="2DTo2D", mesh=inputmesh, space=space, numbering=numberings[0])
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
RunTransfer(inputFEField, data, b)
return "ok"
[docs]def CheckIntegrity1DTo2D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
inputmesh = CreateUniformMeshOfBars(0, 1, 10)
inputmesh.nodes[:, 1] = inputmesh.nodes[:, 0] ** 2
inputmesh.nodes = inputmesh.nodes[:, 0:2].copy(order="C")
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
N = 10
b = CreateSquare(dimensions=[N, N], origin=[-0.5, -0.5], spacing=[2 / (N - 1), 2 / (N - 1)])
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="1DTo2D", mesh=inputmesh, space=space, numbering=numberings[0])
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
RunTransfer(inputFEField, data, b)
return "ok"
[docs]def CheckIntegrity2DTo3D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare, CreateCube
N = 10
inputmesh = CreateSquare(dimensions=[N, N], origin=[0, 0], spacing=[1 / (N - 1), 1 / (N - 1)])
inputmesh = ExtractElementsByElementFilter(inputmesh,ElementFilter(dimensionality=2))
n = inputmesh.nodes
inputmesh.nodes = np.zeros((n.shape[0], 3))
inputmesh.nodes[:, 0:2] = n
inputmesh.nodes[:, 2] = n[:, 0] ** 3
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="2DTo3D", mesh=inputmesh, space=space, numbering=numberings[0])
N = 10
b = CreateCube(dimensions=[N, N, N], origin=[-0.501] * 3, spacing=[2 / (N - 1), 2 / (N - 1), 2 / (N - 1)])
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
RunTransfer(inputFEField, data, b, methods = ["Interp/Nearest", "Nearest/Nearest", "Interp/Clamp"])
return "ok"
[docs]def CheckIntegrity3DTo3DNative(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateCube
N = 10
inputmesh = CreateCube(dimensions=[N, N, N], origin=[-1.0] * 3, spacing=[2 / (N - 1), 2 / (N - 1), 2 / (N - 1)])
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="3DTo3DNative", mesh=inputmesh, space=space, numbering=numberings[0])
N = 6
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
outmesh = CreateConstantRectilinearMesh(dimensions=[N, N, N], origin=[-1.0] * 3, spacing=[2 / (N - 1), 2 / (N - 1), 2 / (N - 1)])
# b = CreateCube(dimensions=[N, N, N], origin=[-1.]*3, spacing=[2/(N-1), 2/(N-1), 2/(N-1)])
from Muscat.LinAlg.Transform import Transform
op = Transform()
op.keepNormalized = False
op.keepOrthogonal = False
op.SetFirst([1.4, 0.56, 0])
op.SetSecond([-1.135, 1.42486, 1.6102])
outmesh.nodes = np.ascontiguousarray(op.ApplyTransform(outmesh.nodes))
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
print("Start transfer timing...")
print(f"Number of elements in the origin mesh = {inputmesh.GetNumberOfElements()}")
print(f"Number of points in the target cloud point = {outmesh.GetNumberOfNodes()}")
import time
from Muscat.Helpers.TextFormatHelper import TFormat
def PrintRow(datarow):
res = "|"
for d in datarow:
if type(d) is str:
res += TFormat.Center(str(d), fill=" ", width=20) + "|"
else:
res += TFormat.Center("%.4f" % d, fill=" ", width=20) + "|"
print(res)
from Muscat.Containers.NativeTransfer import NativeTransfer
nt = NativeTransfer()
nt.SetVerbose(False)
nt.SetUsePointSearch(True)
nt.SetUseElementSearch(True)
nt.SetUseElementSearchFast(True)
nt.SetUseEdgeSearch(True)
setFieldTime = time.time()
nt.SetSourceFEField(inputFEField)
setFieldTime = time.time() - setFieldTime
print(f"Set Up time (SetSourceFEField) {setFieldTime} [ms]")
output = ["method"]
output.extend(f"cpp [ms] {th} thread " for th in [1, 2, 3, 4])
output.extend(f"speedup {th}/1 " for th in [2, 3, 4])
PrintRow(output)
for method in ["Interp/Nearest", "Nearest/Nearest", "Interp/Clamp", "Interp/Extrap"]:
nt.SetTransferMethod(method)
nt.SetTargetPoints(outmesh.nodes)
times = []
for nbThreads in [1, 2, 3, 4]:
nt.SetMaxConcurrency(nbThreads)
cppStartTime = time.time()
nt.Compute()
cppStopTime = time.time()
op = nt.GetOperator()
times.append(cppStopTime - cppStartTime)
output = [method]
output.extend(times)
output.extend([times[0] / (t if t > 0 else 1) for t in times[1:]])
PrintRow(output)
return "ok"
[docs]def CheckIntegrity3DTo3D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateCube
N = 10
inputmesh = CreateCube(dimensions=[N, N, N], origin=[-1.0] * 3, spacing=[2 / (N - 1), 2 / (N - 1), 2 / (N - 1)])
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="3DTo3D", mesh=inputmesh, space=space, numbering=numberings[0])
N = 4
b = CreateCube(dimensions=[N, N, N], origin=[-1.0] * 3, spacing=[2 / (N - 1), 2 / (N - 1), 2 / (N - 1)])
from Muscat.LinAlg.Transform import Transform
op = Transform()
op.keepNormalized = False
op.keepOrthogonal = False
op.SetFirst([1.4, 0.56, 0])
op.SetSecond([-1.135, 1.42486, 1.6102])
b.nodes = np.ascontiguousarray(op.ApplyTransform(b.nodes))
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
RunTransfer(inputFEField, data, b)
return "ok"
[docs]def CheckIntegrity2DTo1DDerivativeOP(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
N = 10
a = CreateSquare(dimensions=[N, N], origin=[-0.5, -0.5], spacing=[2 / (N - 1), 2 / (N - 1)])
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(a)
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="2DTo2DDerivative", mesh=a, space=space, numbering=numberings[0])
inputFEField.data = 2 * a.nodes[:, 0] + a.nodes[:, 1] ** 2
outmesh = CreateSquare(dimensions=[N - 1, N - 1], origin=[-0.5 + 1 / (N - 1), -0.5 + 1 / (N - 1)], spacing=[2 / (N - 1), 2 / (N - 1)])
op, status0, entitiesIds = GetFieldTransferOpCpp(inputFEField, outmesh.nodes, method="Interp/Clamp", options={"DifferentialOperator": 0}, nbThreads=1)
derivative0 = op.dot(inputFEField.data)
from Muscat.Containers.NativeTransfer import NativeTransfer
nt = NativeTransfer()
nt.SetVerbose(False)
#nt.SetMaxConcurrency(1)
nt.SetTargetPoints(outmesh.nodes)
nt.SetTransferMethod("Interp/Clamp")
nt.SetDifferentialOperator(1)
nt.SetSourceFEField(inputFEField, None)
nt.Compute()
op = nt.GetOperator()
status1 = nt.GetStatus()
entities1 = nt.GetEntitiesIds()
derivative1 = op.dot(inputFEField.data)
try:
derivative0_exact = np.full(outmesh.GetNumberOfNodes(), dtype=MuscatFloat, fill_value=2)
np.testing.assert_allclose(derivative0, derivative0_exact)
derivative1_exact = 2 * outmesh.nodes[:, 1]
np.testing.assert_allclose(derivative1, derivative1_exact)
except: # pragma: no cover
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tempdir = TemporaryDirectory.GetTempPath()
WriteMeshToXdmf(tempdir + "GetFieldTransferOp_OriginMesh" + inputFEField.name + ".xdmf", a, PointFields=[inputFEField.data], PointFieldsNames=["data"])
WriteMeshToXdmf(tempdir + "GetFieldTransferOp_TargetMesh" + inputFEField.name + ".xdmf", outmesh, PointFields=[derivative0, derivative1,status0,status1,entities1], PointFieldsNames=["derivative0", "derivative1","status0","status1","entities1"])
print(f"Last file with the data : \n {tempdir}GetFieldTransferOp_TargetMesh{inputFEField.name}.xdmf")
raise
return "OK"
[docs]def CheckIntegrity1DSecondOrderTo2D(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateSquare
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
inputmesh = CreateUniformMeshOfBars(0, 1, 11, secondOrder=True)
inputmesh.nodes[:, 1] = inputmesh.nodes[:, 0] ** 2
inputmesh.nodes = inputmesh.nodes[:, 0:2].copy(order="C")
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, _offset, _NGauss = PrepareFEComputation(inputmesh)
N = 10
# b = CreateSquare(dimensions = [N,N],origin=[-0.5,-0.5],spacing=[2/(N-1),2/(N-1)])
# b = CreateSquare(dimensions = [N,N],origin=[-0.1,0.0],spacing=[1./(N-1),0.7/(N-1)])
b = CreateSquare(dimensions=[N, N], origin=[-0.5, -0.5], spacing=[2 / (N - 1), 2 / (N - 1)])
from Muscat.FE.Fields.FEField import FEField
inputFEField = FEField(name="1DSecondTo2D", mesh=inputmesh, space=space, numbering=numberings[0])
x = inputmesh.nodes[:, 0]
y = inputmesh.nodes[:, 1]
data = (x - 0.5) ** 2 - y * 0.5 + x * y * 0.25
RunTransfer(inputFEField, data, b)
return "ok"
[docs]def CheckIntegrity_PointToCellData(GUI: bool = False) -> str:
myMesh = Mesh()
myMesh.nodes = np.array([[0, 0, 0], [1, 0, 0], [2, 0, 0]], dtype=MuscatFloat)
myMesh.originalIDNodes = np.array([0, 1, 2], dtype=int)
tag = myMesh.GetNodalTag("linPoints")
tag.AddToTag(0)
tag.AddToTag(1)
tag.AddToTag(2)
import Muscat.Containers.ElementsDescription as ED
elements = myMesh.GetElementsOfType(ED.Bar_2)
elements.AddNewElement([0, 1], 3)
elements.AddNewElement([1, 2], 4)
myMesh.AddElementToTagUsingOriginalId(3, "LinElements")
myMesh.AddElementToTagUsingOriginalId(4, "LinElements")
myMesh.PrepareForOutput()
# print(myMesh)
res = PointToCellData(myMesh, np.array([[-2, 2, 4]]).T)
ExactData = np.array([[0, 3]], dtype=float).T
# print(res - ExactData)
if (res - ExactData).any():
return "Error CheckIntegrity_PointToCellData" # pragma: no cover
res = PointToCellData(myMesh, np.array([-2, 2, 4]))
res = PointToCellData(myMesh, np.array([[-2, 2, 4]]).T, dim=2)
return "ok"
[docs]def CheckIntegrityCopyFieldsFromOriginalMeshToTargetMesh(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateMeshOf
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
points = [[-0.5, -0.5, -0.5], [2.5, -0.5, -0.5], [-0.5, 2.5, -0.5], [-0.5, -0.5, 2.5], [2.5, 2.5, 2.5]]
tets = [[0, 1, 2, 3]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
mesh.nodeFields["data"] = mesh.nodes.view()
mesh2 = ExtractElementsByElementFilter(mesh, ElementFilter(elementType=ED.Tetrahedron_4))
CopyFieldsFromOriginalMeshToTargetMesh(mesh, mesh2)
return "ok"
[docs]def CheckIntegrityTransportPos(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateMeshOf
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
points = [[-0.5, -0.5, -0.5], [2.5, -0.5, -0.5], [-0.5, 2.5, -0.5], [-0.5, -0.5, 2.5], [2.5, 2.5, 2.5]]
tets = [[0, 1, 2, 3]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
recMesh = CreateConstantRectilinearMesh(dimensions=[5, 5, 5], origin=[-1, -1, -1], spacing=[1, 1, 1])
from Muscat.FE.FETools import PrepareFEComputation
space, numberings, offset, NGauss = PrepareFEComputation(mesh)
TransportPos(recMesh, mesh, space, tnumbering=numberings)
return "ok"
[docs]def CheckIntegrityQuadFieldToLinField(GUI: bool = False) -> str:
"""Check done in Muscat.Containers.MeshCreationTools.CheckIntegrity_QuadToLin"""
return "ok"
[docs]def CheckIntegrity(GUI: bool = False) -> str:
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
totest = [
# we dont mesure speed up in the CheckIntegrity1D, so we us it
# first to load all the cpp libs (.so)
CheckIntegrity1D,
CheckIntegrity2DTo1DDerivativeOP,
CheckIntegrity1D,
CheckIntegrity_PointToCellData,
CheckIntegrity1DSecondOrderTo2D,
CheckIntegrity1DTo2D,
CheckIntegrity2D,
CheckIntegrity2DTo3D,
CheckIntegrityApplyRotationMatrixTensorField,
CheckIntegrity3DTo3D,
CheckIntegrityCopyFieldsFromOriginalMeshToTargetMesh,
CheckIntegrityTransportPos,
CheckIntegrityQuadFieldToLinField,
CheckIntegrity3DTo3DNative,
]
return RunListOfCheckIntegrities(totest)
if __name__ == "__main__":
print(CheckIntegrity(True)) # pragma: no cover