Choral
Functions/Subroutines
diffusion_Dirichlet_Neumann.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

program diffusion_dirichlet_neumann
  SOLVES THE DIFFUSION PROBLEM with mixed Dirichlet / Neumann boundary conditions More...
 
real(rp) function u (x)
 exact solution 'u' More...
 
real(rp) function, dimension(3) grad_u (x)
 gradient of the exact solution More...
 
real(rp) function f (x)
 right hand side 'f' More...
 
real(rp) function x_le_0 (x)
 To characterise {x=0}. More...
 
real(rp) function x_ge_1 (x)
 To characterise {x=1}. More...
 
real(rp) function y_le_0 (x)
 To characterise {x=0}. More...
 
real(rp) function y_ge_1 (x)
 To characterise {x=1}. More...
 
real(rp) function dx_u_1 (x)
 u on {x=1} More...
 
real(rp) function dy_u_1 (x)
 u on {y=1} More...
 

Function/Subroutine Documentation

◆ diffusion_dirichlet_neumann()

program diffusion_dirichlet_neumann ( )

SOLVES THE DIFFUSION PROBLEM with mixed Dirichlet / Neumann boundary conditions

\(~~~~~~~~~ -\Delta u + u = f ~~~\) on \(~~~ \Omega= [0,1]^2 \)

Dirichlet boundary condition: \( u = g_D \) on \( \Gamma_D\subset \Gamma = \partial \Omega\):

  • \( u = 1+y^2 \) on \(\Gamma \cap \{x=0\} \)
  • \( u = \cos x \) on \(\Gamma \cap \{y=0\} \)

Neumann boundary condition: \( \nabla u\cdot n = g_N \) on \( \Gamma_N\subset \Gamma\):

  • \( \partial_x u = -\sin(1) (1+y^2) \) on \(\Gamma \cap \{x=1\} \)
  • \( \partial_y u = 2 \cos(x) \) on \(\Gamma \cap \{y=0\} \)

Volumic source term and exact solution:

  • \( u = cos(x) (1+y^2)\)
  • \( f = 2u -2 cos x\)



VARIATIONAL FORMULATION

  • The finite element space is \( X_h\subset \Hu\).
  • Let \( X_{h,D} \subset X_h\) the subspace of all functions vanishing on \( \Gamma_D \).

Numerical probelm: find \(u\in X_h\) such that \( ~~~~\forall ~v \in X_{h,D}, ~~~~ \)

\[ \int_\Omega \nabla u\cdot\nabla v \,\dx ~+~ \int_\Omega u\, v \,\dx ~=~ \int_\Omega f \, v \,\dx ~+~ \int_{\Gamma_N} g_N \, v \,\dx \]

and such that \( u(x_i) = g_D(x_i) \) for all finite element nodes \( x_i\in \Gamma_D\)



NUMERICAL RESOLUTION

The basis functions are \( (v_i)_{1\le i\le N} \) (basis of \(X_h\)).

By construction, a subset of the basis functions form a basis of \( X_{h,D}\).

  • Computation of the stiffness matrix

    \[ S =[s_{i,\,j}]_{1\le i,\,j\le N}, \quad \quad s_{i,\,j} = \int_\Omega \nabla v_i\cdot\nabla v_j \,\dx \]

  • Computation of the mass matrix

    \[ M =[s_{i,\,j}]_{1\le i,\,j\le N}, \quad \quad m_{i,\,j} = \int_\Omega v_i\, v_j \,\dx \]

  • Computation of the right hand side for the volumoc source term \( f \)

    \[ F = (f_i)_{1\le i\le N}, \quad \quad f_i = \int_\Omega f \, v_i \,\dx \]

  • Computation of the right hand side for the surfacic source term \( g_N \)

    \[ G_N = (g_i)_{1\le i\le N}, \quad \quad g_i = \int_{\Gamma_N} g_N \, v_i \,\dx \]

  • \( K = M + S \) and \( \text{RHS} = F + G_N \)

Penalisation of the linear system \( K U =\)RHS:

  • let \( 1\le i \le N \) be so that the finite element node \( x_i \in\Gamma_D\),
  • \( \text{RHS}_i =\text{RHS}_i + \rho g_D(x_i)\)
  • \( K_{i,i} =K_{i,i} + \rho \)
  • \(\rho \gg 1\) a penalisation coefficient.

Resolution of the (symmetric positive definite) system

\[ K U_h = \text{RHS}\]



POST TREATMENT

  • Computation of the numerical error between the exact solution \(u\) and the numerical solution \(u_h\),

    \[ \int_\Omega |u-u_h|^2 \dx~,\quad\quad \int_\Omega |\nabla u-\nabla u_h|^2 \dx \]

  • Graphical display of the numerical solution with gmsh

Charles PIERRE, November 2019

Definition at line 95 of file diffusion_Dirichlet_Neumann.f90.

Here is the call graph for this function:

◆ dx_u_1()

real(rp) function diffusion_dirichlet_neumann::dx_u_1 ( real(rp), dimension(3), intent(in)  x)

u on {x=1}

Definition at line 397 of file diffusion_Dirichlet_Neumann.f90.

◆ dy_u_1()

real(rp) function diffusion_dirichlet_neumann::dy_u_1 ( real(rp), dimension(3), intent(in)  x)

u on {y=1}

Definition at line 407 of file diffusion_Dirichlet_Neumann.f90.

◆ f()

real(rp) function diffusion_dirichlet_neumann::f ( real(rp), dimension(3), intent(in)  x)

right hand side 'f'

Definition at line 342 of file diffusion_Dirichlet_Neumann.f90.

Here is the call graph for this function:

◆ grad_u()

real(rp) function, dimension(3) diffusion_dirichlet_neumann::grad_u ( real(rp), dimension(3), intent(in)  x)

gradient of the exact solution

Definition at line 329 of file diffusion_Dirichlet_Neumann.f90.

◆ u()

real(rp) function diffusion_dirichlet_neumann::u ( real(rp), dimension(3), intent(in)  x)

exact solution 'u'

Definition at line 318 of file diffusion_Dirichlet_Neumann.f90.

◆ x_ge_1()

real(rp) function diffusion_dirichlet_neumann::x_ge_1 ( real(rp), dimension(3), intent(in)  x)

To characterise {x=1}.

Definition at line 365 of file diffusion_Dirichlet_Neumann.f90.

◆ x_le_0()

real(rp) function diffusion_dirichlet_neumann::x_le_0 ( real(rp), dimension(3), intent(in)  x)

To characterise {x=0}.

Definition at line 354 of file diffusion_Dirichlet_Neumann.f90.

◆ y_ge_1()

real(rp) function diffusion_dirichlet_neumann::y_ge_1 ( real(rp), dimension(3), intent(in)  x)

To characterise {x=1}.

Definition at line 386 of file diffusion_Dirichlet_Neumann.f90.

◆ y_le_0()

real(rp) function diffusion_dirichlet_neumann::y_le_0 ( real(rp), dimension(3), intent(in)  x)

To characterise {x=0}.

Definition at line 375 of file diffusion_Dirichlet_Neumann.f90.