# -*- 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, Dict, Optional, List
import numpy as np
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.Types import ArrayLike, MuscatFloat
from Muscat.FE.KR.KRBase import KRBase
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.FE.Fields.FEField import FEField
from Muscat.LinAlg.ConstraintsHolder import ConstraintsHolder
[docs]
class KRFromDistribution(KRBase):
"""Class to impose values in the primary unknown of the system
for the moment this only work with iso-parametric fields
"""
def __init__(self):
super().__init__()
self.type: str = "RFromVector",
self.fieldsData: Dict[str,FEField] = {}
[docs]
def FixUsingNodeFields(self, mesh: Mesh, nodalField:Dict[str,ArrayLike] ) -> KRFromDistribution :
"""Set the fields to be used to impose the dof values.
Internally we convert the numpy array into fields with the function :
Muscat.FE.Fields.FieldTools.NodeFieldToFEField()
Parameters
----------
mesh : Mesh
the mesh
nodalField : Dict[str,ArrayLike]
a dict with key:value in the form of "u":numpy(3,NbNodes)
Returns
-------
KRFromDistribution
self
"""
from Muscat.FE.Fields.FieldTools import NodeFieldToFEField
dataFields = NodeFieldToFEField(mesh, nodalField)
return self.FixUsingFEFields(dataFields)
[docs]
def FixUsingFEFields(self, inputFEFields:Dict[str,FEField]) -> KRFromDistribution:
""" Set the fields to be used to impose the dof values
Parameters
----------
inputFEFields : List[FEField]
a list of FEField (the order is not important)
Returns
-------
KRFromDistribution
self
"""
self.fieldsData = inputFEFields
return self
[docs]
def GenerateEquations(self, mesh: Mesh, fields:List[FEField], CH:Optional[ConstraintsHolder]=None, solutionFields:Optional[List[FEField]]=None, reactionFields:Optional[List[FEField]]=None) -> ConstraintsHolder :
"""Generate the equation to impose the constraints defined by this class.
More information about this
Parameters
----------
mesh : Mesh
The mesh
fields : List[FEField]
The list of fields to be used to compute the numbering (in the order on the linear system to be solved)
CH : Optional[ConstraintsHolder], optional
the ConstraintsHolder to store the equations, if None is provided a new object is created, by default None
solutionFields : Optional[List[FEField]], optional
In the case of an iterative solver (i.e. in a newton loop) the current solution, by default None
reactionFields : Optional[List[FEField]], optional
The user can provide the reactions fields ( the integral of the internal forces), by default None
Returns
-------
ConstraintsHolder
the ConstraintsHolder with the new equations
Raises
------
RuntimeError
if no data is provided. the user must call FixUsingFEFields or FixUsingNodeField to provide the fields
LookupError
If the user impose a restriction on a field (using .AddArg()) but not available on the problem
"""
if len(self.fieldsData) == 0:
raise Exception("Please set data to use (FixUsingFEFields,FixUsingNodeField) ")
offsets, fieldOffsets, totalNumberOfDofs = self._ComputeOffsets(fields)
CH = self._GetConstraintHolder(CH)
if CH.nbdof == 0:
CH.nbdof = totalNumberOfDofs
unknownFieldDic = {f.name:f for f in fields }
ef = self._GetOnElementFilter()
dataFieldDic = {f.name:f for f in self.fieldsData.values() }
class FieldData():
def __init__(self, name:str, unknownField, dataField, offset, solutionFields=None, reactionFields=None ):
self.name = name
self.unknownField = unknownField[name]
self.dataField = dataField[name]
self.offset = offset[name]
if solutionFields is None:
self.solutionField = None
else:
solutionFields = {f.name:f for f in solutionFields }
self.solutionField = solutionFields[name[1:]]
if reactionFields is None:
reactionFields= None
else:
reactionFields = {f.name:name for f in reactionFields }
self.reactionField = reactionFields["R"+name[1:]]
fieldsToTreat = []
# sanity check and prep of dofs to impose
for arg in self.args:
if "_" in arg:
if arg in unknownFieldDic.keys() and arg in dataFieldDic.keys() :
fieldsToTreat.append( FieldData(arg, unknownFieldDic, dataFieldDic, fieldOffsets, solutionFields,reactionFields ) )
else:
raise LookupError(f"Field ({arg}) not found.")
else:
if arg in unknownFieldDic.keys() and arg in dataFieldDic.keys() :
fieldsToTreat.append( FieldData(arg, unknownFieldDic, dataFieldDic, fieldOffsets, solutionFields,reactionFields ) )
else:
for i in range(3):
composedName = f"{arg}_{i}"
if composedName in fieldOffsets:
fieldsToTreat.append( FieldData(composedName, unknownFieldDic, dataFieldDic, fieldOffsets, solutionFields,reactionFields ) )
else:
break
for fieldToTreat in fieldsToTreat:
totalNumberOfShapeFuntions = fieldToTreat.unknownField.numbering.size
treated = np.zeros(totalNumberOfShapeFuntions, dtype=bool)
offset = fieldToTreat.offset
for selection in ef(mesh):
numberOfShapeFunctions = fieldToTreat.unknownField.space[selection.elementType].GetNumberOfShapeFunctions()
for elemId in selection.indices:
for i in range(numberOfShapeFunctions):
dofIid = fieldToTreat.unknownField.numbering[selection.elementType][elemId,i]
if treated[dofIid]:
continue
treated[dofIid] = True
dataValue = fieldToTreat.dataField.data[fieldToTreat.dataField.numbering[selection.elementType][elemId,i]]
if fieldToTreat.solutionField is not None:
solValue = fieldToTreat.solutionField.data[fieldToTreat.solutionField.numbering[selection.elementType][elemId,i]]
else:
solValue = 0.
CH.AddFactor(dofIid+offset,1)
CH.AddConstant(dataValue-solValue)
CH.NextEquation()
return CH
def __str__(self) -> str:
res = self.type + " "
if len(self.args) > 1:
res += "("
res += " and ".join(str(x) for x in self.args)
if len(self.args) > 1:
res += ")"
return res
[docs]
def CheckIntegrity(GUI:bool=False):
from Muscat.MeshTools.MeshCreationTools import CreateSquare
mesh = CreateSquare()
print(mesh)
imposedU = mesh.nodes + 3
obj = KRFromDistribution()
obj.AddArg("u").On("X0")
obj.FixUsingNodeFields(mesh, {"u":imposedU})
from Muscat.FE.Fields.FieldTools import NodeFieldToFEField
uFields = list(NodeFieldToFEField(mesh, {"u":imposedU}).values())
ch = obj.GenerateEquations(mesh,uFields)
mat, dofs = ch.ToSparse()
print(mesh.nodes)
print(imposedU.T.ravel()[dofs[:-1]])
print(mat.toarray())
print(dofs)
from Muscat.MeshTools.MeshCreationTools import CreateCube
nx = 5; ny = 5; nz = 5;
mesh = CreateCube(dimensions=[nx,ny,nz], origin=[0,0,0.], spacing=[1/(nx-1),1/(ny-1), 1/(nz-1)], ofTetras=False )
from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym
problem = UnstructuredFeaSym()
from Muscat.FE.SymPhysics import MechPhysics
mechPhysics = MechPhysics()
mechPhysics.SetSpaceToLagrange(1)
mechPhysics.AddBFormulation( "3D",mechPhysics.GetBulkFormulation(1.0,0.3) )
problem.physics.append(mechPhysics)
from Muscat.FE.KR.KRBlock import KRBlockVector
dirichlet = KRBlockVector()
dirichlet.AddArg("u").On('Z0').Fix0().Fix1().Fix2()
problem.solver.constraints.AddConstraint(dirichlet)
obj = KRFromDistribution()
# this to impose only displacement in the z direction
#obj.AddArg("u_2").On("Z1")
obj.AddArg("u").On("Z1")
imposedU = np.full_like(mesh.nodes,0.1)
imposedU[:,2] += np.cos(3*mesh.nodes[:,0])*0.1
obj.FixUsingNodeFields(mesh, {"u":imposedU})
problem.solver.constraints.AddConstraint(obj)
#problem.solver.constraints.SetConstraintsMethod("Penalisation")
problem.solver.SetAlgo("Direct")
problem.SetMesh(mesh)
problem.ComputeDofNumbering()
k,f = problem.GetLinearProblem()
problem.ComputeConstraintsEquations()
problem.Solve(k,f)
problem.PushSolutionVectorToUnknownFields()
from Muscat.FE.Fields.FieldTools import GetPointRepresentation
problem.mesh.nodeFields["sol"] = GetPointRepresentation(problem.unknownFields)
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh,filename="Test_KRFromVector.xmf",run=GUI)
from Muscat.MeshContainers.Filters.FilterObjects import NodeFilter
ids = NodeFilter(eTag=["Z1"]).GetNodesIndices(mesh)
error = np.abs(problem.mesh.nodeFields["sol"][ids,:] - imposedU[ids,:])
if np.linalg.norm(error/imposedU[ids,:]) > 1e-8:
RuntimeError("Error in the CheckIntegrity ")
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(GUI=True))