# -*- 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 __future__ import annotations
from typing import Union, List, Tuple, Callable, Optional, Iterable, Dict
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.Filters.FilterBase import FilterBase
from Muscat.Containers.Filters.FilterObjects import ElementFilter, NodeFilter
from Muscat.Containers.Filters.FilterOperators import IntersectionFilter, FilterOperatorBase
import Muscat.Containers.ElementsDescription as ED
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Fields.IPField import FieldBase, IPField, RestrictedIPField
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo, FESpace
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.IntegrationRules import MeshQuadrature
[docs]
def Maximum(A, B):
return A.binaryOp(B, np.maximum)
[docs]
def Minimum(A, B):
return A.binaryOp(B, np.minimum)
[docs]
def ElementWiseIpToFETransferOp(integrationRule: MeshQuadrature, space: FESpace) -> Dict[str, np.ndarray]:
"""Generate transfer operator (element wise) to pass information from integration points to a FE field. This is done
by solving a least-squares problem.
Parameters
----------
integrationRule : IntegrationRuleType
The integration rule of the original data
space : FESpace
The target space
Returns
-------
Dict[str,np.ndarray]
the operator to pass the data from the integration points to the FEField
"""
res: Dict[str, np.ndarray] = dict()
for name, ir in integrationRule.items():
space_ipValues = space[name].GetSpaceEvaluatedAt(ir.points, ir.weights)
valN = np.asarray(space_ipValues.valN, dtype=MuscatFloat)
sol = np.linalg.lstsq(valN, np.eye(valN.shape[0], valN.shape[0]), rcond=None)[0]
res[name] = sol
return res
[docs]
def ElementWiseFEToFETransferOp(originSpace: FESpace, targetSpace: FESpace) -> Dict[str, np.ndarray]:
"""Generate transfer operator (element wise) to pass information from a FE Field to a different FE field. This is done
by solving a least-squares problem.
Parameters
----------
originSpace : FESpace
The initial space
targetSpace : FESpace
The target space
Returns
-------
Dict[str,np.ndarray]
the operator to pass the data from the initial space to the target FEField
"""
res: Dict[str, np.ndarray] = dict()
for name, ir in targetSpace.items():
ir = (originSpace[name].posN, np.ones(originSpace[name].GetNumberOfShapeFunctions(), dtype=MuscatFloat))
spaceIpValues = targetSpace[name].GetSpaceEvaluatedAt(ir[0], ir[1])
valN = np.asarray(spaceIpValues.valN, dtype=MuscatFloat)
sol = np.linalg.lstsq(valN, np.eye(valN.shape[0], valN.shape[1]), rcond=None)[0]
res[name] = sol
return res
[docs]
def NodeFieldToFEField(mesh: Mesh, nodeFields: Dict[str, ArrayLike] = None, copy: bool = False) -> Dict[str, FEField]:
"""Create FEField(P isoparametric) from the node field data.
if nodesFields is None the mesh.nodeFields is used
Parameters
----------
mesh : Mesh
The support for the FEFields
nodeFields : Dict[str,ArrayLike], optional
the dictionary containing the nodes fields to be converted to FEFields,
if None the mesh.nodeFields is used, by default None
Returns
-------
Dict[str,FEField]
A dictionary the keys are field names and the values are the FEFields
"""
if nodeFields is None:
nodeFields = mesh.nodeFields
numbering = ComputeDofNumbering(mesh, LagrangeSpaceGeo, fromConnectivity=True)
res = {}
for name, values in nodeFields.items():
if len(values.shape) == 2 and values.shape[1] > 1:
for i in range(values.shape[1]):
data = np.asarray(values[:, i], dtype=MuscatFloat, order="C")
if copy:
data = data.copy()
res[name + "_" + str(i)] = FEField(name=name + "_" + str(i), mesh=mesh, space=LagrangeSpaceGeo, numbering=numbering, data=data)
else:
data = np.asarray(values, dtype=MuscatFloat, order="C")
if copy:
data = data.copy()
res[name] = FEField(name=name, mesh=mesh, space=LagrangeSpaceGeo, numbering=numbering, data=data)
return res
[docs]
def ElemFieldsToFEField(mesh: Mesh, elemFields: Dict[str, ArrayLike] = None, copy=False) -> Dict[str, FEField]:
"""Create FEField(P0) from the elements field data.
if elemFields is None the mesh.elemFields is used
Parameters
----------
mesh : Mesh
The support for the FEFields
elemFields : Dict[str,ArrayLike], optional
a dict containing the elements fields to be converted to FEFields
if elemFields is None the mesh.elemFields is used, by default None
Returns
-------
Dict[str,FEField]
A dictionary the keys are field names and the values are the FEFields
"""
if elemFields is None:
elemFields = mesh.elemFields
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
numbering = ComputeDofNumbering(mesh, LagrangeSpaceP0)
res = {}
for name, values in elemFields.items():
if len(values.shape) == 2 and values.shape[1] > 1:
for i in range(values.shape[1]):
data = np.asarray(values[:, i], dtype=MuscatFloat, order="C")
if copy:
data = data.copy()
res[name + "_" + str(i)] = FEField(name=name + "_" + str(i), mesh=mesh, space=LagrangeSpaceP0, numbering=numbering, data=data)
else:
data = np.asarray(values, dtype=MuscatFloat, order="C")
if copy:
data = data.copy()
res[name] = FEField(name=name, mesh=mesh, space=LagrangeSpaceP0, numbering=numbering, data=data)
return res
[docs]
def FEFieldsDataToVector(listOfFields: List[FEField], outVector: Optional[ArrayLike] = None, copy: bool = False) -> np.ndarray:
"""Put FEFields data into a vector format
Parameters
----------
listOfFields : List[FEField]
list of FEFields (the order is important)
outVector : Optional[ArrayLike], optional
if provided the output vector, by default None
Returns
-------
np.ndarray
_description_
"""
if outVector is None:
s = sum([f.numbering.size for f in listOfFields])
outVector = np.zeros(s, dtype=MuscatFloat)
offset = 0
for f in listOfFields:
if copy:
outVector[offset : offset + f.numbering.size] = f.data.copy()
else:
outVector[offset : offset + f.numbering.size] = f.data
offset += f.numbering.size
return outVector
[docs]
def VectorToFEFieldsData(inVector: ArrayLike, listOfFields: List[FEField]):
"""Put vector data in FEFields
Parameters
----------
inVector : ArrayLike
vector with the data
listOfFields : List[FEField]
list of FEFields to store the data
"""
offset = 0
for f in listOfFields:
f.data = inVector[offset : offset + f.numbering.size]
offset += f.numbering.size
[docs]
def GetPointRepresentation(listOfFEFields: List[FEField], fillValue: MuscatFloat = 0.0) -> np.ndarray:
"""Get a np.ndarray compatible with the points of the mesh
for a list of FEFields. Each field in the list is treated
as a scalar component for the output field
Parameters
----------
listOfFEFields : List[FEField]
list of FEFields to extract the information
fillValue : MuscatFloat, optional
value to use in the case some field are not defined
over all the nodes, by default 0.
Returns
-------
np.ndarray
the output and nd array of size (nb points, nb of FEField)
"""
nbFields = len(listOfFEFields)
res = np.empty((listOfFEFields[0].mesh.GetNumberOfNodes(), nbFields))
for i, f in enumerate(listOfFEFields):
res[:, i] = f.GetPointRepresentation(fillValue=fillValue)
return res
[docs]
def GetCellRepresentation(listOfFEFields: List[FEField], fillValue: MuscatFloat = 0.0, method="mean") -> np.ndarray:
"""Get a np.ndarray compatible with the points of the mesh
for a list of FEFields. Each field in the list is treated
as a scalar component for the output field
Parameters
----------
listOfFEFields : list
list of FEFields to extract the information
fillValue: float
value to use in the case some field are not defined
over all the nodes, by default 0.
method: str | mean
Read FEField.GetCellRepresentation for more information, by default "mean".
Returns
-------
np.array
Array of size (number of elements, number of fields)
"""
nbFields = len(listOfFEFields)
res = np.empty((listOfFEFields[0].mesh.GetNumberOfElements(), nbFields))
for i, f in enumerate(listOfFEFields):
res[:, i] = f.GetCellRepresentation(fillValue=fillValue, method=method)
return res
[docs]
class IntegrationPointWrapper(FieldBase):
"""Class to generate a FEField at the integration points
Two important function are available : diff and GetIpField.
Parameters
----------
field : FEField
the field to evaluate
rule : _type_
the integration rule
elementFilter : ElementFilter, optional
in this case a RestrictedIPField is generated, by default None
"""
def __init__(self, field: FEField, rule, elementFilter: Optional[ElementFilter] = None):
if not isinstance(field, FEField): # pragma no cover
raise TypeError("IntegrationPointWrapper work only on FEFields")
self.feField = field
self.rule = rule
self._opipCache = None
self._opdiffCache = {}
self.elementFilter = elementFilter
@property
def name(self):
return self.feField.name
[docs]
def ResetCache(self):
self._ipCache = None
self._diffIpCache = {}
[docs]
def diff(self, compName: str) -> Union[IPField, RestrictedIPField]:
"""Get the IPField representation of this field with a derivation in the direction of compName
Parameters
----------
compName : str
compName are available in : from Muscat.FE.SymWeakForm import space
Returns
-------
(IPField | RestrictedIPField)
Depending on the presence of an elementFilter
"""
from Muscat.FE.SymWeakForm import space3D
from sympy import Symbol
for cm in range(3):
if space3D[cm, 0] == compName:
break
elif isinstance(compName, str) and space3D[cm, 0] == Symbol(compName):
break
else:
cm = compName
from Muscat.FE.Fields.FieldTools import TransferFEFieldToIPField
if cm not in self._opdiffCache:
op = GetTransferOpToIPField(self.feField, rule=self.rule, der=cm, elementFilter=self.elementFilter)
self._opdiffCache[cm] = op
res = TransferFEFieldToIPField(self.feField, der=cm, rule=self.rule, elementFilter=self.elementFilter, op=self._opdiffCache[cm])
res.name = "d" + self.feField.name + "d" + str(compName)
return res
[docs]
def GetIpField(self) -> Union[IPField, RestrictedIPField]:
"""Get the integration point representation of the field
Returns
-------
Union[IPField, RestrictedIPField]
Depending on the presence of an elementFilter
"""
from Muscat.FE.Fields.FieldTools import TransferFEFieldToIPField
if self._opipCache is None:
op = GetTransferOpToIPField(self.feField, rule=self.rule, der=-1, elementFilter=self.elementFilter)
self._opipCache = op
return TransferFEFieldToIPField(self.feField, der=-1, rule=self.rule, elementFilter=self.elementFilter, op=self._opipCache)
[docs]
def unaryOp(self, op):
return self.GetIpField().unaryOp(op)
[docs]
def binaryOp(self, other, op):
return self.GetIpField().binaryOp(other, op)
@property
def data(self):
return self.GetIpField().data
DescriptionValue = Union[float, Callable[[np.ndarray], float]]
DescriptionUnit = Tuple[FilterBase, DescriptionValue]
[docs]
def CreateFieldFromDescription(mesh: Mesh, fieldDefinition: Iterable[DescriptionUnit], fieldType: str = "FE") -> Union[FEField, IPField]:
"""
Create a FEField from a field description
Parameters
----------
mesh : Mesh
fieldDefinition : List[Tuple[Union[ElementFilter,NodeFilter],Union[float,Callable[[np.ndarray],float]]]]
A field definition is a list of Tuples (with 2 element each ).
The first element of the tuple is a ElementFilter or a NodeFilter
The second element of the tuple is a float or float returning callable (argument is a point in the space)
fieldType : str, optional
type of field created FEField or IPField
("FE", "FE-P0", "IP") by default "FE" (isoparametric)
Returns
-------
Union[FEField,IPField]
created Field with values coming from the field description
"""
if fieldType == "FE":
from Muscat.FE.FETools import PrepareFEComputation
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, numberOfComponents=1)
res = FEField(mesh=mesh, space=spaces, numbering=numberings[0])
res.Allocate()
FillFEField(res, fieldDefinition)
elif fieldType == "FE-P0":
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
numbering = ComputeDofNumbering(mesh, LagrangeSpaceP0)
res = FEField(mesh=mesh, space=LagrangeSpaceP0, numbering=numbering)
res.Allocate()
FillFEField(res, fieldDefinition)
elif fieldType == "IP":
res = IPField(mesh=mesh, ruleName="LagrangeIsoParamQuadrature")
res.Allocate()
FillIPField(res, fieldDefinition)
else:
raise RuntimeError("fieldType not valid") # pragma no cover
return res
[docs]
def GetTransferOpToIPField(inputField: FEField, ruleName: Optional[str] = None, rule: Optional[MeshQuadrature] = None, der: int = -1, elementFilter: Optional[ElementFilter] = None) -> FeToIPOp:
"""Compute the transfer operator (.dot()) to construct a IPField (of RestrictedIPField)
IPField = FeToIPOp.dot(FEField)
Parameters
----------
inputField : FEField
the input FEField
ruleName : Optional[str], optional
the ruleName of the integration rule to use, by default None
(see Muscat.FE.IntegrationsRules:GetRule for more info)
rule : Optional[Tuple[np.ndarray,np.ndarray]], optional
the integration rule to use, by default None
(see Muscat.FE.IntegrationsRules:GetRule for more info)
der : int, optional
the coordinate to be derived
-1 no derivations only evaluation at IP
[0,1,2] the coordinate number to compute derivative of the FEField
by default -1
elementFilter : Optional[ElementFilter], optional
An element filter to restrict the operator.
In this case a RestrictedIPField is generated
_description_, by default None
Returns
-------
FeToIPOp
An instance of an object with the .dot operator
"""
from Muscat.FE.IntegrationRules import GetIntegrationRule
from Muscat.FE.Spaces.IPSpaces import GenerateSpaceForIntegrationPointInterpolation
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.DofNumbering import ComputeDofNumbering
iRule = GetIntegrationRule(ruleName=ruleName, rule=rule)
gaussSpace = GenerateSpaceForIntegrationPointInterpolation(iRule)
mesh = inputField.mesh
class FeToIPOp(dict):
"""Helper class to hold the transfer matrices between the FEField and the IPField
Parameters
----------
iRule :
The integration rule to be used
elementFilter : ElementFilter
The ElementFilter to generate the RestrictedIPField
"""
def __init__(self, iRule, elementFilter):
super().__init__()
self.iRule = iRule
self.elementFilter = elementFilter
def dot(self, inputField: FEField) -> Union[IPField, RestrictedIPField]:
"""Implementation of the transfer
Parameters
----------
inputField : FEField
input FEField
Returns
-------
Union[IPField,RestrictedIPField]
the type depends on the presence or not of a self.elementFilter
"""
return TransferFEFieldToIPField(inputField, rule=self.iRule, elementFilter=self.elementFilter, op=self)
res = FeToIPOp(iRule, elementFilter=elementFilter)
for elementType in mesh.elements.keys():
if elementFilter is not None:
eF = IntersectionFilter((ElementFilter(elementType=[elementType]), elementFilter))
else:
eF = ElementFilter(elementType=[elementType])
idsToTreat = eF.GetIdsToTreat(mesh, elementType)
if len(idsToTreat) == 0:
continue
numberingRight = ComputeDofNumbering(mesh, Space=gaussSpace, elementFilter=eF)
rightField = FEField(name="Gauss'", numbering=numberingRight, mesh=mesh, space=gaussSpace)
from Muscat.FE.SymWeakForm import GetField, GetTestField, space3D
LF = GetField(inputField.name, 1)
RF = GetTestField("Gauss", 1)
if der == -1:
symForm = LF.T * RF
else:
symForm = LF.diff(space3D[der]).T * RF
transferMatrixMatrix, _ = IntegrateGeneral(
mesh=mesh, constants={}, fields=[], wform=symForm, unknownFields=[inputField], testFields=[rightField], onlyEvaluation=True, integrationRule=iRule, elementFilter=eF
)
res[elementType] = transferMatrixMatrix
return res
[docs]
def TransferFEFieldToIPField(inFEField: FEField, ruleName: Optional[str] = None, rule=None, der: int = -1, elementFilter: Optional[ElementFilter] = None, op=None) -> Union[IPField, RestrictedIPField]:
"""Transfer FEField to a IPField
Parameters
----------
inFEField : FEField
the input FEField
ruleName : Optional[str], optional
the ruleName of the integration rule to use, by default None
(see Muscat.FE.IntegrationsRules:GetRule for more info)
rule : Optional[Tuple[np.ndarray,np.ndarray]], optional
the integration rule to use, by default None
(see Muscat.FE.IntegrationsRules:GetRule for more info)
der : int, optional
the coordinate to be derived
-1 no derivations only evaluation at IP
[0,1,2] the coordinate number to compute derivative of the FEField
by default -1
elementFilter : Optional[ElementFilter], optional
An element filter to restrict the operator.
In this case a RestrictedIPField is generated
_description_, by default None
op : _type_, optional
an object returned by the GetTransferOpToIPField function
if not given then this op is calculated internally
by default None
Returns
-------
Union[IPField,RestrictedIPField]
The field evaluated at the integration points
"""
if op is None:
op = GetTransferOpToIPField(inputField=inFEField, ruleName=ruleName, rule=rule, der=der, elementFilter=elementFilter)
from Muscat.FE.IntegrationRules import GetIntegrationRule
iRule: MeshQuadrature = GetIntegrationRule(ruleName=ruleName, rule=rule)
if elementFilter is None:
outField = IPField(name=inFEField.name, mesh=inFEField.mesh, rule=iRule)
outField.Allocate()
else:
outField = RestrictedIPField(name=inFEField.name, mesh=inFEField.mesh, rule=iRule, elementFilter=elementFilter)
outField.Allocate()
mesh = inFEField.mesh
for elementType, d in mesh.elements.items():
if elementType not in op:
continue
nbElements = op[elementType].shape[0] // iRule[elementType].GetNumberOfPoints()
data = op[elementType].dot(inFEField.data)
outField.data[elementType] = np.reshape(data, (nbElements, iRule[elementType].GetNumberOfPoints()), "F")
outField.data[elementType] = np.ascontiguousarray(outField.data[elementType])
return outField
[docs]
def TransferPosToIPField(mesh: Mesh, ruleName: Optional[str] = None, rule: Optional[Tuple[np.ndarray, np.ndarray]] = None, elementFilter: Optional[ElementFilter] = None) -> List[IPField]:
"""
Create (2 or 3) IPFields with the values of the space coordinates
Parameters
----------
mesh : Mesh
ruleName : str, optional
The Integration rule name, by default None ("LagrangeIsoParamQuadrature")
rule : Tuple[np.ndarray,np.ndarray], optional
The Integration Rule, by default None ("LagrangeIsoParamQuadrature")
elementFilter : ElementFilter, optional
the zone where the transfer must be applied, by default None (all the elements)
Returns
-------
List[IPField]
The IPFields containing the coordinates
"""
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
numbering = ComputeDofNumbering(mesh, LagrangeSpaceGeo, fromConnectivity=True)
inputField = FEField(mesh=mesh, space=LagrangeSpaceGeo, numbering=numbering, data=None)
op = GetTransferOpToIPField(inputField=inputField, ruleName=ruleName, rule=rule, der=-1, elementFilter=elementFilter)
d = mesh.nodes.shape[1]
outField = []
for c, name in enumerate(["posx", "posy", "posz"][0:d]):
inputField.data = mesh.nodes[:, c]
f = op.dot(inputField)
f.name = name
outField.append(f)
return outField
[docs]
def FillIPField(field: IPField, fieldDefinition) -> None:
"""
Fill a IPField using a definition
Parameters
----------
field : IPField
IPField to treat
fieldDefinition : List[Tuple[ElementFilter,Union[float,Callable[[np.ndarray],float]]]]
A field definition is a list of Tuples (with 2 element each ).
The first element of the tuple is a ElementFilter or a NodeFilter
The second element of the tuple is a float or float returning callable (argument is a point in the space)
"""
for f, val in fieldDefinition:
if callable(val):
fValue = val
else:
fValue = lambda x: val
if isinstance(f, ElementFilter):
for selection in f(field.mesh):
geoSpace = LagrangeSpaceGeo[selection.elementType]
rule = field.GetRuleFor(selection.elementType)
nbIp = rule.GetNumberOfPoints()
geoSpaceIpValues = geoSpace.GetSpaceEvaluatedAt(rule.points, rule.weights)
for elementsId in selection.indices:
for i in range(nbIp):
valN = geoSpaceIpValues.valN[i]
elementNodesCoordinates = field.mesh.nodes[selection.elements.connectivity[elementsId, :], :]
pos = np.dot(valN, elementNodesCoordinates).T
field.data[selection.elementType][elementsId, i] = fValue(pos)
continue
raise ValueError(f"Cant use this type of filter to fill an IPField : {type(f)}") # pragma no cover
[docs]
def FillFEField(field: FEField, fieldDefinition) -> None:
"""Fill a FEField using a definition
Parameters
----------
field : FEField
FEField to treat
fieldDefinition : List[Tuple[Union[ElementFilter, NodeFilter], Union[float, Callable[[np.ndarray], float]]]]
A field definition is a list of Tuples (with 2 element each ).
The first element of the tuple is a ElementFilter or a NodeFilter
The second element of the tuple is a float or float returning callable (argument is a point in the space)
"""
for f, val in fieldDefinition:
needPos = False
if callable(val):
fValue = val
needPos = True
else:
fValue = lambda x: val
if isinstance(f, (ElementFilter, FilterOperatorBase)):
treated = np.zeros(field.numbering.size, dtype=bool)
nodes = field.mesh.nodes
for selection in f(field.mesh):
elementType = selection.elements.elementType
numbering = field.numbering[elementType]
if needPos == False:
idDofs = numbering[selection.indices, :].flatten()
field.data[idDofs] = fValue(None)
else:
sp = field.space[elementType]
nbShapeFunctions = sp.GetNumberOfShapeFunctions()
geoSpace = LagrangeSpaceGeo[elementType]
geoSpaceIpValues = geoSpace.GetSpaceEvaluatedAt(sp.posN, np.ones(nbShapeFunctions))
valN = geoSpaceIpValues.valN
connectivity = selection.elements.connectivity
for elementId in selection.indices:
elementNodesCoordinates = nodes[connectivity[elementId, :], :]
dofs = numbering[elementId, :]
pos = np.dot(valN, elementNodesCoordinates)
for i in range(nbShapeFunctions):
if treated[dofs[i]]:
continue
field.data[dofs[i]] = fValue(pos[i, :])
treated[dofs[i]] = True
# for i in range(nbShapeFunctions):
# dofId = field.numbering[name][elementId,i]
# if treated[dofId]:
# continue
# treated[dofId] =True
# valN = geoSpaceIpValues.valN[i]
# pos = np.dot(valN ,elementNodesCoordinates).T
# field.data[dofId] = fValue(pos)
elif isinstance(f, NodeFilter):
ids = f.GetNodesIndices(field.mesh)
for pid in ids:
dofId = field.numbering.GetDofOfPoint(pid)
pos = field.mesh.nodes[pid, :]
field.data[dofId] = fValue(pos)
else: # pragma no cover
raise ValueError(f"Don't know how to treat this type of field {type(f)}")
[docs]
def FieldsAtIp(listOfFields: List[Union[FEField, IPField, RestrictedIPField]], rule, elementFilter: Optional[ElementFilter] = None) -> List[Union[IntegrationPointWrapper, IPField, RestrictedIPField]]:
"""Convert a list of Field (FEFields,IPField,RestrictedIPField) to a list of IPFields
Parameters
----------
listOfFields : List[Union[IPField,RestrictedIPField]]
the list of fields to be treated
rule : _type_
the integration rule to be used for the evaluation
elementFilter : Optional[ElementFilter], optional
the filter to select the element to treat, by default None
Returns
-------
List[Union[IntegrationPointWrapper,IPField,RestrictedIPField]]
the list of IPFilters
Raises
------
Exception
in the case of a internal error
"""
from Muscat.FE.Fields.FieldTools import IntegrationPointWrapper
res = []
for f in listOfFields:
if isinstance(f, FEField):
res.append(IntegrationPointWrapper(f, rule, elementFilter=elementFilter))
elif isinstance(f, IPField):
if f.rule == rule:
if elementFilter is None:
res.append(f)
else:
res.append(f.GetRestrictedIPField(elementFilter))
else: # pragma no cover
raise ValueError(f"The integration rule of the field {f.name} not compatible with the rule provided")
else: # pragma no cover
raise Exception(f"Don't know how to treat this type of field {type(f)}")
return res
[docs]
class FieldsMeshTransportation:
"""Class to help the transfer of field (FEField and IPField) between old and new meshes"""
def __init__(self, oldMesh, newMesh):
self.cache_numbering = {}
self.oldMesh = oldMesh
self.newMesh = newMesh
[docs]
def ResetCacheData(self):
self.cache_numbering = {}
[docs]
def GetNumbering(self, mesh, space, fromConnectivity=False, discontinuous=False):
meshId = str(id(mesh))
spaceId = str(id(space))
key = (meshId, spaceId, fromConnectivity, discontinuous)
if key not in self.cache_numbering:
self.cache_numbering[key] = ComputeDofNumbering(mesh, space, fromConnectivity=fromConnectivity, discontinuous=discontinuous)
return self.cache_numbering[key]
[docs]
def TransportFieldToNewMesh(self, inputField: Union[FEField, IPField]) -> Union[FEField, IPField]:
"""Function to create a Field on the new mesh, the new mesh must be a
transformation of the mesh in the inputField. This means the new mesh originalIds
(for nodes and elements) must be with respect to the mesh of the inputField
Parameters
----------
inputField : FEField or a IPField
the FEField to be transported to the new mesh
Returns
-------
FEField or IPField
a FEField or a IPField (depending on the inputField Type) on the new mesh filled with the values of the input field
"""
if id(inputField.mesh) == id(self.newMesh):
return inputField
if isinstance(inputField, IPField):
outputData = {}
for elemType, data in self.newMesh.elements.items():
outputData[elemType] = inputField.data[elemType][data.originalIds, :]
res = IPField(name=inputField.name, mesh=self.newMesh, rule=inputField.rule, data=outputData)
return res
if isinstance(inputField, FEField):
name = inputField.name
space = inputField.space
if inputField.numbering.fromConnectivity:
numbering = ComputeDofNumbering(self.newMesh, space, fromConnectivity=True)
res = FEField(name=name, mesh=self.newMesh, space=space, numbering=numbering)
res.data = inputField.data[self.newMesh.originalIDNodes]
else:
numbering = ComputeDofNumbering(self.newMesh, space)
res = FEField(name=name, mesh=self.newMesh, space=space, numbering=numbering)
res.Allocate()
for name, data in self.newMesh.elements.items():
res.data[numbering[name].flatten()] = inputField.data[inputField.numbering[name][data.originalIds, :].flatten()]
return res
raise TypeError(f" input field cant be of type {type()}") # pragma no cover
[docs]
class FieldsEvaluator:
"""Helper to evaluate expression using FEField and IPFields.
The user can add fields and constants.
Then this class can be used to compute an expression (given using an callable object).
The callable can capture some of the argument but the last must be `**args` (so extra argument are ignored).
"""
def __init__(self, fields=None):
self.originals = {}
from Muscat.FE.IntegrationRules import GetIntegrationRuleByName
self.rule = GetIntegrationRuleByName("LagrangeIsoParamQuadrature")
self.ruleAtCenter = GetIntegrationRuleByName("ElementCenter")
self.atIp = {}
self.atCenter = {}
if fields is not None:
for f in fields:
self.AddField(f)
self.constants = {}
self.elementFilter = None
self.modified = True
[docs]
def AddField(self, field: FieldBase):
"""Add a field to the internal almanac
Parameters
----------
field : FieldBase
A field to be added to the almanac
"""
self.originals[field.name] = field
[docs]
def AddConstant(self, name: str, val: MuscatFloat):
"""Add a constant value to the almanac
Parameters
----------
name : str
name of the value
val : MuscatFloat
Value
"""
self.constants[name] = val
[docs]
def Update(self, what: str = "all"):
"""Function to update the field calculated at different support
(a FEField can be internally evaluated at the integration points)
The objective of this function is to update this evaluation
Parameters
----------
what : str, optional
["all", "IPField", "Centroids"], by default "all"
"""
for name, field in self.originals.items():
if what == "all" or what == "IPField":
self.atIp[name] = FieldsAtIp([field], self.rule, elementFilter=self.elementFilter)[0]
if what == "all" or what == "Centroids":
self.atCenter[name] = FieldsAtIp([field], self.ruleAtCenter, elementFilter=self.elementFilter)[0]
[docs]
def GetFieldsAt(self, on: str) -> Dict:
"""Get the internal representations of the user fields ant different support
Parameters
----------
on : str
["IPField", "Centroids", "FEField", "Nodes"]
Returns
-------
Dict
A dictionary containing the different representation of the user fields
Raises
------
ValueError
In the case the "on" parameter is not valid
"""
if on == "IPField":
res = self.atIp
elif on == "Nodes":
res = {}
for f in self.originals.values():
res[f.name] = f.GetPointRepresentation()
elif on == "Centroids":
res = self.atCenter
elif on == "FEField":
res = self.originals
else:
raise ValueError(f"Target support not supported ({on})") # pragma no cover
result = dict(self.constants)
result.update(res)
return result
[docs]
def GetOptimizedFunction(self, func):
import sympy
import inspect
from Muscat.FE.SymWeakForm import space3D
class OptimizedFunction:
"""Internal function to convert a function written in python into a
compiled sympy function.
"""
def __init__(self, func, constants):
self.constants = dict(constants)
# get arguments names
args = inspect.getfullargspec(func).args
# clean self, in the case of a member function
if args[0] == "self":
args.pop(0)
symbolicSymbols = {}
self.args = args
for n in args:
if n in self.constants:
symbolicSymbols[n] = self.constants[n]
else:
symbolicSymbols[n] = sympy.Function(n)(*space3D)
# symbolically evaluation of the function
funcValues = func(**symbolicSymbols)
repl, residual = sympy.cse(funcValues)
restKeys = list(symbolicSymbols.keys())
restKeys.extend(["x", "y", "z"])
self.auxFunctions = []
self.auxNames = []
for i, v in enumerate(repl):
funcLam = sympy.lambdify(restKeys, v[1])
self.auxFunctions.append(funcLam)
funcName = str(v[0])
self.auxNames.append(funcName)
restKeys.append(str(v[0]))
self.mainFunc = sympy.lambdify(restKeys, residual)
def __call__(self, **args):
numericFields = {"x": 0, "y": 0, "z": 0}
class Toto:
def __init__(self, obj):
self.obj = obj
def __call__(self, *args, **kwds):
return self.obj
for n in self.args:
if n not in args:
numericFields[n] = self.constants[n]
else:
numericFields[n] = Toto(args[n])
for name, func in zip(self.auxNames, self.auxFunctions):
numericFields[name] = func(**numericFields)
return self.mainFunc(**numericFields)[0]
return OptimizedFunction(func, self.constants)
[docs]
def Compute(self, func: Callable, on: str, useSympy: bool = False):
"""Compute the function func using the user fields on the support "on"
Parameters
----------
func : Callable
the function to evaluate
on : str
["IPField", "Centroids", "FEField", "Nodes"]
useSympy : bool, optional
if true the function is compiled using sympy , by default False
Returns
-------
Any
The callable evaluated
"""
if useSympy:
func = self.GetOptimizedFunction(func)
fields = self.GetFieldsAt(on)
return func(**fields)
[docs]
def ComputeTransferOp(field1: FEField, field2: FEField, force: bool = False) -> Tuple[np.ndarray, np.ndarray]:
"""
Compute the transfer Operator from field1 to field2
Both fields must be compatibles fields: same mesh, same space
but with different numberings
- example of use:
lIndex,rIndex = ComputeTransferOp(field1,field2):
field2.data[lIndex] = field1.data[rIndex]
Parameters
----------
field1 : FEField
field from where the information is extracted
field2 : FEField
field to receive the information
force : bool, optional
true to bypass the compatibility check, by default False
Returns
-------
Tuple(np.ndarray,np.ndarray)
index vector for field2, index vector for field1
"""
if field1.mesh != field2.mesh and not force: # pragma no cover
raise Exception("The fields must share the same mesh")
if field1.space != field2.space and not force: # pragma no cover
raise Exception("The fields must share the same space")
_extractorLeft = []
_extractorRight = []
for name in field1.mesh.elements.keys():
a = field1.numbering.get(name, None)
b = field2.numbering.get(name, None)
if a is not None and b is not None:
mask = np.logical_and(a >= 0, b >= 0)
_extractorLeft.extend(b[mask])
_extractorRight.extend(a[mask])
extractorLeft = np.array(_extractorLeft, dtype=MuscatIndex)
extractorRight = np.array(_extractorRight, dtype=MuscatIndex)
v, index = np.unique(extractorLeft, return_index=True)
extractorLeft = v
extractorRight = extractorRight[index]
return extractorLeft, extractorRight
# def CheckIntegrity_ElementWiseIpToFETransferOp(GUI:bool=False):
# ElementWiseIpToFETransferOp()
# return "ok"
[docs]
def CheckIntegrity_ElementWiseFEToFETransferOp(GUI: bool = False):
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, LagrangeSpaceP2
ElementWiseFEToFETransferOp(LagrangeSpaceP1, LagrangeSpaceP2)
return "ok"
[docs]
def CheckIntegrity_FieldsAtIp(GUI: bool = False):
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
from Muscat.FE.IntegrationRules import GetIntegrationRuleByName
mesh = CreateUniformMeshOfBars(0, 1, 10)
field1 = CreateFieldFromDescription(mesh, [(ElementFilter(), 1)])
field2 = CreateFieldFromDescription(mesh, [(ElementFilter(), 2)], fieldType="FE-P0")
field3 = CreateFieldFromDescription(mesh, [(ElementFilter(), 2)], fieldType="IP")
FieldsAtIp([field1, field2, field3], GetIntegrationRuleByName("LagrangeIsoParamQuadrature"))
FieldsAtIp([field1, field2, field3], GetIntegrationRuleByName("LagrangeIsoParamQuadrature"), elementFilter=ElementFilter())
return "ok"
[docs]
def CheckIntegrity_Legacy(GUI: bool = False):
CheckIntegrity_ElementWiseFEToFETransferOp(GUI=GUI)
from Muscat.Containers.MeshCreationTools import CreateUniformMeshOfBars
mesh = CreateUniformMeshOfBars(0, 1, 10)
bars = mesh.GetElementsOfType(ED.Bar_2)
bars.tags.CreateTag("first3").SetIds([0, 1, 2])
bars.tags.CreateTag("next3").SetIds([3, 4, 5])
bars.tags.CreateTag("Last").SetIds([8])
mesh.nodesTags.CreateTag("FirstPoint").SetIds([0])
mesh.nodesTags.CreateTag("LastPoint").SetIds([9])
print(mesh)
fieldDefinition = [(ElementFilter(), 5)]
fieldDefinition.append((ElementFilter(eTag=["first3"]), 3))
fieldDefinition.append((ElementFilter(eTag=["Last"]), -1))
fieldDefinition.append((ElementFilter(eTag=["next3"]), lambda x: x[0]))
field = CreateFieldFromDescription(mesh, fieldDefinition, fieldType="IP")
print(field)
print(field.data)
field = CreateFieldFromDescription(mesh, fieldDefinition, fieldType="FE-P0")
print(field)
print(field.data)
fieldDefinition.append((NodeFilter(nTag=["FirstPoint"]), -10))
fieldDefinition.append((NodeFilter(nTag=["LastPoint"]), lambda x: x[0] + 1.2))
field = CreateFieldFromDescription(mesh, fieldDefinition)
field.name = "tempFieldName"
print(field.data)
field.Allocate(val=0)
FillFEField(field, fieldDefinition)
print("IntegrationPointWrapper --*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-")
from Muscat.FE.IntegrationRules import GetIntegrationRuleByName
ipw = IntegrationPointWrapper(field, GetIntegrationRuleByName("LagrangeIsoParamQuadrature"), elementFilter=ElementFilter(dimensionality=1))
np.testing.assert_equal(ipw.name, field.name)
ipw.ResetCache()
from sympy import Symbol
ipw.diff(Symbol("x"))
ipw.diff("x") ## derivative with respect to X (0)
ipw.diff(0) ## derivative with respect to X (0)
ipw.GetIpField()
ipw.GetIpField()
ipw.data
a = ipw + 1
b = -ipw
print("ops on fields --*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-")
fmax = Maximum(field, 1)
fmin = Minimum(field, 1)
print("--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-")
nodalTransferredField = TransferFEFieldToIPField(field, ruleName="LagrangeIsoParamQuadrature", elementFilter=ElementFilter(dimensionality=1, eTag="next3"))
print("--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*--")
print("input")
print(field.data)
print("output")
print(nodalTransferredField.data)
print(nodalTransferredField.GetIpFieldRepresentation(1).data)
res = TransferPosToIPField(mesh, ruleName="LagrangeIsoParamQuadrature", elementFilter=ElementFilter(dimensionality=1, eTag="next3"))
print(res[0].data)
print(res[0].GetIpFieldRepresentation().data)
field.name = "E"
FEval = FieldsEvaluator(fields=[field])
FEval.AddConstant("alpha", 0)
def op(E, alpha, **args):
return (2 * E + 1 + alpha) + (2 * E + 1 + alpha) ** 2 + (2 * E + 1 + alpha) ** 3
print("--------")
print(field.data)
print("--------")
res = FEval.Compute(op, "FEField")
print(res.data)
FEval.Update("IPField")
res = FEval.Compute(op, "IPField")
print(res.data)
FEval.Update("Centroids")
FEval.GetFieldsAt("IPField")
FEval.GetFieldsAt("Nodes")
FEval.GetFieldsAt("Centroids")
FEval.GetFieldsAt("FEField")
FEField
res = FEval.Compute(op, "FEField", useSympy=True)
print(res.data)
mesh.nodeFields["FEField_onPoints"] = GetPointRepresentation(
[
res,
]
)
mesh.nodeFields["FEField_onPoints_Vector"] = GetPointRepresentation(
[
res,
res,
]
)
mesh.elemFields["FEField_onElements"] = GetCellRepresentation(
[
res,
]
)
mesh.elemFields["FEField_onElements_Vector"] = GetCellRepresentation(
[
res,
res,
]
)
print(NodeFieldToFEField(mesh))
print(NodeFieldToFEField(mesh, copy=True))
print(ElemFieldsToFEField(mesh))
print(ElemFieldsToFEField(mesh, copy=True))
vector = FEFieldsDataToVector([res])
vector = FEFieldsDataToVector([res], copy=True)
VectorToFEFieldsData(vector, [res])
res = FEval.GetOptimizedFunction(op)
# ---- TransportFieldToNewMesh -------------
fieldIP = CreateFieldFromDescription(mesh, [(ElementFilter(), 5)], fieldType="IP")
fieldP0 = CreateFieldFromDescription(mesh, [(ElementFilter(), 5)], fieldType="FE-P0")
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
mesh2 = ExtractElementsByElementFilter(field.mesh, ElementFilter(dimensionality=1))
obj = FieldsMeshTransportation(field.mesh, mesh2)
obj.TransportFieldToOldMesh(field)
fieldNew = obj.TransportFieldToNewMesh(field)
almostField = obj.TransportFieldToOldMesh(fieldNew)
obj.TransportFieldToNewMesh(fieldNew)
fieldIPNew = obj.TransportFieldToNewMesh(fieldIP)
almostFieldIP = obj.TransportFieldToOldMesh(fieldIPNew)
fieldP0New = obj.TransportFieldToNewMesh(fieldP0)
almostfieldP0 = obj.TransportFieldToOldMesh(fieldP0New)
obj.ResetCacheData()
extractorLeft, extractorRight = ComputeTransferOp(field, field)
# basic test if the same field, the extrators are equals
np.testing.assert_equal(extractorLeft, extractorRight)
return "ok"
[docs]
def CheckIntegrity(GUI: bool = False):
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
listToCheck = [
# CheckIntegrity_ElementWiseIpToFETransferOp,
CheckIntegrity_ElementWiseFEToFETransferOp,
CheckIntegrity_FieldsAtIp,
CheckIntegrity_Legacy,
]
return RunListOfCheckIntegrities(listToCheck)
if __name__ == "__main__":
print(CheckIntegrity(GUI=True)) # pragma no cover