104 integer ,
parameter :: verb = 2
114 integer ,
parameter :: fe_v =
fe_p2_2d 115 integer ,
parameter :: fe_s =
fe_p2_1d 117 character(LEN=1),
parameter :: mesh_idx =
"2" 127 type(quadmesh) :: qdm
130 character(len=100),
parameter :: mesh_file= &
131 & trim(GMSH_DIR)//
"square/square_"//mesh_idx//
".msh" 143 real(RP),
dimension(:),
allocatable :: rhs, u_h
144 real(RP),
dimension(:),
allocatable :: u_h_1, u_h_2
149 real(RP) :: err_L2, err_H1
151 character(LEN=300) :: shell_com
164 write(*,*)
'elasticity_Dirichlet_2D: start' 166 write(*,*)
" EXAMPLE OF A LINEAR ELASTICITY PROBLEM" 168 write(*,*)
" -div(A(x) e(u)) = f on \Omega" 170 write(*,*)
" WITH A DIRICHLET BOUNDARY CONDITION " 172 write(*,*)
" A(x)e(u)) n = g on \partial\Omega" 174 write(*,*)
" AND WITH NON CONSTANT LAME COEFFICIENTS:" 176 write(*,*)
" A(x) \xi = lambda(x) Tr(\xi) Id + 2 mu(x) \si" 178 write(*,*)
" The problem data 'lambda', 'mu', 'u' and 'f' " 179 write(*,*)
" have been derived in " 180 write(*,*)
" choral/maple/example_elasticity_Dirichlet.mw" 186 write(*,*)
"================================== MESH ASSEMBLING" 187 msh = mesh(mesh_file,
'gmsh')
192 write(*,*)
"================================== FINITE ELEMENT& 199 y = fespacexk(x_h, 2)
203 write(*,*)
"================================== INTEGRATION METHOD& 204 & ON \Omega ASSEMBLING" 211 write(*,*)
"================================== STIFFNESS& 213 call elasticity_stiffmat(stiff,
lambda,
mu, &
218 write(*,*)
"================================== RIGHT HAND SIDE" 219 write(*,*)
" 1- COMPUTATION OF F" 222 call l2_product(rhs, y, qdm,
f_1,
f_2)
226 write(*,*)
"================================== DIRICHLET CONDITION" 227 call elasticity_dirichlet(stiff, rhs, y,
u_1,
u_2,
gam)
231 write(*,*)
"================================== SOLVE LINEAR SYSTEM& 235 kry = krylov(
kry_gmres, tol=1
e-6_rp, itmax=1000, restart=25, verb=2)
239 call interp_vect_func(u_h, y,
u_1,
u_2)
242 call solve(u_h, kry, rhs, stiff)
249 write(*,*)
"================================== POST-TREATMENT" 253 call extract_component(u_h_2, u_h, y, 2)
254 call extract_component(u_h_1, u_h, y, 1)
261 shell_com =
'elasticity_Dirichlet_u.msh' 264 call write(x_h, trim(shell_com),
'gmsh')
268 call gmsh_addview(x_h, u_h_1, &
270 &
'Elasticity problem : numerical solution', 0.0_rp, 1)
271 call gmsh_addview(x_h, u_h_2, &
273 &
'Elasticity problem : numerical solution', 1.0_rp, 2)
276 shell_com =
'gmsh -option '//trim(gmsh_dir)//
'gmsh-options-view'&
277 & //
' '//trim(shell_com)
278 call system(shell_com)
285 write(*,*)
"Numerical errors" 286 err_l2 = l2_dist(u_h, y, qdm,
u_1,
u_2)
288 write(*,*)
" Error | u - u_h|_L2 =", err_l2
289 write(*,*)
" Error |\nabla (u - u_h)|_L2 =", err_h1
299 function lambda(x)
result(r)
301 real(RP),
dimension(3),
intent(in) :: x
303 r = 1.0_rp +
cc4(x) * 0.2_rp
309 function mu(x)
result(r)
311 real(RP),
dimension(3),
intent(in) :: x
313 r = 1.0_rp +
ss4(x) * 0.2_rp
319 function u_1(x)
result(r)
321 real(RP),
dimension(3),
intent(in) :: x
326 function u_2(x)
result(r)
328 real(RP),
dimension(3),
intent(in) :: x
339 real(RP),
dimension(3) :: grad
340 real(RP),
dimension(3),
intent(in) :: x
342 grad(1) = -pi *
ss(x)
348 real(RP),
dimension(3) :: grad
349 real(RP),
dimension(3),
intent(in) :: x
352 grad(2) = -pi *
ss(x)
360 function f_1(x)
result(r)
361 real(RP),
dimension(3),
intent(in) :: x
366 r = r + 0.2_rp*
cc4(x)
367 r = r + 0.4_rp*
ss4(x)
370 r = r +
cs4(x)*
ss(x) * re(8,5)
371 r = r -
sc4(x)*
ss(x) * re(8,5)
372 r = r -
sc4(x)*
cc(x) * re(8,5)
377 function f_2(x)
result(r)
379 real(RP),
dimension(3),
intent(in) :: x
383 r = r + 0.2_rp*
cc4(x)
384 r = r + 0.4_rp*
ss4(x)
387 r = r +
sc4(x)*
ss(x) * re(8,5)
388 r = r -
cs4(x)*
cc(x) * re(8,5)
389 r = r -
cs4(x)*
ss(x) * re(8,5)
397 function gam(x)
result(r)
399 real(RP),
dimension(3),
intent(in) :: x
401 r = - (x(1) - 1.0_rp) * (x(2) - 1.0_rp) * x(1) * x(2)
406 function cc(x)
result(r)
408 real(RP),
dimension(3),
intent(in) :: x
410 r = cos(pi*x(1) ) * cos(pi*x(2))
413 function ss(x)
result(r)
415 real(RP),
dimension(3),
intent(in) :: x
417 r = sin(pi*x(1) ) * sin(pi*x(2))
420 function cs(x)
result(r)
422 real(RP),
dimension(3),
intent(in) :: x
424 r = cos(pi*x(1) ) * sin(pi*x(2))
427 function sc(x)
result(r)
429 real(RP),
dimension(3),
intent(in) :: x
431 r = sin(pi*x(1) ) * cos(pi*x(2))
436 function cc4(x)
result(r)
438 real(RP),
dimension(3),
intent(in) :: x
440 r = cos(4.0_rp*pi*x(1) ) * cos(4.0_rp*pi*x(2))
443 function ss4(x)
result(r)
445 real(RP),
dimension(3),
intent(in) :: x
447 r = sin(4.0_rp*pi*x(1) ) * sin(4.0_rp*pi*x(2))
450 function cs4(x)
result(r)
452 real(RP),
dimension(3),
intent(in) :: x
454 r = cos(4.0_rp*pi*x(1) ) * sin(4.0_rp*pi*x(2))
457 function sc4(x)
result(r)
459 real(RP),
dimension(3),
intent(in) :: x
461 r = sin(4.0_rp*pi*x(1) ) * cos(4.0_rp*pi*x(2))
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
integer, parameter fe_p2_2d
integer, parameter kry_gmres
GmRes linear solver.
program elasticity_dirichlet_2d
SOLVES THE LINEAR ELASTICITY PROBLEM with a Dirichlet boundary conditions
real(rp) function, dimension(3) grad_u_1(x)
gradient of the exact solution
real(rp) function e(x, v1, v2)
real(rp) function, dimension(3) grad_u_2(x)
real(rp) function mu(x)
Lame coefficient 'mu'.
real(rp) function gam(x)
To characterise .
real(rp) function lambda(x)
Lame coefficient 'lambda'.
integer, parameter quad_gauss_trg_12
real(rp) function f_1(x)
right hand side 'f'= =(f_1, f_2)
integer, parameter fe_p2_1d
real(rp) function u_1(x)
exact solution u=(u_1, u_2)