#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 27 11:48:03 2023
@author: nchevaug
Tools for Linear fracture mechanics.
"""
import functools
from numba.cuda import popc
import numpy as np
from scipy.integrate import quad
from pyxfem.core.mesh_simplex_one_level import SMesh
from pyxfem.solidmechanics.elasticity import planeStrainC2222CSR, to_tensor2d, to_voigt2d, \
PLANESTRAIN, LinearIsotropeElasticMaterial, ElasticField, P1FemElasticfield
from pyxfem.solidmechanics.mechanics_solver import ElasticitySolver
[docs]
def cartesian2polar(x,y):
''' return r theta from x and y'''
return np.hypot(x,y), np.arctan2(y,x)
[docs]
def polar2cartesian(r,theta):
''' return x,y from r, theta '''
return r*np.cos(theta), r*np.sin(theta)
[docs]
class crack_coordinate_system():
def __init__(self, crack_tip_xy, crack_angle):
self.crack_tip_xy = crack_tip_xy
self.crack_angle = crack_angle
[docs]
def global_to_crack_xy(self, xy):
c,s = (np.cos(self.crack_angle), np.sin(self.crack_angle))
mat = np.array([[c,-s],[s,c]])
return np.vecmat(xy-self.crack_tip_xy[np.newaxis], mat)
[docs]
def global_to_crack_rtheta(self, xy):
xyc = self.global_to_crack_xy(xy)
r, theta = cartesian2polar(xyc[...,0], xyc[...,1])
return r, theta
[docs]
def crack_xy_to_global(self, xyc):
c,s = (np.cos(self.crack_angle), np.sin(self.crack_angle))
mat = np.array([[c,s],[-s,c]])
return np.vecmat(xyc, mat) + self.crack_tip_xy[np.newaxis]
[docs]
def crack_rtheta_to_global(self, r, theta):
x, y = polar2cartesian(r = r, theta = theta)
return self.crack_xy_to_global(np.stack((x,y), axis =-1))
[docs]
class WilliamsSerie():
''' Tools for crack tip williams serie expansion of the elastic field
Mode 1 williams serie with T and B stress
A0, B0, B1 -> rigid body motion (x, y, theta respectivelly)
ux = A0 * (kappa+1)/2/mu
uy = B0 * (kappa+1)/2/mu
uxr = (kappa+1)/2/mu B2(-rsin(theta))
uyr = (kappa+1)/2/mu B2( rcos(theta))
A1 = KI/sqrt(2*pi)
B1 = -KII/sqrt(2*pi)
A2 = 4T
'''
def __init__(self, material_param : LinearIsotropeElasticMaterial):
self.material_param = material_param
[docs]
@staticmethod
def fI_tenseur2(n, theta, notation = "tensor2D"): # pylint: disable=C0103
''' return the part that depends only on theta and n (the term of the serie)
of the stress tensor for mode 1 (openning)
t2.shape = (3, theta.shape)
using voight notation :
t2[..., 0] : s11
t2[;.., 1] : s22
t2[..., 2] : s12
'''
cn1 = np.cos((n/2-1)*theta)
sn1 = np.sin((n/2-1)*theta)
cn3 = np.cos((n/2-3)*theta)
sn3 = np.sin((n/2-3)*theta)
s11 = n/2*((2+(-1.)**n+n/2)*cn1 - (n/2-1)*cn3)
s22 = n/2*((2-(-1.)**n-n/2)*cn1 + (n/2-1)*cn3)
s12 = n/2*((n/2-1.)*sn3 - (n/2+ (-1.)**n)*sn1)
if notation == "tensor2D" :
return to_tensor2d(s11, s22, s12)
if notation == "voigt2D" :
return to_voigt2d(s11, s22, s12)
raise ValueError
[docs]
@staticmethod
def fII_tenseur2(n, theta, notation = "tensor2D"): # pylint: disable=C0103
''' return the part that depends only on theta and n (the term of the serie)
of the stress tensor for mode II (shear)
t2.shape = (3, theta.shape)
using voight notation :
t2[0] : s11
t2[1] : s22
t2[2] : s12
'''
cn1 = np.cos((n/2-1)*theta)
sn1 = np.sin((n/2-1)*theta)
cn3 = np.cos((n/2-3)*theta)
sn3 = np.sin((n/2-3)*theta)
s11 = n/2*((-2+(-1.)**n-n/2)*sn1 + (n/2-1)*sn3 )
s22 = n/2*((-2-(-1.)**n+n/2)*sn1 - (n/2-1)*sn3 )
s12 = n/2*(( n/2-1.)*cn3 - (n/2- (-1.)**n)*cn1 )
if notation == "voigt2D" :
return to_voigt2d(s11, s22, s12)
if notation == "tensor2D" :
return to_tensor2d(s11, s22, s12)
raise ValueError
[docs]
@staticmethod
def gI_vector2(kappa, n, theta): # pylint: disable=C0103
''' Return the part of the displacement field that depend on theta for symetric mode n '''
u1 = ( kappa + (-1.)**n + n/2)*np.cos(n/2*theta) - n/2*np.cos((n/2-2)*theta)
u2 = ( kappa - (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
return np.stack((u1, u2),-1)
[docs]
@staticmethod
def gII_vector2(kappa, n, theta): # pylint: disable=C0103
''' Return the part of the displacement field that depend on theta for
anti symetric mode n '''
u1 = (-kappa + (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
u2 = ( kappa + (-1)**n - n/2)*np.cos(n/2*theta) + n/2*np.cos((n/2-2)*theta)
return np.stack((u1, u2),-1)
[docs]
@staticmethod
def u_rtheta_mode(mode, n, r, theta, kappa, mu):
''' return the displacement corresponding to mode 1 or 2, r^(n/2)'''
r = np.atleast_1d(r)
fr = r**(n/2)/2/mu
if mode == 1:
ftheta = WilliamsSerie.gI_vector2(kappa, n, theta)
elif mode ==2 :
ftheta = WilliamsSerie.gII_vector2(kappa, n, theta)
else :
raise ValueError
return fr[..., np.newaxis]*ftheta
[docs]
@staticmethod
def u_xy_mode(mode, n, x, y, kappa, mu):
''' return the displacement corresponding to mode 1 or 2, r^(n/2)'''
r, theta = cartesian2polar(x, y)
return WilliamsSerie.u_rtheta_mode(mode, n, r, theta, kappa, mu)
[docs]
@staticmethod
def stress_rtheta_mode(mode, n, r, theta):
''' return the stress corresponding to mode 1 or 2, r^(n/2)'''
r = np.atleast_1d(r)
fr = r**(n/2-1)
if mode == 1:
ftheta = WilliamsSerie.fI_tenseur2(n, theta)
elif mode ==2 :
ftheta = WilliamsSerie.fII_tenseur2(n, theta)
else :
raise ValueError
return fr[..., np.newaxis, np.newaxis]*ftheta
[docs]
@staticmethod
def stress_xy_mode(mode, n, x, y):
''' return the stress corresponding to mode 1 or 2, r^(n/2)'''
r, theta = cartesian2polar(x, y)
return WilliamsSerie.stress_rtheta_mode(mode, n, r, theta)
[docs]
def stress2d_rtheta( self, r=1., theta=np.linspace(-np.pi, np.pi, 10), A=None, B=None): # pylint: disable=C0103
''' return the stress field for williams expansion defined by coefficient A and B
sigm(r, theta) = An*fn(r)*fIn(theta) + Bn*fn(r)*fIIn(theta)
(crack axi - coordinate)
'''
if A is None :
A = dict()
if B is None :
B = dict()
r = np.atleast_1d(r)
stress = np.zeros((*theta.shape, 2, 2))
for (n, a) in A.items() :
stress += a*WilliamsSerie.stress_rtheta_mode(1, n, r, theta)
for (n, b) in B.items() :
stress += b*WilliamsSerie.stress_rtheta_mode(2, n, r, theta)
return stress
[docs]
def u_rtheta(self, r =1., theta=np.linspace(-np.pi, np.pi, 10), A=None, B=None): # pylint: disable=C0103
'''
return the displacement field for williams expansion defined by coefficient A and B
u(r, theta) = An*fn(r)*gIn(theta) + Bn*fn(r)*gIIn(theta)
(crack axi - coordinate)
'''
if A is None :
A = dict()
if B is None :
B = dict()
r = np.atleast_1d(r)
theta = np.atleast_1d(theta)
if np.isscalar(theta) :
u = np.zeros(2)
else :
u = np.zeros(( *theta.shape, 2))
kappa = self.material_param.kappa
mu = self.material_param.mu
for (n, an) in A.items():
u+= an*WilliamsSerie.u_rtheta_mode(1, n, r, theta, kappa, mu)
for (n, bn) in B.items() :
u+= bn*WilliamsSerie.u_rtheta_mode(2, n, r, theta, kappa, mu)
return u
[docs]
def stress_xy(self, xy, A=None, B=None, compute_sigma33 = False): # pylint: disable=C0103
''' compute the stress from xy position and values A and B of the williams coefficients '''
if A is None :
A = dict()
if B is None :
B = dict()
r, theta = cartesian2polar(xy[...,0 ], xy[..., 1])
sigma2d = self.stress2d_rtheta(r, theta, A, B)
if not compute_sigma33 :
return sigma2d
if self.material_param.hypothesis == PLANESTRAIN :
lamb = self.material_param.lamb
mu = self.material_param.mu
sigma33 = lamb/(2*lamb + 2*mu) *(sigma2d[...,0,0] + sigma2d[...,1,1])
else :
sigma33 = np.zeros(xy.shape[:-1])
return sigma2d, sigma33
[docs]
def u_xy(self, xy, A=None, B=None): # pylint: disable=C0103
''' compute the displacement from xy position and values A and B
of the williams coefficients '''
r, theta = cartesian2polar(xy[...,0 ], xy[..., 1])
return self.u_rtheta(r, theta, A, B)
[docs]
class WilliamsField(ElasticField):
''' Exact field from williams expansion (crack tip is at (0,0), x axis aligned with the crack)
^ y
|
|
______o---> x
crack
'''
def __init__(self, mat, coeffs_modeI=None, coeffs_modeII =None): # pylint: disable=C0103
self.will = WilliamsSerie(mat)
if coeffs_modeI is None :
coeffs_modeI = dict()
if coeffs_modeII is None :
coeffs_modeII = dict()
self.a = coeffs_modeI
self.b = coeffs_modeII
[docs]
def displacement(self, xy):
''' Compute the displacement at xy '''
return self.will.u_xy(xy, self.a, self.b)
[docs]
def stress(self, xy):
''' Compute the stress at xy '''
return self.will.stress_xy(xy, self.a, self.b)
[docs]
def strain(self, xy):
raise NotImplementedError
[docs]
def buckner_integrand(u1, t1, u2, t2):
''' return pointwise value of the brickner interaction integrand between two field
on a line u1, t1 and u2, t2 are respectively the displacement and the tension in a direction
normale to the line at the point.
'''
return np.vecdot(t1,u2) - np.vecdot(t2,u1)
[docs]
def construct_face_domain(mesh:SMesh, cracktipxy, r0):
''' Return the list of edges on a contour of a domain of all the 2D triangle that has
at least on node in the dis of radius r0 centered on cracktipxy.
'''
r = np.linalg.norm( mesh.vcoord() - cracktipxy[np.newaxis, :], axis = 1)
vertices = np.nonzero( r < r0)[0]
v2t = mesh.v2t()
f = np.array([], dtype = np.uint64)
for v in vertices :
fs = np.uint64(v2t.index[v])
fe = np.uint64(v2t.index[v+1])
nf =np.array(v2t.data[fs:fe], dtype = np.uint64)
f = np.append(f, nf)
f = np.unique(f)
return f
[docs]
def get_boundary_edges(mesh:SMesh, faces):
''' return the list of boundary edge of the list of faces. Considered Boundary if inside
the domain mesh (edge connected to 2 triangles) and on the boundary of the list of face
(only one of the two faces of each edge is in the list faces)
also, return orient that tell if each edge as the same orientation as the corresponding edge
of it's connected face (usefull to compute the external normal)
'''
istriin = np.zeros(mesh.ntris(), np.bool)
istriin[faces] = True
e2t = mesh.e2t()
nbfaces = e2t.ndata_perentry()
internaledges = np.nonzero(nbfaces==2)[0]
e_face0 = e2t.data[ np.array(e2t.index[internaledges], dtype = np.uint64)]
e_face1 = e2t.data[ np.array(e2t.index[internaledges]+1, dtype = np.uint64)]
e_face0_in = istriin[e_face0]
e_face1_in = istriin[e_face1]
edges = np.nonzero( np.logical_xor(e_face0_in, e_face1_in))[0]
orient = np.zeros(edges.shape, dtype = np.bool)
boundary_edges_id = np.array(internaledges[edges], dtype =np.int64)
nbe = boundary_edges_id.shape[0]
orient = np.zeros(nbe, dtype = np.bool)
e_f0 = e2t.data[np.array(e2t.index[boundary_edges_id], dtype = np.uint64)]
e_f1 = e2t.data[np.array(e2t.index[boundary_edges_id]+1, dtype = np.uint64)]
e_f_in = np.where( np.isin(e_f0, faces), e_f0, e_f1)
t2e, t2eflip = mesh.t2e_flip()
ief = np.nonzero(t2e[e_f_in] == boundary_edges_id[:, np.newaxis])
orient = t2eflip[e_f_in][ief]
return boundary_edges_id, orient
[docs]
def construct_domain_boundary(mesh:SMesh, cracktipxy, r0):
''' return list of faces boundary edges and orientation of the boundary edge of a domain
containing all the elements that have at least one node in the disk of radius r0
'''
domain_faces = construct_face_domain(mesh, cracktipxy, r0)
boundary_edges, boundary_edges_flip = get_boundary_edges(mesh, domain_faces)
return domain_faces, boundary_edges, boundary_edges_flip
#edge_stress = field.stress(edge_cog)
# def KMPotFourrier(imode, id, kappa):
# harm1 = np.array([id, (4-id), (-id)])
# if imode==1 :
# amp1 = np.array([kappa, -id/2, id/2 + -1**id])
# elif imode==2 :
# amp1 = np.array([-1j*kappa, -1j*(id/2), 1j*(id/2-(-1)**id)])
# else :
# print("Error in KMPotFourrier. Mode should be 1 or 2")
# raise ValueError
# return harm1, amp1
# def crack_strain(xs, ys, mode = 1, K = 1., T = 0., B= 0.,
# elasticcoef = elastic_coeffs(), thetap = 0.):
# ''' thetap : angle de la fissure avec les axes '''
# kappa = elasticcoef['kappa']
# mu = elasticcoef['mu']
# scal=2*mu*np.sqrt(2*np.pi)
# rs, ts = cartesian2polar(xs,ys)
# km= np.array([K,T/(4/np.sqrt(2*np.pi)),B])/scal
# epsxx=0
# epsyy=0
# epsxy=0
# for hi in [1,2,3] :
# harm1,amp1=KMPotFourrier(mode,hi,kappa)
# ft=np.zeros(rs.shape)
# fr=np.zeros(rs.shape)
# for amp, harm in zip(amp1, harm1):
# fr=fr+amp*np.exp(harm*ts*1j/2)
# ft=ft+harm*1j/2*amp*np.exp(harm*ts*1j/2)
# ft=ft*(rs**(hi/2-1))
# fr=hi/2*fr*(rs**(hi/2-1))
# dfndth = np.imag(ft)
# dfndr = np.imag(fr)
# dftdth = np.real(ft)
# dftdr = np.real(fr)
# drdt = np.cos(ts)
# drdn = np.sin(ts)
# dthdt = -np.sin(ts)
# dthdn = np.cos(ts)
# dfndt=dfndr*drdt+dfndth*dthdt
# dfndn=dfndr*drdn+dfndth*dthdn
# dftdt=dftdr*drdt+dftdth*dthdt
# dftdn=dftdr*drdn+dftdth*dthdn
# c0 = np.cos(thetap)
# s0 = np.sin(thetap)
# dftdx=dftdt*c0-dftdn*s0
# dftdy=dftdt*s0+dftdn*c0
# dfndx=dfndt*c0-dfndn*s0
# dfndy=dfndt*s0+dfndn*c0
# dfxdx=(dftdx*c0-dfndx*s0)
# dfxdy=(dftdy*c0-dfndy*s0)
# dfydx=(dftdx*s0+dfndx*c0)
# dfydy=(dftdy*s0+dfndy*c0)
# epsxx=epsxx+km[hi-1]*dfxdx
# epsyy=epsyy+km[hi-1]*dfydy
# epsxy=epsxy+km[hi-1]*(dfxdy+dfydx)/2
# return np.array([epsxx, epsyy, epsxy])
# def crack_VM(xs, ys, mode = 1, K = 1., T = 0., B= 0.,
# elasticcoef = elastic_coeffs(), thetap = 0.):
# eps = crack_strain(xs, ys, mode, K, T, B, elasticcoef, thetap)
# sigma2d, sigmazz = stress_from_strain(eps, elasticcoef['lambda'], elasticcoef['mu'])
# return VonMises(sigma2d, sigmazz )
[docs]
def plot_modes(mesh, mat,
modeI:None|dict = None, # pylint: disable=C0103
modeII:None|dict = None): # pylint: disable=C0103
''' plot each of the williams modes from n = 0 to 4, symetric then unsymetric mode '''
default_modes = {0:1., 1:1., 2:1., 3:1., 4:1.}
# pylint: disable=C0415
import pyxfem.core.gmsh_mesh_generator as mg
import pyxfem.core.mesh_simplex_one_level as sm
import pyxfem.core.function_spaces as fs
import pyxfem.postpro.gmsh_post_pro_interface as gp
from pyxfem.core.quadrature import T3_nodes, T3_gauss
import pyxfem.core.sparse_tools as st
fem_space = fs.FEMSpaceP1("DisplacementSpace", mesh, vecdim=2)
will = WilliamsSerie(mat)
xy = mesh.getTris2VerticesCoord()
xycog = mesh.getTris2VerticesCog()
NatVertices = fem_space.operator_dof2val_tri(T3_nodes().uv) # pylint: disable=C0103
B = fem_space.operator_dof2grad_tri(T3_gauss(0).uv) # pylint: disable=C0103
J = fem_space.get_tris_mapping().J # pylint: disable=C0103
nf = len(J)
C = planeStrainC2222CSR(mat = mat) # pylint: disable=C0103
C = C.todense().flatten() # pylint: disable=C0103
wD = np.einsum("i, ij -> ij", T3_gauss(0).w[0]*J, C[np.newaxis, :]).reshape((-1, 4, 4)) # pylint: disable=C0103
wD = st.regular_block_diag_to_csr(wD) # pylint: disable=C0103
D = np.broadcast_to(C, (nf, 16)).reshape((-1, 4, 4)) # pylint: disable=C0103
D = st.regular_block_diag_to_csr(D) # pylint: disable=C0103
stress_op = D@B
kappa = mat.kappa
mu = mat.mu
def plots(mode, n, scal = 1.):
def u_xy(mode, n):
return lambda xy : scal*WilliamsSerie.u_xy_mode(mode, n, xy[...,0], xy[..., 1], kappa, mu)
U = fem_space.interpolate(u_xy(mode, n)).flatten() # pylint: disable=C0103
u = (NatVertices@U).reshape((-1, 3, 2))
stress_calc = (stress_op@U).reshape((-1, 1, 2, 2))
if mode ==1 :
stress_exac = will.stress_xy(xycog, A = {n:scal}, compute_sigma33 = False)
else :
stress_exac = will.stress_xy(xycog, B = {n:scal}, compute_sigma33 = False)
stress_exac = WilliamsSerie.stress_xy_mode(mode, n, xycog[..., 0], xycog[..., 1])
#print()
gp.listPlotFieldTri(xy, u, P0=False,
viewname= f"displacement mode {mode:d}, n = {n:d}", VectorType = "Displacement")
gp.listPlotFieldTri(xy, stress_calc, P0=True,
viewname=f'Stress calc mode {mode:d}, n = {n:d}', Range=[0.,10.])
gp.listPlotFieldTri(xy, stress_exac, P0=True,
viewname=f'Stress exac mode {mode:d}, n = {n:d}', Range=[0.,10.])
if modeI is None:
modeI = default_modes
if modeII is None:
modeII = default_modes
for mode, coefs in zip([1,2], [modeI, modeII]):
for n, scal in coefs.items():
print(n, scal)
plots(mode, n, scal)
[docs]
def test_crack_coordinate_system(gmsh_popup = False):
xy = np.array([[0.,0.],
[1.,0.],
[0.,1.]])
crack_axis = crack_coordinate_system(np.array([0.,0.]), np.pi/6.)
import gmsh
gmsh.initialize()
from pyxfem.core.gmsh_mesh_generator import crack_in_square_contour, crack_in_square, square_with_inner_seg # pylint: disable=W0611
from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh
from pyxfem.postpro.gmsh_post_pro_interface import listPlotFieldTri
gmsh.option.setNumber("PostProcessing.AnimationCycle", 1)
gmsh.option.setNumber('General.Terminal', 0)
#model_name, ____ = crack_in_square_contour(relh = 0.1, r = 0.5)
model_name, ____ = square_with_inner_seg(length=2., relh1=0.1, relh2=0.1, xc0=(-0.5, 0.), xc1=(0.5,0.))
if gmsh_popup:
gmsh.fltk.run()
mesh = gmshmodel2smesh(model_name, edge_group_names= ['Crack'])
#mesh = gmshModel2sMesh(model_name, edgeGroupeNames= ['contour','bottom', 'top', 'left_top', 'left_bottom', 'right'])
#dirichlet_vertices = mesh.get_vertices_on_edges_groups(['bottom', 'top', 'left_top', 'left_bottom', 'right'])
xy = mesh.t2v_coord()
data = np.ones(xy.shape[:-1])
listPlotFieldTri(xy, data)
listPlotFieldTri(crack_axis.global_to_crack_xy(xy), data)
if gmsh_popup:
gmsh.fltk.run()
gmsh.finalize()
# = crack_axis.global_to_crack_xy(xy)
# print(xyc)
# crack_axis = crack_coordinate_system(np.array([1.,1.]), np.pi/6.)
# xyc = crack_axis.global_to_crack_xy(xy)
# print(xyc)
# xy = crack_axis.crack_xy_to_global(xyc)
# print(xy)
[docs]
def check_extract_domain_boundary(showdom=False):
''' check domain construction for interactions integrales '''
# pylint: disable=C0415
import pyxfem.core.gmsh_mesh_generator as mg
from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh
import pyxfem.postpro.gmsh_post_pro_interface as gp
model_name, ____ = mg.crack_in_square_contour(relh= 0.025, r=0.45)
mesh = gmshmodel2smesh(model_name, edge_group_names= ["contour"],
tri_group_names = ["Omega_in", "Omega_ext"])
domain_faces, boundary_edges, boundary_edges_flip = \
construct_domain_boundary(mesh, np.array([0.,0.]), 0.601)
xy_triplot = mesh.t2v_coord(fids = domain_faces)
xy_edgeplot = mesh.e2v_coord(eids = boundary_edges)
edge_normal = mesh.enormal(eids=boundary_edges, flips=boundary_edges_flip)
if showdom:
gp.listPlotFieldTri(xy_triplot, np.ones(domain_faces.shape[0]), P0=True, viewname='Domain')
gp.listPlotFieldLine(xy_edgeplot, np.ones(boundary_edges.shape[0]), P0=True,
viewname='Boundary Domain', LineType = '3D cylinder')
gp.listPlotFieldLine(xy_edgeplot, edge_normal, P0=True,
viewname='Boundary Domain', VectorType = 'Arrow')
E = 1. # pylint: disable=C0103
nu = 0.3
se = 1.
mat = LinearIsotropeElasticMaterial(E = E, nu = nu, elastic_limit = se)
A = {1:1., 2:2., 3:3., 4:4., 5:5.} # pylint: disable=C0103
B = {1:-1., 2:-2., 3:-3., 4:-4., 5:-5.} # pylint: disable=C0103
field = WilliamsField(mat, A, B)
print("input :", A, B)
edge_cog = mesh.ecog(eids = boundary_edges)
edge_lenght = mesh.elength(eids = boundary_edges)
edge_normal = mesh.enormal(eids=boundary_edges, flips=boundary_edges_flip)
edge_disp = field.displacement(edge_cog)
edge_stress = field.stress(edge_cog)
A, B = extract_williams_from_edge_list(edge_cog, edge_lenght, edge_normal, edge_disp, edge_stress, # pylint: disable=C0103
mat, modeIterms=[1,2,3,4,5,6], modeIIterms=[1,2,3])
print("res :", A, B)
[docs]
class ElasticPbWilliamsonBC():
''' class to solve a standard lefm crack problem on a mesh conform to the crack
where the crack goes from (-inf, 0.) to (0,0)
input : - A 2d sMesh
a map that give average traction on each edge with neumann BC
the material stiffness operator
'''
# pylint: disable=C0415
def __init__(self, mesh:SMesh, dirichlet_vertices, mat:LinearIsotropeElasticMaterial,
coef_modeI, coef_modeII): # pylint: disable=C0103
elasticity_tensor = planeStrainC2222CSR(mat = mat)
bcfield = WilliamsField(coeffs_modeI=coef_modeI, coeffs_modeII = coef_modeII, mat = mat)
bc_vxy = mesh.xy[dirichlet_vertices]
bc_values = bcfield.displacement(bc_vxy)
bc_node2val = {"idvx":dirichlet_vertices, "vx":bc_values[:,0], \
"idvy":dirichlet_vertices, "vy":bc_values[:,1]}
elasticity_tensor = planeStrainC2222CSR(mat = mat)
self.mesh = mesh
self.esolver = ElasticitySolver(mesh, elasticity_tensor, node_id_to_imposed_value=bc_node2val)
self.vdisp = self.esolver.solve().reshape((-1,2))
self.fstress = self.esolver.getStress(self.vdisp.flatten())
[docs]
def plot(self):
''' plot displacement and Stress field in gmsh format '''
import pyxfem.postpro.gmsh_post_pro_interface as gp
gp.listPlotFieldTri(*self.esolver.postProDisplacement(self.vdisp.flatten()), viewname='displacement fem', VectorType="displacement")
gp.listPlotFieldTri(*self.esolver.postProStress(self.vdisp.flatten()), P0= True, viewname='stress fem', Range =[0.,10.])
[docs]
def get_elatic_field(self):
mesh = self.esolver.mesh
return P1FemElasticfield(mesh.xy, mesh.t2v, self.vdisp, self.fstress)
[docs]
def convergence():
from pyxfem.core.gmsh_mesh_generator import crack_in_square_contour, crack_in_square
from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh
import gmsh
gmsh.option.setNumber("PostProcessing.AnimationCycle", 1)
gmsh.option.setNumber('General.Terminal', 0)
#check_extract_williams_circle()
# check_extract_domain_boundary(showdom=True)
#model_name, ____ = crack_in_square_contour(relh= 0.05, r=0.45)
for h in [0.5**i for i in range(1,8)] :
model_name, ____ = crack_in_square_contour(relh = h,r = 0.5)
mesh = gmshmodel2smesh(model_name, edge_group_names= ['contour','bottom', 'top', 'left_top', 'left_bottom', 'right'])
dirichlet_vertices = mesh.get_vertices_on_edges_groups(['bottom', 'top', 'left_top', 'left_bottom', 'right'])
mat = LinearIsotropeElasticMaterial(E = 1., nu = 0.3, hyp2d=PLANESTRAIN)
#plot_modes(mesh, mat, modeI={1:1.}, modeII={})
fempb = ElasticPbWilliamsonBC(mesh, dirichlet_vertices, mat= mat,
coef_modeI ={1:1., 2:3.}, coef_modeII={})
#fempb.plot()
exact_field = WilliamsField(mat, coeffs_modeI = {1:1., 2:3.})
fem_field = fempb.get_elatic_field()
# cog_test = mesh.getTris2VerticesCog()
# trixy_test = mesh.getTris2VerticesCoord()
# stress_test = fem_field.stress(cog_test)
# disp_fem_test = fem_field.displacement(trixy_test)
# disp_exact_test = exact_field.displacement(trixy_test)
# print(disp_fem_test.shape)
# print(trixy_test.shape)
# from pyxfem.gmsh_post_pro_interface import listPlotFieldTri
# listPlotFieldTri(trixy_test, disp_fem_test, viewname='displacement fem thru interp', VectorType="Arrow")
# listPlotFieldTri(trixy_test, disp_exact_test, viewname='displacement exact thru interp', VectorType="Arrow")
# listPlotFieldTri(trixy_test, stress_test, P0= True, viewname='stress fem thru interp', Range =[0.,10.])
# A, B, disp1f, tension1f = extract_williams_circle(fem_field, mat, modeIterms=[1], modeIIterms=[], r=0.7, # pylint: disable=C0103
# method='trapeze', npt = 500)
# print("res")
# print(A,B)
# A, B, disp1e, tension1e = extract_williams_circle(exact_field, mat, modeIterms=[1], modeIIterms=[], r=0.7, # pylint: disable=C0103
# method='trapeze', npt = 500)
A, B = extract_williams_from_field_contour(fem_field, cracktipxy=[0.,0.], cracktheta =0., mat = mat, modeIterms=[1,2], modeIIterms=[],# pylint: disable=C0103
r=0.7, method='trapeze', npt = int(10*1/h))
#print("res")
#print(A,B)
print(mesh.nv(), h, A[1]-1.)
# A, B = extract_williams_circle(exact_field, mat, modeIterms=[1,2], modeIIterms=[], r=0.7, # pylint: disable=C0103
# method='trapeze', npt = 500)
# print("res")
# print(A,B)
# print(tension1f)
# print(tension1e)
# print(np.linalg.norm(disp1f-disp1e))
# print(np.linalg.norm(tension1f-tension1e))
# print(np.max(np.abs(tension1f-tension1e)))
[docs]
def test_insert_crack_by_node_duplication():
# pylint: disable=import-outside-toplevel
from pyxfem.core.gmsh_mesh_generator import square_with_inner_seg
from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh
from pyxfem.postpro.gmsh_post_pro_interface import listPlotFieldTri
from pyxfem.core.mesh_modifications import duplicate_nodes
import gmsh
gmsh.initialize()
mat = LinearIsotropeElasticMaterial(E = 1., nu = 0.3, hyp2d=PLANESTRAIN)
model_name, geomdic = square_with_inner_seg(length=2, relh1 = 0.05, relh2 = 0.005, xc0=[-0.5,0.], xc1 = [0.5,0.])
mesh = gmshmodel2smesh(model_name, edge_group_names=geomdic.keys())
def compute(mesh):
elasticity_tensor = planeStrainC2222CSR(mat = mat)
diribc = {'Bottom': {'x':0., 'y':0.}, 'Top': {'y':0.1,'x':0.} }
esolver = ElasticitySolver(mesh, elasticity_tensor, physName2Dirichlet=diribc)
u = esolver.solve() # pylint: disable=C0103
listPlotFieldTri(*esolver.postProStress(u), P0=True, viewname='Stress')
listPlotFieldTri(*esolver.postProDisplacement(u), P0=False, viewname='Displacement',
VectorType="Displacement", DisplacementFactor=1.)
compute(mesh)
duplicate_nodes(mesh, mesh.get_edgegroupe('Crack')[0])
compute(mesh)
gmsh.finalize()
if __name__ == '__main__':
test_crack_coordinate_system(gmsh_popup= True)
# pylint: disable=C0415
#import gmsh
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
#test_extract_from_mesh()
test_insert_crack_by_node_duplication()
#convergence()
gmsh.fltk.run()
gmsh.finalize()
# print("checking Orthogonality")
# for m in rangem :
# for n in rangen :
# A1 = {m:1.}
# B1 = {}
# A2 = {n:1.}
# B2 = {}
# fsigma1 = lambda x, y : will.stress_xy(x, y, A =A1, B=B1 )
# fu1 = lambda x, y : will.u_xy(x, y, A=A1 , B=B1)
# fsigma2 = lambda x, y : will.stress_xy(x, y, A = A2 )
# fu2 = lambda x, y : will.u_xy(x, y, A =A2 )
# theta = np.linspace(-np.pi, np.pi, npt)
# x = r*np.cos(theta)
# y = r*np.sin(theta)
# normals= np.vstack ([ np.cos(theta), np.sin(theta)])
# w = 2.*np.pi*r/npt*np.ones(npt)
# I = buckner(x,y, normals, w, fu1, fsigma1, fu2, fsigma2)
# A = -mat['mu']/(mat['kappa']+1)/(n*np.pi)/(-1.)**n*I
# print(f'm={m:d}, n={n:d}, I={I:.3e}, A= {A:.3e}')
# A = {1:1., 2:1.}
# B = {1:1., 2:1.}
# class elastic_field():
# def stress(xy):
# return will.stress_xy(x, y, A = A , B=B )
# def u(xy) :
# return will.u_xy_(x, y, A = A , B=B )