Coverage for pyxfem/solidmechanics/mechanics_solver.py: 0%
247 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Tue Feb 27 07:43:08 2024
6@author: nchevaug
7"""
8from __future__ import annotations
9import os
11# # Set the maximum number of threads for OpenMP and related backends
12os.environ["OMP_NUM_THREADS"] = "8" # Change '4' to your desired core count
13os.environ["MKL_NUM_THREADS"] = "8"
14os.environ["OPENBLAS_NUM_THREADS"] = "8"
16from warnings import warn
17import time
18import unittest
19import numpy as np
20import scipy as sp
21try:
22 import sksparse.cholmod # pylint: disable=E0401
23 have_cholmod = True
24except ModuleNotFoundError:
25 have_cholmod = False
26 warn('cholmod not installed, using default scipy spsolve')
28import gmsh
30import pyxfem.core.function_spaces as fs
31import pyxfem.core.quadrature as quad
32import pyxfem.core.sparse_tools as st
33from pyxfem.core.mesh_simplex_one_level import SMesh, gmshmodel2smesh
34import pyxfem.postpro.gmsh_post_pro_interface as gp
35from pyxfem.core.exec_tools import Profile
36import pyxfem.core.gmsh_mesh_generator as mg
37from pyxfem.solidmechanics.boundary_conditions import compute_dirichlet, compute_neumann, \
38 get_dofs_and_interpolation_points, compute_dirichlet_fem_p1_vect2d
39from pyxfem.solidmechanics.elasticity import planeStrainC2222CSR, YoungPoisson2LambdaMu
41def getSpaceEntity(edgeGroupeNames, space):
42 ''' get the space's edge id corresponding to the edgeGroupe Name'''
43 edgeGroupeNames2EId = {}
44 for name in edgeGroupeNames:
45 edgeGroupeNames2EId[name] = space.m.getEdgeGroupe(name)[0]
46 return edgeGroupeNames2EId
48def computeVolumeLoad(space, source):
49 """ compute the vector F such as F.U = int_Omega f.u_h dOmega_h """
50 assert isinstance(source, np.ndarray)
51 quadSource = quad.T3_gauss(0)
52 N = space.operator_dof2val_tri(quadSource.uv)
53 J = space.get_tris_mapping().J
54 wSource = quadSource.w[0]*source[np.newaxis, :]*J[:, np.newaxis]
55 Fint = N.T.dot(wSource.flatten())
56 return Fint
58class ElasticitySolver():
59 """ A class to set up and compute a 2d linear elasticity problem with FEM method """
60 def __init__(self, mesh:SMesh, C=planeStrainC2222CSR(lam =12., mu =7.),
61 node_id_to_imposed_value:None|dict=None,
62 physName2Dirichlet:None|dict=None,
63 physName2Neumann:None|dict=None,
64 source=np.array([0., 0.]),
65 space_constructor:type[fs.FEMSpaceP1]|type[fs.FEMSpaceP1NC]=fs.FEMSpaceP1,
66 stab=0.):
67 self.quadBilinForm = quad.T3_gauss(0)
68 self.quadAtCorners = quad.T3_nodes()
69 self.C = C
70 self.stab = stab
71 self.mesh = mesh
72 self.setModel(space_constructor)
73 if node_id_to_imposed_value is not None :
74 if physName2Dirichlet is not None :
75 warn("can't impose dirichlet BC both from physName and list of nodes, ignoring physName")
76 if space_constructor != fs.FEMSpaceP1 :
77 warn("imposing dirichlet at node value make sens only for p1 fem space. Expect errors ...")
78 nv = mesh.nv()
79 idvx = node_id_to_imposed_value.get("idvx", np.array([], dtype = np.int64))
80 idvy = node_id_to_imposed_value.get("idvy", np.array([], dtype = np.int64))
81 vx = node_id_to_imposed_value.get("vx", np.array([], dtype = np.int64))
82 vy = node_id_to_imposed_value.get("vy", np.array([], dtype = np.int64))
83 self.freeDofs, self.fixedDofs, self.fixedValues = \
84 compute_dirichlet_fem_p1_vect2d(nvertex=nv, nodex= idvx, valuex= vx, nodey = idvy, valuey=vy)
85 else:
86 self.physName2Dirichlet:dict = {} if physName2Dirichlet is None else physName2Dirichlet
87 gname2dof, gname2dofxy = get_dofs_and_interpolation_points(mesh,
88 self.physName2Dirichlet.keys(), self.space)
89 self.freeDofs, self.fixedDofs, self.fixedValues = \
90 compute_dirichlet(self.space, self.physName2Dirichlet, gname2dof, gname2dofxy)
91 self.physName2Neumann:dict = {} if physName2Neumann is None else physName2Neumann
92 self.Fneumann = compute_neumann(self.space, self.physName2Neumann)
93 self.Fint = computeVolumeLoad(self.space, source)
94 self.Bfree = self.B[:, self.freeDofs]
95 self.Bfixed = self.B[:, self.fixedDofs]
96 nbfixed = len(self.fixedDofs)
97 if nbfixed < 3:
98 warn(f' Only {nbfixed:d} fixed dofs: not enough to prevent rigid-body motion in {compute_dirichlet.__name__}')
99 if isinstance(self.space, fs.FEMSpaceP1NC):
100 self.setStabilisation(self.stab, quad.Edge_gauss(2))
101 def setD(self):
102 ''' set the sparse stifness operator (K = BtDB)'''
103 J = self.space.get_tris_mapping().J
104 nf = len(J)
105 C = self.C.todense().flatten()
106 wD = np.einsum("i, ij -> ij", self.quadBilinForm.w[0]*J,
107 C[np.newaxis, :]).reshape((-1, 4, 4))
108 self.wD = st.regular_block_diag_to_csr(wD)
109 D = np.broadcast_to(C, (nf, 16)).reshape((-1, 4, 4))
110 self.D = st.regular_block_diag_to_csr(D)
111 def setModel(self, spaceConstructor:type[fs.FEMSpaceP1]|type[fs.FEMSpaceP1NC]):
112 self.space = spaceConstructor("DisplacementSpace", self.mesh, vecdim=2)
113 self.B = self.space.operator_dof2grad_tri(self.quadBilinForm.uv)
114 self.NatVertices = self.space.operator_dof2val_tri(self.quadAtCorners.uv)
115 self.setD()
116 def setStabilisation(self, s, quadEdge):
117 if not isinstance(self.space, fs.FEMSpaceP1NC):
118 raise ValueError("Stabilisation term only implemented for P1NC")
119 space = self.space
120 lE = space.m.elength()
121 sTri2edge = space.operator_trival2edge().dot(0.5*space.get_tris_mapping().J).reshape((-1, 2))
122 hE = np.sum(sTri2edge, axis=1)/2/lE
123 N = space.operator_edgetrace(quadEdge)
124 self.jumpOp = N[1::2] - N[:-1:2]
125 weights = s*((lE/hE)[:, np.newaxis] * quadEdge.w).reshape(-1)
126 DjumpBlock = (weights[:, np.newaxis]*np.eye(2).flatten()[np.newaxis, :]).reshape(-1, 2, 2)
127 self.Djump = st.regular_block_diag_to_csr(DjumpBlock)
128 internalEdges = self.space.m.intern_edges()
129 vecDim = self.space.vecdim
130 nptPerEdge = quadEdge.s.size
131 rowsPerEdge = vecDim*nptPerEdge
132 internalEdgesDofs = (vecDim*internalEdges[:, np.newaxis]
133 +np.arange(vecDim)[np.newaxis, : ]).flatten()
134 stabDofs = np.concatenate([internalEdgesDofs, self.fixedDofs])
135 stabRows = np.sort(((rowsPerEdge*(stabDofs//vecDim)
136 +stabDofs%vecDim)[:, np.newaxis]
137 +nptPerEdge*np.arange(vecDim)[np.newaxis, :]).flatten())
138 Nedges = self.jumpOp[stabRows, :]
139 self.NedgesFree = Nedges[:, self.freeDofs]
140 self.NedgesFixed = Nedges[:, self.fixedDofs]
141 self.Djump = self.Djump[stabRows, :][:, stabRows]
142 self.SFree = self.NedgesFree.T.dot(self.Djump.dot(self.NedgesFree))
144 def assembleK(self):
145 if not hasattr(self, "Kfree"):
146 Bfree = self.Bfree
147 self.Kfree = Bfree.T.dot(self.wD).dot(Bfree)
148 if isinstance(self.space, fs.FEMSpaceP1NC):
149 self.Kfree = self.Kfree + self.SFree
150 return self.Kfree
152 def assembleFU(self):
153 Fneumann = self.Fneumann
154 U = np.zeros(self.space.size())
155 U[self.fixedDofs] = self.fixedValues
156 Bfree = self.Bfree
157 Bfixed = self.Bfixed
158 gradUdiri = Bfixed.dot(self.fixedValues)
159 Fdiri = -Bfree.T.dot(self.wD.dot(gradUdiri))
160 F = Fdiri + self.Fneumann[self.freeDofs] + self.Fint[self.freeDofs]
161 if isinstance(self.space, fs.FEMSpaceP1NC):
162 SDiri = -self.NedgesFree.T.dot(self.Djump.dot(self.NedgesFixed.dot(self.fixedValues)))
163 F = F+SDiri
164 return F,U
166 def assemble(self):
167 Kfree = self.assembleK()
168 F,U = self.assembleFU()
169 return Kfree, F, U
170 # def assembleMass(self):
172 def solve(self, solver="default"):
173 """ solve the problem using the solver to choose between "default" : spsolve or "cholmod" : "sksparse.cholmod.cholesky" """
174 Kfree, F, U = self.assemble()
175 def lin_solve():
176 if solver == "default" or not have_cholmod:
177 U[self.freeDofs] = sp.sparse.linalg.spsolve(Kfree, F, use_umfpack=False)
178 elif solver == "umfpack":
179 U[self.freeDofs] = sp.sparse.linalg.spsolve(Kfree, F)
180 elif solver == "cholmod":
181 print("cholmod solver")
182 if not have_cholmod :
183 raise ValueError
184 chol = sksparse.cholmod.cho_factor(Kfree)
185 U[self.freeDofs] = chol.solve(F)
186 lin_solve()
187 # if compute_rcond:
188 # diag = chol.D()
189 # rcond = (np.min(diag)/np.max(diag))**2
190 return U
191 def postProDisplacement(self, U):
192 return self.mesh.t2v_coord(), self.NatVertices.dot(U).reshape((-1, 3, 2))
193 def postProStress(self, U):
194 return self.mesh.t2v_coord(), self.getStress(U)
195 def getStress(self, U):
196 ''' return per element stress /per gp stress '''
197 return self.D.dot(self.B.dot(U)).reshape((-1, 2, 2))
198 def postProStrain(self, U):
199 return self.mesh.t2v_coord(), self.B.dot(U).reshape((-1, 2, 2))
200 def L2StressError(self, U, exactStressFun, quadrature=quad.T3_gauss(0)):
201 uv = quadrature.uv
202 weight = quadrature.w
203 npt = len(uv)
204 B = self.space.operator_dof2grad_tri(uv)
205 gradNum = B.dot(U).reshape((-1, npt, 2, 2))
206 C = self.C.todense().reshape((1,1,2,2,2,2))
207 stressNum = np.einsum('egijkl, egkl-> egij', C, gradNum)
208 NP1 = fs.FEMSpaceP1.N(uv)
209 trisQuadPoints = np.einsum('ijk, ikl-> ijl', NP1[np.newaxis, :, :], self.mesh.t2v_coord())
210 stressExact = exactStressFun(trisQuadPoints)
211 stressError2 = np.sum((stressNum - stressExact)**2, axis = (2,3))
212 J = self.space.get_tris_mapping().J
213 wStressError2 = J[:, np.newaxis]*weight[np.newaxis, :]*stressError2
214 return np.sqrt(np.sum(wStressError2))
216class testElasticity(unittest.TestCase):
217 def testPatch(self):
218 """ - patch test on a squared load by Dirchlet and then Neumann """
219 gmsh.initialize()
220 gmsh.option.setNumber('General.Terminal', 0)
221 L = 0.7
222 E = 12.2
223 nu = 0.3
224 ux = 0.98
225 lamb, mu = YoungPoisson2LambdaMu(E, nu)
226 C = planeStrainC2222CSR(lam = lamb, mu =mu)
227 phys2Dirichlet = {'left': {'x':0.}, 'bottom': {'y': 0.}, 'right': {'x': ux}}
228 modelName, physical_dict = mg.square(length=L, relh=0.3)
229 mesh = gmshmodel2smesh(modelName, phys2Dirichlet.keys())
230 pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet,
231 space_constructor=fs.FEMSpaceP1)
232 U = pb.solve()
233 ___, stress = pb.postProStress(U)
234 e11 = ux/L
235 s11 = 4*mu*(lamb+mu)/(lamb+2*mu)*e11
236 stressExact = np.array([[s11, 0.], [0.,0.]])
237 error = np.linalg.norm(stress - stressExact[np.newaxis, np.newaxis] )
238 self.assertTrue(error < 1.e-9)
239 phys2Dirichlet = {'left': {'x':0.}, 'bottom': {'y': 0.}}
240 phys2Neumann = {'right': {'x':s11}}
241 pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet, physName2Neumann=phys2Neumann,
242 space_constructor=fs.FEMSpaceP1)
243 e22 = -lamb*e11/(lamb+2*mu)
244 strainExact = np.array([[e11, 0.], [0., e22]])
245 U = pb.solve()
246 ___, strain = pb.postProStrain(U)
247 error = np.linalg.norm(strain - strainExact[np.newaxis, np.newaxis] )
248 self.assertTrue(error < 1.e-9)
249 gmsh.finalize()
251def main(*, element_size=0.05, profile=True, solver = "default", popUpGmsh = True):
252 """ main Program and example : elasticity on L shape domain """
253 E = 1.
254 nu = 0.4999
255 lamb, mu = YoungPoisson2LambdaMu(E, nu)
256 meshgenerator = mg.l_shape
257 #meshgenerator = mg.square
259 phys2Dirichlet = {'dOmega_left': {'x':0., 'y':0.}}
260 phys2Neumann = {'dOmega_right':{'x': 0.01}}
261 #phys2Dirichlet = {'left': {'x':0., 'y':0.}, 'right':{'x': 0.01, 'y':0.}}
262 #phys2Neumann = {}
263 spaceConstructors = [fs.FEMSpaceP1, fs.FEMSpaceP1NC]
264 Profile.doprofile = profile
265 gmsh.initialize()
266 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
267 gmsh.option.setNumber('Mesh.SurfaceEdges', 1 if element_size > 0.1 else 0)
268 gmsh.option.setNumber('General.Terminal', False)
269 modelName, physical_dict = meshgenerator(element_size)
270 edgeGroupeNames = list(phys2Dirichlet.keys())+ list(phys2Neumann.keys())
271 mesh = gmshmodel2smesh(modelName, edgeGroupeNames)
272 nv, ne, nt = mesh.nv(), mesh.nedges(), mesh.ntris()
273 print('nv {:d}, ne {:d}, nf {:d}'.format(nv, ne, nt))
274 C = planeStrainC2222CSR(lam = lamb, mu = mu)
275 for spaceConstructor in spaceConstructors:
276 spaceName = spaceConstructor.__name__+'_'
277 Profile.enable()
278 tic = time.perf_counter()
279 pb = ElasticitySolver(mesh, C, physName2Dirichlet=phys2Dirichlet,
280 physName2Neumann=phys2Neumann,
281 space_constructor=spaceConstructor,
282 stab=2*mu*0.5)
285 U = pb.solve(solver)
286 toc = time.perf_counter()
287 print(f"assemble +solve: {toc - tic:0.4f} seconds")
288 Profile.disable()
289 gp.listPlotFieldTri(*pb.postProStress(U), P0=True, viewname=spaceName+'Stress')
290 gp.listPlotFieldTri(*pb.postProStrain(U), P0=True, viewname=spaceName+'Strain')
291 gp.listPlotFieldTri(*pb.postProDisplacement(U), P0=False,
292 viewname=spaceName+'Displacement', VectorType="Displacement",
293 DisplacementFactor=1.)
295 stress = pb.postProStress(U)[1]
296 neumannEdges, neumannEdgeFlip = mesh.get_edgegroupe('dOmega_other')
297 if neumannEdges.size:
298 e2t = mesh.e2t()
299 stressNeumann = stress[e2t.data[e2t.index[neumannEdges]]]
300 normals = mesh.enormal(neumannEdges, neumannEdgeFlip)
301 tractionNeumann = np.einsum("i...jk,i...k->i...j", stressNeumann, normals)
302 gp.listPlotFieldLine(mesh.e2v_coord()[neumannEdges], tractionNeumann,
303 P0=True, viewname=spaceName+'bcTraction', VectorType="Arrow")
304 print('Max neumann norm ', np.max(np.linalg.norm(tractionNeumann.squeeze(), axis=1)))
306 if popUpGmsh:
307 gmsh.fltk.run()
308 gmsh.finalize()
309 Profile.print_stats()
311if __name__ == '__main__':
312 #example running just one of the unittest :
313 # import unittest
315 # # Create a test suite and add your specific test method
316 # suite = unittest.TestSuite()
317 # suite.addTest(testElasticity('testPatch'))
319 # # Run the suite using the standard text test runner
320 # runner = unittest.TextTestRunner()
321 # runner.run(suite)
323 # for solver in ["default", "cholmod", "umfpack"] :
324 # print(f"### run with Solver {solver}")
325 # main(element_size=0.05, solver = solver, profile = False, popUpGmsh=False)
327 main(element_size=0.008, solver = "cholmod", profile = False, popUpGmsh=True)