Coverage for pyxfem/solidmechanics/exact_solutions.py: 0%
288 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4This file contains exact solutions to some classic pde problems
6Use full to test convergence
7Created on Wed Feb 28 09:49:32 2024
8@author: nchevaug
9"""
11from abc import ABC, abstractmethod
12import unittest
13import numpy as np
14#import numpy.typing as npt
15import gmsh
17import pyxfem.core.gmsh_mesh_generator as mg
18import pyxfem.postpro.gmsh_post_pro_interface as gp
19import pyxfem.core.mesh_simplex_one_level as sm
21from pyxfem.solidmechanics.elasticity import MaterialHypothesis,PLANESTRAIN, PLANESTRESS, LinearIsotropeElasticMaterial
22from pyxfem.solidmechanics.lefm import WilliamsSerie
24class LaplacianGeneral():
25 '''Base class for the definition of Laplacian type problem and its exact solution '''
27 def __init__(self, lsfuncs=None, classifyer=None, filt=0, conductivity=None, source=None,
28 BNDDefinition=None, DirichletBCs=None, NeumannBCs=None,
29 exactSolution=None, exactGradient=None) -> None:
30 super().__init__()
31 self.source = source
32 if lsfuncs is None:
33 self.lsfuncs = []
34 else:
35 self.lsfuncs = lsfuncs
36 self.classifyer = classifyer
37 self.BNDDefinition = {} if BNDDefinition is None else BNDDefinition
38 self.exactSolution = exactSolution
39 self.exactGradient = exactGradient
40 self.DirichletBCs = {} if DirichletBCs is None else DirichletBCs
41 self.NeumannBCs = {} if NeumannBCs is None else NeumannBCs
42 self.conductivity = conductivity
43 self.filter = filt
45 def get_source_term(self):
46 """ return the souce term associated to the problem """
47 return self.source
49 def get_ls(self):
50 """ get The level-set function used to define the problem
51 (when multiple conductiviies are present )"""
52 return self.lsfuncs
54 def get_exact_sol(self):
55 """ return the exact solution of the problem """
56 return self.exactSolution
58 def get_exact_grad(self):
59 """ return the gradient exact solution of the problem """
60 return self.exactGradient
62 def get_dirichlet_bc_val(self):
63 """ return a dictionary of boundary name to dirichlet value """
64 return self.DirichletBCs
66 def get_neumann_bc_val(self):
67 """ return a dictionary of boundary name to neumann value """
68 return self.NeumannBCs
70 def get_conductivity_val(self):
71 """ return a dictionary of conductivity values """
72 return self.conductivity
75class LaplacianCircularInclusionQuad(LaplacianGeneral):
76 """ Circular inclusion, with quadratic solution
78 u equal zero at re, source term 's = 1.' is constant, 2 conductivity : k1 inside, k2 outside
79 """
81 def __init__(self, k1, k2, ri=0.4, re=np.sqrt(2)) -> None:
82 super().__init__()
83 self.ri = ri
84 self.re = re
85 self.lsfuncs = [lambda xyz: np.sqrt(
86 xyz[..., 0]**2 + xyz[..., 1]**2) - ri]
87 self.BNDDefinition = {'up': lambda x: x[1] == +1.0,
88 'down': lambda x: x[1] == -1.0,
89 'right': lambda x: x[0] == +1.0,
90 'left': lambda x: x[0] == -1.0, }
91 self.NeumannBCs = {}
92 self.conductivity = {'lsneg': k1, 'lspos': k2}
94 def r(xy):
95 return np.hypot(xy[..., 0], xy[..., 1])
97 def exact_sol(xy):
98 rr = r(xy)
99 return np.where(rr <= ri,
100 (k1*(re**2 - ri**2) + k2*(-rr**2 + ri**2))/(4*k1*k2),
101 (-rr**2 + re**2)/(4*k2))
103 def exact_grad(xy):
104 return np.where((r(xy) < ri)[..., np.newaxis], -xy/(2*k1), -xy/(2*k2))
105 self.source = 1.
106 self.exactSolution = exact_sol
107 self.exactGradient = exact_grad
108 self.DirichletBCs = {'down': self.exactSolution, 'up': self.exactSolution,
109 'left': self.exactSolution, 'right': self.exactSolution}
112class LaplacianCircularInclusionQuintic(LaplacianGeneral):
113 """ Circular inclusion with order 5 displacement source (r-re)**3 """
115 def __init__(self, k1, k2, ri=0.4, re=np.sqrt(2)) -> None:
116 super().__init__()
117 self.ri = ri
118 self.re = re
119 lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
120 # self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
121 conductivity = {'lsneg': k1, 'lspos': k2}
123 def r(xy):
124 return np.hypot(xy[..., 0], xy[..., 1])
126 def qroverr(rr):
127 return -0.5*re**3 + re**2*rr - 3./4*re*rr**2 + 1./5*rr**3
129 def primitiveqr(rr):
130 return (-48*rr**5 + 225*rr**4*re-400*rr**3*re**2+300*rr**2*re**3)/1200
132 def exact_sol(xy):
133 rr = r(xy)
134 uOut = -primitiveqr(rr)/k2 + primitiveqr(re)/k2
135 uIn = -primitiveqr(ri)/k2 + primitiveqr(re)/k2 - \
136 primitiveqr(rr)/k1 + primitiveqr(ri)/k1
137 return np.where(rr <= ri, uIn, uOut)
139 def exact_grad(xy):
140 rr = r(xy)
141 qrroverr = qroverr(rr)[..., np.newaxis]
142 return np.where(rr[..., np.newaxis] < ri, xy*qrroverr/k1, xy*qrroverr/k2)
144 def source(xy):
145 return -(r(xy) - re)**3
147 dirichlet_bc = {'down': exact_sol, 'up': exact_sol,
148 'left': exact_sol, 'right': exact_sol}
149 super().__init__(lsfuncs=lsfuncs, filt=-1, conductivity=conductivity, source=source,
150 DirichletBCs=dirichlet_bc, NeumannBCs={},
151 exactSolution=exact_sol, exactGradient=exact_grad)
153class LaplacianSquareHole(LaplacianGeneral):
154 """ Manufactured solution : square with a square hole
156 +--------------+
157 | |
158 | +--+ |
159 | | | |
160 | +--+ |
161 | |
162 +--------------+
163 """
165 def __init__(self, a=0.4) -> None:
166 super().__init__()
167 self.a = a
168 lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a,
169 lambda xyz: -xyz[1]-a, lambda xyz: xyz[1]-a]
170 # self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
171 BNDDefinition = {'up': lambda x: x[1] == +1.0,
172 'down': lambda x: x[1] == -1.0,
173 'right': lambda x: x[0] == +1.0,
174 'left': lambda x: x[0] == -1.0, }
175 NeumannBCs = {}
176 conductivity = {'lsneg': 1., 'lspos': 1.}
178 def exact_grad(x, y):
179 gex = np.zeros((2, x.shape[0], x.shape[1]))
180 gex[0] = -np.pi/self.a * \
181 np.sin(np.pi * x / self.a) * np.cos(np.pi * y / self.a)
182 gex[1] = -np.pi/self.a * \
183 np.cos(np.pi * x / self.a) * np.sin(np.pi * y / self.a)
184 return gex
186 def exact_sol(x, y):
187 return np.cos(np.pi * x / self.a) * np.cos(np.pi * y / self.a)
189 def source(x, y):
190 return 2.*np.pi**2/self.a**2*np.cos(np.pi*x/self.a)*np.cos(np.pi*y/self.a)
192 DirichletBCs = {'down': exact_sol, 'up': exact_sol,
193 'left': exact_sol, 'right': exact_sol}
194 super().__init__(lsfuncs=lsfuncs, filt=-1, conductivity=conductivity, source=source,
195 BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs,
196 NeumannBCs=NeumannBCs,
197 exactSolution=exact_sol, exactGradient=exact_grad)
199class LinearElasticityExacSol(ABC):
200 ''' Base class for exact solution of linear elasticity problems '''
201 @abstractmethod
202 def displacement(self, xy):
203 ''' Return the displacement at each point in the xy input array
204 input array of float such as [..., :] contain the x and y coordinate of the point.
205 return an array of the sameshape containing the displacement at the input point
206 '''
207 @abstractmethod
208 def strain(self, xy):
209 ''' Return the displacement at each point in the xy input array
210 input array of float such as [..., :] contain the x and y coordinate of the point.
211 return an array of shape xy.shape[:-1] + (2,2) containing the 2d strain.
212 '''
213 @abstractmethod
214 def stress(self, xy):
215 ''' Return the displacement at each point in the xy input array
216 input array of float such as [..., :] contain the x and y coordinate of the point.
217 return an array of shape xy.shape[:-1] + (2,2) containing the 2d stress.
218 '''
220class WilliamsExpansion():
221 ''' Exact field from williams expansion (crack tip is at (0,0), x axis aligned with the crack)
222 ^ y
223 |
224 |
225 ______o---> x
226 crack
227 '''
228 def __init__(self, mat, A=None, B =None): # pylint: disable=C0103
229 self.will = WilliamsSerie(mat)
230 if A is None :
231 A = dict()
232 if B is None :
233 B = dict()
234 self.a = A
235 self.b = B
236 def displacement(self, xy):
237 ''' Compute the displacement at xy '''
238 return self.will.u_xy(xy, self.a, self.b)
239 def stress(self, xy):
240 ''' Compute the stress at xy '''
241 return self.will.stress_xy(xy, self.a, self.b)
243 def strain(self, xy):
244 ''' Compute the stress at xy '''
245 raise ValueError
249# def cartesian2polar(x,y):
250# ''' convert x, y cartesian coordinate to r theta polar coordinate '''
251# return np.hypot(x,y), np.arctan2(y,x)
252# class williamsExpansion():
253# ''' This is the space for exact solution of a straight crack in 2d elasticity
254# crack tip is supposed to be at (0, 0)
255# crack lip are along y = 0, x<0 half axe.
256# '''
257# # A0, B0, B1 -> rigid body motion (x, y, theta respectivelly)
258# # ux = A0 * (kappa+1)/2/mu
259# # uy = B0 * (kappa+1)/2/mu
260# # uxr = (kappa+1)/2/mu B2(-rsin(theta))
261# # uyr = (kappa+1)/2/mu B2( rcos(theta))
262# # A1 = KI/sqrt(2*pi)
263# # B1 = -KII/sqrt(2*pi)
264# # A2 = 4T
265# def __init__(self, elastic_coeffs):
266# self.elastic_coeffs = elastic_coeffs
268# @staticmethod
269# def set_elastic_coeffs(E:float = 1.,
270# nu:float = 0.3,
271# se:float = 1.,
272# hyp2d:MaterialHypothesis=PLANESTRAIN):
273# ''' return a dictionnary for elastic coefficient usefull for LEFM
274# input :
275# E : Young Modulus
276# nu : Poisson coefficient
277# se : Elastic limite
278# hyp2d : PLANSTRESS or PLANSTRAIN
279# '''
280# lamb = E*nu/((1.+nu)*(1.-2.*nu))
281# mu = 0.5*E/(1+nu)
282# if hyp2d == PLANESTRAIN :
283# kappa = (3-4*nu) # kolossov constant en deformation plane (echantillon massif)
284# elif hyp2d == PLANESTRESS :
285# kappa = (3-nu)/(1+nu)
286# else :
287# print("Error : unknown 2d hypothesis :", hyp2d)
288# raise ValueError
289# elastic_coefs = {'E':E, 'nu':nu, 'se':se,'lambda': lamb, 'mu': mu, 'hyp2d':'plane_strain', 'kappa':kappa}
290# return elastic_coefs
292# @staticmethod
293# def fI_11(n:int, theta:float|np.ndarray) :
294# ''' fI_ij are fII_ij are component ij of tensor for mode I and II (symetric and antisymetric ) respectivelly '''
295# return n/2*( (2+(-1.)**n+n/2)*np.cos((n/2-1)*theta) - (n/2-1)*np.cos((n/2-3)*theta ) )
296# @staticmethod
297# def fI_22(n, theta):
298# return n/2*( (2-(-1.)**n-n/2)*np.cos((n/2-1)*theta) + (n/2-1)*np.cos((n/2-3)*theta ) )
299# @staticmethod
300# def fI_12(n, theta):
301# return n/2*( (n/2-1.)*np.sin((n/2-3)*theta) - (n/2+ (-1.)**n )*np.sin((n/2-1)*theta ) )
302# @staticmethod
303# def fII_11(n, theta) :
304# return n/2*( (-2+(-1)**n-n/2)*np.sin((n/2-1)*theta) + (n/2-1)*np.sin((n/2-3)*theta ) )
305# @staticmethod
306# def fII_22(n, theta) :
307# return n/2*( (-2-(-1)**n+n/2)*np.sin((n/2-1)*theta) - (n/2-1)*np.sin((n/2-3)*theta ) )
308# @staticmethod
309# def fII_12(n, theta) :
310# return n/2*( (n/2-1)*np.cos((n/2-3)*theta) - (n/2 - (-1)**n)*np.cos((n/2-1)*theta ) )
311# @staticmethod
312# def gI_1(kappa, n, theta):
313# return ( kappa + (-1.)**n + n/2)*np.cos(n/2*theta) - n/2*np.cos((n/2-2)*theta)
314# @staticmethod
315# def gI_2(kappa, n, theta):
316# return ( kappa - (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
317# @staticmethod
318# def gII_1(kappa, n, theta):
319# return (-kappa + (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
320# @staticmethod
321# def gII_2(kappa, n, theta):
322# return ( kappa + (-1)**n - n/2)*np.cos(n/2*theta) + n/2*np.cos((n/2-2)*theta)
323# def stress2d_rtheta( self, r=1., theta=np.linspace(-np.pi, np.pi, 10),
324# A=None,
325# B=None):
326# if A is None:
327# A=dict()
328# if B is None:
329# B=dict()
330# stress = np.zeros((3,*theta.shape))
331# for (n, An) in A.items() :
332# stress[0] += r**(n/2-1)*An*self.fI_11(n, theta)
333# stress[1] += r**(n/2-1)*An*self.fI_22(n, theta)
334# stress[2] += r**(n/2-1)*An*self.fI_12(n, theta)
335# for (n, Bn) in B.items() :
336# stress[0] += r**(n/2-1)*Bn*self.fII_11(n, theta)
337# stress[1] += r**(n/2-1)*Bn*self.fII_22(n, theta)
338# stress[2] += r**(n/2-1)*Bn*self.fII_12(n, theta)
339# return stress
341# def u_rtheta(self, r =1., theta=np.linspace(-np.pi, np.pi, 10),
342# A=None, B=None):
343# if A is None:
344# A=dict()
345# if B is None:
346# B=dict()
347# if np.isscalar(theta) :
348# u = np.zeros(2)
349# else :
350# u = np.zeros((2, *theta.shape))
351# kappa = self.elastic_coeffs['kappa']
352# mu = self.elastic_coeffs['mu']
353# for (n, An) in A.items() :
354# u[0] += r**(n/2)/2/mu*An*self.gI_1(kappa, n, theta)
355# u[1] += r**(n/2)/2/mu*An*self.gI_2(kappa, n, theta)
356# for (n, Bn) in B.items() :
357# u[0] += r**(n/2)/2/mu*Bn*self.gII_1(kappa, n, theta)
358# u[1] += r**(n/2)/2/mu*Bn*self.gII_2(kappa, n, theta)
359# return u
361# def stress_xy(self, x, y, A=None, B=None, compute_sigma33 = False):
362# if A is None:
363# A=dict()
364# if B is None:
365# B=dict()
366# r, theta = cartesian2polar(x, y)
367# sigma2d = self.stress2d_rtheta(r, theta, A, B)
368# if not compute_sigma33 :
369# return sigma2d
370# if self.elastic_coeffs['hyp2d'] == 'plane_strain':
371# lamb = self.elastic_coeffs['lambda']
372# mu = self.elastic_coeffs['mu']
373# sigma33 = lamb/(2*lamb + 2*mu) *(sigma2d[0] + sigma2d[1])
374# else :
375# sigma33 = np.zeros(x.shape)
376# return sigma2d, sigma33
378# def u_xy(self, x, y, A=None, B=None):
379# if A is None:
380# A=dict()
381# if B is None:S
382# B=dict()
383# r, theta = cartesian2polar(x, y)
384# return self.u_rtheta(r, theta, A, B)
388class LinearElasticity2dCircularHoleInfiniteSpace(LinearElasticityExacSol):
389 """ linear Elasticity exact solution : hole in a plate """
391 def __init__(self, mat:LinearIsotropeElasticMaterial, r=0.5, p=1.):
392 nu = mat.nu
393 if mat.hypothesis == PLANESTRAIN:
394 k2 = (1.-2.*nu) / (2.*(1.-nu))
395 elif mat.hypothesis == PLANESTRESS:
396 k2 = (1.-nu) / 2.
397 else:
398 raise ValueError("hyp must be either planestrain or planestress")
399 self.young = mat.E
400 self.nu = nu
401 self.mu = mat.mu
402 self.k2 = k2
403 self.hyp = mat.hypothesis
404 self.r = r
405 self.p = p
407 def f(self, z):
408 """ complex potential f evaluation """
409 assert z.dtype == np.complex128
410 p, r = (self.p, self.r)
411 return p*r/4.*(z/r + 2*r/z)
413 def df(self, z):
414 """ fist derivative of the complex potential f evaluation """
415 assert z.dtype == np.complex128
416 p, r = (self.p, self.r)
417 return p*r/4.*(1/r - 2*r/z/z)
419 def ddf(self, z):
420 """ fist derivative of the complex potential f evaluation """
421 assert z.dtype == np.complex128
422 p, r = (self.p, self.r)
423 return p*r*r/z/z/z
425 def dg(self, z):
426 """ fist derivative of the complex potential g evaluation """
427 assert z.dtype == np.complex128
428 p, r = (self.p, self.r)
429 return -p*r/2.*(z/r + r/z - np.power(r, 3)*np.power(z, -3))
431 def ddg(self, z):
432 """ second derivative of the complex potential g evaluation """
433 assert z.dtype == np.complex128
434 p, r = (self.p, self.r)
435 return -p*r/2.*(1/r - r/z/z + 3 * pow(r, 3)*pow(z, -4))
437 def displacement(self, xy):
438 """ return the displacement value at each position xy """
439 if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
440 shape = xy.shape
441 z = xy.view(dtype=np.complex128).squeeze()
442 else:
443 z = xy
444 shape = xy.shape + (2,)
445 k2 = self.k2
446 mu = self.mu
447 mu2upiv = (1+k2)/(1-k2) * self.f(z) - z * \
448 np.conj(self.df(z)) - np.conj(self.dg(z))
449 return mu2upiv.view(dtype=np.float64).reshape(shape)/2./mu
451 def stress(self, xy, dim=2):
452 """ return the stress field
454 if dim = 2 return array shape is (..,2,2)
455 if dim = 3 return array shape is (..,3,3)
456 """
457 if dim not in [2, 3]:
458 raise ValueError("dim not in [2,3]")
459 if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
460 shape = xy.shape
461 z = xy.view(dtype=np.complex128).reshape(shape[:-1])
462 else:
463 z = xy
464 shape = xy.shape + (2,)
465 # sigmaxx + sigmayy
466 sxxpsyy = 4*np.real(self.df(z))
467 # sigmayy - sigmaxx + 2isigmaxy
468 syymsxxp2isxy = 2.*(np.conj(z)*self.ddf(z) + self.ddg(z))
469 sxx = (sxxpsyy - np.real(syymsxxp2isxy))*0.5
470 syy = (sxxpsyy + np.real(syymsxxp2isxy))*0.5
471 sxy = 0.5*np.imag(syymsxxp2isxy)
472 stress_ = np.zeros(shape[:-1]+(dim, dim))
473 stress_[..., 0, 0] = sxx
474 stress_[..., 0, 1] = sxy
475 stress_[..., 1, 0] = sxy
476 stress_[..., 1, 1] = syy
477 if (self.hyp == PLANESTRAIN and dim == 3):
478 stress_[..., 2, 2] = self.nu*(sxx+syy)
479 return stress_
481 def strain(self, xy, dim=2):
482 """ return the displacement field
484 if dim = 2 return array shape is (..,2,2)
485 if dim = 3 return array shape is (..,3,3)
486 """
487 if dim not in [2, 3]:
488 raise ValueError("dim not in [2,3]")
489 if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
490 z = xy.view(dtype=np.complex128).squeeze()
491 else:
492 z = xy
493 stress_ = self.stress(z, dim)
494 E, nu = (self.young, self.nu)
495 strain_ = np.zeros(stress_.shape)
496 if self.hyp == PLANESTRAIN:
497 strain_[..., 0, 0] = (1.+nu)/E * ((1-nu) *
498 stress_[..., 0, 0] - nu*stress_[..., 1, 1])
499 strain_[..., 1, 1] = (1.+nu)/E * ((1-nu) *
500 stress_[..., 1, 1] - nu*stress_[..., 0, 0])
501 strain_[..., 0, 1] = (1.+nu)/E * stress_[..., 0, 1]
502 elif self.hyp == PLANESTRESS:
503 strain_[..., 0, 0] = (
504 1./E) * (stress_[..., 0, 0] - nu*stress_[..., 1, 1])
505 strain_[..., 1, 1] = (
506 1./E) * (stress_[..., 1, 1] - nu*stress_[..., 0, 0])
507 strain_[..., 0, 1] = (1./E)*(1.-nu) * stress_[..., 0, 1]
508 if dim == 3:
509 strain_[..., 2, 2] = -nu / \
510 (1.-nu) * (strain_[..., 0, 0]+strain_[..., 1, 1])
511 else:
512 raise ValueError("mat.hypothesis should be PLANESTRAIN or PLANESTRESS")
513 strain_[..., 1, 0] = strain_[..., 0, 1]
514 return strain_
518def laplacian_inclusion(exactsolbuilder, popUpGmsh=False):
519 ''' plot exact solution for laplacien with 2 conductivities on a mesh '''
520 gmsh.initialize()
521 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
522 gmsh.option.setNumber('General.Terminal', 0)
523 # Conductivities :
524 k1 = 1.
525 k2 = 2.*k1
526 # length of the square
527 L = 2.
528 # radius of the disk separting conductivities.
529 ri = 0.4
530 re = np.sqrt(2.)
531 # mesh density
532 lc = 0.05
533 # generate gmshmesh
534 modelName, groupeName2PhysDimTag = mg.disk_in_ref_square(
535 length=L, radius=ri, relh=lc)
536 # convert gmshmesh to smesh
537 mesh = sm.gmshmodel2smesh(modelName, groupeName2PhysDimTag.keys())
538 trisxy = mesh.t2v_coord()
539 triscogxy = mesh.t2v_cog()
540 exact = exactsolbuilder(k1=k1, k2=k2, ri=ri, re=re)
541 if callable(exact.source):
542 source = exact.source(trisxy)
543 else:
544 source = exact.source*np.ones(trisxy.shape[:-1])
545 temp = exact.exactSolution(trisxy)
546 grad = exact.exactGradient(triscogxy)
547 # plotting the fields in gmsh
548 gp.listPlotFieldTri(trisxy, source, viewname="source")
549 gp.listPlotFieldTri(trisxy, temp, viewname="temperature")
550 gp.listPlotFieldTri(trisxy, grad, P0=True, viewname="gradT")
551 # run gmsh ide
552 if popUpGmsh:
553 gmsh.fltk.run()
554 gmsh.finalize()
555 return temp
558def elasticity_hole_in_plate(popUpGmsh=False):
559 """ test exactSolutions.linearElasticity2dCircularHoleInfiniteSpace and plot it"""
560 gmsh.initialize()
561 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
562 gmsh.option.setNumber('General.Terminal', 0)
563 L = 2.
564 r0 = 0.2
565 mat = LinearIsotropeElasticMaterial(E = 1., nu = 0.3)
566 lc = 0.05
567 P = 1.
568 # generate gmshmesh
569 modelName, groupeName2PhysDimTag = mg.hole_in_square(
570 length=L, radius=r0, relh=lc)
571 # convert gmshmesh to smesh
572 mesh = sm.gmshmodel2smesh(modelName, groupeName2PhysDimTag.keys())
573 trisxy = mesh.t2v_coord()
574 # getting all the boundary edges and their normal
575 eids, eflips = (np.array([], dtype=np.int64), np.array([], dtype=np.int64))
576 for name in mesh.get_edgegroupe_names():
577 eids_n, eflips_n = mesh.get_edgegroupe(name)
578 eids = np.append(eids, eids_n)
579 eflips = np.append(eflips, eflips_n)
580 normals = mesh.enormal(eids, eflips)
581 boundaryedgexy = mesh.e2v_coord(eids)
582 # constructing the exact solution
583 exact = LinearElasticity2dCircularHoleInfiniteSpace(mat, r=r0, p=P)
584 disp = exact.displacement(trisxy)
585 stress = exact.stress(trisxy)
586 strain = exact.strain(trisxy)
587 boundarystress = exact.stress(boundaryedgexy)
588 boundarytraction = np.einsum(
589 "ijkl,ijl->ijk", boundarystress, normals[:, np.newaxis, :])
590 # plotting the fields in gmsh
591 gp.listPlotFieldTri(trisxy, disp, viewname="displacement",
592 VectorType="Displacement", DisplacementFactor=0.1)
593 gp.listPlotFieldTri(trisxy, stress, viewname="stress")
594 gp.listPlotFieldTri(trisxy, strain, viewname="strain")
595 gp.listPlotFieldLine(boundaryedgexy, boundarytraction,
596 viewname='boundary tractions')
598 # run gmsh ide
599 if popUpGmsh:
600 gmsh.fltk.run()
601 gmsh.finalize()
602 return stress
604# OTHER EXACT SOLUTIONS COURTESIE OF G. LEGRAIN. Will updated ans corrected when needed
605# class laplacianUniaxialDirichlet(laplacianGeneral):
606# def __init__(self, source=0) -> None:
607# super().__init__()
608# # self.diri = {'up': 1., 'down': 0.}
609# self.source = source
610# self.lsfuncs = [lambda xyz: xyz[1]]
611# self.classifyer = classifyers.classifyerAllIn
612# self.BNDDefinition = { 'up': lambda x: x[1] == 1.0,
613# 'down': lambda x: x[1] == -1.0,}
614# self.DirichletBCs = {'up': 1., 'down': 0.}
615# self.NeumannBCs = {}
616# self.conductivity = {}
617# self.exactSolution = lambda x,y: 0.5*(y+1.)
618# self.exactGradient = lambda x,y: 0.5*np.array([x*0.,y*0.+1])
620# class laplacianUniaxialGREG(laplacianGeneral):
621# def __init__(self, source=0) -> None:
622# super().__init__()
623# # self.diri = {'up': 1., 'down': 0.}
624# self.source = source
625# self.lsfuncs = [lambda xyz: xyz[1]]
626# self.classifyer = classifyers.classifyerAllIn
627# self.BNDDefinition = { 'up': lambda x: x[1] == 1.0,
628# 'down': lambda x: x[1] == -1.0,}
629# self.DirichletBCs = {'down': 0.}
630# self.NeumannBCs = {'up': 1.}
631# self.conductivity = {}
632# self.exactSolution = lambda x,y: 0.5*(y+1.)
633# self.exactGradient = lambda x,y: 0.5*np.array([x*0.,y*0.+1])
635# # s s
636# # |-------------|--------------|
637# # 0 k1 k2 1
638# class laplacianUniaxialSourceMultimat(laplacianGeneral):
639# def __init__(self, source=1, k1=1., k2=1.) -> None:
640# super().__init__()
641# lsfuncs = [lambda xyz: xyz[1]]
642# classifyer = classifyers.classifyerAllIn
643# BNDDefinition = { 'up': lambda x: x[1] == 1.0,
644# 'down': lambda x: x[1] == -1.0,}
645# DirichletBCs = {'down': 0., 'up': 1.}
646# NeumannBCs = {}
647# conductivity = {'lsneg': k1, 'lspos': k2}
648# s = source
649# def exactSol(x,y):
650# uex = np.zeros(x.shape)
651# uex[y < 0.] = -(y[y < 0.] + 1)*(-2*k1*k2 + k1*s*y[y < 0.] - 2*k1*s + k2*s*y[y < 0.])
652# /(2*k1*(k1 + k2))
653# uex[y > 0.] = -(-2*k1*k2*y[y > 0.] + k1*s*y[y > 0.]**2 - k1*s*y[y > 0.] -
654# 2*k2**2 + k2*s*y[y > 0.]**2 + k2*s*y[y > 0.] - 2*k2*s)/(2*k2*(k1 + k2))
655# return uex
656# def exactGrad(x,y):
657# gex = np.zeros((2,x.shape[0], x.shape[1]))
658# gex[1][y<0.] = (k1*k2-k1*s*y[y<0.] + k1*s/2 - k2*s*y[y<0.] - k2*s/2)/(k1*(k1+k2))
659# gex[1][y>0.] = (k1*k2-k1*s*y[y>0.] + k1*s/2 - k2*s*y[y>0.] - k2*s/2)/(k2*(k1+k2))
660# return gex
661# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer, filter=-1,
662# conductivity=conductivity, source=source,
663# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
664# exactSolution=exactSol, exactGradient=exactGrad)
666# # Circular hole in cylindrical coordinates
667# class laplacianCircularHoleDirichletRi(laplacianGeneral):
668# def __init__(self, ri=0.4, re=2, k=6) -> None:
669# super().__init__()
670# NOT_WORKING
671# self.ri = ri
672# self.re = re
673# self.k = k
674# lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
675# # self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
676# classifyer = classifyers.classifyerAllIn
677# BNDDefinition = {'all': lambda x:
678# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
679# np.logical_or(x[0] == 1.0, x[0] == -1.0)),
680# 'hole': lambda x: np.abs(np.sqrt(x[0]**2 + x[1]**2) - ri) < 1.e-5}
681# NeumannBCs = {}
682# conductivity = {'lsneg': 1., 'lspos': 1.}
683# def r(x,y):
684# return np.hypot(x,y)
685# def exactGrad(x,y):
686# gex = np.zeros((2,x.shape[0], x.shape[1]))
687# rr = r(x,y)
688# dr = self.re - self.ri
689# gex[0] = np.pi*self.k*x*np.cos(np.pi*self.k*(self.ri - rr)/(dr))/((dr)*rr)
690# gex[1] = np.pi*self.k*y*np.cos(np.pi*self.k*(self.ri - rr)/(dr))/((dr)*rr)
691# return gex
692# def exactSol(x,y):
693# return np.sin(np.pi*self.k*(r(x,y) - self.ri)/(self.re - self.ri))
694# source = lambda x,y: -np.pi*self.k*(-np.pi*self.k*r(x,y)*np.sin(np.pi*self.k
695# *(r(x,y) - self.ri)/(self.re - self.ri))
696# + (self.re - self.ri)*np.cos(np.pi*self.k*(r(x,y) - self.ri)
697# /(self.re - self.ri)))/(r(x,y)*(self.re - self.ri)**2)
698# DirichletBCs = {'all': exactSol, 'hole': exactSol}
699# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer,
700# filter=-1, conductivity=conductivity, source=source,
701# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
702# exactSolution=exactSol, exactGradient=exactGrad)
703# class laplacianCircularHoleNeumannRi(laplacianGeneral):
704# def __init__(self, ri=0.4, re=4, k=3) -> None:
705# super().__init__()
706# self.ri = ri
707# self.re = re
708# self.k = k
709# lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
710# classifyer = classifyers.classifyerAllIn
711# BNDDefinition = {'all': lambda x:
712# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
713# np.logical_or(x[0] == 1.0, x[0] == -1.0)),}
714# NeumannBCs = {}
715# conductivity = {'lsneg': 1., 'lspos': 1.}
716# def r(x,y):
717# return np.sqrt(x**2 + y**2)
718# def exactGrad(x,y):
719# gex = np.zeros((2,x.shape[0], x.shape[1]))
720# tmp = (ri - np.sqrt(x**2 + y**2))/(2*(re - ri)))/(2*(re - ri)*np.sqrt(x**2 + y**2))
721# gex[0] = np.pi*x*(2*k + 1)*np.sin(np.pi*(2*k + 1)*tmp
722# gex[1] = np.pi*y*(2*k + 1)*np.sin(np.pi*(2*k + 1)*tmp
723# return gex
724# def exactSol(x,y):
725# # return 0.*x
726# return np.cos(np.pi*(k + 1/2)*(r(x,y) - ri)/(re - ri))
727# source = lambda x,y: -np.pi*(2*k + 1)*(-np.pi*r(x,y)*(2*k + 1)
728# *np.cos(np.pi*(2*k + 1)*(r(x,y) - ri)/(2*(re - ri)))
729# + 2*(-re + ri)*np.sin(np.pi*(2*k + 1)*(r(x,y) - ri)
730# /(2*(re - ri))))/(4*r(x,y)*(re - ri)**2)
731# DirichletBCs = {'all': exactSol}
732# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer,
733# filter=-1, conductivity=conductivity, source=source,
734# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
735# exactSolution=exactSol, exactGradient=exactGrad)
738# class laplacianCorner(laplacianGeneral):
739# def __init__(self) -> None:
740# lsfuncs = [lambda xyz: xyz[0], lambda xyz: xyz[1]]
741# BNDDefinition = {'all': lambda x:
742# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
743# np.logical_or(x[0] == 1.0, x[0] == -1.0)),}
744# conductivity = {'lsneg': 1., 'lspos': 1.}
745# # Coordinate change to define alpha angle: theta = 3pi/4 + alpha
746# # Used so that tha discontinuity in the arctan2 function is not inside the domain
747# P=np.array([[np.cos(5*np.pi/4), np.cos(-np.pi/4)],
748# [np.sin(5*np.pi/4), np.sin(-np.pi/4)]]).T
749# def r(x,y):
750# return np.sqrt(x**2 + y**2)
751# def computeTheta(x,y):
752# initShape = x.shape
753# xy = np.vstack([x.flatten(), y.flatten()])
754# xyt = P@xy
755# xt = xyt[0].reshape(initShape)
756# yt = xyt[1].reshape(initShape)
757# alpha = np.arctan2(yt,xt)
758# theta = 3.*np.pi/4. + alpha
759# # aaaa
760# return theta
761# def exactSol(x,y):
762# uex = np.zeros(x.shape)
763# rr = r(x,y)
764# theta = computeTheta(x,y)
765# uex = 2./5. * rr**(2./3.)*np.cos(2./3. * theta)
766# return uex
767# def exactGrad(x,y):
768# gex = np.zeros((2,x.shape[0], x.shape[1]))
769# theta = computeTheta(x,y)
770# rr = r(x,y)
771# gex_r = 4./15. * rr**(-1./3.) * np.cos(2./3. * theta)
772# gex_t = -4./15. * rr**(-1./3.) * np.sin(2./3. * theta)
773# # Project back in cartesian coordinates
774# gex_x = gex_r * np.cos(theta + np.pi/2.) - gex_t * np.sin(theta + np.pi/2.)
775# gex_y = gex_r * np.sin(theta + np.pi/2.) + gex_t * np.cos(theta + np.pi/2.)
776# gex[0] = gex_x
777# gex[1] = gex_y
778# return gex
779# DirichletBCs = {'all': exactSol}
780# super().__init__(lsfuncs=lsfuncs, classifyer=classifyers.classifyerAnyIn,
781# filter=1, conductivity=conductivity, source=0.,
782# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs={},
783# exactSolution=exactSol, exactGradient=exactGrad)
784#
785# class laplacianStraightPatchTest(laplacianGeneral):
786# def __init__(self, type='lin') -> None:
787# # super().__init__()
788# assert(type == 'lin' or type == 'quad')
789# lsfuncs = [lambda xyz: xyz[1]]
790# self.type = type
791# classifyer = classifyers.classifyerAllIn
792# BNDDefinition = { 'down': lambda x: x[1] == -1.0,
793# 'right': lambda x: x[0] == 1.0,
794# 'left': lambda x: x[0] == -1.0,}
795# NeumannBCs = {}
796# conductivity = {'lsneg': 1., 'lspos': 1.}
797# def exactSol(x,y):
798# if self.type == 'lin':
799# uex = 0.5 * (1-x)
800# else:
801# uex = x**2 - 0.5 * (1+x)
802# return uex
803# def exactGrad(x,y):
804# gex = np.zeros((2,x.shape[0], x.shape[1]))
805# if self.type == 'lin':
806# gex[0] -= 0.5
807# else:
808# gex[0] = 2. * x - 0.5
809# return gex
810# if self.type == 'lin':
811# source = lambda x,y: 0.
812# else:
813# source = lambda x,y: -2.
814# DirichletBCs = {'left': exactSol, 'right': exactSol}
815# super().__init__(lsfuncs=lsfuncs,
816# classifyer=classifyers.classifyerAnyIn, filter=1,
817# conductivity=conductivity, source=source,
818# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs={},
819# exactSolution=exactSol, exactGradient=exactGrad)
823class testAnaly(unittest.TestCase):
824 """ testing Analytical solutions """
826 def testElasticityHoleInPlate(self, popUpGmsh=False):
827 """ - Compute max exact stress on mesh nodes and compare to ref value """
828 stress = elasticity_hole_in_plate(popUpGmsh)
829 max_stressxx = np.max(stress[:, :, 0, 0])
830 self.assertTrue((max_stressxx - 3) < 1.e-9)
832 def testLaplacianCircularInclusionQuad(self, popUpGmsh=False):
833 """ - Compute max Temperature on mesh nodes and compare to ref value """
834 temp = laplacian_inclusion(
835 LaplacianCircularInclusionQuad, popUpGmsh=popUpGmsh)
836 self.assertTrue(np.abs(np.max(temp) - 0.2698578163107752) < 1.e-6)
838 def testLaplacianCircularInclusionQuintic(self, popUpGmsh=False):
839 """ - Compute max Temperature on mesh nodes and compare to ref value """
840 temp = laplacian_inclusion(
841 LaplacianCircularInclusionQuintic, popUpGmsh=popUpGmsh)
842 self.assertTrue(np.abs(np.max(temp) - 0.21952206273616834) < 1.e-6)
845if __name__ == '__main__':
846 laplacian_inclusion(LaplacianCircularInclusionQuad, popUpGmsh=True)
847 laplacian_inclusion(LaplacianCircularInclusionQuintic, popUpGmsh=True)
848 elasticity_hole_in_plate(popUpGmsh=True)