Commit a27dbfc0 authored by rabgra's avatar rabgra
Browse files

Merge branch 'master' of git.math.uzh.ch:remi.abgrall/RD_public

parents 4f7fe35f 52a6e305
File added
\documentclass{article}
\usepackage{fullpage,url}
\usepackage[dvipsnames]{xcolor}
\definecolor{airforceblue}{rgb}{0.36, 0.54, 0.66}
\begin{document}
\title{Euler code, using Residual distribution scheme, unsteady, order 2,3,4 in time/space and hybrid meshes. Uses Mood technology, and compatible with additional conservation laws (Kinetic momentum for example)}
\author{R. Abgrall, P. Bacigaluppi, P. Oeffner, S. Tokareva, D. Torlo
......@@ -10,22 +12,24 @@ remi.abgrall@math.uzh.ch}
Warning: works for Euler, the other models need to be checked.
\section{Compilation}
Makefile in Make/Makefile\_2D.gfortran. Uses gfortran
Please find the Makefile in Test1D and Test2D accordingly.\\\\
In 1D the file is called \texttt{Makefile\_1D.gfortran}, while in 2D it is \texttt{Makefile\_2D.gfortran}. To compile go in the Test1D or Test2D folder and type:
\begin{verbatim}
\make -f Make//Makefile_2D.gfortran dec
make -f Make/Makefile_2D.gfortran dec
\end{verbatim}
\begin{verbatim}
\make -f Make//Makefile_2D.gfortran clean
make -f Make/Makefile_2D.gfortran clean
\end{verbatim}
Note:LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib for mac
Note \#1: The Makefile in Test\# D uses \textbf{gfortran}. To compile with different settings, check the folder Make on the main folder.\\
Note \#2: LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib for mac.\\\\
Note \#3: \textcolor{airforceblue}{According to the desired system set-up, uncomment/comment the required modules in the Makefile. The only modules to be changed are the \texttt{variable\_def\_*} and \texttt{init\_bc\_*}.}
\section{Meshes}
gmsh format msh2: gmsh -format msh2 geofile.geo
\section{Run}
test case: see Boundary\_euler and init_bc_euler (they should be consistant).
test case: see Boundary\_euler and init\_bc\_euler (they should be consistant).
\section{Reference:}
\begin{itemize}
......
\documentclass{article}
\usepackage{fullpage,url}
\begin{document}
\title{Euler code, using Residual distribution scheme, unsteady, order 2,3,4 in time/space and hybrid meshes. Uses Mood technology, and compatible with additional conservation laws (Kinetic momentum for example)}
\author{R. Abgrall, P. Bacigaluppi, P. Oeffner, S. Tokareva, D. Torlo
F. Mojarrad\\
remi.abgrall@math.uzh.ch}
\date{}
\maketitle
Warning: works for Euler, the other models need to be checked.
\section{Compilation}
Please find the Makefile in Test1D and Test2D accordingly.\\\\
In 1D the file is called \texttt{Makefile\_2D.gfortran}, while in 2D it is \texttt{Makefile\_2D.gfortran}. To compile go in the Test1D or Test2D folder and type:
\begin{verbatim}
make -f Make/Makefile_2D.gfortran dec
\end{verbatim}
\begin{verbatim}
make -f Make/Makefile_2D.gfortran clean
\end{verbatim}
Note \#1: The Makefile in Test\# D uses \textbf{gfortran}. To compile with different settings, check the folder Make on the main folder.\\
Note \#2: LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib for mac.\\\\
Note \#3: \alert{According to the desired system set-up, uncomment/comment the required modules in the Makefile. The only modules to be changed are the \texttt{variable\_def\_*} and \texttt{init\_bc\_*}.}
\section{Meshes}
gmsh format msh2: gmsh -format msh2 geofile.geo
\section{Run}
test case: see Boundary\_euler and init\_bc\_euler (they should be consistant).
\section{Reference:}
\begin{itemize}
\item R. Abgrall, J. Nordstr\"om, P. \"Offner and S. Tokareva, Analysis of the SBP-SAT stabilisation for finite element methods, part II: the non-linear case, Communications in Applied Mathematics and Computation, 2021, DOI: \url{10.1007/s42967-020-00086-2}
\item R. Abgrall, J. Nordstr\"om, P. \"Offner and S. Tokareva, Analysis of the SBP-SAT stabilisation for finite element methods, part I: the linear case, R. Abgrall, J. Nordstr\"om, P. \"Offner and S. Tokareva, Analysis of the SBP-SAT stabilisation for finite element methods, part I: the linear case, J. Sci. Comput. 85 (2020), no. 2, Paper No. 43.
\item R. Abgrall and D. Torlo, High Order Asymptotic Preserving Deferred Correction Implicit-Explicit Schemes for Kinetic Models, SIAM SISC, 2020, v42(3), pp B816-845, \url{https://arxiv.org/abs/1811.09284}
\item R. Abgrall, A general framework to construct schemes satisfying additional conservation relations, application to entropy conservative and entropy dissipative schemes, J. Comput. Phys, vol 372(1), 2018
\item R. Abgrall, P. Bacigaluppi and S. Tokareva, High-order residual distribution scheme for the time-dependent Euler equations of fluid dynamics, Computer \& Mathematics with Applications, 2019, vol 78 (2), pages 274-297
\item R. Abgrall, Some remarks about conservation for residual distribution schemes, Computational Methods in Applied Mathematics, v18(3), pp 327-350, 2018, doi:\url{https://doi.org/10.1515/cmam-2017-0056}
\item R. Abgrall, P. Bacigaluppi and S. Tokareva, A high-order nonconservative approach for hyperbolic equations in fluid dynamics, Computers and Fluids, vol 169, pages 10-22, 2018
doi: \url{https://doi.org/10.1016/j.compfluid.2017.08.019}
\item R. Abgrall, High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices., Journal of Scientific Computing, 73(2-3), pp 461-494, 2017
\end{itemize}
\end{document}
\ No newline at end of file
......@@ -29,7 +29,7 @@ CONTAINS
u_cons=u
SELECT CASE(e%itype)
CASE(1,3,4) ! Lagrange
CASE(1,3,4,7,11,12,13,14) ! Lagrange
CASE(2,5,6) ! Bezier
......@@ -54,7 +54,7 @@ CONTAINS
u=u_cons
SELECT CASE(e%itype)
CASE(1,3,4) ! Lagrange
CASE(1,3,4,7,11,12,13,14) ! Lagrange
CASE(2,5,6)! cubic bezier
......
......@@ -37,14 +37,20 @@ MODULE element_class
TYPE, PUBLIC:: element
INTEGER:: type_flux=5
INTEGER:: diag=-1
INTEGER:: type_flux=-10
INTEGER:: diag=-1, diag2=-1
INTEGER:: nsommets, itype, nvertex ! nbre de dofs, type element:
!1->P1
!2->B2
!3->P2
!4->P3
!5->B2
!5->B3
!6->B4
!7->P4
!11->PGL1
!12->PGL2
!13->PGL3
!14->PGL4
REAL(dp), DIMENSION(:), POINTER :: coor =>NULL() !
! *-----*-----* for quadratic, *---*---*---* for cubic elements
! 1 3 2 1 3 4 2
......@@ -77,6 +83,7 @@ MODULE element_class
PROCEDURE, PUBLIC:: eval_der2=>eval_der2_element
PROCEDURE, PUBLIC:: der_sec=>der_sec_element
PROCEDURE, PUBLIC:: eval_coeff=>eval_coeff_element
PROCEDURE, PUBLIC:: der2_poly=>der2_poly_element
FINAL:: clean
END TYPE element
......@@ -104,9 +111,9 @@ CONTAINS
CLASS(element), INTENT(in):: e
INTEGER, INTENT(in):: k ! index of basis function
REAL(dp), DIMENSION(2), INTENT(in):: x ! barycentric coordinate
REAL(dp) :: alpha, beta, s
SELECT CASE(e%itype)
CASE(1) ! P1
CASE(1,11) ! P1
SELECT CASE(k)
CASE(2)
base_element=x(1)
......@@ -134,7 +141,7 @@ CONTAINS
STOP
END SELECT
CASE(3)! P2
CASE(3,12)! P2
SELECT CASE(k)
CASE(1)
base_element=x(2)*(x(2)-x(1) )!(1.0-x(1))*(1.00-2.0*x(1))
......@@ -198,8 +205,61 @@ CONTAINS
PRINT*, "B4, numero base ", k
STOP
END SELECT
CASE(7) ! P4
s= x(1)
SELECT CASE(k)
CASE(1)
base_element = (32._dp/3._dp*(s - 1._dp)*(s - 1._dp/2._dp)*(s - 1._dp/4._dp)*(s - 3._dp/4._dp))
CASE(2)
base_element = (32._dp*s*(s - 1._dp/2._dp)*(s - 1._dp/4._dp)*(s - 3._dp/4._dp))/3._dp
CASE(3)
base_element = -(128._dp*s*(s - 1._dp)*(s - 1._dp/2._dp)*(s - 3._dp/4._dp))/3._dp
CASE(4)
base_element = 64._dp*s*(s - 1._dp)*(s - 1._dp/4._dp)*(s - 3._dp/4._dp)
CASE(5)
base_element = -(128._dp*s*(s - 1._dp)*(s - 1._dp/2._dp)*(s - 1._dp/4._dp))/3._dp
CASE default
PRINT*, "P4, numero base ", k
STOP
END SELECT
CASE(13) ! P3 Gauss Lobatto
alpha=0.5_dp-SQRT(5._dp)/10._dp
beta = 1._dp -alpha
SELECT CASE(k)
CASE(1)
base_element = (-x(1)+alpha)*(x(1)-beta)*(x(1)-1.0)/(alpha*beta)
CASE(2)
base_element = x(1)*(x(1)-alpha)*(x(1)-beta)/alpha/beta
CASE(3)
base_element = x(1)*(x(1)-beta)*(x(1)-1.0_dp)/alpha/(-2._dp*alpha+1._dp)/beta
CASE(4)
base_element = x(1)*(x(1)-alpha)*(x(1)-1.0_dp)/alpha/(2._dp*alpha-1._dp)/beta
CASE default
PRINT*, "P3 Gauss Lobatto, numero base ", k
STOP
END SELECT
CASE(14) ! P4 Gauss Lobatto
alpha=0.5_dp-SQRT(21._dp)/14._dp
beta = 1._dp -alpha
SELECT CASE(k)
CASE(1)
base_element = (2._dp*(alpha-x(1))*(x(1)-1._dp)*(x(1)-0.5)*(x(1)-beta))/(-alpha*beta)
CASE(2)
base_element = (2._dp*x(1)*(-alpha + x(1))*(x(1) - 0.5)*(x(1) - beta))/(alpha*beta)
CASE(3)
base_element = (x(1)*(x(1)-1._dp)*(x(1)-0.5_dp)*(x(1)-beta))/(-alpha*(2._dp*alpha-1)*beta*(alpha-0.5_dp))
CASE(4)
base_element = -(4._dp*x(1)*(alpha -x(1))*(x(1)-1._dp)*(x(1)-beta))/(alpha-0.5)**2._dp
CASE(5)
base_element =-(x(1)*(alpha-x(1))*(x(1)-1._dp)*(x(1)-0.5_dp))/(-alpha*(2._dp*alpha-1._dp)*beta*(alpha-0.5))
CASE default
PRINT*, "P4 Gauss Lobatto, numero base ", k
STOP
END SELECT
CASE default
PRINT*, mod_name
PRINT*, "Type non existant", e%itype
STOP
END SELECT
......@@ -276,22 +336,77 @@ CONTAINS
END FUNCTION der_sec_element
FUNCTION der2_poly_element(e,u) RESULT(Uder)
CLASS(element), INTENT(in)::e
REAL(8),DIMENSION(e%nsommets), INTENT(in)::u
REAL(8),DIMENSION(MAX(1,e%nsommets-2)):: coeff
REAL(8),DIMENSION(e%nsommets):: Uder
REAL(8):: u1,u2,u3,u4,x0,x1,x2,x3
SELECT CASE(e%itype)
CASE(1)!P1
coeff(1)=0.
Uder(1:2)=0.
CASE(2)
x0=0.
x1=0.5
x2=1.
u1=u(1)
u2=u(3)
u3=u(2)
coeff(1)=(2*(u1*x1 - u2*x0 - u1*x2 + u3*x0 + u2*x2 - u3*x1))&
&/((x0 - x1)*(x0 - x2)*(x1 - x2))
Uder(1:3)=coeff(1)/e%volume**2
CASE(5)
x0=0.
x1=1./3.
x2=2./3.
x3=1.
u1=u(1)
u2=u(3)
u3=u(4)
u4=u(2)
coeff(1)=(2*(u1*x1*(x2**3) - u1*(x1**3)*x2 - u2*x0*(x2**3) &
&+ u2*(x0**3)*x2 + u3*x0*(x1**3) - u3*(x0**3)*x1 - u1*x1*(x3**3 )&
& + u1*(x1**3)*x3 + u2*x0*(x3**3) - u2*(x0**3)*x3 - u4*x0*(x1**3) &
&+ u4*(x0**3)*x1 + u1*x2*(x3**3) - u1*(x2**3)*x3 - u3*x0*(x3**3) &
&+ u3*(x0**3)*x3 + u4*x0*(x2**3) - u4*(x0**3)*x2 - u2*x2*(x3**3) &
&+ u2*(x2**3)*x3 + u3*x1*(x3**3) - u3*(x1**3)*x3 - u4*x1*(x2**3) + &
& u4*(x1**3)*x2))/((x0 - x1)*(x0 - x2)*(x0 - x3)*(x1 - x2)*(x1 - x3)*(x2 - x3))
coeff(2)= -(6*(u1*x1*(x2**2) - u1*(x1**2)*x2 - u2*x0*(x2**2) + u2*(x0**2)*x2 &
& + u3*x0*(x1**2) - u3*(x0**2)*x1 - u1*x1*(x3**2) + u1*(x1**2)*x3 + u2*x0*(x3**2) &
& - u2*(x0**2)*x3 - u4*x0*(x1**2) + u4*(x0**2)*x1 + u1*x2*(x3**2) - u1*(x2**2)*x3 &
& - u3*x0*(x3**2) + u3*(x0**2)*x3 + u4*x0*(x2**2) - u4*(x0**2)*x2 - u2*x2*(x3**2) &
& + u2*(x2**2)*x3 + u3*x1*(x3**2) - u3*(x1**2)*x3 - u4*x1*(x2**2) + u4*(x1**2)*x2)) &
&/((x0 - x1)*(x0 - x2)*(x0 - x3)*(x1 - x2)*(x1 - x3)*(x2 - x3))
Uder(1)=(coeff(1)+coeff(2)*x0)/e%volume**2
Uder(3)=(coeff(1)+coeff(2)*x1)/e%volume**2
Uder(4)=(coeff(1)+coeff(2)*x2)/e%volume**2
Uder(2)=(coeff(1)+coeff(2)*x3)/e%volume**2
CASE default
PRINT*, 'der2_poly, element 1d'
STOP
END SELECT
END FUNCTION der2_poly_element
FUNCTION gradient_element(e,k,x) RESULT (grad) ! gradient in reference element
CLASS(element), INTENT(in):: e
INTEGER, INTENT(in):: k ! numero de la fonction de base
REAL(dp), DIMENSION(2), INTENT(in):: x ! coordonnees barycentriques
REAL(dp),DIMENSION(n_dim):: grad
REAL(dp):: fx,fy,fz
REAL(dp):: fx,fy,fz, alpha, s
SELECT CASE(e%itype)
CASE(1)! P1
CASE(1,11)! P1
SELECT CASE(k)
CASE(1)
fx=-1.0_dp
CASE(2)
fx= 1.0_dp
END SELECT
CASE(3) !P2
CASE(3,12) !P2
SELECT CASE(k)
CASE(1)
fx=-3.0_dp+4.0_dp*x(1)
......@@ -323,36 +438,88 @@ CONTAINS
CASE(5)
SELECT CASE(k)
CASE(1)
fx=x(1)*( -3.0_dp*x(1)+6.0_dp) -3.0_dp
fx=-3.0_dp*x(2)*x(2)
!fx=-3.0*x(2)*x(2)
CASE(2)
fx=3.0*x(1)*x(1)
fx=3.0_dp*x(1)*x(1)
!fx=3.0*x(1)*x(1)
CASE(3)
fx=(9.0_dp*x(1)-12.0_dp)*x(1)+3.0_dp
fx= 3.0_dp*x(2)*( x(2)-2._dp*x(1) )
!fx=3.0*x(2)*( x(2)-2.*x(1) )
CASE(4)
fx= (-9.0_dp*x(1)+6.0_dp)*x(1)
fx= 3.0_dp*x(1)*(2._dp*x(2)-x(1))
!fx=6*x(1)*x(2)-3*x(1)*x(1)
END SELECT
CASE(6)
CASE(6) !B4
SELECT CASE(k)
CASE(1)
!fx=-4.0+12.0*x(1)-12.0*x(1)*x(1)+4.0*x(1)*x(1)*x(1)
fx=-4.0_dp*x(2)*x(2)*x(2)
CASE(2)
!fx=4.0*x(1)*x(1)*x(1)
fx= 4.0_dp* x(1)*x(1)*x(1)
CASE(3)
!fx=4.0-24.0*x(1)+36.0*x(1)*x(1)-16.0*x(1)*x(1)*x(1)
fx= 4.0_dp*( x(2)-3.0_dp*x(1) )*x(2)*x(2)
CASE(4)
!fx=12.0*x(1)-36.0*x(1)*x(1)+24.0*x(1)*x(1)*x(1)
fx=12.0_dp * x(1)*x(2) * ( x(2)-x(1) )
CASE(5)
!fx=12.0*x(1)*x(1)-16.0*x(1)*x(1)*x(1)
fx=4.0_dp*( 3._dp*x(2)-x(1) )*x(1)*x(1)
END SELECT
CASE(7) !P4
s=x(1)
SELECT CASE(k)
CASE(1)
!fx=-4.0+12.0*x(1)-12.0*x(1)*x(1)+4.0*x(1)*x(1)*x(1)
fx= (((128._dp*s)/3._dp - 80._dp)*s + 140._dp/3._dp)*s - 25._dp/3._dp
CASE(2)
!fx=4.0*x(1)*x(1)*x(1)
fx=(((128._dp*s)/3._dp - 48._dp)*s + (44._dp)/3._dp)*s - 1._dp
CASE(3)
!fx=4.0-24.0*x(1)+36.0*x(1)*x(1)-16.0*x(1)*x(1)*x(1)
fx= ((- (512._dp*s)/3._dp + 288._dp)*s - (416._dp)/3._dp)*s + 16._dp
CASE(4)
!fx=12.0*x(1)-36.0*x(1)*x(1)+24.0*x(1)*x(1)*x(1)
fx=((256._dp*s - 384._dp)*s + 152._dp)*s - 12._dp
CASE(5)
!fx=12.0*x(1)*x(1)-16.0*x(1)*x(1)*x(1)
fx= s*(s*(-(512._dp*s)/3._dp + 224._dp) - 224._dp/3._dp) + 16._dp/3._dp
END SELECT
CASE(13) !P3 Gauss Lobatto
alpha=0.5_dp-SQRT(5._dp)/10._dp
SELECT CASE(k)
CASE(1)
fx=(- alpha**2._dp + alpha + 3._dp*x(1)**2._dp - 4._dp*x(1) + 1._dp)/(alpha*(alpha - 1._dp))
!fx=-3.0*x(2)*x(2)
CASE(2)
fx=-(- alpha**2._dp + alpha + 3._dp*x(1)**2._dp - 2._dp*x(1))/(alpha*(alpha - 1._dp))
!fx=3.0*x(1)*x(1)
CASE(3)
fx=(2._dp*alpha*x(1) - 4._dp*x(1) - alpha + 3._dp*x(1)**2._dp + 1._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
!fx=3.0*x(2)*( x(2)-2.*x(1) )
CASE(4)
fx= -(alpha - 2._dp*x(1) - 2._dp*alpha*x(1) + 3._dp*x(1)**2._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
!fx=6*x(1)*x(2)-3*x(1)*x(1)
END SELECT
CASE(14) !P4 Gauss Lobatto
alpha=0.5_dp-SQRT(21._dp)/14._dp
s=x(1)
SELECT CASE(k)
CASE(1)
fx=(4*alpha**2*s - 3*alpha**2 - 4*alpha*s + 3*alpha - 8*s**3 + 15*s**2 - 8*s + 1)/(alpha*(alpha - 1))
!fx=-3.0*x(2)*x(2)
CASE(2)
fx=-(- 4*alpha**2*s + alpha**2 + 4*alpha*s - alpha + 8*s**3 - 9*s**2 + 2*s)/(alpha*(alpha - 1))
!fx=3.0*x(1)*x(1)
CASE(3)
fx=(alpha + 8*s - 6*alpha*s + 6*alpha*s**2 - 15*s**2 + 8*s**3 - 1)/(alpha*(2*alpha - 1)**2*(alpha - 1))
!fx=3.0*x(2)*( x(2)-2.*x(1) )
CASE(4)
fx= (16*(2*s - 1)*(- alpha**2 + alpha + 2*s**2 - 2*s))/(2*alpha - 1)**2
!fx=6*x(1)*x(2)-3*x(1)*x(1)
CASE(5)
fx = -(alpha - 2*s - 6*alpha*s + 6*alpha*s**2 + 9*s**2 - 8*s**3)/(alpha*(2*alpha - 1)**2*(alpha - 1))
END SELECT
CASE default
PRINT*, "Type non existant", e%itype
STOP
......@@ -365,17 +532,17 @@ CONTAINS
INTEGER, INTENT(in):: k ! numero de la fonction de base
REAL(dp), DIMENSION(2), INTENT(in):: x ! coordonnees barycentriques
REAL(dp),DIMENSION(n_dim):: grad2
REAL(dp):: fxx
REAL(dp):: fxx, alpha,s
SELECT CASE(e%itype)
CASE(1)! P1
CASE(1,11)! P1
SELECT CASE(k)
CASE(1)
fxx=0.0_dp
CASE(2)
fxx=0.0_dp
END SELECT
CASE(3) !P2
CASE(3,12) !P2
SELECT CASE(k)
CASE(1)
fxx=4.0_dp
......@@ -394,15 +561,20 @@ CONTAINS
fxx=-4.0_dp
END SELECT
CASE(4) ! P3
alpha=1._dp/3._dp
SELECT CASE(k)
CASE(1)
fxx=-27._dp*x(1)
!fxx=-6.0*x(1)+6.0
fxx=(6._dp*x(1) - 4._dp)/(alpha*(alpha - 1._dp))
CASE(2)
fxx=27._dp*x(1)
!fxx=6.0*x(1)
fxx=-(6._dp*x(1) - 2._dp)/(alpha*(alpha - 1._dp))
CASE(3)
fxx= 81._dp*x(1)
!fxx=18.0*x(1)-12.0
fxx= (2._dp*alpha + 6._dp*x(1) - 4._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
CASE(4)
fxx=-81._dp*x(1)
!fxx=-18.0*x(1)+6.0
fxx=(2._dp*alpha - 6._dp*x(1) + 2._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
END SELECT
CASE(5)
SELECT CASE(k)
......@@ -430,6 +602,59 @@ CONTAINS
CASE(5)
fxx=24.0_dp*x(1)-48.0_dp*x(1)*x(1)
END SELECT
CASE(7) !P4
s=x(1)
SELECT CASE(k)
CASE(1)
fxx=(128._dp*s - 160._dp)*s + 140._dp/3._dp
CASE(2)
fxx= (128._dp*s - 96._dp)*s + 44._dp/3._dp
CASE(3)
fxx= (-512._dp*s + 576._dp)*s - 416._dp/3._dp
CASE(4)
fxx= 768._dp*s*(s - 1._dp) + 152._dp
CASE(5)
fxx= (- 512._dp*s + 448._dp)*s - 224._dp/3._dp
END SELECT
CASE(13)
alpha=0.5_dp-SQRT(5._dp)/10._dp
SELECT CASE(k)
CASE(1)
!fxx=-6.0*x(1)+6.0
fxx=(6._dp*x(1) - 4._dp)/(alpha*(alpha - 1._dp))
CASE(2)
!fxx=6.0*x(1)
fxx=-(6._dp*x(1) - 2._dp)/(alpha*(alpha - 1._dp))
CASE(3)
!fxx=18.0*x(1)-12.0
fxx= (2._dp*alpha + 6._dp*x(1) - 4._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
CASE(4)
!fxx=-18.0*x(1)+6.0
fxx=(2._dp*alpha - 6._dp*x(1) + 2._dp)/(alpha*(2._dp*alpha**2._dp - 3._dp*alpha + 1._dp))
END SELECT
CASE(14)
alpha=0.5_dp-SQRT(21._dp)/14._dp
s=x(1)
SELECT CASE(k)
CASE(1)
!fxx=-6.0*x(1)+6.0
fxx=-(- 4*alpha**2 + 4*alpha + 24*s**2 - 30*s + 8)/(alpha*(alpha - 1))
CASE(2)
!fxx=6.0*x(1)
fxx=-(- 4*alpha**2 + 4*alpha + 24*s**2 - 18*s + 2)/(alpha*(alpha - 1))
CASE(3)
!fxx=18.0*x(1)-12.0
fxx= (12*alpha*s - 30*s - 6*alpha + 24*s**2 + 8)/(alpha*(2*alpha - 1)**2*(alpha - 1))
CASE(4)
!fxx=-18.0*x(1)+6.0
fxx=(32*(- alpha**2 + alpha + 6*s**2 - 6*s + 1))/(2*alpha - 1)**2
CASE(5)
fxx = (6*alpha - 18*s - 12*alpha*s + 24*s**2 + 2)/(alpha*(2*alpha - 1)**2*(alpha - 1))
END SELECT
CASE default
PRINT*, "Type non existant", e%itype
STOP
......@@ -488,7 +713,7 @@ CONTAINS
REAL(dp),PARAMETER:: s4_9=1._dp/3._dp*SQRT(5._dp+2._dp*SQRT(10.0_dp/7.0_dp))
SELECT CASE(e%itype)
CASE(1)
CASE(1,11)
e%nquad=2
nquad=e%nquad
ALLOCATE(e%quad(2,e%nquad),e%weight(e%nquad))
......@@ -514,42 +739,83 @@ CONTAINS
! s=sqrt(0.6_dp) !0.7745966692414834_dp !SQRT(0.6d0)
! s=sqrt(0.6_dp) !0.7745966692414834_dp !SQRT(0.6_dp)
e%quad(1,1)=s2_3 !0.5_dp*(1.0_dp - s)
e%quad(2,1)=s3_3 !0.5_dp*(1.0_dp + s)
e%weight(1) =weight1_3 !5._dp/18._dp! 0.2777777777777778 !5.0d0/18.0d0
e%weight(1) =weight1_3 !5._dp/18._dp! 0.2777777777777778 !5.0_dp/18.0_dp
e%quad(1,2)=s3_3 !0.5d0*(1.0d0 + s)
e%quad(2,2)=s2_3 !0.5d0*(1.0d0 - s)
e%weight(2) =weight1_3 !5._dp/18._dp! 0.2777777777777778 !5.0d0/18.0d0
e%quad(1,2)=s3_3 !0.5_dp*(1.0_dp + s)
e%quad(2,2)=s2_3 !0.5_dp*(1.0_dp - s)
e%weight(2) =weight1_3 !5._dp/18._dp! 0.2777777777777778 !5.0_dp/18.0_dp
e%quad(1,3)=0.5d0
e%quad(2,3)=0.5d0
e%weight(3)= weight2_3 !8._dp/18._dp!0.4444444444444444 ! 8.0d0/18.0d0
e%quad(1,3)=0.5_dp
e%quad(2,3)=0.5_dp
e%weight(3)= weight2_3 !8._dp/18._dp!0.4444444444444444 ! 8.0_dp/18.0_dp
CASE(12)!! exact for degree 3, underintegration
e%nquad = 3
nquad = e%nquad
ALLOCATE(e%quad(2,e%nquad),e%weight(e%nquad))
e%quad(1,1)=0._dp
e%quad(2,1)=1.0_dp
e%weight(1)=1._dp/6._dp
e%quad(1,2)=0.5_dp
e%quad(2,2)=0.5_dp
e%weight(2)=2._dp/3._dp
e%quad(1,3)=1._dp
e%quad(2,3)=0._dp
e%weight(3)=1._dp/6._dp
CASE(4,5) !! exact for degree 7
e%nquad=4 ! ordre 7
nquad=e%nquad
ALLOCATE(e%quad(2,e%nquad),e%weight(e%nquad) )
! s=(18._dp-sqrt(30._dp))/72._dp !0.1739274225687269_dp !(18.- SQRT(30.))/72. !SQRT(30.d0)*(-5.0d0+3.d0*SQRT(30.d0))/360.d0
! s=(18._dp-sqrt(30._dp))/72._dp !0.1739274225687269_dp !(18.- SQRT(30.))/72. !SQRT(30._dp)*(-5.0_dp+3._dp*SQRT(30._dp))/360._dp
e%weight(1:2)=s1_7
! s=SQRT(525._dp+70._dp*SQRT(30._dp))/35._dp !0.8611363115940526_dp!SQRT(525.d0+70.d0*SQRT(30.d0))/35.d0
! s=SQRT(525._dp+70._dp*SQRT(30._dp))/35._dp !0.8611363115940526_dp!SQRT(525._dp+70._dp*SQRT(30._dp))/35._dp
e%quad(1,1)=s3_7! 0.50_dp*(1.00_dp-s)
e%quad(1,2)=s4_7 !0.50_dp*(1.00_dp+s)