# -*- 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 Dict, List
import numpy as np
from Muscat.Types import MuscatFloat
from Muscat.FE.FeaBase import FeaBase
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
import Muscat.FE.SymWeakForm as SWF
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.Fields.FEField import FEField
from Muscat.Helpers.Logger import Debug
from Muscat.MeshContainers.Mesh import Mesh
[docs]
class UnstructuredFeaSym(FeaBase):
def __init__(self,spaceDim:int=3 ) -> None:
super().__init__(spaceDim=spaceDim)
self.constants:Dict[str,MuscatFloat] = {}
self.fields: Dict = {}
self.physics:List = []
self.spaces:List = []
self.numberings:List = []
[docs]
def SetMesh(self, support: Mesh) -> None:
super().SetMesh(support)
for phy in self.physics:
phy.Reset()
[docs]
def ComputeDofNumbering(self,tagsToKeep=None,fromConnectivity=False):
self.spaces = []
self.numberings = []
for phy in self.physics:
self.spaces.extend(phy.spaces)
phy.ComputeDofNumbering(self.mesh,tagsToKeep,fromConnectivity=fromConnectivity)
self.numberings.extend(phy.numberings)
self.totalNumberOfDof = 0
for num in self.numberings:
self.totalNumberOfDof += num.size
self.unknownFields = []
for phy in self.physics:
for dim in range(phy.GetNumberOfUnknownFields()):
field = FEField()
field.numbering = phy.numberings[dim]
field.name = phy.GetPrimalNames()[dim]
field.mesh = self.mesh
field.space = phy.spaces[dim]
field.Allocate()
self.unknownFields.append(field)
self.solver.constraints.SetNumberOfDofs(self.totalNumberOfDof)
[docs]
def GetLinearProblem(self,computeK=True, computeF=True, lform = None ,bform = None,unknownFields=None):
rhsRes = None
lhsRes = None
if lform is not None:
if unknownFields is None:
unknownFields = self.unknownFields
for ff,form in lform:
if form is None:
continue
Debug("integration of lform "+ str(ff) )
#Debug("integration of lform "+ str(form) )
_,f = IntegrateGeneral(mesh=self.mesh,
wform=form,
constants=self.constants,
fields=list(self.fields.values()),
unknownFields= unknownFields,
elementFilter=ff)
if rhsRes is None:
rhsRes = f
else:
rhsRes += f
return (lhsRes,rhsRes)
if bform is not None:
if unknownFields is None:
unknownFields = self.unknownFields
for ff,form in bform:
if form is None:
continue
Debug("integration of bform " + str(ff) )
#Debug("integration of bform " + str(form) )
k,f = IntegrateGeneral(mesh=self.mesh,
wform=form,
constants=self.constants,
fields=list(self.fields.values()),
unknownFields= unknownFields,
elementFilter=ff)
if rhsRes is None:
rhsRes = f
else:
rhsRes += f
if lhsRes is None:
lhsRes = k
else:
lhsRes += k
return (lhsRes,rhsRes)
for phy in self.physics:
for nf,data in phy.extraRHSTerms:
nids = nf.GetNodesIndices(self.mesh)
size = 0
offset = []
for n in phy.numberings:
offset.append(size)
size += n.size
for i,val in enumerate(data):
inumbering = phy.numberings[i]
if val == 0 :
continue
if rhsRes is None:
rhsRes = np.zeros(size,float)
for nid in nids:
dof = inumbering.GetDofOfPoint(nid)+offset[i]
rhsRes[dof] += val
if (computeF):
Debug("In Integration F")
for phy in self.physics:
linearWeakFormulations = phy.linearWeakFormulations
for ff,form in linearWeakFormulations:
if form is None:
continue
Debug("In Integration of f "+ str(ff) )
_,f = IntegrateGeneral(mesh=self.mesh,wform=form, constants=self.constants, fields=list(self.fields.values()),unknownFields= self.unknownFields,
integrationRuleName=phy.integrationRule,elementFilter=ff)
if rhsRes is None:
rhsRes = f
else:
rhsRes += f
if (computeK):
Debug("In Integration K")
for phy in self.physics:
bilinearWeakFormulations = phy.bilinearWeakFormulations
for ff,form in bilinearWeakFormulations:
if form is None:
continue
k,f = IntegrateGeneral(mesh=self.mesh,wform=form, constants=self.constants, fields=list(self.fields.values()), unknownFields= self.unknownFields,
integrationRuleName=phy.integrationRule,elementFilter=ff)
if not (f is None):
if rhsRes is None:
rhsRes = f
else:
rhsRes += f
if lhsRes is None:
lhsRes = k
else:
lhsRes += k
return (lhsRes,rhsRes)
[docs]
def CheckIntegrity(GUI:bool=False):
for P in [1,2]:
for tetra in [False,True]:
print("in CheckIntegrityFlexion P="+str(P)+" tetra="+str(tetra))
res = CheckIntegrityFlexion( P = P,tetra = tetra,GUI=GUI);
if res.lower()!="ok": return "ERROR: "+res + " " + str(P) + " " + str(tetra)
return "ok"
[docs]
def CheckIntegrityFlexion(P,tetra,GUI=False):
# the main class
problem = UnstructuredFeaSym()
# the mechanical problem
from Muscat.FE.SymPhysics import MechPhysics
mechPhysics = MechPhysics()
mechPhysics.SetSpaceToLagrange(P=P)
mechPhysics.AddBFormulation( "3D",mechPhysics.GetBulkFormulation(1.0,0.3) )
mechPhysics.AddLFormulation( "Z1", mechPhysics.GetForceFormulation([1,0,0], 0.002) )
mechPhysics.AddLFormulation( "Z0", None )
if P == 1:
mechPhysics.integrationRule = "LagrangeP1Quadrature"
else:
mechPhysics.integrationRule = "LagrangeP2Quadrature"
problem.physics.append(mechPhysics)
# the boundary conditions
from Muscat.FE.KR.KRBlock import KRBlockVector
dirichlet = KRBlockVector()
dirichlet.AddArg("u").On('Z0').Fix0().Fix1().Fix2().To(offset=[1,2,3],first=[1,0,1] )
dirichlet.constraintDirections= "Global"
problem.solver.constraints.AddConstraint(dirichlet)
# the mesh
from Muscat.MeshTools.MeshCreationTools import CreateCube
nx = 11; ny = 12; nz = 13
nx = 3
ny = 3
nz = 4
mesh = CreateCube(dimensions=[nx,ny,nz],origin=[0,0,0.], spacing=[1./(nx-1),1./(ny-1), 10./(nz-1)], ofTetras=tetra )
mesh.ConvertDataForNativeTreatment()
problem.SetMesh(mesh)
print(mesh)
# we compute the numbering
problem.ComputeDofNumbering()
from Muscat.Helpers.Timer import Timer
with Timer("Assembly "):
k,f = problem.GetLinearProblem()
#problem.solver = LinearProblem()
#problem.solver.SetAlgo("EigenCG")
#problem.solver.SetAlgo("EigenLU")
problem.solver.SetAlgo("Direct")
problem.ComputeConstraintsEquations()
print("k.shape " + str(k.shape) )
print("f.shape "+ str(f.shape))
with Timer("Solve"):
problem.Solve(k,f)
problem.PushSolutionVectorToUnknownFields()
from Muscat.FE.Fields.FieldTools import GetPointRepresentation
problem.mesh.nodeFields["sol"] = GetPointRepresentation(problem.unknownFields)
#Creation of a fake fields to export the rhs member
rhsFields = [ FEField(mesh=mesh,space=problem.spaces[0],numbering=problem.numberings[i]) for i in range(3) ]
from Muscat.FE.Fields.FieldTools import VectorToFEFieldsData
VectorToFEFieldsData(f,rhsFields)
problem.mesh.nodeFields["RHS"] = GetPointRepresentation(rhsFields)
print("done solve")
symdep = SWF.GetField("u",3,sdim=3)
from Muscat.FE.MaterialHelp import HookeIso
K = HookeIso(1,0.3,dim=3)
symCellData = SWF.GetField("cellData",1,sdim=3)
symCellDataT = SWF.GetTestField("cellData",1,sdim=3)
print("Post process")
EnerForm = SWF.ToVoigtEpsilon(SWF.Strain(symdep)).T*K*SWF.ToVoigtEpsilon(SWF.Strain(symdep))*symCellDataT + symCellData.T*symCellDataT
print("Post process Eval")
ff = ElementFilter(eTag="3D")
energyDensityField = FEField(name="cellData",mesh=problem.mesh,numbering=problem.numberings[0],space=problem.spaces[0])
m,energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm, constants={},
fields=problem.unknownFields, unknownFields = [ energyDensityField ], integrationRuleName="NodalEvaluationP"+str(P),
onlyEvaluation=True,
elementFilter=ff)
print("energyDensity",energyDensity)
energyDensity /= m.diagonal()
energyDensityField.data = energyDensity
problem.mesh.nodeFields["PEnergy"] = energyDensityField.GetPointRepresentation()
problem.mesh.elemFields["PEnergy_fromP"+str(P) ] = energyDensityField.GetCellRepresentation()
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
from Muscat.FE.DofNumbering import ComputeDofNumbering
P0Numbering = ComputeDofNumbering(mesh,LagrangeSpaceP0,elementFilter=ElementFilter(dimensionality=mesh.GetPointsDimensionality()))
P0energyDensityField = FEField(name="cellData",mesh=problem.mesh,numbering=P0Numbering,space=LagrangeSpaceP0)
m,P0energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm, constants={},
fields=problem.unknownFields, unknownFields = [ P0energyDensityField ],
integrationRuleName="ElementCenter",
onlyEvaluation=True,
elementFilter=ff)
P0energyDensity /= m.diagonal()
P0energyDensityField.data = P0energyDensity
problem.mesh.elemFields["CEnergy"] = P0energyDensityField.GetCellRepresentation()
if GUI :
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh,filename="UnstructuredFeaSym_Sols_P"+str(P)+("Tetra"if tetra else "Hexa")+".xmf",run=True)
print(Timer())
return("ok")
if __name__ == '__main__':
import time
starttime = time.time()
print(CheckIntegrity(True))#pragma: no cover
stoptime = time.time()
print("Total Time {0}s".format(stoptime-starttime))
print("Done")