111 integer ,
parameter :: verb = 2
116 real(RP),
parameter :: penal = 1e10_rp
117 real(RP),
parameter :: tol = 1
e-15_rp
128 integer ,
parameter :: fe_v =
fe_p2_2d 129 integer ,
parameter :: fe_s =
fe_p2_1d 132 character(LEN=1),
parameter :: mesh_idx =
"2" 140 type(quadmesh) :: qdm
143 character(len=100),
parameter :: mesh_file= &
144 & trim(GMSH_DIR)//
"square/square_"//mesh_idx//
".msh" 158 type(csr) :: mass, stiff, K
159 real(RP),
dimension(:),
allocatable :: rhs, u_h, aux
164 real(RP) :: err_L2, err_H1
166 character(LEN=300) :: shell_com
179 write(*,*)
'diffusion_Dirichlet_Neumann: start' 181 write(*,*)
" EXAMPLE OF AN ELLIPTIC PROBLEM" 183 write(*,*)
" -\Delta(u) + u = f on \Omega" 185 write(*,*)
" WITH A MIXED DIRICHLET / NEUMANN BOUNDARY CONDITION " 190 write(*,*)
"================================== MESH ASSEMBLING" 191 msh = mesh(mesh_file,
'gmsh')
196 write(*,*)
"================================== FINITE ELEMENT& 205 write(*,*)
"================================== INTEGRATION METHOD& 206 & ON \Omega ASSEMBLING" 213 write(*,*)
"================================== MASS AND STIFFNESS& 215 call diffusion_massmat(mass , one_r3, x_h, qdm)
216 call diffusion_stiffmat(stiff, emetric, x_h, qdm)
220 write(*,*)
"================================== LINEAR SYSTEM MATRIX" 223 call add(k, stiff, 1._rp, mass, 1._rp)
230 write(*,*)
"================================== RIGHT HAND SIDE" 231 write(*,*)
" 1- COMPUTATION OF & 232 & THE VOLUMIC SOURCE TERM F" 235 call l2_product(rhs,
f, x_h, qdm)
239 write(*,*)
"================================== RIGHT HAND SIDE" 240 write(*,*)
" 2- COMPUTATION OF & 241 & THE SURFACIC SOURCE TERM G (NEUMANN CONDITION)" 244 call diffusion_neumann_rhs(aux,
dx_u_1, x_h, qd_s,
x_ge_1)
248 call diffusion_neumann_rhs(aux,
dy_u_1, x_h, qd_s,
y_ge_1)
252 write(*,*)
"================================== RIGHT HAND SIDE" 253 write(*,*)
" 3- PENALISATION OF & 254 & THE SYSTEM FOR THE DIRICHLET CONDITION" 257 call diffusion_dirichlet(k, rhs,
u, x_h, penal,
x_le_0)
260 call diffusion_dirichlet(k, rhs,
u, x_h, penal,
y_le_0)
263 write(*,*)
"================================== SOLVE LINEAR SYSTEM& 267 kry = krylov(
kry_cg, tol=tol, itmax=1000, verb=2)
271 call interp_scal_func(u_h,
u, x_h)
275 call solve(u_h, kry, rhs, k)
282 write(*,*)
"================================== POST-TREATMENT" 287 shell_com =
'diff_Neumann_u.msh' 290 call write(x_h, trim(shell_com),
'gmsh')
293 call gmsh_addview(x_h, u_h, &
295 &
'Diffusion problem: numerical solution', 0.0_rp, 1)
298 shell_com =
'gmsh -option '//trim(gmsh_dir)//
'gmsh-options-view'&
299 & //
' '//trim(shell_com)
300 call system(shell_com)
304 write(*,*)
"Numerical errors" 305 err_l2 = l2_dist(
u, u_h, x_h, qdm)
306 err_h1 = l2_dist_grad(
grad_u, u_h, x_h, qdm)
307 write(*,*)
" Error | u - u_h|_L2 =", err_l2
308 write(*,*)
" Error |\nabla (u - u_h)|_L2 =", err_h1
319 real(RP),
dimension(3),
intent(in) :: x
321 u = cos( x(1) ) * (1._rp + x(2)**2)
329 real(RP),
dimension(3) :: grad_u
330 real(RP),
dimension(3),
intent(in) :: x
332 grad_u(1) =-sin(x(1)) * (1._rp + x(2)**2)
333 grad_u(2) = cos(x(1)) * 2._rp * x(2)
343 real(RP),
dimension(3),
intent(in) :: x
345 f = 2._rp *
u(x) - 2._rp * cos( x(1) )
353 function x_le_0(x)
result(r)
355 real(RP),
dimension(3),
intent(in) :: x
364 function x_ge_1(x)
result(r)
366 real(RP),
dimension(3),
intent(in) :: x
374 function y_le_0(x)
result(r)
376 real(RP),
dimension(3),
intent(in) :: x
385 function y_ge_1(x)
result(r)
387 real(RP),
dimension(3),
intent(in) :: x
396 function dx_u_1(x)
result(r)
398 real(RP),
dimension(3),
intent(in) :: x
400 r = -sin( 1.0_rp ) * ( 1.0_rp + x(2)**2 )
406 function dy_u_1(x)
result(r)
408 real(RP),
dimension(3),
intent(in) :: x
410 r = cos( x(1) ) * 2.0_rp
real(rp) function x_ge_1(x)
To characterise {x=1}.
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
integer, parameter fe_p2_2d
integer, parameter quad_gauss_edg_4
real(rp) function x_le_0(x)
To characterise {x=0}.
real(rp) function e(x, v1, v2)
real(rp) function y_le_0(x)
To characterise {x=0}.
real(rp) function, dimension(3) grad_u(x)
gradient of the exact solution
real(rp) function f(x)
right hand side 'f'
real(rp) function u(x)
exact solution 'u'
real(rp) function y_ge_1(x)
To characterise {x=1}.
integer, parameter quad_gauss_trg_12
real(rp) function dx_u_1(x)
u on {x=1}
program diffusion_dirichlet_neumann
SOLVES THE DIFFUSION PROBLEM with mixed Dirichlet / Neumann boundary conditions ...
integer, parameter fe_p2_1d
integer, parameter kry_cg
CG linear solver.
real(rp) function dy_u_1(x)
u on {y=1}