Skip to content
GitLab
Commit f9d52e4d authored by Davide Torlo's avatar Davide Torlo
Browse files

Updated 1D code. Models available: linear transport, lineal transport damped,... 

Updated 1D code. Models available: linear transport, lineal transport damped, Burgers, SW, Euler, waves.  There are also cubature elements (Lagrange in gauss lobatto points) There is also a source that must be defined in variable_def, and utils can help defining it. See SW for an example. Tested convergence on lin Scalar, Burges,Euler, more or less correct. Scheme 5 for shocks is working on Sod Euler, Dam Break (wet) SW.
parent 80c342c5
Hide whitespace changes
Inline Side-by-side
Make/Makefile.2D.gfortran
View file @ f9d52e4d
......@@ -4,10 +4,10 @@
# VERSION DEBBUGGEUR -db
# VERSION OPTIMISEUR
F90=gfortran
OBJDIR = obj2D.2
MODDIR = mod2D.2
BINDIR = bin2D.2
SRC = Src2D2.2
OBJDIR = obj2D
MODDIR = mod2D
BINDIR = bin2D
SRC = Src2D2
OPT=-O3
LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib
......
ReadMe.aux 0 → 100644
View file @ f9d52e4d
\relax
\@writefile{toc}{\contentsline {section}{\numberline {1}Compilation}{1}}
\@writefile{toc}{\contentsline {section}{\numberline {2}Meshes}{1}}
\@writefile{toc}{\contentsline {section}{\numberline {3}Run}{1}}
\@writefile{toc}{\contentsline {section}{\numberline {4}Reference:}{1}}
ReadMe.log 0 → 100644
View file @ f9d52e4d
This is pdfTeX, Version 3.14159265-2.6-1.40.19 (TeX Live 2018) (preloaded format=pdflatex 2020.6.26) 19 MAR 2021 15:46
entering extended mode
restricted \write18 enabled.
%&-line parsing enabled.
**ReadMe.tex
(./ReadMe.tex
LaTeX2e <2018-04-01> patch level 5
(/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/tex/latex/base/arti
cle.cls
Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
(/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/tex/latex/base/size
10.clo
File: size10.clo 2014/09/29 v1.4h Standard LaTeX file (size option)
)
\c@part=\count80
\c@section=\count81
\c@subsection=\count82
\c@subsubsection=\count83
\c@paragraph=\count84
\c@subparagraph=\count85
\c@figure=\count86
\c@table=\count87
\abovecaptionskip=\skip41
\belowcaptionskip=\skip42
\bibindent=\dimen102
)
(/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/tex/latex/preprint/
fullpage.sty
Package: fullpage 1999/02/23 1.1 (PWD)
\FP@margin=\skip43
)
(/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/tex/latex/url/url.s
ty
\Urlmuskip=\muskip10
Package: url 2013/09/16 ver 3.4 Verb mode for urls, etc.
) (./ReadMe.aux)
\openout1 = `ReadMe.aux'.
LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 3.
LaTeX Font Info: ... okay on input line 3.
LaTeX Font Info: External font `cmex10' loaded for size
(Font) <12> on input line 9.
LaTeX Font Info: External font `cmex10' loaded for size
(Font) <8> on input line 9.
LaTeX Font Info: External font `cmex10' loaded for size
(Font) <6> on input line 9.
LaTeX Font Info: Try loading font information for OMS+cmr on input line 32.
(/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/tex/latex/base/omsc
mr.fd
File: omscmr.fd 2014/09/29 v2.5h Standard LaTeX font definitions
)
LaTeX Font Info: Font shape `OMS/cmr/m/n' in size <10> not available
(Font) Font shape `OMS/cmsy/m/n' tried instead on input line 32.
LaTeX Font Info: External font `cmex10' loaded for size
(Font) <7> on input line 32.
LaTeX Font Info: External font `cmex10' loaded for size
(Font) <5> on input line 32.
[1
{/usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-var/fonts/map/pdftex/upd
map/pdftex.map}]
Overfull \hbox (7.85497pt too wide) in paragraph at lines 39--40
\OT1/cmr/m/n/10 ods in Ap-plied Math-e-mat-ics, v18(3), pp 327-350, 2018, doi:$
\OT1/cmtt/m/n/10 https : / / doi . org / 10 . 1515 / cmam-[]2017-[]0056$
[]
[2] (./ReadMe.aux) )
Here is how much of TeX's memory you used:
345 strings out of 492639
4380 string characters out of 6117887
64287 words of memory out of 5000000
4297 multiletter control sequences out of 15000+600000
7448 words of font info for 27 fonts, out of 8000000 for 9000
1141 hyphenation exceptions out of 8191
23i,6n,19p,378b,263s stack positions out of 5000i,500n,10000p,200000b,80000s
</usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist
/fonts/type1/public/amsfonts/cm/cmbx12.pfb></usr/sepp1.5.1/drwho/pack-a/texlive
-20181114-cr/texmf-dist/fonts/type1/public/amsfonts/cm/cmr10.pfb></usr/sepp1.5.
1/drwho/pack-a/texlive-20181114-cr/texmf-dist/fonts/type1/public/amsfonts/cm/cm
r12.pfb></usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/texmf-dist/fonts/type1
/public/amsfonts/cm/cmr17.pfb></usr/sepp1.5.1/drwho/pack-a/texlive-20181114-cr/
texmf-dist/fonts/type1/public/amsfonts/cm/cmsy10.pfb></usr/sepp1.5.1/drwho/pack
-a/texlive-20181114-cr/texmf-dist/fonts/type1/public/amsfonts/cm/cmtt10.pfb>
Output written on ReadMe.pdf (2 pages, 87432 bytes).
PDF statistics:
35 PDF objects out of 1000 (max. 8388607)
24 compressed objects within 1 object stream
0 named destinations out of 1000 (max. 500000)
1 words of extra memory for PDF output out of 10000 (max. 10000000)
ReadMe.pdf 0 → 100644
View file @ f9d52e4d
File added
ReadMe.tex
View file @ f9d52e4d
......@@ -16,7 +16,7 @@ Makefile in Make/Makefile\_2D.gfortran. Uses gfortran
\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
......@@ -25,7 +25,7 @@ Note:LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib for mac
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}
......
ReadMe.tex.bak 0 → 100644
View file @ f9d52e4d
\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}
Makefile in Make/Makefile\_2D.gfortran. Uses gfortran
\begin{verbatim}
\make -f Make//Makefile_2D.gfortran dec
\end{verbatim}
\begin{verbatim}
\make -f Make/Makefile_2D.gfortran clean
\end{verbatim}
Note:LIBS=-L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib for mac
\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
Src1D/Model.f90
View file @ f9d52e4d
......@@ -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
......
Src1D/elements_1D.f90
View file @ f9d52e4d
......@@ -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)