#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This file contains exact solutions to some classic pde problems
Use full to test convergence
Created on Wed Feb 28 09:49:32 2024
@author: nchevaug
"""
from abc import ABC, abstractmethod
import unittest
import numpy as np
#import numpy.typing as npt
import gmsh
import pyxfem.core.gmsh_mesh_generator as mg
import pyxfem.postpro.gmsh_post_pro_interface as gp
import pyxfem.core.mesh_simplex_one_level as sm
from pyxfem.solidmechanics.elasticity import MaterialHypothesis,PLANESTRAIN, PLANESTRESS, LinearIsotropeElasticMaterial
from pyxfem.solidmechanics.lefm import WilliamsSerie
[docs]
class LaplacianGeneral():
'''Base class for the definition of Laplacian type problem and its exact solution '''
def __init__(self, lsfuncs=None, classifyer=None, filt=0, conductivity=None, source=None,
BNDDefinition=None, DirichletBCs=None, NeumannBCs=None,
exactSolution=None, exactGradient=None) -> None:
super().__init__()
self.source = source
if lsfuncs is None:
self.lsfuncs = []
else:
self.lsfuncs = lsfuncs
self.classifyer = classifyer
self.BNDDefinition = {} if BNDDefinition is None else BNDDefinition
self.exactSolution = exactSolution
self.exactGradient = exactGradient
self.DirichletBCs = {} if DirichletBCs is None else DirichletBCs
self.NeumannBCs = {} if NeumannBCs is None else NeumannBCs
self.conductivity = conductivity
self.filter = filt
[docs]
def get_source_term(self):
""" return the souce term associated to the problem """
return self.source
[docs]
def get_ls(self):
""" get The level-set function used to define the problem
(when multiple conductiviies are present )"""
return self.lsfuncs
[docs]
def get_exact_sol(self):
""" return the exact solution of the problem """
return self.exactSolution
[docs]
def get_exact_grad(self):
""" return the gradient exact solution of the problem """
return self.exactGradient
[docs]
def get_dirichlet_bc_val(self):
""" return a dictionary of boundary name to dirichlet value """
return self.DirichletBCs
[docs]
def get_neumann_bc_val(self):
""" return a dictionary of boundary name to neumann value """
return self.NeumannBCs
[docs]
def get_conductivity_val(self):
""" return a dictionary of conductivity values """
return self.conductivity
[docs]
class LaplacianCircularInclusionQuad(LaplacianGeneral):
""" Circular inclusion, with quadratic solution
u equal zero at re, source term 's = 1.' is constant, 2 conductivity : k1 inside, k2 outside
"""
def __init__(self, k1, k2, ri=0.4, re=np.sqrt(2)) -> None:
super().__init__()
self.ri = ri
self.re = re
self.lsfuncs = [lambda xyz: np.sqrt(
xyz[..., 0]**2 + xyz[..., 1]**2) - ri]
self.BNDDefinition = {'up': lambda x: x[1] == +1.0,
'down': lambda x: x[1] == -1.0,
'right': lambda x: x[0] == +1.0,
'left': lambda x: x[0] == -1.0, }
self.NeumannBCs = {}
self.conductivity = {'lsneg': k1, 'lspos': k2}
def r(xy):
return np.hypot(xy[..., 0], xy[..., 1])
def exact_sol(xy):
rr = r(xy)
return np.where(rr <= ri,
(k1*(re**2 - ri**2) + k2*(-rr**2 + ri**2))/(4*k1*k2),
(-rr**2 + re**2)/(4*k2))
def exact_grad(xy):
return np.where((r(xy) < ri)[..., np.newaxis], -xy/(2*k1), -xy/(2*k2))
self.source = 1.
self.exactSolution = exact_sol
self.exactGradient = exact_grad
self.DirichletBCs = {'down': self.exactSolution, 'up': self.exactSolution,
'left': self.exactSolution, 'right': self.exactSolution}
[docs]
class LaplacianCircularInclusionQuintic(LaplacianGeneral):
""" Circular inclusion with order 5 displacement source (r-re)**3 """
def __init__(self, k1, k2, ri=0.4, re=np.sqrt(2)) -> None:
super().__init__()
self.ri = ri
self.re = re
lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
# self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
conductivity = {'lsneg': k1, 'lspos': k2}
def r(xy):
return np.hypot(xy[..., 0], xy[..., 1])
def qroverr(rr):
return -0.5*re**3 + re**2*rr - 3./4*re*rr**2 + 1./5*rr**3
def primitiveqr(rr):
return (-48*rr**5 + 225*rr**4*re-400*rr**3*re**2+300*rr**2*re**3)/1200
def exact_sol(xy):
rr = r(xy)
uOut = -primitiveqr(rr)/k2 + primitiveqr(re)/k2
uIn = -primitiveqr(ri)/k2 + primitiveqr(re)/k2 - \
primitiveqr(rr)/k1 + primitiveqr(ri)/k1
return np.where(rr <= ri, uIn, uOut)
def exact_grad(xy):
rr = r(xy)
qrroverr = qroverr(rr)[..., np.newaxis]
return np.where(rr[..., np.newaxis] < ri, xy*qrroverr/k1, xy*qrroverr/k2)
def source(xy):
return -(r(xy) - re)**3
dirichlet_bc = {'down': exact_sol, 'up': exact_sol,
'left': exact_sol, 'right': exact_sol}
super().__init__(lsfuncs=lsfuncs, filt=-1, conductivity=conductivity, source=source,
DirichletBCs=dirichlet_bc, NeumannBCs={},
exactSolution=exact_sol, exactGradient=exact_grad)
[docs]
class LaplacianSquareHole(LaplacianGeneral):
""" Manufactured solution : square with a square hole
+--------------+
| |
| +--+ |
| | | |
| +--+ |
| |
+--------------+
"""
def __init__(self, a=0.4) -> None:
super().__init__()
self.a = a
lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a,
lambda xyz: -xyz[1]-a, lambda xyz: xyz[1]-a]
# self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
BNDDefinition = {'up': lambda x: x[1] == +1.0,
'down': lambda x: x[1] == -1.0,
'right': lambda x: x[0] == +1.0,
'left': lambda x: x[0] == -1.0, }
NeumannBCs = {}
conductivity = {'lsneg': 1., 'lspos': 1.}
def exact_grad(x, y):
gex = np.zeros((2, x.shape[0], x.shape[1]))
gex[0] = -np.pi/self.a * \
np.sin(np.pi * x / self.a) * np.cos(np.pi * y / self.a)
gex[1] = -np.pi/self.a * \
np.cos(np.pi * x / self.a) * np.sin(np.pi * y / self.a)
return gex
def exact_sol(x, y):
return np.cos(np.pi * x / self.a) * np.cos(np.pi * y / self.a)
def source(x, y):
return 2.*np.pi**2/self.a**2*np.cos(np.pi*x/self.a)*np.cos(np.pi*y/self.a)
DirichletBCs = {'down': exact_sol, 'up': exact_sol,
'left': exact_sol, 'right': exact_sol}
super().__init__(lsfuncs=lsfuncs, filt=-1, conductivity=conductivity, source=source,
BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs,
NeumannBCs=NeumannBCs,
exactSolution=exact_sol, exactGradient=exact_grad)
[docs]
class LinearElasticityExacSol(ABC):
''' Base class for exact solution of linear elasticity problems '''
[docs]
@abstractmethod
def displacement(self, xy):
''' Return the displacement at each point in the xy input array
input array of float such as [..., :] contain the x and y coordinate of the point.
return an array of the sameshape containing the displacement at the input point
'''
[docs]
@abstractmethod
def strain(self, xy):
''' Return the displacement at each point in the xy input array
input array of float such as [..., :] contain the x and y coordinate of the point.
return an array of shape xy.shape[:-1] + (2,2) containing the 2d strain.
'''
[docs]
@abstractmethod
def stress(self, xy):
''' Return the displacement at each point in the xy input array
input array of float such as [..., :] contain the x and y coordinate of the point.
return an array of shape xy.shape[:-1] + (2,2) containing the 2d stress.
'''
[docs]
class WilliamsExpansion():
''' 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, A=None, B =None): # pylint: disable=C0103
self.will = WilliamsSerie(mat)
if A is None :
A = dict()
if B is None :
B = dict()
self.a = A
self.b = B
[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):
''' Compute the stress at xy '''
raise ValueError
# def cartesian2polar(x,y):
# ''' convert x, y cartesian coordinate to r theta polar coordinate '''
# return np.hypot(x,y), np.arctan2(y,x)
# class williamsExpansion():
# ''' This is the space for exact solution of a straight crack in 2d elasticity
# crack tip is supposed to be at (0, 0)
# crack lip are along y = 0, x<0 half axe.
# '''
# # 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, elastic_coeffs):
# self.elastic_coeffs = elastic_coeffs
# @staticmethod
# def set_elastic_coeffs(E:float = 1.,
# nu:float = 0.3,
# se:float = 1.,
# hyp2d:MaterialHypothesis=PLANESTRAIN):
# ''' return a dictionnary for elastic coefficient usefull for LEFM
# input :
# E : Young Modulus
# nu : Poisson coefficient
# se : Elastic limite
# hyp2d : PLANSTRESS or PLANSTRAIN
# '''
# lamb = E*nu/((1.+nu)*(1.-2.*nu))
# mu = 0.5*E/(1+nu)
# if hyp2d == PLANESTRAIN :
# kappa = (3-4*nu) # kolossov constant en deformation plane (echantillon massif)
# elif hyp2d == PLANESTRESS :
# kappa = (3-nu)/(1+nu)
# else :
# print("Error : unknown 2d hypothesis :", hyp2d)
# raise ValueError
# elastic_coefs = {'E':E, 'nu':nu, 'se':se,'lambda': lamb, 'mu': mu, 'hyp2d':'plane_strain', 'kappa':kappa}
# return elastic_coefs
# @staticmethod
# def fI_11(n:int, theta:float|np.ndarray) :
# ''' fI_ij are fII_ij are component ij of tensor for mode I and II (symetric and antisymetric ) respectivelly '''
# return n/2*( (2+(-1.)**n+n/2)*np.cos((n/2-1)*theta) - (n/2-1)*np.cos((n/2-3)*theta ) )
# @staticmethod
# def fI_22(n, theta):
# return n/2*( (2-(-1.)**n-n/2)*np.cos((n/2-1)*theta) + (n/2-1)*np.cos((n/2-3)*theta ) )
# @staticmethod
# def fI_12(n, theta):
# return n/2*( (n/2-1.)*np.sin((n/2-3)*theta) - (n/2+ (-1.)**n )*np.sin((n/2-1)*theta ) )
# @staticmethod
# def fII_11(n, theta) :
# return n/2*( (-2+(-1)**n-n/2)*np.sin((n/2-1)*theta) + (n/2-1)*np.sin((n/2-3)*theta ) )
# @staticmethod
# def fII_22(n, theta) :
# return n/2*( (-2-(-1)**n+n/2)*np.sin((n/2-1)*theta) - (n/2-1)*np.sin((n/2-3)*theta ) )
# @staticmethod
# def fII_12(n, theta) :
# return n/2*( (n/2-1)*np.cos((n/2-3)*theta) - (n/2 - (-1)**n)*np.cos((n/2-1)*theta ) )
# @staticmethod
# def gI_1(kappa, n, theta):
# return ( kappa + (-1.)**n + n/2)*np.cos(n/2*theta) - n/2*np.cos((n/2-2)*theta)
# @staticmethod
# def gI_2(kappa, n, theta):
# return ( kappa - (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
# @staticmethod
# def gII_1(kappa, n, theta):
# return (-kappa + (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta)
# @staticmethod
# def gII_2(kappa, n, theta):
# return ( kappa + (-1)**n - n/2)*np.cos(n/2*theta) + n/2*np.cos((n/2-2)*theta)
# def stress2d_rtheta( self, r=1., theta=np.linspace(-np.pi, np.pi, 10),
# A=None,
# B=None):
# if A is None:
# A=dict()
# if B is None:
# B=dict()
# stress = np.zeros((3,*theta.shape))
# for (n, An) in A.items() :
# stress[0] += r**(n/2-1)*An*self.fI_11(n, theta)
# stress[1] += r**(n/2-1)*An*self.fI_22(n, theta)
# stress[2] += r**(n/2-1)*An*self.fI_12(n, theta)
# for (n, Bn) in B.items() :
# stress[0] += r**(n/2-1)*Bn*self.fII_11(n, theta)
# stress[1] += r**(n/2-1)*Bn*self.fII_22(n, theta)
# stress[2] += r**(n/2-1)*Bn*self.fII_12(n, theta)
# return stress
# def u_rtheta(self, r =1., theta=np.linspace(-np.pi, np.pi, 10),
# A=None, B=None):
# if A is None:
# A=dict()
# if B is None:
# B=dict()
# if np.isscalar(theta) :
# u = np.zeros(2)
# else :
# u = np.zeros((2, *theta.shape))
# kappa = self.elastic_coeffs['kappa']
# mu = self.elastic_coeffs['mu']
# for (n, An) in A.items() :
# u[0] += r**(n/2)/2/mu*An*self.gI_1(kappa, n, theta)
# u[1] += r**(n/2)/2/mu*An*self.gI_2(kappa, n, theta)
# for (n, Bn) in B.items() :
# u[0] += r**(n/2)/2/mu*Bn*self.gII_1(kappa, n, theta)
# u[1] += r**(n/2)/2/mu*Bn*self.gII_2(kappa, n, theta)
# return u
# def stress_xy(self, x, y, A=None, B=None, compute_sigma33 = False):
# if A is None:
# A=dict()
# if B is None:
# B=dict()
# r, theta = cartesian2polar(x, y)
# sigma2d = self.stress2d_rtheta(r, theta, A, B)
# if not compute_sigma33 :
# return sigma2d
# if self.elastic_coeffs['hyp2d'] == 'plane_strain':
# lamb = self.elastic_coeffs['lambda']
# mu = self.elastic_coeffs['mu']
# sigma33 = lamb/(2*lamb + 2*mu) *(sigma2d[0] + sigma2d[1])
# else :
# sigma33 = np.zeros(x.shape)
# return sigma2d, sigma33
# def u_xy(self, x, y, A=None, B=None):
# if A is None:
# A=dict()
# if B is None:S
# B=dict()
# r, theta = cartesian2polar(x, y)
# return self.u_rtheta(r, theta, A, B)
[docs]
class LinearElasticity2dCircularHoleInfiniteSpace(LinearElasticityExacSol):
""" linear Elasticity exact solution : hole in a plate """
def __init__(self, mat:LinearIsotropeElasticMaterial, r=0.5, p=1.):
nu = mat.nu
if mat.hypothesis == PLANESTRAIN:
k2 = (1.-2.*nu) / (2.*(1.-nu))
elif mat.hypothesis == PLANESTRESS:
k2 = (1.-nu) / 2.
else:
raise ValueError("hyp must be either planestrain or planestress")
self.young = mat.E
self.nu = nu
self.mu = mat.mu
self.k2 = k2
self.hyp = mat.hypothesis
self.r = r
self.p = p
[docs]
def f(self, z):
""" complex potential f evaluation """
assert z.dtype == np.complex128
p, r = (self.p, self.r)
return p*r/4.*(z/r + 2*r/z)
[docs]
def df(self, z):
""" fist derivative of the complex potential f evaluation """
assert z.dtype == np.complex128
p, r = (self.p, self.r)
return p*r/4.*(1/r - 2*r/z/z)
[docs]
def ddf(self, z):
""" fist derivative of the complex potential f evaluation """
assert z.dtype == np.complex128
p, r = (self.p, self.r)
return p*r*r/z/z/z
[docs]
def dg(self, z):
""" fist derivative of the complex potential g evaluation """
assert z.dtype == np.complex128
p, r = (self.p, self.r)
return -p*r/2.*(z/r + r/z - np.power(r, 3)*np.power(z, -3))
[docs]
def ddg(self, z):
""" second derivative of the complex potential g evaluation """
assert z.dtype == np.complex128
p, r = (self.p, self.r)
return -p*r/2.*(1/r - r/z/z + 3 * pow(r, 3)*pow(z, -4))
[docs]
def displacement(self, xy):
""" return the displacement value at each position xy """
if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
shape = xy.shape
z = xy.view(dtype=np.complex128).squeeze()
else:
z = xy
shape = xy.shape + (2,)
k2 = self.k2
mu = self.mu
mu2upiv = (1+k2)/(1-k2) * self.f(z) - z * \
np.conj(self.df(z)) - np.conj(self.dg(z))
return mu2upiv.view(dtype=np.float64).reshape(shape)/2./mu
[docs]
def stress(self, xy, dim=2):
""" return the stress field
if dim = 2 return array shape is (..,2,2)
if dim = 3 return array shape is (..,3,3)
"""
if dim not in [2, 3]:
raise ValueError("dim not in [2,3]")
if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
shape = xy.shape
z = xy.view(dtype=np.complex128).reshape(shape[:-1])
else:
z = xy
shape = xy.shape + (2,)
# sigmaxx + sigmayy
sxxpsyy = 4*np.real(self.df(z))
# sigmayy - sigmaxx + 2isigmaxy
syymsxxp2isxy = 2.*(np.conj(z)*self.ddf(z) + self.ddg(z))
sxx = (sxxpsyy - np.real(syymsxxp2isxy))*0.5
syy = (sxxpsyy + np.real(syymsxxp2isxy))*0.5
sxy = 0.5*np.imag(syymsxxp2isxy)
stress_ = np.zeros(shape[:-1]+(dim, dim))
stress_[..., 0, 0] = sxx
stress_[..., 0, 1] = sxy
stress_[..., 1, 0] = sxy
stress_[..., 1, 1] = syy
if (self.hyp == PLANESTRAIN and dim == 3):
stress_[..., 2, 2] = self.nu*(sxx+syy)
return stress_
[docs]
def strain(self, xy, dim=2):
""" return the displacement field
if dim = 2 return array shape is (..,2,2)
if dim = 3 return array shape is (..,3,3)
"""
if dim not in [2, 3]:
raise ValueError("dim not in [2,3]")
if (xy.dtype != np.complex128) and xy.shape[-1] == 2:
z = xy.view(dtype=np.complex128).squeeze()
else:
z = xy
stress_ = self.stress(z, dim)
E, nu = (self.young, self.nu)
strain_ = np.zeros(stress_.shape)
if self.hyp == PLANESTRAIN:
strain_[..., 0, 0] = (1.+nu)/E * ((1-nu) *
stress_[..., 0, 0] - nu*stress_[..., 1, 1])
strain_[..., 1, 1] = (1.+nu)/E * ((1-nu) *
stress_[..., 1, 1] - nu*stress_[..., 0, 0])
strain_[..., 0, 1] = (1.+nu)/E * stress_[..., 0, 1]
elif self.hyp == PLANESTRESS:
strain_[..., 0, 0] = (
1./E) * (stress_[..., 0, 0] - nu*stress_[..., 1, 1])
strain_[..., 1, 1] = (
1./E) * (stress_[..., 1, 1] - nu*stress_[..., 0, 0])
strain_[..., 0, 1] = (1./E)*(1.-nu) * stress_[..., 0, 1]
if dim == 3:
strain_[..., 2, 2] = -nu / \
(1.-nu) * (strain_[..., 0, 0]+strain_[..., 1, 1])
else:
raise ValueError("mat.hypothesis should be PLANESTRAIN or PLANESTRESS")
strain_[..., 1, 0] = strain_[..., 0, 1]
return strain_
[docs]
def laplacian_inclusion(exactsolbuilder, popUpGmsh=False):
''' plot exact solution for laplacien with 2 conductivities on a mesh '''
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
gmsh.option.setNumber('General.Terminal', 0)
# Conductivities :
k1 = 1.
k2 = 2.*k1
# length of the square
L = 2.
# radius of the disk separting conductivities.
ri = 0.4
re = np.sqrt(2.)
# mesh density
lc = 0.05
# generate gmshmesh
modelName, groupeName2PhysDimTag = mg.disk_in_ref_square(
length=L, radius=ri, relh=lc)
# convert gmshmesh to smesh
mesh = sm.gmshmodel2smesh(modelName, groupeName2PhysDimTag.keys())
trisxy = mesh.t2v_coord()
triscogxy = mesh.t2v_cog()
exact = exactsolbuilder(k1=k1, k2=k2, ri=ri, re=re)
if callable(exact.source):
source = exact.source(trisxy)
else:
source = exact.source*np.ones(trisxy.shape[:-1])
temp = exact.exactSolution(trisxy)
grad = exact.exactGradient(triscogxy)
# plotting the fields in gmsh
gp.listPlotFieldTri(trisxy, source, viewname="source")
gp.listPlotFieldTri(trisxy, temp, viewname="temperature")
gp.listPlotFieldTri(trisxy, grad, P0=True, viewname="gradT")
# run gmsh ide
if popUpGmsh:
gmsh.fltk.run()
gmsh.finalize()
return temp
[docs]
def elasticity_hole_in_plate(popUpGmsh=False):
""" test exactSolutions.linearElasticity2dCircularHoleInfiniteSpace and plot it"""
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
gmsh.option.setNumber('General.Terminal', 0)
L = 2.
r0 = 0.2
mat = LinearIsotropeElasticMaterial(E = 1., nu = 0.3)
lc = 0.05
P = 1.
# generate gmshmesh
modelName, groupeName2PhysDimTag = mg.hole_in_square(
length=L, radius=r0, relh=lc)
# convert gmshmesh to smesh
mesh = sm.gmshmodel2smesh(modelName, groupeName2PhysDimTag.keys())
trisxy = mesh.t2v_coord()
# getting all the boundary edges and their normal
eids, eflips = (np.array([], dtype=np.int64), np.array([], dtype=np.int64))
for name in mesh.get_edgegroupe_names():
eids_n, eflips_n = mesh.get_edgegroupe(name)
eids = np.append(eids, eids_n)
eflips = np.append(eflips, eflips_n)
normals = mesh.enormal(eids, eflips)
boundaryedgexy = mesh.e2v_coord(eids)
# constructing the exact solution
exact = LinearElasticity2dCircularHoleInfiniteSpace(mat, r=r0, p=P)
disp = exact.displacement(trisxy)
stress = exact.stress(trisxy)
strain = exact.strain(trisxy)
boundarystress = exact.stress(boundaryedgexy)
boundarytraction = np.einsum(
"ijkl,ijl->ijk", boundarystress, normals[:, np.newaxis, :])
# plotting the fields in gmsh
gp.listPlotFieldTri(trisxy, disp, viewname="displacement",
VectorType="Displacement", DisplacementFactor=0.1)
gp.listPlotFieldTri(trisxy, stress, viewname="stress")
gp.listPlotFieldTri(trisxy, strain, viewname="strain")
gp.listPlotFieldLine(boundaryedgexy, boundarytraction,
viewname='boundary tractions')
# run gmsh ide
if popUpGmsh:
gmsh.fltk.run()
gmsh.finalize()
return stress
# OTHER EXACT SOLUTIONS COURTESIE OF G. LEGRAIN. Will updated ans corrected when needed
# class laplacianUniaxialDirichlet(laplacianGeneral):
# def __init__(self, source=0) -> None:
# super().__init__()
# # self.diri = {'up': 1., 'down': 0.}
# self.source = source
# self.lsfuncs = [lambda xyz: xyz[1]]
# self.classifyer = classifyers.classifyerAllIn
# self.BNDDefinition = { 'up': lambda x: x[1] == 1.0,
# 'down': lambda x: x[1] == -1.0,}
# self.DirichletBCs = {'up': 1., 'down': 0.}
# self.NeumannBCs = {}
# self.conductivity = {}
# self.exactSolution = lambda x,y: 0.5*(y+1.)
# self.exactGradient = lambda x,y: 0.5*np.array([x*0.,y*0.+1])
# class laplacianUniaxialGREG(laplacianGeneral):
# def __init__(self, source=0) -> None:
# super().__init__()
# # self.diri = {'up': 1., 'down': 0.}
# self.source = source
# self.lsfuncs = [lambda xyz: xyz[1]]
# self.classifyer = classifyers.classifyerAllIn
# self.BNDDefinition = { 'up': lambda x: x[1] == 1.0,
# 'down': lambda x: x[1] == -1.0,}
# self.DirichletBCs = {'down': 0.}
# self.NeumannBCs = {'up': 1.}
# self.conductivity = {}
# self.exactSolution = lambda x,y: 0.5*(y+1.)
# self.exactGradient = lambda x,y: 0.5*np.array([x*0.,y*0.+1])
# # s s
# # |-------------|--------------|
# # 0 k1 k2 1
# class laplacianUniaxialSourceMultimat(laplacianGeneral):
# def __init__(self, source=1, k1=1., k2=1.) -> None:
# super().__init__()
# lsfuncs = [lambda xyz: xyz[1]]
# classifyer = classifyers.classifyerAllIn
# BNDDefinition = { 'up': lambda x: x[1] == 1.0,
# 'down': lambda x: x[1] == -1.0,}
# DirichletBCs = {'down': 0., 'up': 1.}
# NeumannBCs = {}
# conductivity = {'lsneg': k1, 'lspos': k2}
# s = source
# def exactSol(x,y):
# uex = np.zeros(x.shape)
# uex[y < 0.] = -(y[y < 0.] + 1)*(-2*k1*k2 + k1*s*y[y < 0.] - 2*k1*s + k2*s*y[y < 0.])
# /(2*k1*(k1 + k2))
# uex[y > 0.] = -(-2*k1*k2*y[y > 0.] + k1*s*y[y > 0.]**2 - k1*s*y[y > 0.] -
# 2*k2**2 + k2*s*y[y > 0.]**2 + k2*s*y[y > 0.] - 2*k2*s)/(2*k2*(k1 + k2))
# return uex
# def exactGrad(x,y):
# gex = np.zeros((2,x.shape[0], x.shape[1]))
# 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))
# 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))
# return gex
# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer, filter=-1,
# conductivity=conductivity, source=source,
# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
# exactSolution=exactSol, exactGradient=exactGrad)
# # Circular hole in cylindrical coordinates
# class laplacianCircularHoleDirichletRi(laplacianGeneral):
# def __init__(self, ri=0.4, re=2, k=6) -> None:
# super().__init__()
# NOT_WORKING
# self.ri = ri
# self.re = re
# self.k = k
# lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
# # self.lsfuncs = [lambda xyz: -xyz[0]-a, lambda xyz: xyz[0]-a]
# classifyer = classifyers.classifyerAllIn
# BNDDefinition = {'all': lambda x:
# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
# np.logical_or(x[0] == 1.0, x[0] == -1.0)),
# 'hole': lambda x: np.abs(np.sqrt(x[0]**2 + x[1]**2) - ri) < 1.e-5}
# NeumannBCs = {}
# conductivity = {'lsneg': 1., 'lspos': 1.}
# def r(x,y):
# return np.hypot(x,y)
# def exactGrad(x,y):
# gex = np.zeros((2,x.shape[0], x.shape[1]))
# rr = r(x,y)
# dr = self.re - self.ri
# gex[0] = np.pi*self.k*x*np.cos(np.pi*self.k*(self.ri - rr)/(dr))/((dr)*rr)
# gex[1] = np.pi*self.k*y*np.cos(np.pi*self.k*(self.ri - rr)/(dr))/((dr)*rr)
# return gex
# def exactSol(x,y):
# return np.sin(np.pi*self.k*(r(x,y) - self.ri)/(self.re - self.ri))
# source = lambda x,y: -np.pi*self.k*(-np.pi*self.k*r(x,y)*np.sin(np.pi*self.k
# *(r(x,y) - self.ri)/(self.re - self.ri))
# + (self.re - self.ri)*np.cos(np.pi*self.k*(r(x,y) - self.ri)
# /(self.re - self.ri)))/(r(x,y)*(self.re - self.ri)**2)
# DirichletBCs = {'all': exactSol, 'hole': exactSol}
# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer,
# filter=-1, conductivity=conductivity, source=source,
# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
# exactSolution=exactSol, exactGradient=exactGrad)
# class laplacianCircularHoleNeumannRi(laplacianGeneral):
# def __init__(self, ri=0.4, re=4, k=3) -> None:
# super().__init__()
# self.ri = ri
# self.re = re
# self.k = k
# lsfuncs = [lambda xyz: np.sqrt(xyz[0]**2 + xyz[1]**2) - ri]
# classifyer = classifyers.classifyerAllIn
# BNDDefinition = {'all': lambda x:
# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
# np.logical_or(x[0] == 1.0, x[0] == -1.0)),}
# NeumannBCs = {}
# conductivity = {'lsneg': 1., 'lspos': 1.}
# def r(x,y):
# return np.sqrt(x**2 + y**2)
# def exactGrad(x,y):
# gex = np.zeros((2,x.shape[0], x.shape[1]))
# tmp = (ri - np.sqrt(x**2 + y**2))/(2*(re - ri)))/(2*(re - ri)*np.sqrt(x**2 + y**2))
# gex[0] = np.pi*x*(2*k + 1)*np.sin(np.pi*(2*k + 1)*tmp
# gex[1] = np.pi*y*(2*k + 1)*np.sin(np.pi*(2*k + 1)*tmp
# return gex
# def exactSol(x,y):
# # return 0.*x
# return np.cos(np.pi*(k + 1/2)*(r(x,y) - ri)/(re - ri))
# source = lambda x,y: -np.pi*(2*k + 1)*(-np.pi*r(x,y)*(2*k + 1)
# *np.cos(np.pi*(2*k + 1)*(r(x,y) - ri)/(2*(re - ri)))
# + 2*(-re + ri)*np.sin(np.pi*(2*k + 1)*(r(x,y) - ri)
# /(2*(re - ri))))/(4*r(x,y)*(re - ri)**2)
# DirichletBCs = {'all': exactSol}
# super().__init__(lsfuncs=lsfuncs, classifyer=classifyer,
# filter=-1, conductivity=conductivity, source=source,
# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs=NeumannBCs,
# exactSolution=exactSol, exactGradient=exactGrad)
# class laplacianCorner(laplacianGeneral):
# def __init__(self) -> None:
# lsfuncs = [lambda xyz: xyz[0], lambda xyz: xyz[1]]
# BNDDefinition = {'all': lambda x:
# np.logical_or(np.logical_or(x[1] == 1.0, x[1] == -1.0),
# np.logical_or(x[0] == 1.0, x[0] == -1.0)),}
# conductivity = {'lsneg': 1., 'lspos': 1.}
# # Coordinate change to define alpha angle: theta = 3pi/4 + alpha
# # Used so that tha discontinuity in the arctan2 function is not inside the domain
# P=np.array([[np.cos(5*np.pi/4), np.cos(-np.pi/4)],
# [np.sin(5*np.pi/4), np.sin(-np.pi/4)]]).T
# def r(x,y):
# return np.sqrt(x**2 + y**2)
# def computeTheta(x,y):
# initShape = x.shape
# xy = np.vstack([x.flatten(), y.flatten()])
# xyt = P@xy
# xt = xyt[0].reshape(initShape)
# yt = xyt[1].reshape(initShape)
# alpha = np.arctan2(yt,xt)
# theta = 3.*np.pi/4. + alpha
# # aaaa
# return theta
# def exactSol(x,y):
# uex = np.zeros(x.shape)
# rr = r(x,y)
# theta = computeTheta(x,y)
# uex = 2./5. * rr**(2./3.)*np.cos(2./3. * theta)
# return uex
# def exactGrad(x,y):
# gex = np.zeros((2,x.shape[0], x.shape[1]))
# theta = computeTheta(x,y)
# rr = r(x,y)
# gex_r = 4./15. * rr**(-1./3.) * np.cos(2./3. * theta)
# gex_t = -4./15. * rr**(-1./3.) * np.sin(2./3. * theta)
# # Project back in cartesian coordinates
# gex_x = gex_r * np.cos(theta + np.pi/2.) - gex_t * np.sin(theta + np.pi/2.)
# gex_y = gex_r * np.sin(theta + np.pi/2.) + gex_t * np.cos(theta + np.pi/2.)
# gex[0] = gex_x
# gex[1] = gex_y
# return gex
# DirichletBCs = {'all': exactSol}
# super().__init__(lsfuncs=lsfuncs, classifyer=classifyers.classifyerAnyIn,
# filter=1, conductivity=conductivity, source=0.,
# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs={},
# exactSolution=exactSol, exactGradient=exactGrad)
#
# class laplacianStraightPatchTest(laplacianGeneral):
# def __init__(self, type='lin') -> None:
# # super().__init__()
# assert(type == 'lin' or type == 'quad')
# lsfuncs = [lambda xyz: xyz[1]]
# self.type = type
# classifyer = classifyers.classifyerAllIn
# BNDDefinition = { 'down': lambda x: x[1] == -1.0,
# 'right': lambda x: x[0] == 1.0,
# 'left': lambda x: x[0] == -1.0,}
# NeumannBCs = {}
# conductivity = {'lsneg': 1., 'lspos': 1.}
# def exactSol(x,y):
# if self.type == 'lin':
# uex = 0.5 * (1-x)
# else:
# uex = x**2 - 0.5 * (1+x)
# return uex
# def exactGrad(x,y):
# gex = np.zeros((2,x.shape[0], x.shape[1]))
# if self.type == 'lin':
# gex[0] -= 0.5
# else:
# gex[0] = 2. * x - 0.5
# return gex
# if self.type == 'lin':
# source = lambda x,y: 0.
# else:
# source = lambda x,y: -2.
# DirichletBCs = {'left': exactSol, 'right': exactSol}
# super().__init__(lsfuncs=lsfuncs,
# classifyer=classifyers.classifyerAnyIn, filter=1,
# conductivity=conductivity, source=source,
# BNDDefinition=BNDDefinition, DirichletBCs=DirichletBCs, NeumannBCs={},
# exactSolution=exactSol, exactGradient=exactGrad)
[docs]
class testAnaly(unittest.TestCase):
""" testing Analytical solutions """
[docs]
def testElasticityHoleInPlate(self, popUpGmsh=False):
""" - Compute max exact stress on mesh nodes and compare to ref value """
stress = elasticity_hole_in_plate(popUpGmsh)
max_stressxx = np.max(stress[:, :, 0, 0])
self.assertTrue((max_stressxx - 3) < 1.e-9)
[docs]
def testLaplacianCircularInclusionQuad(self, popUpGmsh=False):
""" - Compute max Temperature on mesh nodes and compare to ref value """
temp = laplacian_inclusion(
LaplacianCircularInclusionQuad, popUpGmsh=popUpGmsh)
self.assertTrue(np.abs(np.max(temp) - 0.2698578163107752) < 1.e-6)
[docs]
def testLaplacianCircularInclusionQuintic(self, popUpGmsh=False):
""" - Compute max Temperature on mesh nodes and compare to ref value """
temp = laplacian_inclusion(
LaplacianCircularInclusionQuintic, popUpGmsh=popUpGmsh)
self.assertTrue(np.abs(np.max(temp) - 0.21952206273616834) < 1.e-6)
if __name__ == '__main__':
laplacian_inclusion(LaplacianCircularInclusionQuad, popUpGmsh=True)
laplacian_inclusion(LaplacianCircularInclusionQuintic, popUpGmsh=True)
elasticity_hole_in_plate(popUpGmsh=True)