Coverage for pyxfem/solidmechanics/membrane.py: 44%
510 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1""" Module to test ideas for hyper elastic membrane. For now splitting F in F1@iF0 """
2# pylint disable=invalid-name
3from __future__ import annotations
4import os
5import sys
6from typing import Literal
7import numba
8import numpy as np
9import scipy as sp
10import gmsh
11from pyxfem.core.config import fhalf, IndexType, IndexArray
12import pyxfem.core.mesh_simplex_one_level as sm
13import pyxfem.solidmechanics.hyperelasticity as hyp
14import pyxfem.core.gmsh_mesh_generator as mg
15from pyxfem.core.smallsize_tensor import FASTMATH, npfloat, FloatArray, metric, dmetric, d2metric,\
16 inv, det, vecprodmat
18from pyxfem.core.nonlinear_solver import newton, arc_lenght_predictor_expansion, arc_lenght_predictor_a, arc_lenght_predictor_eval
19import pyxfem.core.sparse_tools as spt
20from pyxfem.core.sparse_tools import regular_block_diag_to_csc, vblock, mblock, split_vblock
21import pyxfem.core.function_spaces as fs
22from pyxfem.core.exec_tools import Profile
23import pyxfem.solidmechanics.boundary_conditions as bc
24import pyxfem.postpro.gmsh_post_pro_interface as gp
26# # Set the maximum number of threads for OpenMP and related backends
27os.environ["OMP_NUM_THREADS"] = "8" # Change '4' to your desired core count
28os.environ["MKL_NUM_THREADS"] = "8"
29os.environ["OPENBLAS_NUM_THREADS"] = "8"
34class AssemblerVector():
35 def __init__(self, Op, dofs:IndexArray|slice = slice(None,None, None)):
36 # pylint: disable=invalid-name
37 self.Op_f = Op[:, dofs]
39 def assemble(self, D):
40 """ Assemble on the defined dofs """
41 return D.reshape(-1)@self.Op_f
43class AssemblerMatrix():
44 def __init__(self, OpLeft, OpRight, dofs:IndexArray|slice = slice(None)):
45 self.OpL_f = OpLeft[dofs,:]
46 self.OpR_f = OpRight[:, dofs]
47 def assemble(self, D):
48 Dcsr = spt.regular_block_diag_to_csc(D)
49 return self.OpL_f@Dcsr@self.OpR_f
51class Membrane():
52 """ Class to represent membrane problem in 3d """
53 # pylint: disable=C0103
54 def __init__(self, law:hyp.HyperElasticMaterial, mesh:sm.SMesh):
55 assert law.hyp == hyp.MaterialHypothesis.PLANESTRESS
56 self._B = None
57 self._Natv = None
58 self._N = None
59 self.law = law
60 self.mesh = mesh
61 self.dofm = fs.VertexSpaceDofManager(mesh)
62 self.free_dofs, self.fixed_dofs, self.fixed_values = [slice(None, None, None), None, None]
63 self.initial_zstretch = 2.
64 XYZ = mesh.vcoord()
65 self.X = XYZ[self.dofm.vid]
66 if self.X.shape[-1] == 2 :
67 print("Assuming a plannar initial membrane")
68 self.X = np.column_stack([self.X, np.zeros(self.X.shape[:-1])])
69 elif self.X.shape[-1] == 3 :
70 print("The initial membrane is assumed to be a curved surface")
71 else :
72 raise ValueError("The mesh vertices must have 2d or 3d coordinates")
73 # x = self.X.copy()
74 # x[...,2] = initial_zstretch*x[...,2]
75 # self.x = x
76 self.F0 = self.getGrad(self.X)
77 self.iF0 = inv(self.F0)
78 self.M = metric(self.F0)
79 self.J0 = np.sqrt(det(self.M))
80 self.iM = inv(self.M)
81 self.setAssembler()
83 def getGrad(self, A:FloatArray, algo:Literal[0,1] = 0)-> FloatArray:
84 """ return the gradient of the nodal vector field stored in A
85 input A[i,j] coordinate j of field A at node j
86 return G(nelem, 3, 2) matrix, such as :
87 G[ie, j, 0] is the derivative of A in elem ie with regard to reference coordinate u
88 G[ie, j, 1] is the derivative of A in elem ie with regard to reference coordinate v
89 """
90 if algo == 0:
91 return (self.getB()@A.flatten()).reshape((self.dofm.nelem, 3, 2))
92 if algo == 1:
93 t2vf = self.dofm.a2dofs
94 Atris = A[t2vf]
95 BelemsT = np.array([[[-1., 1. ,0.],[-1.,0.,1.]]]).swapaxes(-1,-2)
96 return Atris.swapaxes(-1,-2)@BelemsT
98 def getB(self)-> sp.sparse.spmatrix:
99 """ return the sparse matrix B such as B@x return the values of the (FE) local gradient
100 of the vector field x xith regard to the parameters uv of the point in the space
101 of the referent element, for each element
102 """
103 if self._B is None:
104 nelem = self.dofm.nelem
105 dN = np.array([[[-1., 1. ,0.],[-1.,0.,1.]]])
106 dN = np.broadcast_to(dN, (nelem,) + dN.shape)
107 self._B = fs.op2sparse(dN, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 1, 2 )
108 return self._B
110 def getNatVertices(self)-> sp.sparse.spmatrix:
111 """ return the sparse matrix B such as B@x return the values of the (FE) local gradient
112 of the vector field x xith regard to the parameters uv of the point in the space
113 of the referent element, for each element
114 """
115 if self._Natv is None:
116 nelem = self.dofm.nelem
117 N = np.array([[[1., 0. ,0.]], [[0.,1.,0.]], [[0.,0.,1.]] ])
118 N = np.broadcast_to(N, (nelem,) + N.shape)
119 self._Natv = fs.op2sparse(N, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 3, 1 )
120 return self._Natv
122 def getN(self)-> sp.sparse.spmatrix:
123 """ return the sparse matrix N such as N@x return the values of the (FE) at gauss point
124 (u,v) = (1/3, 1/3) (Gauss point of the triangle)
125 """
126 if self._N is None:
127 nelem = self.dofm.nelem
128 N = np.array([[[1./3., 1./3. , 1./3.]]])
129 N = np.broadcast_to(N, (nelem,) + N.shape)
130 self._N = fs.op2sparse(N, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 1, 1 )
131 return self._N
133 def getF(self, x:FloatArray, algo:Literal[0,1] = 0)-> FloatArray:
134 """
135 return an array of size nelem, 3, 3
136 such as F[k].dX = dx where DX is is a tangent vector included
137 in element k in the reference configuration
138 """
139 if algo == 0 :
140 return self.getGrad(x, algo = 0)@self.iF0
141 if algo == 1 :
142 return self.getGrad(x, algo = 1)@self.iF0
144 def surf0(self):
145 """ return the initial surface of the membrane """
146 return 0.5*np.sum(self.J0)
148 def surf(self, x:FloatArray)-> npfloat:
149 """ return the deformed surface of the membrane """
150 J1 = det(self.getGrad(x))
151 return fhalf*np.sum(J1)
153 def getC(self, F1:FloatArray)-> FloatArray:
154 """ compute the C=FT@F component (2x2) for each element
155 Where F1 contain the gradient with regard to local coordinates of the deformed position
156 dx= F1@du
157 """
158 m = metric(F1)
159 return self.iM@m
161 def getdCdF1(self, F1:FloatArray, usenumba = False)-> FloatArray:
162 """ return the derivative of C with regard to F1"""
163 dm = dmetric(F1)
164 if not usenumba:
165 return np.einsum("...ak, ...kbic -> ...abic", self.iM, dm)
166 else :
167 N=F1.shape[0]
168 dCdF1 = np.zeros((N, 2, 2, 3, 2), dtype = npfloat)
169 Membrane._dCdF1_njit(self.iM, dm, dCdF1)
170 return dCdF1
173 # @staticmethod
174 # @numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
175 # def _dCdF1_njit( iM, dm, dCdF1_out):
176 # N = iM.shape[0]
177 # for elem in numba.prange(N): # type ignore # pylint: disable=E1133
178 # iMe, dme, dCdF1e = iM[elem], dm[elem], dCdF1_out[elem]
179 # for a in range(2):
180 # for c in range(2):
181 # for i in range(3):
182 # for d in range(2):
183 # dCdF1e[a,c,i,d] = 0.
184 # for b in range(2):
185 # dCdF1e[a,c,i,d] += iMe[a,b] * dme[b,c,i,d]
186 @staticmethod
187 @numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
188 def _dCdF1_njit(iM, dm, dCdF1_out):
189 N = iM.shape[0]
190 for elem in numba.prange(N): # type ignore # pylint: disable=E1133
191 iMe, dme, dCdF1e = iM[elem], dm[elem], dCdF1_out[elem]
192 for a in range(2):
193 for c in range(2):
194 for i in range(3):
195 for d in range(2):
196 # Use a local accumulator for register speed
197 acc = 0.0
198 for b in range(2):
199 acc += iMe[a,b] * dme[b,c,i,d]
200 # Write to the element array exactly once
201 dCdF1e[a,c,i,d] = acc
203 def getd2Cd2F1(self, F1:FloatArray, usenumba = False)-> FloatArray:
204 """ return the second derivative of C with regard to F1"""
205 d2m = d2metric(F1)
206 #dm = 0.5*(dm+np.swapaxes(dm, -3,-4))
207 return np.einsum("...ak, ...kbicjd -> ...abicjd", self.iM, d2m)
209 def strain_energy(self, *, x:FloatArray|None=None, F1:FloatArray|None = None)-> npfloat:
210 """ return the integrates strain energy """
211 if F1 is None:
212 if x is not None:
213 F1 = self.getGrad(x, algo = 0)
214 else :
215 raise ValueError("in put either x or F1")
216 C = self.getC(F1)
217 J0 = self.J0
218 return fhalf*np.sum(J0*self.law.potentialC(C))
220 def dpotentialdF1(self, F1: FloatArray, usenumba = False) -> FloatArray:
221 """ return the first derivative of the strain enegy with regard to F1 fo each element
222 (a stress measure)
223 """
224 C = self.getC(F1)
225 dCdF1 = self.getdCdF1(F1, usenumba=usenumba)
226 dphidC = 0.5*self.law.PKII(C)
227 def fint_terms():
228 if usenumba :
229 N = F1.shape[0]
230 dphidF1 = np.zeros((N, 3, 2), dtype = npfloat)
231 Membrane._dpotentialdF1_njit(dphidC, dCdF1, self.J0, dphidF1)
232 return dphidF1
233 else:
234 print("###")
235 return np.einsum("...ef, ...efia, ... -> ...ia", dphidC, dCdF1, 0.5*self.J0 )
236 return fint_terms()
238 @staticmethod
239 @numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
240 def _dpotentialdF1_njit( dphidC, dCdF1, J0, d2phid2F1_out):
241 N = dphidC.shape[0]
242 for elem in numba.prange(N): # type ignore # pylint: disable=E1133
243 dphidCe, dCdF1e, J0e, d2phid2F1e = dphidC[elem], dCdF1[elem], J0[elem], d2phid2F1_out[elem]
244 for i in range(3):
245 for a in range(2):
246 d2phid2F1e[i, a] = 0.
247 for e in range(2):
248 for f in range(2):
249 d2phid2F1e[i, a] += 0.5*J0e*dphidCe[e,f]*dCdF1e[e,f, i,a]
251 def d2potentiald2F1(self, F1: FloatArray, usenumba = False) -> FloatArray:
252 """ return the second derivative of the strain enegy with regard to F1 fo each element
253 (a tangent stifness)
254 """
255 C = self.getC(F1)
256 dCdF1 = self.getdCdF1(F1, usenumba = usenumba)
257 d2CdF1_2 = self.getd2Cd2F1(F1)
258 dphidC = 0.5*self.law.PKII(C)
259 d2phid2C = 0.5*self.law.dPKIIdC(C)
261 def kint_terms():
262 if usenumba:
263 N = F1.shape[0]
264 d2phid2CF1 = np.zeros((N,3,2,3,2), dtype=npfloat)
265 Membrane._d2phid2CF1_njit( dCdF1, d2phid2C, dphidC, d2CdF1_2, self.J0, d2phid2CF1)
266 else :
267 print("ggg")
268 dCdF1_d2phid2C_dCdF1 = np.einsum(" ...efia, ...efgh, ...ghjb -> ...iajb",
269 dCdF1, d2phid2C, dCdF1)
270 dphidC_d2Cd2F1 = np.einsum("...ef, ...efiajb , ...-> ...iajb", dphidC, d2CdF1_2, 0.5*self.J0)
271 d2phid2CF1 = dCdF1_d2phid2C_dCdF1 + dphidC_d2Cd2F1
272 return d2phid2CF1
273 return kint_terms()
276 @staticmethod
277 @numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
278 def _d2phid2CF1_njit(dCdF1, d2phid2C, dphidC, d2CdF1_2, J0, d2phid2F1_out):
279 N = dCdF1.shape[0]
280 for elem in numba.prange(N): # type ignore # pylint: disable=E1133
281 d2phid2F1e = d2phid2F1_out[elem]
282 dCdF1e = dCdF1[elem]
283 d2phid2Ce = d2phid2C[elem]
284 dphidCe = dphidC[elem]
285 d2CdF1_2e = d2CdF1_2[elem]
286 # 1. Hoist the scalar weight calculation outside the tensor loops
287 weight = 0.5 * J0[elem]
288 for i in range(3):
289 for a in range(2):
290 for j in range(3):
291 for b in range(2):
292 # 2. Use a local scalar accumulator to stay in CPU registers
293 acc = 0.0
294 for e in range(2):
295 for f in range(2):
296 # Term 1 contribution
297 acc += dphidCe[e,f] * d2CdF1_2e[e,f,i,a,j,b]
298 # 3. Cache the array lookup for Term 2
299 dCdF1e_efia = dCdF1e[e,f,i,a]
300 for g in range(2):
301 for h in range(2):
302 # Term 2 contribution
303 acc += dCdF1e_efia * d2phid2Ce[e,f,g,h] * dCdF1e[g,h,j,b]
304 # Apply the element weight exactly once per final tensor component
305 d2phid2F1e[i,a,j,b] = acc * weight
307 def internalForce(self, *, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
308 """ return the finite element internal force vector """
309 if (x is None) and (F1 is None) :
310 raise ValueError("input either x or F1")
311 if F1 is None and x is not None:
312 F1 = self.getGrad(x, algo = 0)
313 if F1 is not None:
314 dphidF1 = self.dpotentialdF1(F1, usenumba=True)
315 return self.assemblerFint.assemble(dphidF1)
316 #return dphidF1.flatten()@self.getB()
317 raise ValueError("input either x or F1")
319 def internalStiffness(self, *, x:FloatArray|None =None, F1:FloatArray|None=None) \
320 -> sp.sparse.csr_array:
321 """ return the finite element internal stifness as a asparse matrix """
322 if F1 is None:
323 if x is not None :
324 F1 = self.getGrad(x, algo = 0)
325 else :
326 raise ValueError("input either x or F1")
327 #w = 0.5*self.J0.reshape(self.J0.shape + (1,)*4)
328 d2phid2F1 = (self.d2potentiald2F1(F1, usenumba=True)).reshape(F1.shape[:-2]+(6,6))
329 return self.assemblerKint.assemble(d2phid2F1)
331 # def internalStiffnessSensibility(self, *, x:FloatArray, dx:FloatArray) \
332 # -> sp.sparse.csr_array:
333 # """ return the finite element internal stifness sensibility in direction dx, as a asparse matrix """
334 # eps = 1.e-6
335 # xp = x +eps*dx
336 # F1p = self.getGrad(xp)
337 # d2phid2F1p = self.d2potentiald2F1(F1p, usenumba=True)
339 # xm = x - eps*dx
340 # F1m = self.getGrad(xm)
341 # d2phid2F1m = self.d2potentiald2F1(F1m, usenumba=True)
343 # Dd2phidF1_dx = (d2phid2F1p -d2phid2F1m)/2/eps
345 # w = 0.5*self.J0.reshape(self.J0.shape + (1,)*4)
346 # d2phid2F1 = (w*self.d2potentiald2F1(F1, usenumba=True)).reshape(F1.shape[:-2]+(6,6))
348 # d2phid2F1 = spt.regular_block_diag_to_csr(d2phid2F1)
349 # Ke = self.getB().T@d2phid2F1@self.getB()
350 # return Ke
352 def normalPerElem(self, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
353 """ return the normal to each triangle of the deformed membrane """
354 if (x is None) and (F1 is None) :
355 raise ValueError("input either x or F1")
356 if F1 is None and x is not None:
357 F1 = self.getGrad(x, algo = 0)
358 if F1 is not None :
359 return np.cross(F1[...,0], F1[...,1])
360 raise ValueError("input either x or F1")
362 def dnormalPerElem(self, x:FloatArray|None =None, F1:FloatArray|None=None):
363 """ return the derivatives of the normal per element as a function of x,
364 the nodals positions"""
365 if (x is None) and (F1 is None) :
366 raise ValueError("input either x or F1")
367 if F1 is None and x is not None:
368 F1 = self.getGrad(x, algo = 0)
369 if F1 is not None:
370 Q0 = vecprodmat(F1[...,0])
371 Q1 = vecprodmat(F1[...,1])
372 dn_elem = np.zeros(F1.shape[:-2] + (3,6))
373 dn_elem[...,[0, 2, 4]] = -Q1
374 dn_elem[...,[1, 3, 5]] = Q0
375 dn_elem_sparse = regular_block_diag_to_csc(dn_elem)
376 opdnormal = dn_elem_sparse@self.getB()
377 opdnormal.sort_indices()
378 return opdnormal
379 else:
380 raise ValueError("input either x or F1")
382 def pressureForcePerElem(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
383 -> FloatArray:
384 """ return the pressure force per element, as a function of p and x (or p and F1) """
385 return 0.5*p* self.normalPerElem(x = x, F1 = F1)
387 def pressureForce(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
388 """ return the "nodal" pressure force per node, as a function of p and x (or p and F1) """
389 Fe = self.pressureForcePerElem(p=p, x=x, F1=F1)
390 return self.assemblerFext.assemble(Fe)
391 #N = self.getN()
392 #return N.T@Fe.flatten()
394 def pressureStiffnessPerElem(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
395 -> sp.sparse.csr_array:
396 """ return the derivative of the pressure force per element with reagardto x
397 (the deformed nodal position), as a function of p, x (or p F1)"""
398 return 0.5*p*self.dnormalPerElem(x=x, F1 =F1)
400 def pressureStiffness(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
401 -> sp.sparse.csr_array:
402 """ return the finite element internal stiffness as a sparse matrix """
403 if (x is None) and (F1 is None) :
404 raise ValueError("input either x or F1")
405 if F1 is None and x is not None:
406 F1 = self.getGrad(x, algo = 0)
407 if F1 is not None:
408 Q0 = vecprodmat(F1[...,0])
409 Q1 = vecprodmat(F1[...,1])
410 dn_elem = np.zeros(F1.shape[:-2] + (3,6))
411 dn_elem[...,[0, 2, 4]] = -Q1
412 dn_elem[...,[1, 3, 5]] = Q0
413 return self.assemblerKext.assemble(0.5*p*dn_elem)
415 def postProDofVector(self, U:FloatArray, pos:FloatArray|None = None):
416 """ return X and U in A 3d array : U[i,j,k] is the displacement in direction k of node j of
417 element i. Usefull for gmsh post-pro"""
418 if pos is None :
419 assert U.shape == (self.dofm.nscaldof,3)
420 elemvertexVals = U[self.dofm.a2dofs]
421 return self.mesh.t2v_coord(), elemvertexVals
422 else:
423 raise NotImplementedError()
425 def setDirichlet(self, physname2dirichlet:dict):
426 ''' Compute and return the set of free dofs, fixed dofs and values for fixed dofs '''
427 gname2dofid, gname2dofxy = bc.get_dofs_and_interpolation_points(
428 self.mesh, physname2dirichlet.keys(), self.dofm)
429 self.free_dofs, self.fixed_dofs, self.fixed_values = bc.compute_dirichlet(
430 self.dofm, physname2dirichlet, gname2dofid, gname2dofxy)
431 #self.free_dofs = np.asarray(free_dofs, dtype = IndexType)
432 self.setAssembler()
434 def setAssembler(self):
435 free_dofs = self.free_dofs
436 B = self.getB()
437 N = self.getN()
438 self.assemblerFint = AssemblerVector(B, dofs = free_dofs)
439 self.assemblerFext = AssemblerVector(N, dofs = free_dofs)
440 self.assemblerKint = AssemblerMatrix(B.T, B, dofs = free_dofs)
441 self.assemblerKext = AssemblerMatrix(N.T, B, dofs = free_dofs)
443class membrane_operator_dirichlet():
444 def __init__(self, mem:Membrane, free_dofs):
445 self.mem = mem
446 self.free_dofs = free_dofs
447 self.start_position :None | FloatArray = None
448 self.ndofs = len(free_dofs)
450 def set_start(self, start_position: FloatArray):
451 self.start_position = start_position
453 def x2position(self, x):
454 if self.start_position is None:
455 raise ValueError("Start Position not setup")
456 position = self.start_position.flatten()
457 position[self.free_dofs] = x
458 return position.reshape((-1,3))
460 def position2x(self, position):
461 return position.flatten()[self.free_dofs]
463 def compute_equilibrium_res(self, x, p):
464 position = self.x2position(x)
465 mem = self.mem
466 Fint = mem.internalForce(x = position)
467 Fext = mem.pressureForce(p = p, x = position)
468 return Fint-Fext
470 def compute_fext(self, x, p=1.):
471 position = self.x2position(x)
472 mem = self.mem
473 return mem.pressureForce(p = p, x = position)
475 def compute_Kext(self, x, p=1.):
476 position = self.x2position(x)
477 mem = self.mem
478 return mem.pressureStiffness(p= p, x=position)
480 def compute_equilibrium_jac(self, x, p):
481 ''' compute the jacobian : derivative of the residual with regard to x'''
482 position = self.x2position(x)
483 mem = self.mem
484 Kint = mem.internalStiffness(x=position)
485 Kext = mem.pressureStiffness(p= p, x=position)
486 K = Kint-Kext
487 # rcond = spt.rcond_estimate(K)
488 # print("rcond", rcond)
489 return K
491 def compute_djacdx_u(self, x , p, u):
492 ''' directional derivative of the jacobian in the direction u'''
493 eps = 1.e-5
494 Kt_p_u = self.compute_equilibrium_jac(x = x + eps*u, p = p)@u
495 Kt_m_u = self.compute_equilibrium_jac(x = x - eps*u, p = p)@u
496 return (Kt_p_u-Kt_m_u)/2./eps
497 ## Alternative implementation using second order derivative of the residual, assuming the residual is zero at x
498 # eps = 1.e-5
499 # r_p = self.compute_equilibrium_res(x = x + eps*u, p = p)
500 # r_m = self.compute_equilibrium_res(x = x - eps*u, p = p)
501 # return (r_p + r_m)/eps/eps
503class membrane_solver_imposed_displacement(membrane_operator_dirichlet):
504 def __init__(self, mem, free_dofs):
505 super().__init__(mem, free_dofs)
506 #self.memop = membrane_operator_dirichlet(mem, free_dofs)
507 self.start_position = None
508 self.p = None
510 def get_fun(self):
511 return lambda x: self.compute_equilibrium_res(x, self.p)
513 def get_jac(self):
514 return lambda x: self.compute_equilibrium_jac(x, self.p)
516 def get_linsolver(self):
517 return sp.sparse.linalg.spsolve
519 def solve(self, start_position, p, *newton_args, **newton_kwargs):
520 self.set_start(start_position)
521 self.p = p
522 x0 = self.position2x(start_position)
523 newton_kwargs['jac'] = self.get_jac()
524 newton_kwargs['linsolver'] = self.get_linsolver()
525 # tol = 1.e-6, atol = 1.e-6, nitermax =10,
526 # res = newton(self.get_fun(), x0, jac = self.get_jac(),
527 # linsolver = self.get_linsolver(), rtol = 1.e-6, atol = 1.e-6, nitermax =10, *newton_args, **newton_kwargs)
528 res = newton(self.get_fun(), x0, *newton_args, **newton_kwargs)
529 if not res.success:
530 print(f"newton failed after {res.nit}")
531 return res.success, self.x2position(res.x)
533class membrane_solver_pressureloading_newton(membrane_operator_dirichlet):
534 def __init__(self, mem, free_dofs):
535 super().__init__(mem, free_dofs)
536 self.start_position = None
537 self.p = None
539 def get_fun(self):
540 return lambda x: self.compute_equilibrium_res(x, self.p)
542 def get_jac(self):
543 return lambda x: self.compute_equilibrium_jac(x, self.p)
545 def get_linsolver(self):
546 return sp.sparse.linalg.spsolve
548 def solve(self, start_position, p, *newton_args, **newton_kwargs):
549 self.set_start(start_position)
550 self.p = p
551 x0 = self.position2x(start_position)
552 newton_kwargs['jac'] = self.get_jac()
553 newton_kwargs['linsolver'] = self.get_linsolver()
554 res = newton(self.get_fun(), x0, *newton_args, **newton_kwargs)
555 if not res.success:
556 print(f"newton failed after {res.nit}")
557 return res.success, self.x2position(res.x)
559class membrane_solver_pressureloading_arclenght(membrane_operator_dirichlet):
560 def __init__(self, mem, free_dofs):
561 super().__init__(mem, free_dofs)
562 # self.p0 = None
563 # self.x0 = None
564 self.alpha = None
565 self.du_prev = None
566 self.dp_prev = None
568 def xp2_block_x(self, x, p):
569 return vblock(x, np.array([p]))
571 def block_x2xp(self, block_x):
572 x = block_x[0]
573 p = block_x[1][0]
574 du = x - self.x0
575 dp = p - self.p0
576 return x, p, du, dp
578 def get_fun(self, a, alpha, scal_equ, scal_arc):
579 def fun(block_x):
580 x, p, du, dp = self.block_x2xp(block_x)
581 #plot_disp(self.mem, self.x2position(x), "bl")
582 res_equ = self.compute_equilibrium_res(x, p)
583 n = len(x)
584 res_arc = du.dot(du)/n + alpha*dp**2 - a**2
585 return vblock(scal_equ*res_equ , np.array([scal_arc*res_arc]))
586 return fun
588 def get_jac(self, alpha, scal_arc, scal_equ):
589 def jac(block_x):
590 n = self.ndofs
591 x, p, du, dp = self.block_x2xp(block_x)
592 k00 = scal_equ*self.compute_equilibrium_jac(x, p)
593 k01 = scal_equ*self.compute_fext(x, p = -1.).reshape((-1,1))
594 k10 = scal_arc*2.*du.reshape((1, -1))/n
595 k11 = scal_arc*2.*alpha*dp.reshape((1,1))
596 k = mblock(k00, k01, k10, k11)
597 return k
598 return jac
600 def get_linsolver(self):
601 def solve(block_k, block_x):
602 return spt.block_linsolve(block_k, block_x)
603 return solve
605 def solve(self, start_position, p0, a = None, delta = None, afac=None, amax = 0.1, alpha =None, pred_order = 2, *newton_args, **newton_kwargs):
606 du_prev = self.du_prev
607 dp_prev = self.dp_prev
608 if alpha is None:
609 alpha = self.alpha # after the first step, we always keep the same alpha
610 self.start_position = start_position.copy()
611 x0 = self.position2x(start_position)
612 self.x0 = x0
613 self.p0 = p0
614 K = self.compute_equilibrium_jac(x0, p0)
615 dus, dls, alpha = arc_lenght_predictor_expansion( drdx = K,
616 drdl = -self.compute_fext(x0, 1.),
617 fd2rdx2_u = lambda u : self.compute_djacdx_u(x = x0, p=p0, u=u),
618 d2rdxdl = -self.compute_Kext(x0, 1.),
619 du_prev = du_prev,
620 dl_prev = dp_prev,
621 alpha = alpha,
622 pred_order= pred_order
623 )
624 if a is None:
625 a = arc_lenght_predictor_a(dus, dls, alpha, delta, afac, amax)
626 dx, dp = arc_lenght_predictor_eval(dus, dls, a)
628 x_pred = x0 + dx
629 p_pred = p0 + dp
630 self.alpha = alpha
631 scal_arc = 1./a/a
632 scal_equ = 1./np.max(K.diagonal())/a
633 block_x_pred = self.xp2_block_x(x_pred ,p_pred)
634 fun = self.get_fun(a=a, alpha=alpha, scal_equ=scal_equ, scal_arc= scal_arc)
635 newton_kwargs['jac'] = self.get_jac(alpha, scal_arc, scal_equ)
636 newton_kwargs['linsolver'] = self.get_linsolver()
637 newton_kwargs['inplace_sub'] = spt.block_isub
638 newton_kwargs['norm'] = spt.block_norm
639 newton_kwargs.setdefault('nitermax', 4)
640 newton_kwargs.setdefault('fail_on_increasing_res', True)
641 res = newton(fun, block_x_pred, *newton_args, **newton_kwargs)
642 if not res.success:
643 print(f"newton failed after {res.nit}")
644 x, p = split_vblock(res.x)
645 p = p[0]
646 self.du_prev = x - self.x0
647 self.dp_prev = p - self.p0
648 position = self.x2position(x)
649 return res.success, res.nit, position, p, a
651def plot_disp(mem, x, name):
652 gmsh.option.setNumber("PostProcessing.AnimationCycle", 1)
653 gmsh.option.setNumber("Mesh.SurfaceEdges", 0)
654 trisXYZ, trisDISP = mem.postProDofVector(x - mem.X)
655 gp.listPlotFieldTri(trisXYZ, trisDISP, P0=False,
656 viewname=f"{name}_disp", VectorType="Displacement",
657 DisplacementFactor=1.)
659def cylindrical_membrane_under_pressure_load():
660 # Build hyperelastic law
661 E = 1.
662 nu = 0.2
663 law1 = hyp.CompressibleNeoHookeen(E = E, nu = nu, hyp = hyp.PLANESTRESS)
664 #initialize gmsh
665 argv = sys.argv
666 gmsh.initialize(argv)
667 # build the mesh
668 radius = 1.
669 height = 1.
670 ref = 2
671 n_theta_quarter=10*ref
672 n_z = 20*ref
673 #modelname, physdic = mg.cylinder(radius, height, hscale)
674 modelname, physdic = mg.structured_cylinder(radius, height, n_theta_quarter=n_theta_quarter, n_z = n_z)
675 mesh1 = sm.gmshmodel2smesh(modelname, physdic.keys(), dim=3)
676 nv, nt = (mesh1.nv(), mesh1.ntris())
677 print(f"nv: {nv:d}, nt {nt:d}")
678 # print(physdic) show groups name
679 # Create Membrane model
680 mem = Membrane(law1, mesh1)
681 diriBC = {'bottom': {'x':0., 'y':0., 'z':0.}, 'top': {'y':0.,'x':0., 'z':1.} }
682 gname2dofid, gname2dofxy = bc.get_dofs_and_interpolation_points(
683 mesh1, diriBC.keys(), mem.dofm)
684 free_dofs, fixed_dofs, fixed_values = bc.compute_dirichlet(mem.dofm, diriBC, gname2dofid, gname2dofxy)
685 x = mem.X.copy()
686 p = 0 #1e3
687 # plot initial position
688 plot_disp(mem, x, "init")
689 # set up the solver
690 solver = membrane_solver_pressureloading_newton(mem, free_dofs=free_dofs)
691 solver.mem.setDirichlet(diriBC)
692 solver2 = membrane_solver_pressureloading_arclenght(mem, free_dofs=free_dofs)
693 # First load increment : initial pull in the z direction
694 lz = 1.1 # initial stretch in z direction
695 inc = f"initial z pull {lz}"
696 x[:,2] *= lz
697 Profile.enable()
698 converged, x = solver.solve(start_position = x, p = 0., nitermax = 20, rtol = 1.e-8, atol = 1.e-8)
699 Profile.disable()
700 if not converged :
701 print(f"inc {inc} not converged. Stop")
702 raise ValueError("solver.solve did not converge")
703 print(f"inc {inc} converged ")
704 plot_disp(mem, x, "zpull")
705 inc = 0
706 p = 0.
707 delta0 = 1.e-4
708 delta = delta0
709 a0 = 0.05
710 a = a0
711 amax = 1.
712 afac = 1.
713 max_inc = 20
715 def solve(x, p, afac):
716 return solver2.solve(start_position = x, p0 = p, a = None, delta = delta, afac = afac, amax = amax,
717 pred_order = 2, nitermax = 10, verbose = True, atol = 1.e-10, rtol = 1.e-7)
719 while a > 1.e-6*a0 and delta > 1.e-6 * delta0 and inc < max_inc:
720 xstart = x.copy()
721 pstart = p
722 while 1:
723 Profile.enable()
724 converged, nit, x, p, a = solve(x, p, afac)
725 Profile.disable()
726 if converged :
727 inc += 1
728 print(f"inc {inc:03d} a {a:.3e}, p {p:.3e} converged in {nit} iterations")
729 #a = a/( nit/3)
730 plot_disp(mem, x, f"inc_{inc:03d}_p_{p:.3e}")
731 afac =1.
732 break
733 if not converged :
734 x = xstart.copy()
735 p = pstart
736 afac = afac/2.
737 print(f"inc {inc:03d} a {a:.3e}, p {p:.3e} not converged after {nit} iterations")
739if __name__ == "__main__":
740 Profile.doprofile = True
741 cylindrical_membrane_under_pressure_load()
742 Profile.print_stats()
744 # Launch gmsh API
745 gmsh.fltk.run()
749# inc 001 a 2.000e-02, p 3.021e+04 converged in 3 iterations
750# inc 002 a 2.000e-02, p 5.870e+04 converged in 2 iterations
751# inc 003 a 3.000e-02, p 9.887e+04 converged in 2 iterations
752# inc 004 a 4.500e-02, p 1.541e+05 converged in 3 iterations
753# inc 005 a 4.500e-02, p 2.038e+05 converged in 2 iterations
754# inc 006 a 6.750e-02, p 2.705e+05 converged in 3 iterations
755# inc 007 a 6.750e-02, p 3.276e+05 converged in 3 iterations
756# inc 008 a 6.750e-02, p 3.755e+05 converged in 3 iterations
759# ncalls tottime percall cumtime percall filename:lineno(function)
760# 1 0.000 0.000 0.385 0.385 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:462(solve)
761# 1 0.001 0.001 0.385 0.385 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/nonlinear_solver.py:13(newton)
762# 2 0.000 0.000 0.213 0.107 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:134(spsolve)
763# 2 0.210 0.105 0.210 0.105 {built-in method scipy.sparse.linalg._dsolve._superlu.gssv}
764# 2 0.000 0.000 0.137 0.069 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:457(<lambda>)
765# 2 0.000 0.000 0.137 0.069 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:437(compute_equilibrium_jac)
766# 2 0.001 0.000 0.114 0.057 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:288(internalStiffness)
767# 2 0.000 0.000 0.071 0.035 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:230(d2potentiald2F1)
768# 24 0.000 0.000 0.044 0.002 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_base.py:788(_matmul_dispatch)
769# 21 0.000 0.000 0.043 0.002 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_base.py:919(__matmul__)
770# 8 0.000 0.000 0.042 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_compressed.py:415(_matmul_sparse)
771# 29 0.000 0.000 0.031 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numpy/_core/einsumfunc.py:1244(einsum)
772# 29 0.031 0.001 0.031 0.001 {built-in method numpy._core._multiarray_umath.c_einsum}
773# 3 0.000 0.000 0.030 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:454(<lambda>)
774# 3 0.000 0.000 0.030 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:414(compute_equilibrium_res)
775# 5 0.000 0.000 0.029 0.006 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:180(getdCdF1)
776# 3 0.000 0.000 0.028 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:276(internalForce)
777# 3 0.000 0.000 0.026 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:200(dpotentialdF1)
778# 8 0.025 0.003 0.025 0.003 {built-in method scipy.sparse._sparsetools.csr_matmat}
779# 5 0.011 0.002 0.023 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/core/smallsize_tensor.py:153(dmetric)
780# 2 0.010 0.005 0.023 0.012 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:240(kint_terms)
781# 2 0.000 0.000 0.021 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/dispatcher.py:344(_compile_for_args)
782# 2 0.000 0.000 0.020 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/dispatcher.py:862(compile)
783# 2 0.000 0.000 0.020 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/caching.py:713(load_overload)
784# 47/41 0.000 0.000 0.018 0.000 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_compressed.py:30(__init__)
785# 2 0.000 0.000 0.018 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:372(pressureStiffness)
786# 2 0.000 0.000 0.017 0.008 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/caching.py:724(_load_overload)
787# 2 0.000 0.000 0.016 0.008 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/caching.py:477(rebuild)
788# 2 0.000 0.000 0.016 0.008 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/compiler.py:227(_rebuild)
789# 2 0.004 0.002 0.013 0.007 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:500(dPKIIdC)
790# 8 0.012 0.001 0.012 0.001 {built-in method scipy.sparse._sparsetools.csr_matmat_maxnnz}
791# 2 0.000 0.000 0.012 0.006 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/codegen.py:1168(unserialize_library)
792# 2 0.000 0.000 0.012 0.006 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/codegen.py:936(_unserialize)
793# 5 0.004 0.001 0.011 0.002 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:474(PKII)
794# 357/355 0.009 0.000 0.011 0.000 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/llvmlite/binding/ffi.py:210(__call__)
795# 4 0.004 0.001 0.011 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_coo.py:383(_coo_to_compressed)
796# 2 0.000 0.000 0.011 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:335(dnormalPerElem)
797# 2 0.000 0.000 0.010 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/core/sparse_tools.py:54(regular_block_diag_to_csr)
798# 7 0.000 0.000 0.009 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:430(C33_planestress)
799# 7 0.009 0.001 0.009 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/special/_lambertw.py:6(lambertw)
800# 3 0.001 0.000 0.008 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:207(fint_terms)
801# 2 0.000 0.000 0.008 0.004 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:186(getd2Cd2F1)
802# 4 0.007 0.002 0.007 0.002 {built-in method scipy.sparse._sparsetools.coo_tocsr}
803# 2 0.000 0.000 0.006 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/codegen.py:792(_finalize_final_module)
804# 2 0.000 0.000 0.006 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/codegen.py:1007(_finalize_specific)
805# 8 0.000 0.000 0.006 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/numba/core/codegen.py:1069(wrapper)
806# 2 0.000 0.000 0.006 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/llvmlite/binding/executionengine.py:94(finalize_object)
807# 2 0.000 0.000 0.005 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/core/sparse_tools.py:72(regular_block_diag_to_csr_v2)
808# 5 0.001 0.000 0.005 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:172(getC)
809# 10 0.000 0.000 0.004 0.000 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/.venv/lib/python3.12/site-packages/scipy/sparse/_base.py:451(asformat)
812# The initial membrane is assumed to be a curved surface
813# 0 0.00440687981439875
814# 1 0.0007090203268567412
815# 2 7.995410736685281e-07
816# inc initial z pull 1.1 converged
817# 32.22576608864403 40.067667319006084
818# 0 3.4310598963475734
819# 1 0.2935674403395363
820# 2 0.01827820455432498
821# 3 0.0002816147637323649
822# 4 9.036886193951496e-08
823# inc 001 a 5.000e-02, 3.426e-01, p 8.403e-01 converged in 4 iterations
824# 32.225766088644036 42.96510692768051
825# 0 0.36539286272311966
826# 1 0.029378488740486756
827# 2 0.00046313210328144475
828# 3 5.788205975849805e-07
829# 4 2.2703713941695908e-13
830# inc 002 a 5.000e-02, 3.347e-01, p 1.071e+00 converged in 4 iterations
831# 32.22576608864404 34.1846510483814
832# 0 0.31496387608195736
833# 1 2.1415147031750412
834# newton failed after 1
835# inc 002 p 1.071e+00 not converged after 1 iterations
836# 32.22576608864404 34.1846510483814
837# 0 0.19127906246138682
838# 1 2.221491558000836
839# newton failed after 1
840# inc 002 p 1.071e+00 not converged after 1 iterations
841# 32.22576608864404 34.1846510483814
842# 0 0.11732610750436605
843# 1 41.41408859196228
844# newton failed after 1
845# inc 002 p 1.071e+00 not converged after 1 iterations
846# 32.22576608864404 34.1846510483814
847# 0 0.07214008537315857
848# 1 0.883440686809872
849# newton failed after 1
850# inc 002 p 1.071e+00 not converged after 1 iterations
851# 32.22576608864404 34.1846510483814
852# 0 0.044518527641326175
853# 1 0.08708125448554915
854# newton failed after 1
855# inc 002 p 1.071e+00 not converged after 1 iterations
856# 32.22576608864404 34.1846510483814
857# 0 0.0275742096485467
858# 1 0.017107224332357146
859# 2 8.750321381946716e-05
860# 3 3.631399800953893e-09
861# 4 1.6966004494157222e-14
862# inc 003 a 1.563e-03, 1.138e-01, p 1.091e+00 converged in 4 iterations
863# 32.22576608864403 2584.6899989463413
864# 0 4.660028981796859
865# 1 1.0883871245650745
866# 2 0.24051205442431087
867# 3 0.07826732787845271
868# 4 0.015068006598989521
869# 5 0.0016632268498018973
870# 6 1.168588040514311e-05
871# 7 3.1948361642820114e-10
872# inc 004 a 1.563e-03, 2.691e-02, p 1.090e+00 converged in 7 iterations
873# 32.22576608864404 921.3652756624772
874# 0 0.9095356315019443
875# 1 0.13097402901504196
876# 2 0.006575055412094633
877# 3 8.462309488545935e-05
878# 4 4.5710282161511506e-07
879# 5 1.958417172507177e-12