""" Module to test ideas for hyper elastic membrane. For now splitting F in F1@iF0 """
# pylint disable=invalid-name
from __future__ import annotations
import os
import sys
from typing import Literal
import numba
import numpy as np
import scipy as sp
import gmsh
from pyxfem.core.config import fhalf, IndexType, IndexArray
import pyxfem.core.mesh_simplex_one_level as sm
import pyxfem.solidmechanics.hyperelasticity as hyp
import pyxfem.core.gmsh_mesh_generator as mg
from pyxfem.core.smallsize_tensor import FASTMATH, npfloat, FloatArray, metric, dmetric, d2metric,\
inv, det, vecprodmat
from pyxfem.core.nonlinear_solver import newton, arc_lenght_predictor_expansion, arc_lenght_predictor_a, arc_lenght_predictor_eval
import pyxfem.core.sparse_tools as spt
from pyxfem.core.sparse_tools import regular_block_diag_to_csc, vblock, mblock, split_vblock
import pyxfem.core.function_spaces as fs
from pyxfem.core.exec_tools import Profile
import pyxfem.solidmechanics.boundary_conditions as bc
import pyxfem.postpro.gmsh_post_pro_interface as gp
# # 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"
[docs]
class AssemblerVector():
def __init__(self, Op, dofs:IndexArray|slice = slice(None,None, None)):
# pylint: disable=invalid-name
self.Op_f = Op[:, dofs]
[docs]
def assemble(self, D):
""" Assemble on the defined dofs """
return D.reshape(-1)@self.Op_f
[docs]
class AssemblerMatrix():
def __init__(self, OpLeft, OpRight, dofs:IndexArray|slice = slice(None)):
self.OpL_f = OpLeft[dofs,:]
self.OpR_f = OpRight[:, dofs]
[docs]
def assemble(self, D):
Dcsr = spt.regular_block_diag_to_csc(D)
return self.OpL_f@Dcsr@self.OpR_f
[docs]
class Membrane():
""" Class to represent membrane problem in 3d """
# pylint: disable=C0103
def __init__(self, law:hyp.HyperElasticMaterial, mesh:sm.SMesh):
assert law.hyp == hyp.MaterialHypothesis.PLANESTRESS
self._B = None
self._Natv = None
self._N = None
self.law = law
self.mesh = mesh
self.dofm = fs.VertexSpaceDofManager(mesh)
self.free_dofs, self.fixed_dofs, self.fixed_values = [slice(None, None, None), None, None]
self.initial_zstretch = 2.
XYZ = mesh.vcoord()
self.X = XYZ[self.dofm.vid]
if self.X.shape[-1] == 2 :
print("Assuming a plannar initial membrane")
self.X = np.column_stack([self.X, np.zeros(self.X.shape[:-1])])
elif self.X.shape[-1] == 3 :
print("The initial membrane is assumed to be a curved surface")
else :
raise ValueError("The mesh vertices must have 2d or 3d coordinates")
# x = self.X.copy()
# x[...,2] = initial_zstretch*x[...,2]
# self.x = x
self.F0 = self.getGrad(self.X)
self.iF0 = inv(self.F0)
self.M = metric(self.F0)
self.J0 = np.sqrt(det(self.M))
self.iM = inv(self.M)
self.setAssembler()
[docs]
def getGrad(self, A:FloatArray, algo:Literal[0,1] = 0)-> FloatArray:
""" return the gradient of the nodal vector field stored in A
input A[i,j] coordinate j of field A at node j
return G(nelem, 3, 2) matrix, such as :
G[ie, j, 0] is the derivative of A in elem ie with regard to reference coordinate u
G[ie, j, 1] is the derivative of A in elem ie with regard to reference coordinate v
"""
if algo == 0:
return (self.getB()@A.flatten()).reshape((self.dofm.nelem, 3, 2))
if algo == 1:
t2vf = self.dofm.a2dofs
Atris = A[t2vf]
BelemsT = np.array([[[-1., 1. ,0.],[-1.,0.,1.]]]).swapaxes(-1,-2)
return Atris.swapaxes(-1,-2)@BelemsT
[docs]
def getB(self)-> sp.sparse.spmatrix:
""" return the sparse matrix B such as B@x return the values of the (FE) local gradient
of the vector field x xith regard to the parameters uv of the point in the space
of the referent element, for each element
"""
if self._B is None:
nelem = self.dofm.nelem
dN = np.array([[[-1., 1. ,0.],[-1.,0.,1.]]])
dN = np.broadcast_to(dN, (nelem,) + dN.shape)
self._B = fs.op2sparse(dN, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 1, 2 )
return self._B
[docs]
def getNatVertices(self)-> sp.sparse.spmatrix:
""" return the sparse matrix B such as B@x return the values of the (FE) local gradient
of the vector field x xith regard to the parameters uv of the point in the space
of the referent element, for each element
"""
if self._Natv is None:
nelem = self.dofm.nelem
N = np.array([[[1., 0. ,0.]], [[0.,1.,0.]], [[0.,0.,1.]] ])
N = np.broadcast_to(N, (nelem,) + N.shape)
self._Natv = fs.op2sparse(N, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 3, 1 )
return self._Natv
[docs]
def getN(self)-> sp.sparse.spmatrix:
""" return the sparse matrix N such as N@x return the values of the (FE) at gauss point
(u,v) = (1/3, 1/3) (Gauss point of the triangle)
"""
if self._N is None:
nelem = self.dofm.nelem
N = np.array([[[1./3., 1./3. , 1./3.]]])
N = np.broadcast_to(N, (nelem,) + N.shape)
self._N = fs.op2sparse(N, self.dofm.nscaldof, self.dofm.a2dofs, 3, 3, 1, 1 )
return self._N
[docs]
def getF(self, x:FloatArray, algo:Literal[0,1] = 0)-> FloatArray:
"""
return an array of size nelem, 3, 3
such as F[k].dX = dx where DX is is a tangent vector included
in element k in the reference configuration
"""
if algo == 0 :
return self.getGrad(x, algo = 0)@self.iF0
if algo == 1 :
return self.getGrad(x, algo = 1)@self.iF0
[docs]
def surf0(self):
""" return the initial surface of the membrane """
return 0.5*np.sum(self.J0)
[docs]
def surf(self, x:FloatArray)-> npfloat:
""" return the deformed surface of the membrane """
J1 = det(self.getGrad(x))
return fhalf*np.sum(J1)
[docs]
def getC(self, F1:FloatArray)-> FloatArray:
""" compute the C=FT@F component (2x2) for each element
Where F1 contain the gradient with regard to local coordinates of the deformed position
dx= F1@du
"""
m = metric(F1)
return self.iM@m
[docs]
def getdCdF1(self, F1:FloatArray, usenumba = False)-> FloatArray:
""" return the derivative of C with regard to F1"""
dm = dmetric(F1)
if not usenumba:
return np.einsum("...ak, ...kbic -> ...abic", self.iM, dm)
else :
N=F1.shape[0]
dCdF1 = np.zeros((N, 2, 2, 3, 2), dtype = npfloat)
Membrane._dCdF1_njit(self.iM, dm, dCdF1)
return dCdF1
# @staticmethod
# @numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
# def _dCdF1_njit( iM, dm, dCdF1_out):
# N = iM.shape[0]
# for elem in numba.prange(N): # type ignore # pylint: disable=E1133
# iMe, dme, dCdF1e = iM[elem], dm[elem], dCdF1_out[elem]
# for a in range(2):
# for c in range(2):
# for i in range(3):
# for d in range(2):
# dCdF1e[a,c,i,d] = 0.
# for b in range(2):
# dCdF1e[a,c,i,d] += iMe[a,b] * dme[b,c,i,d]
@staticmethod
@numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
def _dCdF1_njit(iM, dm, dCdF1_out):
N = iM.shape[0]
for elem in numba.prange(N): # type ignore # pylint: disable=E1133
iMe, dme, dCdF1e = iM[elem], dm[elem], dCdF1_out[elem]
for a in range(2):
for c in range(2):
for i in range(3):
for d in range(2):
# Use a local accumulator for register speed
acc = 0.0
for b in range(2):
acc += iMe[a,b] * dme[b,c,i,d]
# Write to the element array exactly once
dCdF1e[a,c,i,d] = acc
[docs]
def getd2Cd2F1(self, F1:FloatArray, usenumba = False)-> FloatArray:
""" return the second derivative of C with regard to F1"""
d2m = d2metric(F1)
#dm = 0.5*(dm+np.swapaxes(dm, -3,-4))
return np.einsum("...ak, ...kbicjd -> ...abicjd", self.iM, d2m)
[docs]
def strain_energy(self, *, x:FloatArray|None=None, F1:FloatArray|None = None)-> npfloat:
""" return the integrates strain energy """
if F1 is None:
if x is not None:
F1 = self.getGrad(x, algo = 0)
else :
raise ValueError("in put either x or F1")
C = self.getC(F1)
J0 = self.J0
return fhalf*np.sum(J0*self.law.potentialC(C))
[docs]
def dpotentialdF1(self, F1: FloatArray, usenumba = False) -> FloatArray:
""" return the first derivative of the strain enegy with regard to F1 fo each element
(a stress measure)
"""
C = self.getC(F1)
dCdF1 = self.getdCdF1(F1, usenumba=usenumba)
dphidC = 0.5*self.law.PKII(C)
def fint_terms():
if usenumba :
N = F1.shape[0]
dphidF1 = np.zeros((N, 3, 2), dtype = npfloat)
Membrane._dpotentialdF1_njit(dphidC, dCdF1, self.J0, dphidF1)
return dphidF1
else:
print("###")
return np.einsum("...ef, ...efia, ... -> ...ia", dphidC, dCdF1, 0.5*self.J0 )
return fint_terms()
@staticmethod
@numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
def _dpotentialdF1_njit( dphidC, dCdF1, J0, d2phid2F1_out):
N = dphidC.shape[0]
for elem in numba.prange(N): # type ignore # pylint: disable=E1133
dphidCe, dCdF1e, J0e, d2phid2F1e = dphidC[elem], dCdF1[elem], J0[elem], d2phid2F1_out[elem]
for i in range(3):
for a in range(2):
d2phid2F1e[i, a] = 0.
for e in range(2):
for f in range(2):
d2phid2F1e[i, a] += 0.5*J0e*dphidCe[e,f]*dCdF1e[e,f, i,a]
[docs]
def d2potentiald2F1(self, F1: FloatArray, usenumba = False) -> FloatArray:
""" return the second derivative of the strain enegy with regard to F1 fo each element
(a tangent stifness)
"""
C = self.getC(F1)
dCdF1 = self.getdCdF1(F1, usenumba = usenumba)
d2CdF1_2 = self.getd2Cd2F1(F1)
dphidC = 0.5*self.law.PKII(C)
d2phid2C = 0.5*self.law.dPKIIdC(C)
def kint_terms():
if usenumba:
N = F1.shape[0]
d2phid2CF1 = np.zeros((N,3,2,3,2), dtype=npfloat)
Membrane._d2phid2CF1_njit( dCdF1, d2phid2C, dphidC, d2CdF1_2, self.J0, d2phid2CF1)
else :
print("ggg")
dCdF1_d2phid2C_dCdF1 = np.einsum(" ...efia, ...efgh, ...ghjb -> ...iajb",
dCdF1, d2phid2C, dCdF1)
dphidC_d2Cd2F1 = np.einsum("...ef, ...efiajb , ...-> ...iajb", dphidC, d2CdF1_2, 0.5*self.J0)
d2phid2CF1 = dCdF1_d2phid2C_dCdF1 + dphidC_d2Cd2F1
return d2phid2CF1
return kint_terms()
@staticmethod
@numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
def _d2phid2CF1_njit(dCdF1, d2phid2C, dphidC, d2CdF1_2, J0, d2phid2F1_out):
N = dCdF1.shape[0]
for elem in numba.prange(N): # type ignore # pylint: disable=E1133
d2phid2F1e = d2phid2F1_out[elem]
dCdF1e = dCdF1[elem]
d2phid2Ce = d2phid2C[elem]
dphidCe = dphidC[elem]
d2CdF1_2e = d2CdF1_2[elem]
# 1. Hoist the scalar weight calculation outside the tensor loops
weight = 0.5 * J0[elem]
for i in range(3):
for a in range(2):
for j in range(3):
for b in range(2):
# 2. Use a local scalar accumulator to stay in CPU registers
acc = 0.0
for e in range(2):
for f in range(2):
# Term 1 contribution
acc += dphidCe[e,f] * d2CdF1_2e[e,f,i,a,j,b]
# 3. Cache the array lookup for Term 2
dCdF1e_efia = dCdF1e[e,f,i,a]
for g in range(2):
for h in range(2):
# Term 2 contribution
acc += dCdF1e_efia * d2phid2Ce[e,f,g,h] * dCdF1e[g,h,j,b]
# Apply the element weight exactly once per final tensor component
d2phid2F1e[i,a,j,b] = acc * weight
[docs]
def internalForce(self, *, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
""" return the finite element internal force vector """
if (x is None) and (F1 is None) :
raise ValueError("input either x or F1")
if F1 is None and x is not None:
F1 = self.getGrad(x, algo = 0)
if F1 is not None:
dphidF1 = self.dpotentialdF1(F1, usenumba=True)
return self.assemblerFint.assemble(dphidF1)
#return dphidF1.flatten()@self.getB()
raise ValueError("input either x or F1")
[docs]
def internalStiffness(self, *, x:FloatArray|None =None, F1:FloatArray|None=None) \
-> sp.sparse.csr_array:
""" return the finite element internal stifness as a asparse matrix """
if F1 is None:
if x is not None :
F1 = self.getGrad(x, algo = 0)
else :
raise ValueError("input either x or F1")
#w = 0.5*self.J0.reshape(self.J0.shape + (1,)*4)
d2phid2F1 = (self.d2potentiald2F1(F1, usenumba=True)).reshape(F1.shape[:-2]+(6,6))
return self.assemblerKint.assemble(d2phid2F1)
# def internalStiffnessSensibility(self, *, x:FloatArray, dx:FloatArray) \
# -> sp.sparse.csr_array:
# """ return the finite element internal stifness sensibility in direction dx, as a asparse matrix """
# eps = 1.e-6
# xp = x +eps*dx
# F1p = self.getGrad(xp)
# d2phid2F1p = self.d2potentiald2F1(F1p, usenumba=True)
# xm = x - eps*dx
# F1m = self.getGrad(xm)
# d2phid2F1m = self.d2potentiald2F1(F1m, usenumba=True)
# Dd2phidF1_dx = (d2phid2F1p -d2phid2F1m)/2/eps
# w = 0.5*self.J0.reshape(self.J0.shape + (1,)*4)
# d2phid2F1 = (w*self.d2potentiald2F1(F1, usenumba=True)).reshape(F1.shape[:-2]+(6,6))
# d2phid2F1 = spt.regular_block_diag_to_csr(d2phid2F1)
# Ke = self.getB().T@d2phid2F1@self.getB()
# return Ke
[docs]
def normalPerElem(self, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
""" return the normal to each triangle of the deformed membrane """
if (x is None) and (F1 is None) :
raise ValueError("input either x or F1")
if F1 is None and x is not None:
F1 = self.getGrad(x, algo = 0)
if F1 is not None :
return np.cross(F1[...,0], F1[...,1])
raise ValueError("input either x or F1")
[docs]
def dnormalPerElem(self, x:FloatArray|None =None, F1:FloatArray|None=None):
""" return the derivatives of the normal per element as a function of x,
the nodals positions"""
if (x is None) and (F1 is None) :
raise ValueError("input either x or F1")
if F1 is None and x is not None:
F1 = self.getGrad(x, algo = 0)
if F1 is not None:
Q0 = vecprodmat(F1[...,0])
Q1 = vecprodmat(F1[...,1])
dn_elem = np.zeros(F1.shape[:-2] + (3,6))
dn_elem[...,[0, 2, 4]] = -Q1
dn_elem[...,[1, 3, 5]] = Q0
dn_elem_sparse = regular_block_diag_to_csc(dn_elem)
opdnormal = dn_elem_sparse@self.getB()
opdnormal.sort_indices()
return opdnormal
else:
raise ValueError("input either x or F1")
[docs]
def pressureForcePerElem(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
-> FloatArray:
""" return the pressure force per element, as a function of p and x (or p and F1) """
return 0.5*p* self.normalPerElem(x = x, F1 = F1)
[docs]
def pressureForce(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)-> FloatArray:
""" return the "nodal" pressure force per node, as a function of p and x (or p and F1) """
Fe = self.pressureForcePerElem(p=p, x=x, F1=F1)
return self.assemblerFext.assemble(Fe)
#N = self.getN()
#return N.T@Fe.flatten()
[docs]
def pressureStiffnessPerElem(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
-> sp.sparse.csr_array:
""" return the derivative of the pressure force per element with reagardto x
(the deformed nodal position), as a function of p, x (or p F1)"""
return 0.5*p*self.dnormalPerElem(x=x, F1 =F1)
[docs]
def pressureStiffness(self, *, p:float, x:FloatArray|None =None, F1:FloatArray|None=None)\
-> sp.sparse.csr_array:
""" return the finite element internal stiffness as a sparse matrix """
if (x is None) and (F1 is None) :
raise ValueError("input either x or F1")
if F1 is None and x is not None:
F1 = self.getGrad(x, algo = 0)
if F1 is not None:
Q0 = vecprodmat(F1[...,0])
Q1 = vecprodmat(F1[...,1])
dn_elem = np.zeros(F1.shape[:-2] + (3,6))
dn_elem[...,[0, 2, 4]] = -Q1
dn_elem[...,[1, 3, 5]] = Q0
return self.assemblerKext.assemble(0.5*p*dn_elem)
[docs]
def postProDofVector(self, U:FloatArray, pos:FloatArray|None = None):
""" return X and U in A 3d array : U[i,j,k] is the displacement in direction k of node j of
element i. Usefull for gmsh post-pro"""
if pos is None :
assert U.shape == (self.dofm.nscaldof,3)
elemvertexVals = U[self.dofm.a2dofs]
return self.mesh.t2v_coord(), elemvertexVals
else:
raise NotImplementedError()
[docs]
def setDirichlet(self, physname2dirichlet:dict):
''' Compute and return the set of free dofs, fixed dofs and values for fixed dofs '''
gname2dofid, gname2dofxy = bc.get_dofs_and_interpolation_points(
self.mesh, physname2dirichlet.keys(), self.dofm)
self.free_dofs, self.fixed_dofs, self.fixed_values = bc.compute_dirichlet(
self.dofm, physname2dirichlet, gname2dofid, gname2dofxy)
#self.free_dofs = np.asarray(free_dofs, dtype = IndexType)
self.setAssembler()
[docs]
def setAssembler(self):
free_dofs = self.free_dofs
B = self.getB()
N = self.getN()
self.assemblerFint = AssemblerVector(B, dofs = free_dofs)
self.assemblerFext = AssemblerVector(N, dofs = free_dofs)
self.assemblerKint = AssemblerMatrix(B.T, B, dofs = free_dofs)
self.assemblerKext = AssemblerMatrix(N.T, B, dofs = free_dofs)
[docs]
class membrane_operator_dirichlet():
def __init__(self, mem:Membrane, free_dofs):
self.mem = mem
self.free_dofs = free_dofs
self.start_position :None | FloatArray = None
self.ndofs = len(free_dofs)
[docs]
def set_start(self, start_position: FloatArray):
self.start_position = start_position
[docs]
def x2position(self, x):
if self.start_position is None:
raise ValueError("Start Position not setup")
position = self.start_position.flatten()
position[self.free_dofs] = x
return position.reshape((-1,3))
[docs]
def position2x(self, position):
return position.flatten()[self.free_dofs]
[docs]
def compute_equilibrium_res(self, x, p):
position = self.x2position(x)
mem = self.mem
Fint = mem.internalForce(x = position)
Fext = mem.pressureForce(p = p, x = position)
return Fint-Fext
[docs]
def compute_fext(self, x, p=1.):
position = self.x2position(x)
mem = self.mem
return mem.pressureForce(p = p, x = position)
[docs]
def compute_Kext(self, x, p=1.):
position = self.x2position(x)
mem = self.mem
return mem.pressureStiffness(p= p, x=position)
[docs]
def compute_equilibrium_jac(self, x, p):
''' compute the jacobian : derivative of the residual with regard to x'''
position = self.x2position(x)
mem = self.mem
Kint = mem.internalStiffness(x=position)
Kext = mem.pressureStiffness(p= p, x=position)
K = Kint-Kext
# rcond = spt.rcond_estimate(K)
# print("rcond", rcond)
return K
[docs]
def compute_djacdx_u(self, x , p, u):
''' directional derivative of the jacobian in the direction u'''
eps = 1.e-5
Kt_p_u = self.compute_equilibrium_jac(x = x + eps*u, p = p)@u
Kt_m_u = self.compute_equilibrium_jac(x = x - eps*u, p = p)@u
return (Kt_p_u-Kt_m_u)/2./eps
## Alternative implementation using second order derivative of the residual, assuming the residual is zero at x
# eps = 1.e-5
# r_p = self.compute_equilibrium_res(x = x + eps*u, p = p)
# r_m = self.compute_equilibrium_res(x = x - eps*u, p = p)
# return (r_p + r_m)/eps/eps
[docs]
class membrane_solver_imposed_displacement(membrane_operator_dirichlet):
def __init__(self, mem, free_dofs):
super().__init__(mem, free_dofs)
#self.memop = membrane_operator_dirichlet(mem, free_dofs)
self.start_position = None
self.p = None
[docs]
def get_fun(self):
return lambda x: self.compute_equilibrium_res(x, self.p)
[docs]
def get_jac(self):
return lambda x: self.compute_equilibrium_jac(x, self.p)
[docs]
def get_linsolver(self):
return sp.sparse.linalg.spsolve
[docs]
def solve(self, start_position, p, *newton_args, **newton_kwargs):
self.set_start(start_position)
self.p = p
x0 = self.position2x(start_position)
newton_kwargs['jac'] = self.get_jac()
newton_kwargs['linsolver'] = self.get_linsolver()
# tol = 1.e-6, atol = 1.e-6, nitermax =10,
# res = newton(self.get_fun(), x0, jac = self.get_jac(),
# linsolver = self.get_linsolver(), rtol = 1.e-6, atol = 1.e-6, nitermax =10, *newton_args, **newton_kwargs)
res = newton(self.get_fun(), x0, *newton_args, **newton_kwargs)
if not res.success:
print(f"newton failed after {res.nit}")
return res.success, self.x2position(res.x)
[docs]
class membrane_solver_pressureloading_newton(membrane_operator_dirichlet):
def __init__(self, mem, free_dofs):
super().__init__(mem, free_dofs)
self.start_position = None
self.p = None
[docs]
def get_fun(self):
return lambda x: self.compute_equilibrium_res(x, self.p)
[docs]
def get_jac(self):
return lambda x: self.compute_equilibrium_jac(x, self.p)
[docs]
def get_linsolver(self):
return sp.sparse.linalg.spsolve
[docs]
def solve(self, start_position, p, *newton_args, **newton_kwargs):
self.set_start(start_position)
self.p = p
x0 = self.position2x(start_position)
newton_kwargs['jac'] = self.get_jac()
newton_kwargs['linsolver'] = self.get_linsolver()
res = newton(self.get_fun(), x0, *newton_args, **newton_kwargs)
if not res.success:
print(f"newton failed after {res.nit}")
return res.success, self.x2position(res.x)
[docs]
class membrane_solver_pressureloading_arclenght(membrane_operator_dirichlet):
def __init__(self, mem, free_dofs):
super().__init__(mem, free_dofs)
# self.p0 = None
# self.x0 = None
self.alpha = None
self.du_prev = None
self.dp_prev = None
[docs]
def xp2_block_x(self, x, p):
return vblock(x, np.array([p]))
[docs]
def block_x2xp(self, block_x):
x = block_x[0]
p = block_x[1][0]
du = x - self.x0
dp = p - self.p0
return x, p, du, dp
[docs]
def get_fun(self, a, alpha, scal_equ, scal_arc):
def fun(block_x):
x, p, du, dp = self.block_x2xp(block_x)
#plot_disp(self.mem, self.x2position(x), "bl")
res_equ = self.compute_equilibrium_res(x, p)
n = len(x)
res_arc = du.dot(du)/n + alpha*dp**2 - a**2
return vblock(scal_equ*res_equ , np.array([scal_arc*res_arc]))
return fun
[docs]
def get_jac(self, alpha, scal_arc, scal_equ):
def jac(block_x):
n = self.ndofs
x, p, du, dp = self.block_x2xp(block_x)
k00 = scal_equ*self.compute_equilibrium_jac(x, p)
k01 = scal_equ*self.compute_fext(x, p = -1.).reshape((-1,1))
k10 = scal_arc*2.*du.reshape((1, -1))/n
k11 = scal_arc*2.*alpha*dp.reshape((1,1))
k = mblock(k00, k01, k10, k11)
return k
return jac
[docs]
def get_linsolver(self):
def solve(block_k, block_x):
return spt.block_linsolve(block_k, block_x)
return solve
[docs]
def solve(self, start_position, p0, a = None, delta = None, afac=None, amax = 0.1, alpha =None, pred_order = 2, *newton_args, **newton_kwargs):
du_prev = self.du_prev
dp_prev = self.dp_prev
if alpha is None:
alpha = self.alpha # after the first step, we always keep the same alpha
self.start_position = start_position.copy()
x0 = self.position2x(start_position)
self.x0 = x0
self.p0 = p0
K = self.compute_equilibrium_jac(x0, p0)
dus, dls, alpha = arc_lenght_predictor_expansion( drdx = K,
drdl = -self.compute_fext(x0, 1.),
fd2rdx2_u = lambda u : self.compute_djacdx_u(x = x0, p=p0, u=u),
d2rdxdl = -self.compute_Kext(x0, 1.),
du_prev = du_prev,
dl_prev = dp_prev,
alpha = alpha,
pred_order= pred_order
)
if a is None:
a = arc_lenght_predictor_a(dus, dls, alpha, delta, afac, amax)
dx, dp = arc_lenght_predictor_eval(dus, dls, a)
x_pred = x0 + dx
p_pred = p0 + dp
self.alpha = alpha
scal_arc = 1./a/a
scal_equ = 1./np.max(K.diagonal())/a
block_x_pred = self.xp2_block_x(x_pred ,p_pred)
fun = self.get_fun(a=a, alpha=alpha, scal_equ=scal_equ, scal_arc= scal_arc)
newton_kwargs['jac'] = self.get_jac(alpha, scal_arc, scal_equ)
newton_kwargs['linsolver'] = self.get_linsolver()
newton_kwargs['inplace_sub'] = spt.block_isub
newton_kwargs['norm'] = spt.block_norm
newton_kwargs.setdefault('nitermax', 4)
newton_kwargs.setdefault('fail_on_increasing_res', True)
res = newton(fun, block_x_pred, *newton_args, **newton_kwargs)
if not res.success:
print(f"newton failed after {res.nit}")
x, p = split_vblock(res.x)
p = p[0]
self.du_prev = x - self.x0
self.dp_prev = p - self.p0
position = self.x2position(x)
return res.success, res.nit, position, p, a
[docs]
def plot_disp(mem, x, name):
gmsh.option.setNumber("PostProcessing.AnimationCycle", 1)
gmsh.option.setNumber("Mesh.SurfaceEdges", 0)
trisXYZ, trisDISP = mem.postProDofVector(x - mem.X)
gp.listPlotFieldTri(trisXYZ, trisDISP, P0=False,
viewname=f"{name}_disp", VectorType="Displacement",
DisplacementFactor=1.)
[docs]
def cylindrical_membrane_under_pressure_load():
# Build hyperelastic law
E = 1.
nu = 0.2
law1 = hyp.CompressibleNeoHookeen(E = E, nu = nu, hyp = hyp.PLANESTRESS)
#initialize gmsh
argv = sys.argv
gmsh.initialize(argv)
# build the mesh
radius = 1.
height = 1.
ref = 2
n_theta_quarter=10*ref
n_z = 20*ref
#modelname, physdic = mg.cylinder(radius, height, hscale)
modelname, physdic = mg.structured_cylinder(radius, height, n_theta_quarter=n_theta_quarter, n_z = n_z)
mesh1 = sm.gmshmodel2smesh(modelname, physdic.keys(), dim=3)
nv, nt = (mesh1.nv(), mesh1.ntris())
print(f"nv: {nv:d}, nt {nt:d}")
# print(physdic) show groups name
# Create Membrane model
mem = Membrane(law1, mesh1)
diriBC = {'bottom': {'x':0., 'y':0., 'z':0.}, 'top': {'y':0.,'x':0., 'z':1.} }
gname2dofid, gname2dofxy = bc.get_dofs_and_interpolation_points(
mesh1, diriBC.keys(), mem.dofm)
free_dofs, fixed_dofs, fixed_values = bc.compute_dirichlet(mem.dofm, diriBC, gname2dofid, gname2dofxy)
x = mem.X.copy()
p = 0 #1e3
# plot initial position
plot_disp(mem, x, "init")
# set up the solver
solver = membrane_solver_pressureloading_newton(mem, free_dofs=free_dofs)
solver.mem.setDirichlet(diriBC)
solver2 = membrane_solver_pressureloading_arclenght(mem, free_dofs=free_dofs)
# First load increment : initial pull in the z direction
lz = 1.1 # initial stretch in z direction
inc = f"initial z pull {lz}"
x[:,2] *= lz
Profile.enable()
converged, x = solver.solve(start_position = x, p = 0., nitermax = 20, rtol = 1.e-8, atol = 1.e-8)
Profile.disable()
if not converged :
print(f"inc {inc} not converged. Stop")
raise ValueError("solver.solve did not converge")
print(f"inc {inc} converged ")
plot_disp(mem, x, "zpull")
inc = 0
p = 0.
delta0 = 1.e-4
delta = delta0
a0 = 0.05
a = a0
amax = 1.
afac = 1.
max_inc = 20
def solve(x, p, afac):
return solver2.solve(start_position = x, p0 = p, a = None, delta = delta, afac = afac, amax = amax,
pred_order = 2, nitermax = 10, verbose = True, atol = 1.e-10, rtol = 1.e-7)
while a > 1.e-6*a0 and delta > 1.e-6 * delta0 and inc < max_inc:
xstart = x.copy()
pstart = p
while 1:
Profile.enable()
converged, nit, x, p, a = solve(x, p, afac)
Profile.disable()
if converged :
inc += 1
print(f"inc {inc:03d} a {a:.3e}, p {p:.3e} converged in {nit} iterations")
#a = a/( nit/3)
plot_disp(mem, x, f"inc_{inc:03d}_p_{p:.3e}")
afac =1.
break
if not converged :
x = xstart.copy()
p = pstart
afac = afac/2.
print(f"inc {inc:03d} a {a:.3e}, p {p:.3e} not converged after {nit} iterations")
if __name__ == "__main__":
Profile.doprofile = True
cylindrical_membrane_under_pressure_load()
Profile.print_stats()
# Launch gmsh API
gmsh.fltk.run()
# inc 001 a 2.000e-02, p 3.021e+04 converged in 3 iterations
# inc 002 a 2.000e-02, p 5.870e+04 converged in 2 iterations
# inc 003 a 3.000e-02, p 9.887e+04 converged in 2 iterations
# inc 004 a 4.500e-02, p 1.541e+05 converged in 3 iterations
# inc 005 a 4.500e-02, p 2.038e+05 converged in 2 iterations
# inc 006 a 6.750e-02, p 2.705e+05 converged in 3 iterations
# inc 007 a 6.750e-02, p 3.276e+05 converged in 3 iterations
# inc 008 a 6.750e-02, p 3.755e+05 converged in 3 iterations
# ncalls tottime percall cumtime percall filename:lineno(function)
# 1 0.000 0.000 0.385 0.385 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:462(solve)
# 1 0.001 0.001 0.385 0.385 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/nonlinear_solver.py:13(newton)
# 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)
# 2 0.210 0.105 0.210 0.105 {built-in method scipy.sparse.linalg._dsolve._superlu.gssv}
# 2 0.000 0.000 0.137 0.069 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:457(<lambda>)
# 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)
# 2 0.001 0.000 0.114 0.057 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:288(internalStiffness)
# 2 0.000 0.000 0.071 0.035 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:230(d2potentiald2F1)
# 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)
# 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__)
# 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)
# 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)
# 29 0.031 0.001 0.031 0.001 {built-in method numpy._core._multiarray_umath.c_einsum}
# 3 0.000 0.000 0.030 0.010 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:454(<lambda>)
# 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)
# 5 0.000 0.000 0.029 0.006 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:180(getdCdF1)
# 3 0.000 0.000 0.028 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:276(internalForce)
# 3 0.000 0.000 0.026 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:200(dpotentialdF1)
# 8 0.025 0.003 0.025 0.003 {built-in method scipy.sparse._sparsetools.csr_matmat}
# 5 0.011 0.002 0.023 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/core/smallsize_tensor.py:153(dmetric)
# 2 0.010 0.005 0.023 0.012 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:240(kint_terms)
# 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)
# 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)
# 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)
# 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__)
# 2 0.000 0.000 0.018 0.009 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:372(pressureStiffness)
# 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)
# 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)
# 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)
# 2 0.004 0.002 0.013 0.007 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:500(dPKIIdC)
# 8 0.012 0.001 0.012 0.001 {built-in method scipy.sparse._sparsetools.csr_matmat_maxnnz}
# 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)
# 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)
# 5 0.004 0.001 0.011 0.002 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:474(PKII)
# 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__)
# 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)
# 2 0.000 0.000 0.011 0.005 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:335(dnormalPerElem)
# 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)
# 7 0.000 0.000 0.009 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/hyperelasticity.py:430(C33_planestress)
# 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)
# 3 0.001 0.000 0.008 0.003 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:207(fint_terms)
# 2 0.000 0.000 0.008 0.004 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:186(getd2Cd2F1)
# 4 0.007 0.002 0.007 0.002 {built-in method scipy.sparse._sparsetools.coo_tocsr}
# 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)
# 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)
# 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)
# 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)
# 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)
# 5 0.001 0.000 0.005 0.001 /home/chevaugeon-n/Documents/RechercheUN/develop/PY_XFEM/pyxfem/solidmechanics/membrane.py:172(getC)
# 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)
# The initial membrane is assumed to be a curved surface
# 0 0.00440687981439875
# 1 0.0007090203268567412
# 2 7.995410736685281e-07
# inc initial z pull 1.1 converged
# 32.22576608864403 40.067667319006084
# 0 3.4310598963475734
# 1 0.2935674403395363
# 2 0.01827820455432498
# 3 0.0002816147637323649
# 4 9.036886193951496e-08
# inc 001 a 5.000e-02, 3.426e-01, p 8.403e-01 converged in 4 iterations
# 32.225766088644036 42.96510692768051
# 0 0.36539286272311966
# 1 0.029378488740486756
# 2 0.00046313210328144475
# 3 5.788205975849805e-07
# 4 2.2703713941695908e-13
# inc 002 a 5.000e-02, 3.347e-01, p 1.071e+00 converged in 4 iterations
# 32.22576608864404 34.1846510483814
# 0 0.31496387608195736
# 1 2.1415147031750412
# newton failed after 1
# inc 002 p 1.071e+00 not converged after 1 iterations
# 32.22576608864404 34.1846510483814
# 0 0.19127906246138682
# 1 2.221491558000836
# newton failed after 1
# inc 002 p 1.071e+00 not converged after 1 iterations
# 32.22576608864404 34.1846510483814
# 0 0.11732610750436605
# 1 41.41408859196228
# newton failed after 1
# inc 002 p 1.071e+00 not converged after 1 iterations
# 32.22576608864404 34.1846510483814
# 0 0.07214008537315857
# 1 0.883440686809872
# newton failed after 1
# inc 002 p 1.071e+00 not converged after 1 iterations
# 32.22576608864404 34.1846510483814
# 0 0.044518527641326175
# 1 0.08708125448554915
# newton failed after 1
# inc 002 p 1.071e+00 not converged after 1 iterations
# 32.22576608864404 34.1846510483814
# 0 0.0275742096485467
# 1 0.017107224332357146
# 2 8.750321381946716e-05
# 3 3.631399800953893e-09
# 4 1.6966004494157222e-14
# inc 003 a 1.563e-03, 1.138e-01, p 1.091e+00 converged in 4 iterations
# 32.22576608864403 2584.6899989463413
# 0 4.660028981796859
# 1 1.0883871245650745
# 2 0.24051205442431087
# 3 0.07826732787845271
# 4 0.015068006598989521
# 5 0.0016632268498018973
# 6 1.168588040514311e-05
# 7 3.1948361642820114e-10
# inc 004 a 1.563e-03, 2.691e-02, p 1.090e+00 converged in 7 iterations
# 32.22576608864404 921.3652756624772
# 0 0.9095356315019443
# 1 0.13097402901504196
# 2 0.006575055412094633
# 3 8.462309488545935e-05
# 4 4.5710282161511506e-07
# 5 1.958417172507177e-12