# -*- 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 Union
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex
from Muscat.Helpers.TextFormatHelper import TFormat
from Muscat.FE.Fields.FieldBase import FieldBase
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
[docs]
class FEField(FieldBase):
def __init__(self,name=None,mesh=None,space=None,numbering=None,data=None):
super().__init__(name=name,mesh = mesh)
self.data = data
if space is None:
space = LagrangeSpaceGeo
self.space = space
if (numbering is None) and (mesh is not None) and (space is not None):
from Muscat.FE.DofNumbering import ComputeDofNumbering
if space is LagrangeSpaceGeo:
numbering = ComputeDofNumbering(mesh, space, fromConnectivity=True)
else:
numbering = ComputeDofNumbering(mesh, space)
self.numbering = numbering
[docs]
def Allocate(self,val=0):
if val == 0:
self.data = np.zeros(self.numbering.size,dtype=MuscatFloat)
else:
self.data = np.full(self.numbering.size,val,dtype=MuscatFloat)
[docs]
def copy(self):
return FEField(name=self.name,mesh=self.mesh,space=self.space,numbering=self.numbering, data=self.data.copy())
[docs]
def GetPointRepresentation(self,fillValue=0):
"""Generates a vector representation of the field data corresponding to the mesh.
This function transfers the field's data into a vector that is compatible in size
with the mesh, which can be useful for visualization purposes. If an entity in the
mesh does not have degrees of freedom (dofs), it is filled with a specified fill value,
which defaults to 0 if not provided.
Parameters
----------
fillValue : float, optional
The value used to fill entities with no degrees of freedom. Default is 0.
Returns
-------
numpy.ndarray
A numpy array representing the field data across the entire mesh. Entities
that originally had no dofs are filled with the specified fill value.
"""
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
if self.space is LagrangeSpaceGeo and self.numbering.fromConnectivity:
return self.data
res = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatFloat)
restemp = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatFloat)
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
## loop over the elements in descending dimensionality order to fill the data and
#mask to track updated points
mask = np.zeros(self.mesh.GetNumberOfNodes(), dtype=bool)
#counter to tack the number of contribution to a point
cpt = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatIndex)
cpttmp = np.zeros(self.mesh.GetNumberOfNodes(), dtype=MuscatIndex)
for dim in [3,2,1,0]:
localMask = np.zeros(self.mesh.GetNumberOfNodes(), dtype=bool)
restemp.fill(0.)
cpttmp.fill(0)
for selection in ElementFilter(dimensionality=dim).IterOnTypesOnly(self.mesh, force=True):
elementType = selection.elements.elementType
selfNumbering = self.numbering[elementType]
if selfNumbering is None:
continue
geoElSpace = LagrangeSpaceGeo[elementType]
selfElSpace = self.space[elementType]
#evaluation of the shape function of this space at the LagrangeSpaceGeo pos
ir = (geoElSpace.posN, np.ones(geoElSpace.GetNumberOfShapeFunctions(), dtype=MuscatFloat))
saip = selfElSpace.GetSpaceEvaluatedAt(ir[0], ir[1])
# mask des element with a numbering
emask = np.all(selfNumbering >= 0,axis=1)
dataOnElementShapeFunctions = self.data[selfNumbering[emask,:]]
localConnectivity = selection.elements.connectivity[emask,:]
for sfn in range(geoElSpace.GetNumberOfShapeFunctions()):
# evaluation of the field at the geo space points
restemp[localConnectivity[:,sfn]] += saip.valN[sfn].dot(dataOnElementShapeFunctions.T)
cpttmp[localConnectivity[:,sfn]] += 1
localMask[localConnectivity[:,sfn]] = True
notmask = np.logical_not(mask)
res[notmask] += restemp[notmask]
cpt[notmask] += cpttmp[notmask]
np.logical_or(mask, localMask, out=mask)
res[mask] = res[mask]/cpt[mask]
if fillValue != 0.:
# last time we use mask, so we do an inplace operation
np.logical_not(mask,out=mask)
res[mask] = fillValue
return res
[docs]
def SetDataFromPointRepresentation(self,userdata, fillValue=0.):
if fillValue==0.:
self.data = np.zeros(self.numbering.size)
else:
self.data = np.full(self.numbering.size,fillValue)
self.data[self.numbering.doftopointRight] = userdata[self.numbering.doftopointLeft]
[docs]
def GetCellRepresentation(self,fillValue:MuscatFloat=0, method:Union[str,int]='mean') -> np.ndarray:
"""
Function to push the data from the field into a vector homogeneous to
the mesh cell (for visualisation for example). Entities with no dofs are filled
with the fillValue (default 0)
the method controls to transfer function
"""
if fillValue==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():
nbelems = data.GetNumberOfElements()
numbering = self.numbering[name]
if name is None:
cpt += nbelems
continue
if method == 'mean':
data = np.mean(self.data[numbering],axis=1)
elif method == 'max':
data = np.max(self.data[numbering],axis=1)
elif method == 'min':
data = np.min(self.data[numbering],axis=1)
elif method == 'maxdiff' or method == "maxdifffraction":
cols = self.data[numbering].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[numbering].dot(op)),axis=1)
if method == "maxdifffraction":
data /= np.mean(self.data[numbering],axis=1)
else:
col = min(int(method),self.data[numbering].shape[1])
data = self.data[numbering][:,col]
res[cpt:cpt+nbelems] = data
cpt += nbelems
return res
[docs]
def CheckCompatibility(self,B):
if isinstance(B,type(self)):
if id(self.mesh) != id(B.mesh):
raise (Exception("The support of the fields are not the same"))
if id(self.space) != id(B.space):
raise (Exception("The space of the fields are not the same"))
if self.numbering != B.numbering:
raise (Exception("The numbering of the fields are not the same"))
[docs]
def unaryOp(self,op):
res = type(self)(name = None,mesh=self.mesh,space=self.space, numbering = self.numbering )
res.data = op(self.data)
return res
[docs]
def binaryOp(self,other,op):
self.CheckCompatibility(other)
res = type(self)(name = None,mesh=self.mesh,space=self.space, numbering = self.numbering )
if isinstance(other,type(self)):
res.data = op(self.data,other.data)
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 = op(self.data,other)
return res
else:
raise Exception(f"operator {op} not valid for types :{type(self)} and {type(other)} ")
[docs]
def GetTestField(self):
from Muscat.FE.SymWeakForm import GetTestSufixChar
tc = GetTestSufixChar()
if len(self.name) == 0:
raise Exception("FEField must have a name")
elif self.name[-1] == tc:
raise Exception("this FEField is already a test field")
else:
return FEField(name=self.name+tc,mesh=self.mesh,space=self.space,numbering=self.numbering)
def __str__(self):
return f"<FEField object '{self.name}' ({id(self)})>"
def __repr__(self):
return str(self)
[docs]
def CheckIntegrity(GUI:bool=False):
from Muscat.MeshTools.MeshCreationTools import CreateCube
mesh = CreateCube([2.,3.,4.],[-1.0,-1.0,-1.0],[2./10, 2./10,2./10])
sig11 = FEField(name = "temp",mesh=mesh)
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
sig11 = FEField(name = "temp",mesh=mesh, space= LagrangeSpaceGeo)
sig11 = FEField(name = "temp",mesh=mesh)
print(sig11.GetTestField())
sig11.Allocate()
print(sig11)
sig11 = sig11.GetPointRepresentation()
sig22 = sig11+0.707107
sig12 = 2*(-sig22)*5
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(np.linalg.norm([sig22, sig22 ] ).data )
print(A/C)
dummyField = FEField()
dummyField.data = np.arange(3)+1
print("dummyField")
print(dummyField*dummyField-dummyField**2/dummyField)
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(True))# pragma: no cover