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.

Includes

from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym
from Muscat.FE.KR.KRBlock import KRBlockVector
from Muscat.IO.GmshReader import ReadGmsh
from Muscat.FE.Fields.FieldTools import GetPointRepresentation
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.SymPhysics import MecaPhysics
import Muscat.FE.SymWeakForm as SWF
from Muscat.FE.MaterialHelp import HookeIso
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.FE.DofNumbering import  ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.IO.XdmfWriter import XdmfWriter
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.FE.Fields.FieldTools import VectorToFEFieldsData
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory

import numpy as np
import pyvista
pyvista.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots

import pathlib

Unstructured FEA with mechanical

problem = UnstructuredFeaSym()

# the mecanical problem
mecaPhysics = MecaPhysics()

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

# add weak form terms to the tanget matrix
youngModulusInclusionBase = 1.0
mecaPhysics.AddBFormulation( "Bulk",mecaPhysics.GetBulkFormulation(youngModulusInclusionBase,0.3)  )
mecaPhysics.AddBFormulation( "Inclusion1",mecaPhysics.GetBulkFormulation(5*youngModulusInclusionBase,0.3)  )
youngModulusInclusionII = youngModulusInclusionBase/2
mecaPhysics.AddBFormulation( "Inclusion2",mecaPhysics.GetBulkFormulation(youngModulusInclusionII,0.3)  )

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

Boundary conditions

# the boundary conditions
dirichlet = KRBlockVector()
dirichlet.AddArg("u").On('Left').Fix0().Fix1().Fix2().To()

problem.solver.constraints.AddConstraint(dirichlet)

Numerical part (mesh)

# Read The mesh

mesh = ReadGmsh(pathlib.Path("").resolve() / pathlib.Path("TwoInclusions.msh"))
problem.SetMesh(mesh)

print(mesh)

Mesh
  Number Of Nodes    : 6750
    Tags :
  Number Of Elements : 29129
    ElementsContainer,   Type : (ElementType.Tetrahedron_4,28889),   Tags : (Bulk:17170) (GeoTag1:17170) (Inclusion1:4348) (GeoTag2:4348) (Inclusion2:7371) (GeoTag3:7371)
    ElementsContainer,   Type : (ElementType.Triangle_3,240),   Tags : (Right:120) (GeoTag33:120) (Left:120) (GeoTag41:120)

Resolution

# 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()

# Recover a point representation of the displacement
mesh.nodeFields["sol"] = GetPointRepresentation(problem.unknownFields)
print(mesh)

Mesh
  Number Of Nodes    : 6750
    Tags :
  Number Of Elements : 29129
    ElementsContainer,   Type : (ElementType.Tetrahedron_4,28889),   Tags : (Bulk:17170) (GeoTag1:17170) (Inclusion1:4348) (GeoTag2:4348) (Inclusion2:7371) (GeoTag3:7371)
    ElementsContainer,   Type : (ElementType.Triangle_3,240),   Tags : (Right:120) (GeoTag33:120) (Left:120) (GeoTag41:120)
  nodeFields         : ['sol']

Post Processing

#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) ]
VectorToFEFieldsData(f,rhsFields)
problem.mesh.nodeFields["RHS"] = GetPointRepresentation(rhsFields)

Compute of the strain energy only on the second inclusion (integral in each element)

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")

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
print("Strain energy on the second inclusion:", np.sum(energyDensity) )

Post process
Post process Eval
energyDensity [0. 0. 0. ... 0. 0. 0.]
Strain energy on the second inclusion: 0.0777788948445724

Graphical Output

# generate a element field with un index per material
mesh.elemFields["inclussions"]  = np.zeros(mesh.GetNumberOfElements())
for i, name in enumerate(("Inclusion1","Inclusion2")):
    for selection in ElementFilter(eTag=name)(mesh):
        mesh.elemFields["inclussions"][selection.GetMeshSlice()]= 1 + i

import pyvista as pv
from Muscat.Bridges.PyVistaBridge import MeshToPyVista

plotter = pv.Plotter(shape=(2, 2))

camera = [(1.4464393507963245, 1.760895123462273, 2.9110300316336843),
 (-0.22723828652575934, -0.08454810909441325, 0.21480857439407788),
 (-0.25495941474725625, 0.8643967380870753, -0.43337509851877826)]

plotter.subplot(0, 0)
plotter.add_text("Two Inclusions", font_size=10)
plotter.add_mesh(MeshToPyVista(mesh), scalars="inclussions", show_edges=True, show_scalar_bar=False, edge_color="grey", cmap="rainbow_r")
plotter.camera_position = camera

plotter.subplot(0, 1)
plotter.add_text("Deformation", font_size=10)
deformedMesh = mesh.View()
deformedMesh.nodes = mesh.nodes + problem.mesh.nodeFields["sol"]*.3
deformedMesh.nodeFields["sol"] = problem.mesh.nodeFields["sol"]
plotter.add_mesh(MeshToPyVista(deformedMesh), scalars="sol", show_edges=False, show_scalar_bar=False, cmap="coolwarm_r")
plotter.camera_position = camera
rightMesh = ExtractElementsByElementFilter(deformedMesh,ElementFilter(eTag="Right"))
rightMesh.nodeFields = deformedMesh.nodeFields
glyphs = MeshToPyVista(rightMesh).glyph(orient="RHS", scale="RHS", factor=10, geom=pv.Arrow())
plotter.add_mesh(glyphs, show_scalar_bar=False)

plotter.subplot(1, 0)
plotter.add_text("Inclusion 2 energy dencity", font_size=10)
plotter.add_mesh(MeshToPyVista(ExtractElementsByElementFilter(mesh,ElementFilter(eTag="Bulk")) ).extract_feature_edges(45) )
in1Mesh = ExtractElementsByElementFilter(mesh,ElementFilter(eTag="Inclusion2"))
in1Mesh.elemFields["Energy"] = mesh.elemFields["Energy"][in1Mesh.GetElementsOriginalIDs()]
plotter.add_mesh(MeshToPyVista(in1Mesh), scalars="Energy", style="surface" , cmap="coolwarm" )
plotter.camera_position = camera

plotter.subplot(1, 1)
plotter.add_text("Inclusion deformation", font_size=10)
deformedMesh = mesh.View()
deformedMesh.nodes = mesh.nodes + problem.mesh.nodeFields["sol"]*.3
buldkDeformedMesh = ExtractElementsByElementFilter(deformedMesh,ElementFilter(eTag="Bulk"))
buldkDeformedMesh.nodeFields["sol"] = problem.mesh.nodeFields["sol"]
plotter.add_mesh(MeshToPyVista(buldkDeformedMesh).extract_feature_edges(45),scalars="sol", show_edges=False, show_scalar_bar=False)
in1Mesh.nodes = mesh.nodes + problem.mesh.nodeFields["sol"]*.3
plotter.add_mesh(MeshToPyVista(in1Mesh), scalars="Energy", style="surface" , cmap="coolwarm" )
plotter.add_mesh(glyphs)
plotter.camera_position = camera

plotter.show()

Export data to file

# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
writer = XdmfWriter(TemporaryDirectory.GetTempPath() + '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()