#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 27 07:43:08 2024
@author: nchevaug
"""
from __future__ import annotations
import os
# # Set the maximum number of threads for OpenMP and related backends
os.environ["OMP_NUM_THREADS"] = "8" # Change '4' to your desired core count
os.environ["MKL_NUM_THREADS"] = "8"
os.environ["OPENBLAS_NUM_THREADS"] = "8"
from warnings import warn
import time
import unittest
import numpy as np
import scipy as sp
try:
import sksparse.cholmod # pylint: disable=E0401
have_cholmod = True
except ModuleNotFoundError:
have_cholmod = False
warn('cholmod not installed, using default scipy spsolve')
import gmsh
import pyxfem.core.function_spaces as fs
import pyxfem.core.quadrature as quad
import pyxfem.core.sparse_tools as st
from pyxfem.core.mesh_simplex_one_level import SMesh, gmshmodel2smesh
import pyxfem.postpro.gmsh_post_pro_interface as gp
from pyxfem.core.exec_tools import Profile
import pyxfem.core.gmsh_mesh_generator as mg
from pyxfem.solidmechanics.boundary_conditions import compute_dirichlet, compute_neumann, \
get_dofs_and_interpolation_points, compute_dirichlet_fem_p1_vect2d
from pyxfem.solidmechanics.elasticity import planeStrainC2222CSR, YoungPoisson2LambdaMu
[docs]
def getSpaceEntity(edgeGroupeNames, space):
''' get the space's edge id corresponding to the edgeGroupe Name'''
edgeGroupeNames2EId = {}
for name in edgeGroupeNames:
edgeGroupeNames2EId[name] = space.m.getEdgeGroupe(name)[0]
return edgeGroupeNames2EId
[docs]
def computeVolumeLoad(space, source):
""" compute the vector F such as F.U = int_Omega f.u_h dOmega_h """
assert isinstance(source, np.ndarray)
quadSource = quad.T3_gauss(0)
N = space.operator_dof2val_tri(quadSource.uv)
J = space.get_tris_mapping().J
wSource = quadSource.w[0]*source[np.newaxis, :]*J[:, np.newaxis]
Fint = N.T.dot(wSource.flatten())
return Fint
[docs]
class ElasticitySolver():
""" A class to set up and compute a 2d linear elasticity problem with FEM method """
def __init__(self, mesh:SMesh, C=planeStrainC2222CSR(lam =12., mu =7.),
node_id_to_imposed_value:None|dict=None,
physName2Dirichlet:None|dict=None,
physName2Neumann:None|dict=None,
source=np.array([0., 0.]),
space_constructor:type[fs.FEMSpaceP1]|type[fs.FEMSpaceP1NC]=fs.FEMSpaceP1,
stab=0.):
self.quadBilinForm = quad.T3_gauss(0)
self.quadAtCorners = quad.T3_nodes()
self.C = C
self.stab = stab
self.mesh = mesh
self.setModel(space_constructor)
if node_id_to_imposed_value is not None :
if physName2Dirichlet is not None :
warn("can't impose dirichlet BC both from physName and list of nodes, ignoring physName")
if space_constructor != fs.FEMSpaceP1 :
warn("imposing dirichlet at node value make sens only for p1 fem space. Expect errors ...")
nv = mesh.nv()
idvx = node_id_to_imposed_value.get("idvx", np.array([], dtype = np.int64))
idvy = node_id_to_imposed_value.get("idvy", np.array([], dtype = np.int64))
vx = node_id_to_imposed_value.get("vx", np.array([], dtype = np.int64))
vy = node_id_to_imposed_value.get("vy", np.array([], dtype = np.int64))
self.freeDofs, self.fixedDofs, self.fixedValues = \
compute_dirichlet_fem_p1_vect2d(nvertex=nv, nodex= idvx, valuex= vx, nodey = idvy, valuey=vy)
else:
self.physName2Dirichlet:dict = {} if physName2Dirichlet is None else physName2Dirichlet
gname2dof, gname2dofxy = get_dofs_and_interpolation_points(mesh,
self.physName2Dirichlet.keys(), self.space)
self.freeDofs, self.fixedDofs, self.fixedValues = \
compute_dirichlet(self.space, self.physName2Dirichlet, gname2dof, gname2dofxy)
self.physName2Neumann:dict = {} if physName2Neumann is None else physName2Neumann
self.Fneumann = compute_neumann(self.space, self.physName2Neumann)
self.Fint = computeVolumeLoad(self.space, source)
self.Bfree = self.B[:, self.freeDofs]
self.Bfixed = self.B[:, self.fixedDofs]
nbfixed = len(self.fixedDofs)
if nbfixed < 3:
warn(f' Only {nbfixed:d} fixed dofs: not enough to prevent rigid-body motion in {compute_dirichlet.__name__}')
if isinstance(self.space, fs.FEMSpaceP1NC):
self.setStabilisation(self.stab, quad.Edge_gauss(2))
[docs]
def setD(self):
''' set the sparse stifness operator (K = BtDB)'''
J = self.space.get_tris_mapping().J
nf = len(J)
C = self.C.todense().flatten()
wD = np.einsum("i, ij -> ij", self.quadBilinForm.w[0]*J,
C[np.newaxis, :]).reshape((-1, 4, 4))
self.wD = st.regular_block_diag_to_csr(wD)
D = np.broadcast_to(C, (nf, 16)).reshape((-1, 4, 4))
self.D = st.regular_block_diag_to_csr(D)
[docs]
def setModel(self, spaceConstructor:type[fs.FEMSpaceP1]|type[fs.FEMSpaceP1NC]):
self.space = spaceConstructor("DisplacementSpace", self.mesh, vecdim=2)
self.B = self.space.operator_dof2grad_tri(self.quadBilinForm.uv)
self.NatVertices = self.space.operator_dof2val_tri(self.quadAtCorners.uv)
self.setD()
[docs]
def setStabilisation(self, s, quadEdge):
if not isinstance(self.space, fs.FEMSpaceP1NC):
raise ValueError("Stabilisation term only implemented for P1NC")
space = self.space
lE = space.m.elength()
sTri2edge = space.operator_trival2edge().dot(0.5*space.get_tris_mapping().J).reshape((-1, 2))
hE = np.sum(sTri2edge, axis=1)/2/lE
N = space.operator_edgetrace(quadEdge)
self.jumpOp = N[1::2] - N[:-1:2]
weights = s*((lE/hE)[:, np.newaxis] * quadEdge.w).reshape(-1)
DjumpBlock = (weights[:, np.newaxis]*np.eye(2).flatten()[np.newaxis, :]).reshape(-1, 2, 2)
self.Djump = st.regular_block_diag_to_csr(DjumpBlock)
internalEdges = self.space.m.intern_edges()
vecDim = self.space.vecdim
nptPerEdge = quadEdge.s.size
rowsPerEdge = vecDim*nptPerEdge
internalEdgesDofs = (vecDim*internalEdges[:, np.newaxis]
+np.arange(vecDim)[np.newaxis, : ]).flatten()
stabDofs = np.concatenate([internalEdgesDofs, self.fixedDofs])
stabRows = np.sort(((rowsPerEdge*(stabDofs//vecDim)
+stabDofs%vecDim)[:, np.newaxis]
+nptPerEdge*np.arange(vecDim)[np.newaxis, :]).flatten())
Nedges = self.jumpOp[stabRows, :]
self.NedgesFree = Nedges[:, self.freeDofs]
self.NedgesFixed = Nedges[:, self.fixedDofs]
self.Djump = self.Djump[stabRows, :][:, stabRows]
self.SFree = self.NedgesFree.T.dot(self.Djump.dot(self.NedgesFree))
[docs]
def assembleK(self):
if not hasattr(self, "Kfree"):
Bfree = self.Bfree
self.Kfree = Bfree.T.dot(self.wD).dot(Bfree)
if isinstance(self.space, fs.FEMSpaceP1NC):
self.Kfree = self.Kfree + self.SFree
return self.Kfree
[docs]
def assembleFU(self):
Fneumann = self.Fneumann
U = np.zeros(self.space.size())
U[self.fixedDofs] = self.fixedValues
Bfree = self.Bfree
Bfixed = self.Bfixed
gradUdiri = Bfixed.dot(self.fixedValues)
Fdiri = -Bfree.T.dot(self.wD.dot(gradUdiri))
F = Fdiri + self.Fneumann[self.freeDofs] + self.Fint[self.freeDofs]
if isinstance(self.space, fs.FEMSpaceP1NC):
SDiri = -self.NedgesFree.T.dot(self.Djump.dot(self.NedgesFixed.dot(self.fixedValues)))
F = F+SDiri
return F,U
[docs]
def assemble(self):
Kfree = self.assembleK()
F,U = self.assembleFU()
return Kfree, F, U
# def assembleMass(self):
[docs]
def solve(self, solver="default"):
""" solve the problem using the solver to choose between "default" : spsolve or "cholmod" : "sksparse.cholmod.cholesky" """
Kfree, F, U = self.assemble()
def lin_solve():
if solver == "default" or not have_cholmod:
U[self.freeDofs] = sp.sparse.linalg.spsolve(Kfree, F, use_umfpack=False)
elif solver == "umfpack":
U[self.freeDofs] = sp.sparse.linalg.spsolve(Kfree, F)
elif solver == "cholmod":
print("cholmod solver")
if not have_cholmod :
raise ValueError
chol = sksparse.cholmod.cho_factor(Kfree)
U[self.freeDofs] = chol.solve(F)
lin_solve()
# if compute_rcond:
# diag = chol.D()
# rcond = (np.min(diag)/np.max(diag))**2
return U
[docs]
def postProDisplacement(self, U):
return self.mesh.t2v_coord(), self.NatVertices.dot(U).reshape((-1, 3, 2))
[docs]
def postProStress(self, U):
return self.mesh.t2v_coord(), self.getStress(U)
[docs]
def getStress(self, U):
''' return per element stress /per gp stress '''
return self.D.dot(self.B.dot(U)).reshape((-1, 2, 2))
[docs]
def postProStrain(self, U):
return self.mesh.t2v_coord(), self.B.dot(U).reshape((-1, 2, 2))
[docs]
def L2StressError(self, U, exactStressFun, quadrature=quad.T3_gauss(0)):
uv = quadrature.uv
weight = quadrature.w
npt = len(uv)
B = self.space.operator_dof2grad_tri(uv)
gradNum = B.dot(U).reshape((-1, npt, 2, 2))
C = self.C.todense().reshape((1,1,2,2,2,2))
stressNum = np.einsum('egijkl, egkl-> egij', C, gradNum)
NP1 = fs.FEMSpaceP1.N(uv)
trisQuadPoints = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :, :], self.mesh.t2v_coord())
stressExact = exactStressFun(trisQuadPoints)
stressError2 = np.sum((stressNum - stressExact)**2, axis = (2,3))
J = self.space.get_tris_mapping().J
wStressError2 = J[:, np.newaxis]*weight[np.newaxis, :]*stressError2
return np.sqrt(np.sum(wStressError2))
[docs]
class testElasticity(unittest.TestCase):
[docs]
def testPatch(self):
""" - patch test on a squared load by Dirchlet and then Neumann """
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 0)
L = 0.7
E = 12.2
nu = 0.3
ux = 0.98
lamb, mu = YoungPoisson2LambdaMu(E, nu)
C = planeStrainC2222CSR(lam = lamb, mu =mu)
phys2Dirichlet = {'left': {'x':0.}, 'bottom': {'y': 0.}, 'right': {'x': ux}}
modelName, physical_dict = mg.square(length=L, relh=0.3)
mesh = gmshmodel2smesh(modelName, phys2Dirichlet.keys())
pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet,
space_constructor=fs.FEMSpaceP1)
U = pb.solve()
___, stress = pb.postProStress(U)
e11 = ux/L
s11 = 4*mu*(lamb+mu)/(lamb+2*mu)*e11
stressExact = np.array([[s11, 0.], [0.,0.]])
error = np.linalg.norm(stress - stressExact[np.newaxis, np.newaxis] )
self.assertTrue(error < 1.e-9)
phys2Dirichlet = {'left': {'x':0.}, 'bottom': {'y': 0.}}
phys2Neumann = {'right': {'x':s11}}
pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet, physName2Neumann=phys2Neumann,
space_constructor=fs.FEMSpaceP1)
e22 = -lamb*e11/(lamb+2*mu)
strainExact = np.array([[e11, 0.], [0., e22]])
U = pb.solve()
___, strain = pb.postProStrain(U)
error = np.linalg.norm(strain - strainExact[np.newaxis, np.newaxis] )
self.assertTrue(error < 1.e-9)
gmsh.finalize()
[docs]
def main(*, element_size=0.05, profile=True, solver = "default", popUpGmsh = True):
""" main Program and example : elasticity on L shape domain """
E = 1.
nu = 0.4999
lamb, mu = YoungPoisson2LambdaMu(E, nu)
meshgenerator = mg.l_shape
#meshgenerator = mg.square
phys2Dirichlet = {'dOmega_left': {'x':0., 'y':0.}}
phys2Neumann = {'dOmega_right':{'x': 0.01}}
#phys2Dirichlet = {'left': {'x':0., 'y':0.}, 'right':{'x': 0.01, 'y':0.}}
#phys2Neumann = {}
spaceConstructors = [fs.FEMSpaceP1, fs.FEMSpaceP1NC]
Profile.doprofile = profile
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
gmsh.option.setNumber('Mesh.SurfaceEdges', 1 if element_size > 0.1 else 0)
gmsh.option.setNumber('General.Terminal', False)
modelName, physical_dict = meshgenerator(element_size)
edgeGroupeNames = list(phys2Dirichlet.keys())+ list(phys2Neumann.keys())
mesh = gmshmodel2smesh(modelName, edgeGroupeNames)
nv, ne, nt = mesh.nv(), mesh.nedges(), mesh.ntris()
print('nv {:d}, ne {:d}, nf {:d}'.format(nv, ne, nt))
C = planeStrainC2222CSR(lam = lamb, mu = mu)
for spaceConstructor in spaceConstructors:
spaceName = spaceConstructor.__name__+'_'
Profile.enable()
tic = time.perf_counter()
pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet,
physName2Neumann=phys2Neumann,
space_constructor=spaceConstructor,
stab=2*mu*0.5)
U = pb.solve(solver)
toc = time.perf_counter()
print(f"assemble +solve: {toc - tic:0.4f} seconds")
Profile.disable()
gp.listPlotFieldTri(*pb.postProStress(U), P0=True, viewname=spaceName+'Stress')
gp.listPlotFieldTri(*pb.postProStrain(U), P0=True, viewname=spaceName+'Strain')
gp.listPlotFieldTri(*pb.postProDisplacement(U), P0=False,
viewname=spaceName+'Displacement', VectorType="Displacement",
DisplacementFactor=1.)
stress = pb.postProStress(U)[1]
neumannEdges, neumannEdgeFlip = mesh.get_edgegroupe('dOmega_other')
if neumannEdges.size:
e2t = mesh.e2t()
stressNeumann = stress[e2t.data[e2t.index[neumannEdges]]]
normals = mesh.enormal(neumannEdges, neumannEdgeFlip)
tractionNeumann = np.einsum("i...jk,i...k->i...j", stressNeumann, normals)
gp.listPlotFieldLine(mesh.e2v_coord()[neumannEdges], tractionNeumann,
P0=True, viewname=spaceName+'bcTraction', VectorType="Arrow")
print('Max neumann norm ', np.max(np.linalg.norm(tractionNeumann.squeeze(), axis=1)))
if popUpGmsh:
gmsh.fltk.run()
gmsh.finalize()
Profile.print_stats()
if __name__ == '__main__':
#example running just one of the unittest :
# import unittest
# # Create a test suite and add your specific test method
# suite = unittest.TestSuite()
# suite.addTest(testElasticity('testPatch'))
# # Run the suite using the standard text test runner
# runner = unittest.TextTestRunner()
# runner.run(suite)
# for solver in ["default", "cholmod", "umfpack"] :
# print(f"### run with Solver {solver}")
# main(element_size=0.05, solver = solver, profile = False, popUpGmsh=False)
main(element_size=0.008, solver = "cholmod", profile = False, popUpGmsh=True)