Two inclusions mechanical analysis

Here we present a study case of a thick plate with 2 inclusions. One softer and the other stiffer than the base material. Then we compute the energy deformation on one inclusion.

Fig. 1 Example of mechanical simulation and post processing

from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym

#Main class to
problem = UnstructuredFeaSym()

# the mecanical problem
from Muscat.FE.SymPhysics import MecaPhysics
mecaPhysics = MecaPhysics()

# Definition of the degree of the spaces [1 or 2]
mecaPhysics.SetSpaceToLagrange(P=1)

# add weak form terms to the tanget matrix
mecaPhysics.AddBFormulation( "Bulk",mecaPhysics.GetBulkFormulation(1.0,0.3)  )
mecaPhysics.AddBFormulation( "Inclusion1",mecaPhysics.GetBulkFormulation(5.0,0.3)  )
youngModulusInclusionII = 0.5
mecaPhysics.AddBFormulation( "Inclusion2",mecaPhysics.GetBulkFormulation(0.5,0.3)  )

# add weak form term to the rhs
mecaPhysics.AddLFormulation( "Right", mecaPhysics.GetForceFormulation([1,0,0],1)  )

problem.physics.append(mecaPhysics)

# the boundary conditions
from Muscat.FE.KR.KRBlock import KRBlockVector
dirichlet = KRBlockVector()
dirichlet.AddArg("u").On('Left').Fix0().Fix1().Fix2().To()

problem.solver.constraints.AddConstraint(dirichlet)

# Read The mesh
from Muscat.IO.GmshReader import ReadGmsh
import pathlib

mesh = ReadGmsh(pathlib.Path(__file__).parent.resolve() / pathlib.Path("TwoInclusions.msh"))

print(mesh)

problem.SetMesh(mesh)

# we compute the numbering (maping from the mesh to the linear system)
problem.ComputeDofNumbering()

from Muscat.Helpers.Timer import Timer
with Timer("Assembly "):
    k,f = problem.GetLinearProblem()

#compute the constraints to add to the system
problem.ComputeConstraintsEquations()


with Timer("Solve"):
    problem.Solve(k,f)

problem.PushSolutionVectorToUnknownFields()

from Muscat.FE.Fields.FieldTools import GetPointRepresentation
# Recover a point representation of the displacement
problem.mesh.nodeFields["sol"] = GetPointRepresentation(problem.unknownFields)

from Muscat.FE.Fields.FEField import FEField
#Creation of a fake fields to export the rhs member
rhsFields = [ FEField(mesh=mesh,space=None,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")
print("Compute of the strain energy only on the second inclusion (integral in each element) ")

import Muscat.FE.SymWeakForm as SWF
from Muscat.FE.MaterialHelp import HookeIso
from Muscat.Containers.Filters.FilterObjects import ElementFilter

symdep = SWF.GetField("u",3)
K = HookeIso(youngModulusInclusionII,0.3,dim=3)
symCellData = SWF.GetField("cellData",1)
symCellDataT = SWF.GetTestField("cellData",1)

print("Post process")

EnerForm = SWF.ToVoigtEpsilon(SWF.Strain(symdep)).T*K*SWF.ToVoigtEpsilon(SWF.Strain(symdep))*symCellDataT

print("Post process Eval")
ff = ElementFilter(eTag="Inclusion2")
from Muscat.FE.DofNumbering import  ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
from Muscat.FE.Integration import IntegrateGeneral
p0Numbering = ComputeDofNumbering(mesh,LagrangeSpaceP0)

energyDensityField = FEField(name="cellData",mesh=problem.mesh,numbering=p0Numbering,space=LagrangeSpaceP0)

m,energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm,  constants={},
                        fields=problem.unknownFields, unknownFields = [ energyDensityField ], elementFilter=ff)
print("energyDensity",energyDensity)
energyDensityField.data = energyDensity

problem.mesh.elemFields["Energy"] = energyDensity

import numpy as np
print("Strain energy on the second inclusion:", np.sum(energyDensity) )
# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
from Muscat.IO import XdmfWriter as XW
writer = XW.XdmfWriter('TwoInclusions_Output.xdmf')
writer.SetHdf5(False)
writer.Open()
writer.Write(mesh,PointFields=list(mesh.nodeFields.values()), PointFieldsNames=list(mesh.nodeFields.keys()),
                CellFields=list(mesh.elemFields.values()), CellFieldsNames=list(mesh.elemFields.keys()))
writer.Close()

print("Done")