Choral
elasticity_Dirichlet_2D.f90
Go to the documentation of this file.
1 
87 
89 
91  use choral
92 
93  implicit none
94 
95  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97  !!
98  !! VARAIBLE DEFINITION : BEGIN
99  !!
100  !!
101 
102  !! verb = verbosity level
103  !!
104  integer , parameter :: verb = 2
105 
106 
107  !! SPACE DISCRETISATION
108  !!
109  !! fe_v = finite element method (volumic)
110  !! fe_s = finite element method (surfacic)
111  !! qd_v = quadrature method (volumic)
112  !! mesh_idx = index for the mesh
113  !!
114  integer , parameter :: fe_v = fe_p2_2d
115  integer , parameter :: fe_s = fe_p2_1d
116  integer , parameter :: qd_v = quad_gauss_trg_12
117  character(LEN=1), parameter :: mesh_idx = "2"
118  !!
119  !! m = mesh
120  !! X_h = finite element space
121  !! Y = [X_h]^dim
122  !! qdm = integration method
123  !!
124  type(mesh) :: msh
125  type(fespace) :: X_h
126  type(fespacexk) :: Y
127  type(quadmesh) :: qdm
128  !!
129  !! msh_file
130  character(len=100), parameter :: mesh_file= &
131  & trim(GMSH_DIR)//"square/square_"//mesh_idx//".msh"
132 
133 
134  !! LINEAR SYSTEM
135  !!
136  !! kry = krylov method def.
137  !! stiff = stiffness matrix
138  !! rhs = right hand side
139  !! u_h = numerical solution
140  !!
141  type(krylov) :: kry
142  type(csr) :: stiff
143  real(RP), dimension(:), allocatable :: rhs, u_h
144  real(RP), dimension(:), allocatable :: u_h_1, u_h_2
145 
146 
147  !! NUMERICAL ERRORS
148  !!
149  real(RP) :: err_L2, err_H1
150 
151  character(LEN=300) :: shell_com
152 
153  !!
154  !! VARAIBLE DEFINITION : END
155  !!
156  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
157  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158 
159 
160  !! !!!!!!!!!!!!!!!!!!!!! INITIALISATION
161  !!
162  call choral_init(verb=1)
163  write(*,*)
164  write(*,*)'elasticity_Dirichlet_2D: start'
165  write(*,*)""
166  write(*,*)" EXAMPLE OF A LINEAR ELASTICITY PROBLEM"
167  write(*,*)""
168  write(*,*)" -div(A(x) e(u)) = f on \Omega"
169  write(*,*)""
170  write(*,*)" WITH A DIRICHLET BOUNDARY CONDITION "
171  write(*,*)""
172  write(*,*)" A(x)e(u)) n = g on \partial\Omega"
173  write(*,*)""
174  write(*,*)" AND WITH NON CONSTANT LAME COEFFICIENTS:"
175  write(*,*)""
176  write(*,*)" A(x) \xi = lambda(x) Tr(\xi) Id + 2 mu(x) \si"
177  write(*,*)""
178  write(*,*)" The problem data 'lambda', 'mu', 'u' and 'f' "
179  write(*,*)" have been derived in "
180  write(*,*)" choral/maple/example_elasticity_Dirichlet.mw"
181  write(*,*)""
182  write(*,*)""
183  write(*,*)""
184 
185  write(*,*) ""
186  write(*,*)"================================== MESH ASSEMBLING"
187  msh = mesh(mesh_file, 'gmsh')
188 
189  call print(msh)
190 
191  write(*,*) ""
192  write(*,*)"================================== FINITE ELEMENT&
193  & SPACE ASSEMBLING"
194  x_h =fespace(msh)
195  call set(x_h, fe_v)
196  call set(x_h, fe_s)
197  call assemble(x_h)
198  call print(x_h)
199  y = fespacexk(x_h, 2)
200  call print(y)
201 
202  write(*,*) ""
203  write(*,*)"================================== INTEGRATION METHOD&
204  & ON \Omega ASSEMBLING"
205  qdm = quadmesh(msh)
206  call set(qdm, qd_v)
207  call assemble(qdm)
208  call print(qdm)
209 
210  write(*,*) ""
211  write(*,*)"================================== STIFFNESS&
212  & MATRIX ASSEMBLING"
213  call elasticity_stiffmat(stiff, lambda, mu, &
214  & y, qdm)
215 
216 
217  write(*,*) ""
218  write(*,*)"================================== RIGHT HAND SIDE"
219  write(*,*)" 1- COMPUTATION OF F"
220  !!
221  !! \int_\Omega f u_i dx
222  call l2_product(rhs, y, qdm, f_1, f_2)
223 
224 
225  write(*,*) ""
226  write(*,*)"================================== DIRICHLET CONDITION"
227  call elasticity_dirichlet(stiff, rhs, y, u_1, u_2, gam)
228 
229 
230  write(*,*) ""
231  write(*,*)"================================== SOLVE LINEAR SYSTEM&
232  & Stiff U_h = RHS "
233  !!
234  !! Solver setting
235  kry = krylov(kry_gmres, tol=1e-6_rp, itmax=1000, restart=25, verb=2)
236  call print(kry)
237  !!
238  !! initial guess
239  call interp_vect_func(u_h, y, u_1, u_2)
240  !!
241  !! linear system inversion
242  call solve(u_h, kry, rhs, stiff)
243  !!
244  !! cleaning
245  call clear(stiff)
246  call freemem(rhs)
247 
248  write(*,*) ""
249  write(*,*)"================================== POST-TREATMENT"
250  !!
251  !! extract component
252  !!
253  call extract_component(u_h_2, u_h, y, 2)
254  call extract_component(u_h_1, u_h, y, 1)
255 
256  !!
257  !! graphical display
258  !!
259  if (verb>0) then
260  !! output file def.
261  shell_com = 'elasticity_Dirichlet_u.msh'
262  !!
263  !! first write the 'finite element mesh'
264  call write(x_h, trim(shell_com), 'gmsh')
265  !!
266  !! second save the finite element solution u_h
267  !!
268  call gmsh_addview(x_h, u_h_1, &
269  & trim(shell_com), &
270  &'Elasticity problem : numerical solution', 0.0_rp, 1)
271  call gmsh_addview(x_h, u_h_2, &
272  & trim(shell_com), &
273  &'Elasticity problem : numerical solution', 1.0_rp, 2)
274  !!
275  !! shell command to visualise with gmsh
276  shell_com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'&
277  & //' '//trim(shell_com)
278  call system(shell_com)
279 
280  end if
281 
282  !!
283  !! L2 and H1_0 numerical errors
284  write(*,*) ""
285  write(*,*) "Numerical errors"
286  err_l2 = l2_dist(u_h, y, qdm, u_1, u_2)
287  err_h1 = l2_dist_grad(u_h, y, qdm, grad_u_1, grad_u_2)
288  write(*,*)" Error | u - u_h|_L2 =", err_l2
289  write(*,*)" Error |\nabla (u - u_h)|_L2 =", err_h1
290 
291  call freemem(u_h)
292  call freemem(u_h_1)
293  call freemem(u_h_2)
294 
295 contains
296 
299  function lambda(x) result(r)
300  real(RP) :: r
301  real(RP), dimension(3), intent(in) :: x
302 
303  r = 1.0_rp + cc4(x) * 0.2_rp
304 
305  end function lambda
306 
309  function mu(x) result(r)
310  real(RP) :: r
311  real(RP), dimension(3), intent(in) :: x
312 
313  r = 1.0_rp + ss4(x) * 0.2_rp
314 
315  end function mu
316 
319  function u_1(x) result(r)
320  real(RP) :: r
321  real(RP), dimension(3), intent(in) :: x
322 
323  r = cs(x)
324 
325  end function u_1
326  function u_2(x) result(r)
327  real(RP) :: r
328  real(RP), dimension(3), intent(in) :: x
329 
330  r = sc(x)
331 
332  end function u_2
333 
334 
335 
338  function grad_u_1(x) result(grad)
339  real(RP), dimension(3) :: grad
340  real(RP), dimension(3), intent(in) :: x
341 
342  grad(1) = -pi * ss(x)
343  grad(2) = pi * cc(x)
344  grad(3) = 0._rp
345 
346  end function grad_u_1
347  function grad_u_2(x) result(grad)
348  real(RP), dimension(3) :: grad
349  real(RP), dimension(3), intent(in) :: x
350 
351  grad(1) = pi * cc(x)
352  grad(2) = -pi * ss(x)
353  grad(3) = 0._rp
354 
355  end function grad_u_2
356 
357 
360  function f_1(x) result(r)
361  real(RP), dimension(3), intent(in) :: x
362  real(RP) :: r
363 
364  r = 6.0_rp
365 
366  r = r + 0.2_rp*cc4(x)
367  r = r + 0.4_rp*ss4(x)
368  r = r * cs(x)
369 
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)
373 
374  r = r * pi**2
375 
376  end function f_1
377  function f_2(x) result(r)
378  real(RP) :: r
379  real(RP), dimension(3), intent(in) :: x
380 
381  r = 6.0_rp
382 
383  r = r + 0.2_rp*cc4(x)
384  r = r + 0.4_rp*ss4(x)
385  r = r * sc(x)
386 
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)
390 
391  r = r * pi**2
392 
393  end function f_2
394 
397  function gam(x) result(r)
398  real(RP) :: r
399  real(RP), dimension(3), intent(in) :: x
400 
401  r = - (x(1) - 1.0_rp) * (x(2) - 1.0_rp) * x(1) * x(2)
402 
403  end function gam
404 
405 
406  function cc(x) result(r)
407  real(RP) :: r
408  real(RP), dimension(3), intent(in) :: x
409 
410  r = cos(pi*x(1) ) * cos(pi*x(2))
411 
412  end function cc
413  function ss(x) result(r)
414  real(RP) :: r
415  real(RP), dimension(3), intent(in) :: x
416 
417  r = sin(pi*x(1) ) * sin(pi*x(2))
418 
419  end function ss
420  function cs(x) result(r)
421  real(RP) :: r
422  real(RP), dimension(3), intent(in) :: x
423 
424  r = cos(pi*x(1) ) * sin(pi*x(2))
425 
426  end function cs
427  function sc(x) result(r)
428  real(RP) :: r
429  real(RP), dimension(3), intent(in) :: x
430 
431  r = sin(pi*x(1) ) * cos(pi*x(2))
432 
433  end function sc
434 
435 
436  function cc4(x) result(r)
437  real(RP) :: r
438  real(RP), dimension(3), intent(in) :: x
439 
440  r = cos(4.0_rp*pi*x(1) ) * cos(4.0_rp*pi*x(2))
441 
442  end function cc4
443  function ss4(x) result(r)
444  real(RP) :: r
445  real(RP), dimension(3), intent(in) :: x
446 
447  r = sin(4.0_rp*pi*x(1) ) * sin(4.0_rp*pi*x(2))
448 
449  end function ss4
450  function cs4(x) result(r)
451  real(RP) :: r
452  real(RP), dimension(3), intent(in) :: x
453 
454  r = cos(4.0_rp*pi*x(1) ) * sin(4.0_rp*pi*x(2))
455 
456  end function cs4
457  function sc4(x) result(r)
458  real(RP) :: r
459  real(RP), dimension(3), intent(in) :: x
460 
461  r = sin(4.0_rp*pi*x(1) ) * cos(4.0_rp*pi*x(2))
462 
463  end function sc4
464 
465 
466 
467 end program elasticity_dirichlet_2d
real(rp) function sc4(x)
real(rp) function cc(x)
real(rp) function ss(x)
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
Definition: choral.F90:166
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
Definition: choral.F90:27
integer, parameter fe_p2_2d
integer, parameter kry_gmres
GmRes linear solver.
real(rp) function f_2(x)
program elasticity_dirichlet_2d
SOLVES THE LINEAR ELASTICITY PROBLEM with a Dirichlet boundary conditions
real(rp) function cs(x)
real(rp) function, dimension(3) grad_u_1(x)
gradient of the exact solution
real(rp) function cc4(x)
real(rp) function ss4(x)
real(rp) function e(x, v1, v2)
CHORAL CONSTANTS
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 cs4(x)
real(rp) function lambda(x)
Lame coefficient 'lambda'.
integer, parameter quad_gauss_trg_12
real(rp) function sc(x)
real(rp) function f_1(x)
right hand side 'f'= =(f_1, f_2)
integer, parameter fe_p2_1d
real(rp) function u_2(x)
real(rp) function u_1(x)
exact solution u=(u_1, u_2)