# -*- 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 Dict, Any, Callable, Union, Optional
import numpy as np
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
from Muscat.FE.Fields.FieldBase import FieldBase
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Filters.FilterOperators import IntersectionFilter
from Muscat.MeshContainers.Filters.FilterTools import GetFrozenFilter
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.FE.IntegrationRules import MeshQuadrature, ElementaryQuadrature, GetIntegrationRule
from Muscat.MeshContainers.Mesh import Mesh
[docs]
class IPField(FieldBase):
"""Class to hold information at the integration points"""
def __init__(
self, name: Optional[str] = None, mesh: Optional[Mesh] = None, rule: Optional[MeshQuadrature] = None, ruleName: Optional[str] = None, data: Optional[Dict[ED.ElementType, ArrayLike]] = None
):
"""Constructor
Parameters
----------
name : str, optional
The name of the field, by default None
mesh : Mesh, optional
the mesh , by default None
rule : MeshQuadrature, optional
the integration rule to use, by default None
ruleName : str, optional
the name of the integration rule to use, by default None
data : Dict[str,ArrayLike], optional
the data per element type, the ArrayLike must be of size (number of elements, number of integration points), by default None
"""
super().__init__(name=name, mesh=mesh)
if data is None:
self.data: Dict[ED.ElementType, np.ndarray] = {}
else:
self.data = data
self.rule: MeshQuadrature = GetIntegrationRule(ruleName=ruleName, rule=rule)
[docs]
def SetRule(self, ruleName: Optional[str] = None, rule: Optional[MeshQuadrature] = None):
"""Set the integration rule to be used by the field
Please read the constructor documentation for more information about the arguments
"""
self.rule = GetIntegrationRule(ruleName=ruleName, rule=rule)
[docs]
def GetRuleFor(self, elementType: ED.ElementType) -> ElementaryQuadrature:
"""Helper function to get the integration rule for one type of elements
Parameters
----------
elementType : str
the element name
Returns
-------
ElementIntegrationRuleType
The integration rule, this is a tuple with the positions and the weight of the integration rule
"""
return self.rule[elementType]
[docs]
def GetDataFor(self, elementType: ED.ElementType) -> np.ndarray:
"""Get the numpy array storing the data for one type of element
Parameters
----------
elementType : str
the element name
Returns
-------
np.ndarray
the data of size (number of elements, number of integration points)
"""
return self.data[elementType]
[docs]
def Allocate(self, val: MuscatFloat = 0.0):
"""Allocate the memory to store the data
Parameters
----------
val : MuscatFloat, optional
initial fill value, by default 0.
"""
self.data = dict()
for name, data in self.mesh.elements.items():
nbIntegrationPoints = self.GetRuleFor(name).GetNumberOfPoints()
nbElements = data.GetNumberOfElements()
self.data[name] = np.zeros((nbElements, nbIntegrationPoints), dtype=MuscatFloat) + val
[docs]
def CheckCompatibility(self, B: Any) -> None:
"""Check of two field are compatible to make arithmetic operations
Parameters
----------
B : Any
the other IpField
Raises
------
Exception
In the case the fields are not compatible
"""
if isinstance(B, type(self)):
if id(self.mesh) != id(B.mesh): # pragma: no cover
raise Exception("The support of the fields are not the same")
if id(self.rule) != id(B.rule): # pragma: no cover
raise Exception("The rules of the fields are not the same")
[docs]
def unaryOp(self, op: Callable, out: Optional[IPField] = None) -> IPField:
"""Internal function to apply a unary operation
Parameters
----------
op : Callable
a function to apply to the data
out : IPField, optional
output IPFiled if non is provided a new one is created, by default None
Returns
-------
IPField
The output IPField
"""
if out is None:
res = type(self)(name=None, mesh=self.mesh, rule=self.rule)
else:
res = out
res.data = {key: op(self.data[key]) for key in self.data.keys()}
return res
[docs]
def binaryOp(self, other: Any, op: Callable, out: Optional[IPField] = None) -> IPField:
"""Internal function to apply a binary operator. A + B for example
Parameters
----------
other : Any
a second IPField
op : Callable
the binary operator to apply
out : IPField, optional
output IPFiled if Non is provided a new one is created, by default None, by default None
Returns
-------
IPField
The output IPField
Raises
------
Exception
if the field are not compatible
Exception
if the binary operator cant be apply correctly
"""
self.CheckCompatibility(other)
if out is None:
res = type(self)(name=None, mesh=self.mesh, rule=self.rule)
else:
res = out
if isinstance(other, (type(self), RestrictedIPField)):
res.data = {key: op(self.data[key], other.data[key]) for key in set(self.data.keys()).union(other.data.keys())}
return res
elif type(other).__module__ == np.__name__ and np.ndim(other) != 0:
res = np.empty(other.shape, dtype=object)
for res_data, other_data in np.nditer([res, other], flags=["refs_ok"], op_flags=["readwrite"]):
res_data[...] = op(self, other_data)
return res
elif np.isscalar(other):
res.data = {key: op(data, other) for key, data in self.data.items()}
return res
else:
raise Exception(f"operator {op} not valid for types :{type(self)} and {type(other)} ") # pragma: no cover
[docs]
def GetCellRepresentation(self, fillValue: MuscatFloat = 0.0, method: Union[str, int] = "mean") -> np.ndarray:
"""Function to push the data from the field into a vector homogeneous to
the mesh (for visualization for example).
Parameters
----------
fillValue : MuscatFloat, optional
field value for cell without data, by default 0.
method : Union[str,int], optional
the method to pass the information from the integration points
to the cells:
'mean' : compute the mean value in every cell
'max' : extract the max value for every cell
'min' : extract the min value for every cell
'maxDiff': compute the maximal difference for every cell max(abs(vi-vj)) for i!=j in range(number of integration point)
'maxDiffFraction': same as before but divided by the mean value
'int' or, int: a int select a specific integration point the value is clipped to [0,number of integration point[
, by default 'mean'
Returns
-------
np.ndarray
np.ndarray of size (number of elements in the mesh, 1) with the extracted information
"""
if fillValue == 0.0:
res = np.zeros(self.mesh.GetNumberOfElements(), dtype=MuscatFloat)
else:
res = np.full(self.mesh.GetNumberOfElements(), fillValue, dtype=MuscatFloat)
cpt = 0
for name, data in self.mesh.elements.items():
nbElements = data.GetNumberOfElements()
if name not in self.data:
cpt += nbElements
continue
if method == "mean":
data = np.mean(self.data[name], axis=1)
elif method == "max":
data = np.max(self.data[name], axis=1)
elif method == "min":
data = np.min(self.data[name], axis=1)
elif method == "maxDiff" or method == "maxDiffFraction":
cols = self.data[name].shape[1]
op = np.zeros((cols, (cols * (cols - 1)) // 2))
iCpt = 0
for i in range(0, cols - 1):
for j in range(i + 1, cols):
op[i, iCpt] = 1
op[j, iCpt] = -1
iCpt += 1
data = np.max(abs(self.data[name].dot(op)), axis=1)
if method == "maxDiffFraction":
data /= np.mean(self.data[name], axis=1)
else:
col = min(int(method), self.data[name].shape[1])
data = self.data[name][:, col]
res[cpt : cpt + nbElements] = data
cpt += nbElements
return res
[docs]
def GetPointRepresentation(self, fillValue: MuscatFloat = 0.0, elementFilter: ElementFilter = None, weightScheme="Multiplicity") -> np.ndarray:
"""Generate a point representation of this field, only information from the highest element dimensionality is used.
Use a custom element filter to override this behavior.
Parameters
----------
fillValue : float, optional
value to be used on nodes with no data available, by default 0.
elementFilter : ElementFilter, optional
elements to construct the point representation, by default None
weightScheme: "Multiplicity" of "Area". Choose the strategy to pass information to points.
"Multiplicity", the weight is one for every element contributing to a point
"Area", the area of the elements is used as weight
Returns
-------
np.ndarray
a numpy vector of size mesh.GetNumberOfNodes with the computed values
"""
if elementFilter is None:
elementFilter = ElementFilter(dimensionality=self.mesh.GetElementsDimensionality())
res = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatFloat)
unusedPoints = np.ones(self.mesh.GetNumberOfNodes(), dtype=bool)
totalWeight = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatFloat)
if weightScheme not in ["Multiplicity","Area"]:
raise RuntimeError("weightScheme must be 'Multiplicity' or 'Area' ")
if weightScheme == "Area":
from Muscat.MeshTools.MeshInspectionTools import GetVolumePerElement
volumes = GetVolumePerElement(self.mesh,elementFilter)
else:
localWeight = 1.
from Muscat.FE.Fields.FieldTools import ElementWiseIpToFETransferOp
interpolator = ElementWiseIpToFETransferOp(self.rule, LagrangeSpaceGeo)
for selection in elementFilter(self.mesh):
if weightScheme == "Area":
localWeight = volumes[selection.GetSelectionSlice()]
if selection.elementType not in interpolator or selection.elementType not in self.data:
continue
data = self.data[selection.elementType][selection.indices, :].T
for i in range(selection.elements.GetNumberOfNodesPerElement()):
point_ids = selection.elements.connectivity[selection.indices, i]
res[point_ids] += interpolator[selection.elementType][i, :].dot(data)* localWeight
unusedPoints[point_ids] = False
totalWeight[point_ids] += localWeight
totalWeight[unusedPoints] = 1.0
res[unusedPoints] = fillValue
res /= totalWeight
return res
[docs]
def Flatten(self, elementFilter: ElementFilter = None) -> np.ndarray:
"""Flatten all the data into a single ndarray
see also: self.SetDataFromNumpy
Parameters
----------
elementFilter : ElementFilter, optional
ElementFilter to operate, by default (None) all element are treated
Returns
-------
np.ndarray
the extracted data
"""
if elementFilter is None:
elementFilter = ElementFilter()
nbValues = 0
for selection in elementFilter(self.mesh):
nbValues += np.prod(self.data[selection.elementType].shape)
res = np.empty(nbValues, dtype=MuscatFloat)
cpt = 0
for selection in elementFilter(self.mesh):
localSize = np.prod(self.data[selection.elementType].shape)
res[cpt : cpt + localSize] = self.data[selection.elementType].flatten()
cpt += localSize
return res
[docs]
def SetDataFromNumpy(self, inData: ArrayLike, elementFilter: ElementFilter = None) -> None:
"""Set the internal data from a one dimensional array
The user must allocate the memory before calling this function
see also: self.Flatten
Parameters
----------
inData : ArrayLike
the incoming data
elementFilter : ElementFilter, optional
ElementFilter to operate, by default (None) all element are treated
"""
if elementFilter is None:
elementFilter = ElementFilter()
nbValues = 0
for selection in elementFilter(self.mesh):
nbValues += np.prod(self.data[selection.elementType].shape)
if len(inData.shape) != 1: # pragma: no cover
raise RuntimeError("inData must be a vector")
if inData.size != nbValues: # pragma: no cover
raise (Exception("incompatible sizes"))
cpt = 0
for selection in elementFilter(self.mesh):
localSize = np.prod(self.data[selection.elementType].shape)
self.data[selection.elementType][:, :] = inData[cpt : cpt + localSize].reshape(self.data[selection.elementType].shape)
cpt += localSize
[docs]
def GetRestrictedIPField(self, elementFilter: ElementFilter) -> RestrictedIPField:
"""Create a RestrictedIPField only on element selected by the elementFilter
Parameters
----------
elementFilter : ElementFilter
ElementFilter to work on
Returns
-------
RestrictedIPField
a new RestrictedIPField
"""
res = RestrictedIPField(name=self.name, mesh=self.mesh, rule=self.rule, elementFilter=elementFilter)
res.AllocateFromIpField(self)
return res
def __str__(self) -> str:
"""Nice string representations
Returns
-------
str
a string in the form of "<IPField object 'name' (id(self))>"
"""
return f"<IPField object '{self.name}' ({id(self)})>"
[docs]
class RestrictedIPField(IPField):
"""Class to hold information at integration point on a subdomain of the mesh (ElementFilter)
This class is experimental and for the moment not compatible with the integration module
"""
def __init__(self, name: str = None, mesh: Mesh = None, rule: MeshQuadrature = None, ruleName: str = None, data: Dict[str, ArrayLike] = None, elementFilter: ElementFilter = None) -> None:
"""Constructor
Parameters
----------
name : str, optional
The name of the field, by default None
mesh : Mesh, optional
the mesh , by default None
rule : MeshQuadratureonal
the integration rule to use, by default None
ruleName : str, optional
the name of the integration rule to use, by default None
data : Dict[str,ArrayLike], optional
the data per element type, the ArrayLike must be of size (number of elements, number of integration points), by default None
elementFilter : ElementFilter, optional
the filter to select the element with data
"""
super().__init__(name=name, mesh=mesh, rule=rule, ruleName=ruleName, data=data)
if elementFilter == None:
self.elementFilter = ElementFilter()
else:
self.elementFilter = GetFrozenFilter(elementFilter, mesh)
[docs]
def Allocate(self, val: MuscatFloat = 0.0) -> None:
"""Allocate the memory to store the data
Parameters
----------
val : MuscatFloat, optional
initial fill value, by default 0.
"""
self.data = dict()
for selection in self.elementFilter(self.mesh):
nbIntegrationPoints = self.GetRuleFor(selection.elementType).GetNumberOfPoints()
self.data[selection.elementType] = np.zeros((selection.Size(), nbIntegrationPoints), dtype=MuscatFloat) + val
[docs]
def GetRestrictedIPField(self, elementFilter: ElementFilter) -> RestrictedIPField:
"""Create a RestrictedIPField only on element intersection between elementFilter and the internal elementFilter
Parameters
----------
elementFilter : ElementFilter
ElementFilter to work on
Returns
-------
RestrictedIPField
a new RestrictedIPField
"""
res = RestrictedIPField(name=self.name, mesh=self.mesh, rule=self.rule, elementFilter=IntersectionFilter(filters=[elementFilter, self.elementFilter]))
res.AllocateFromIpField(self)
return res
[docs]
def AllocateFromIpField(self, inputIPField: Union[IPField, RestrictedIPField]) -> None:
"""Fill internal data from a external ipField (data is copied)
Parameters
----------
ipField : Union[IPField,RestrictedIPField]
The IPField to extract the data
Raises
------
Exception
In the case the inputIPField is not of the type IPField or RestrictedIPField
"""
self.name = inputIPField.name
self.mesh = inputIPField.mesh
self.rule = inputIPField.rule
self.data = dict()
if isinstance(inputIPField, RestrictedIPField):
for selection in self.elementFilter(self.mesh):
idsII = inputIPField.elementFilter.GetIdsToTreat(self.mesh, selection.elementType)
tempIds = np.empty(selection.elements.GetNumberOfElements(), dtype=MuscatIndex)
tempIds[idsII] = np.arange(len(idsII))
self.data[selection.elementType] = inputIPField.data[selection.elementType][tempIds[selection.indices], :]
elif isinstance(inputIPField, IPField):
for selection in self.elementFilter(self.mesh):
self.data[selection.elementType] = inputIPField.data[selection.elementType][selection.indices, :]
else: # pragma: no cover
raise Exception(f"don't know how to treat ipField of type {type(inputIPField)}")
[docs]
def SetDataFromNumpy(self, inData: ArrayLike) -> None:
"""Set the internal data from a one dimensional array
The user must allocate the memory before calling this function
Parameters
----------
inData : ArrayLike
the incoming data
elementFilter : ElementFilter, optional
ElementFilter to operate, by default (None) all element are treated
"""
nbValues = 0
for selection in self.elementFilter(self.mesh):
nbValues += selection.Size() * self.GetRuleFor(selection.elementType).GetNumberOfPoints()
if inData.size != nbValues: # pragma: no cover
raise RuntimeError("incompatible sizes")
if len(inData.shape) != 1: # pragma: no cover
raise Exception("inData must be a vector")
cpt = 0
for selection in self.elementFilter(self.mesh):
localSize = selection.Size() * self.GetRuleFor(selection.elementType).GetNumberOfPoints()
self.data[selection.elementType] = inData[cpt : cpt + localSize].reshape((selection.Size(), self.GetRuleFor(selection.elementType).GetNumberOfPoints()))
cpt += localSize
[docs]
def GetIpFieldRepresentation(self, fillValue: MuscatFloat = MuscatFloat(0.0)) -> IPField:
"""Get a IPFieldRepresentation. Function to expand the representation to he hold domain
Parameters
----------
fillValue : MuscatFloat, optional
field value for cell without data, by default 0.
Returns
-------
IPField
IPField over the all the element with fillValue on element not covered by the RestrictedIPField
"""
res = IPField(name=self.name, mesh=self.mesh, rule=self.rule)
res.Allocate(fillValue)
for selection in self.elementFilter(self.mesh):
if selection.elementType not in self.data:
continue
res.data[selection.elementType][selection.indices, :] = self.data[selection.elementType]
return res
[docs]
def CheckCompatibility(self, B: Any) -> None:
"""Check of two field are compatible to make arithmetic operations
Parameters
----------
B : Any
the other RestrictedIPField
Raises
------
Exception
In the case the fields are not compatible
"""
super().CheckCompatibility(B)
if isinstance(B, type(self)):
if not self.elementFilter.IsEquivalent(B.elementFilter):
raise Exception("The elementFilter of the fields are not the same")
[docs]
def unaryOp(self, op: Callable, out: RestrictedIPField = None) -> RestrictedIPField:
"""Internal function to apply a unary operation
Parameters
----------
op : Callable
a function to apply to the data
out : RestrictedIPField, optional
output IPFiled if non is provided a new one is created, by default None
Returns
-------
RestrictedIPField
The output RestrictedIPField
"""
res = type(self)(name=None, mesh=self.mesh, rule=self.rule, elementFilter=self.elementFilter)
return super().unaryOp(op, out=res)
[docs]
def binaryOp(self, other: Any, op: Callable, out: RestrictedIPField = None) -> RestrictedIPField:
"""Internal function to apply a binary operator. A + B for example
Parameters
----------
other : Any
a second RestrictedIPField
op : Callable
the binary operator to apply
out : RestrictedIPField, optional
output RestrictedIPFiled, if None is provided a new one is created, by default None, by default None, by default None
Returns
-------
RestrictedIPField
The output RestrictedIPField
"""
res = type(self)(name=None, mesh=self.mesh, rule=self.rule, elementFilter=self.elementFilter)
return super().binaryOp(other, op, out=res)
[docs]
def GetCellRepresentation(self, fillValue: MuscatFloat = 0.0, method: Union[str, int] = "mean") -> np.ndarray:
if fillValue == 0.0:
res = np.zeros(self.mesh.GetNumberOfElements(), dtype=MuscatFloat)
else:
res = np.full(self.mesh.GetNumberOfElements(), fillValue, dtype=MuscatFloat)
cpt = 0
for elementType, elementData in self.mesh.elements.items():
ids = self.elementFilter.GetIdsToTreat(self.mesh, elementType)
nbElements = elementData.GetNumberOfElements()
if elementType not in self.data:
cpt += nbElements
continue
if method == "mean":
data = np.mean(self.data[elementType], axis=1)
elif method == "max":
data = np.max(self.data[elementType], axis=1)
elif method == "min":
data = np.min(self.data[elementType], axis=1)
elif method == "maxDiff" or method == "maxDiffFraction":
cols = self.data[elementType].shape[1]
op = np.zeros((cols, (cols * (cols - 1)) // 2))
iCpt = 0
for i in range(0, cols - 1):
for j in range(i + 1, cols):
op[i, iCpt] = 1
op[j, iCpt] = -1
iCpt += 1
data = np.max(abs(self.data[elementType].dot(op)), axis=1)
if method == "maxDiffFraction":
data /= np.mean(self.data[elementType], axis=1)
else:
col = min(int(method), self.data[elementType].shape[1])
data = self.data[elementType][:, col]
res[cpt + np.asarray(ids, dtype=MuscatIndex)] = data
cpt += nbElements
return res
[docs]
def CheckIntegrity(GUI: bool = False):
from Muscat.FE.IntegrationRules import LagrangeP1Quadrature
from Muscat.MeshTools.MeshCreationTools import CreateCube
mesh = CreateCube([2.0, 3.0, 4.0], [-1.0, -1.0, -1.0], [2.0 / 10, 2.0 / 10, 2.0 / 10])
print(mesh)
sig11 = IPField("Sig_11", mesh=mesh, rule=LagrangeP1Quadrature)
print(sig11.GetCellRepresentation())
sig11.Allocate()
print(sig11)
print(sig11.GetPointRepresentation())
sig22 = sig11 + 0.707107
sig12 = 2 * (-sig22) * 5 / sig22
A = sig11**2
B = sig11 * sig22
C = sig22**2
D = 1.5 * sig12 * 2
E = A - B + C + (D) ** 2
vonMises = np.sqrt(E)
print(vonMises.data)
print("545454")
print(np.linalg.norm([sig22, sig22]).data)
data = sig22.GetCellRepresentation()
methods = ["min", "max", "mean", "maxDiff", "maxDiffFraction", "-1"]
for method in methods:
data = sig22.GetCellRepresentation(method=method)
data = sig22.GetCellRepresentation(1)
sig22.SetDataFromNumpy(sig22.Flatten())
data2 = sig22.GetCellRepresentation()
if np.linalg.norm(data - data2) > 0:
raise () # pragma: no cover
data2 = sig22.GetPointRepresentation(0)
data2 = sig22.GetPointRepresentation(0,weightScheme="Area")
dummyField = IPField()
print(str(dummyField))
dummyField.data[None] = np.arange(3) + 1
print("dummyField")
print(dummyField * dummyField - dummyField**2 / dummyField)
restrictedIPField = sig22.GetRestrictedIPField(ElementFilter(eTag="Skin"))
restrictedIPField = -(2 * restrictedIPField)
print(restrictedIPField)
data2 = restrictedIPField.GetPointRepresentation(0)
restrictedIPField.GetIpFieldRepresentation()
# restrictedIPField2 = restrictedIPField.GetRestrictedIPField(ElementFilter(dimensionality=3))
# restrictedIPField2.Allocate()
restrictedIPField2 = restrictedIPField.GetRestrictedIPField(ElementFilter(eTag="X0"))
# restrictedIPField2.Allocate()
print(restrictedIPField2.data)
restrictedIPField2 = -(2 * sig22.GetRestrictedIPField(ElementFilter(eTag=["X0"])))
print(restrictedIPField2.data)
restrictedIPField2.Allocate(1)
print(restrictedIPField2 * np.array([0.1, 0.3]))
import Muscat.MeshContainers.ElementsDescription as ED
restrictedIPField2.GetDataFor(ED.Quadrangle_4)
obj = RestrictedIPField(mesh=mesh, data={})
obj.SetRule()
print(obj.GetIpFieldRepresentation())
data = obj.GetCellRepresentation(fillValue=1.0)
NumberOfValues = 0
for selection in obj.elementFilter(mesh):
NumberOfValues += selection.Size() * obj.GetRuleFor(selection.elementType).GetNumberOfPoints()
obj.SetDataFromNumpy(np.ones(NumberOfValues, dtype=MuscatFloat))
methods = ["min", "max", "mean", "maxDiff", "maxDiffFraction", "-1"]
data = obj.GetCellRepresentation(fillValue=1.0)
for method in methods:
data = obj.GetCellRepresentation(method=method)
# must fail
error = False
try:
restrictedIPField + restrictedIPField2
error = True # pragma: no cover
except:
pass
if error: # pragma: no cover
raise
return "ok"
if __name__ == "__main__":
print(CheckIntegrity(True)) # pragma: no cover