add 1D

Src1D/Model.f90

+!!!  HIGH ORDER IN SPACE AND TIME DEFERRED CORRECTION (EXPLICIT) 
+!!!     RESIDUAL DISTRIBUTION METHOD 
+!!!  DESIGNED FOR THE SYSTEM GIVEN BY THE EULER EQUATIONS in 1D and 2D
+!!!
+!!!  Authors:
+!!!  Remi Abgrall (University of Zurich),
+!!!  Paola Bacigaluppi (University of Zurich),
+!!!  Svetlana Tokareva (University of Zurich)
+!!!  Institute of Mathematics and Institute of Computational Sciences
+!!!  University of Zurich
+!!!  July 10, 2018
+!!!  Correspondance:	remi.abgrall@math.uzh.ch
+!!!  ------------------------------------------
+MODULE Model
+
+  ! This module allows to pass from control to physical variables and viceversa
+
+  USE param2d
+  USE PRECISION
+  IMPLICIT NONE
+CONTAINS
+
+  FUNCTION Control_to_Cons(u,e) RESULT (u_cons)
+    ! crom control point, compute the values a physical dofs
+    TYPE(Pvar), DIMENSION(:), INTENT(in):: u
+    TYPE(element), INTENT(in):: e
+    TYPE(Pvar),DIMENSION(SIZE(u,dim=1)):: u_cons
+    INTEGER:: k,l
+    u_cons=u
+
+    SELECT CASE(e%itype)
+    CASE(1,3,4) ! Lagrange
+
+    CASE(2,5,6) ! Bezier
+
+       DO l=1, e%nsommets
+          u_cons(l) = e%eval_func(u(:),e%x(:,l))
+       ENDDO
+
+    CASE default
+       PRINT*, "erreur dans Model/Control_to_Cons"
+       STOP
+    END SELECT
+  END FUNCTION Control_to_cons
+
+  FUNCTION Cons_to_Control(u_cons,e) RESULT (u)
+    ! from values at dofs, compute control points
+    TYPE(Pvar), DIMENSION(:), INTENT(in):: u_cons
+    TYPE(element), INTENT(in):: e
+    TYPE(Pvar),DIMENSION(SIZE(u_cons,dim=1)):: u
+    INTEGER:: l, k
+
+
+    u=u_cons
+
+    SELECT CASE(e%itype)
+    CASE(1,3,4) ! Lagrange
+
+
+    CASE(2,5,6)! cubic bezier
+
+
+       DO k=1, n_vars
+
+          u(:)%u(k) = MATMUL(e%base1,u_cons%u(k))
+
+
+       ENDDO
+    CASE default
+       PRINT*, "erreur dans Model/Control_to_Cons"
+       STOP
+    END SELECT
+    !enddo
+  END FUNCTION Cons_to_Control
+END MODULE Model

Src1D/algebra.f90

+!!!  HIGH ORDER IN SPACE AND TIME DEFERRED CORRECTION (EXPLICIT) 
+!!!     RESIDUAL DISTRIBUTION METHOD 
+!!!  DESIGNED FOR THE SYSTEM GIVEN BY THE EULER EQUATIONS in 1D and 2D
+!!!
+!!!  Authors:
+!!!  Remi Abgrall (University of Zurich),
+!!!  Paola Bacigaluppi (University of Zurich),
+!!!  Svetlana Tokareva (University of Zurich)
+!!!  Institute of Mathematics and Institute of Computational Sciences
+!!!  University of Zurich
+!!!  July 10, 2018
+!!!  Correspondance:	remi.abgrall@math.uzh.ch
+!!!  ------------------------------------------
+MODULE algebra
+  USE PRECISION
+  IMPLICIT NONE
+
+  LOGICAL, SAVE, PRIVATE :: initialise = .FALSE.
+  LOGICAL, SAVE, PRIVATE :: module_debug = .FALSE.
+
+  LOGICAL, SAVE, PUBLIC               :: optim1_plain = .FALSE.
+  LOGICAL, SAVE, PUBLIC               :: optim2_plain = .FALSE.
+  LOGICAL, SAVE, PUBLIC               :: optim3_plain = .FALSE.
+  LOGICAL, SAVE, PUBLIC               :: optim4_plain = .FALSE.
+  LOGICAL, SAVE, PUBLIC               :: optim5_plain = .FALSE.
+CONTAINS
+
+  !------------------------------------------
+  !!--  Inversion d'une matrice carree par LU et pivot partiel par ligne
+  !!--  Ref : Numerical recipes in C
+  !------------------------------------------
+  FUNCTION Inverse( Mat)  RESULT (Mat1)
+    CHARACTER(LEN = *), PARAMETER :: mod_name = "InverseLu"
+    REAL(DP), DIMENSION(:, : ), INTENT(IN) :: Mat
+    REAL(DP), DIMENSION(SIZE(Mat, dim = 1), SIZE( Mat, dim = 1)) :: Mat1, mat2
+    REAL(DP), DIMENSION(SIZE(Mat, dim = 1))   :: col
+    REAL(DP)                        :: d
+    INTEGER    :: j, Nsize, l, zz
+    INTEGER, DIMENSION(SIZE(Mat, dim = 1)) :: indx !-- CCE R.B. 2008/10/21
+
+
+
+    IF (.NOT. initialise) THEN
+       optim1_plain = .FALSE.
+       optim2_plain = .TRUE.
+       optim3_plain = .TRUE. !-- celle là est particulièrement efficace
+       optim4_plain = .TRUE.
+
+       optim2_plain = .FALSE. !-- ces 2 optimisations là mettre MHD en l'air, et pourtant elles sont meilleures en précusion, ce qui ne saute pas aux yeux au niveau des résidus %%%%%%
+       optim4_plain = .FALSE.
+
+       optim3_plain = .FALSE. !-- cette optimisation met MHD en l'air uniquement en -O5, pas en -O2
+    END IF
+
+    Nsize = SIZE(Mat, dim=1)
+    mat2  = Mat
+
+    CALL ludcmp(mat2, Nsize, indx, d) !-- indx OUT
+
+    DO j = 1, Nsize
+       col = 0.0
+       col(j) = 1.0
+       CALL luksb(mat2, Nsize, indx, col) !-- col IN OUT; indx IN
+       Mat1(:, j) = col(: )
+    END DO
+
+
+  CONTAINS
+    SUBROUTINE ludcmp(mat3, Nsize, indx, d)
+      CHARACTER(LEN = *), PARAMETER :: mod_name = "ludcmp"
+
+      INTEGER, INTENT(IN) :: Nsize
+      REAL(DP), DIMENSION(Nsize, Nsize ), INTENT(INOUT) :: mat3 !-- Passer en :,: pose des pbs avec le -O5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      INTEGER, DIMENSION(Nsize ), INTENT(OUT) :: indx !-- Passer en :,: pose des pbs avec le -O5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      REAL(DP), INTENT(OUT) :: d
+
+      REAL(DP), DIMENSION(Nsize) :: vv
+      INTEGER :: i, j, k, imax, zz
+      REAL(DP) :: Big, dum, somme, temp
+
+      d = 1.0_dp
+      DO i = 1, Nsize
+         IF (optim1_plain) THEN
+            Big = MAXVAL(ABS(mat3(i, :)))
+         ELSE
+            Big = 0.0_dp
+            DO j = 1, Nsize
+               temp = ABS(mat3(i, j))
+               IF (temp > Big) THEN
+                  Big = temp
+               END IF
+            END DO
+         END IF
+         IF (Big == 0.0_dp) THEN
+            WRITE( *, *) mod_name, " ERREUR : Matrice singuliere, i==", i
+            STOP
+         END IF
+         vv(i) = 1.0 / Big
+      END DO
+
+      DO j = 1, Nsize
+
+         DO i = 1, j -1
+            IF (optim2_plain) THEN
+               mat3(i, j) = mat3(i, j) - SUM(mat3(i, 1: i - 1) * mat3(1: i - 1, j))
+            ELSE
+               somme = mat3(i, j)
+               DO k = 1, i -1
+                  somme = somme - mat3(i, k) * mat3(k, j)
+               END DO
+               mat3(i, j) = somme
+            END IF
+         END DO !-- boucle sur i
+
+
+         Big = 0.0_dp
+         imax = -1
+         DO i = j, Nsize
+            somme = mat3(i, j)
+            DO k = 1, j -1
+               somme = somme - mat3(i, k) * mat3(k, j)
+
+            END DO
+
+            mat3(i, j) = somme
+            dum = vv(i) * ABS(somme)
+            IF (dum >= Big) THEN
+               Big = dum
+               imax = i
+            END IF
+         END DO !-- boucle sur i
+
+         IF (j /= imax) THEN
+            DO k = 1, Nsize
+               dum = mat3(imax, k)
+               mat3(imax, k) = mat3(j, k)
+               mat3(j, k) = dum
+            END DO !-- boucle sur k
+            d = - d
+            vv(imax ) = vv(j)
+         END IF
+
+
+
+         indx(j) = imax
+
+         IF (ABS(mat3(j, j)) <= 1.0e-20_dp) THEN
+            mat3(j, j) = SIGN(1.0e-20_dp, mat3(j, j)) !-- CCE 2007/04/24 !-- CCE 2008/12/17
+         END IF
+
+         IF (j /= Nsize) THEN
+
+            DO i = j + 1, Nsize
+               mat3(i, j) = mat3(i, j) / mat3(j, j)
+            END DO !-- boucle sur i
+         END IF
+
+
+
+      END DO !-- boucle sur j
+
+    END SUBROUTINE ludcmp
+
+
+    SUBROUTINE luksb(mat2, Nsize, indx, col)
+      CHARACTER(LEN = *), PARAMETER :: mod_name = "luksb"
+
+      INTEGER, INTENT(IN) :: Nsize
+      REAL(DP),  DIMENSION(Nsize, Nsize), INTENT(IN) :: mat2 !-- Passer en :,: pose des pbs avec le -O5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      INTEGER,  DIMENSION(Nsize), INTENT(IN) :: indx !-- Passer en :,: pose des pbs avec le -O5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      REAL(DP), DIMENSION(Nsize), INTENT(INOUT) :: col !-- Passer en :,: pose des pbs avec le -O5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+      INTEGER ::  i, ii, ip, j
+      REAL(DP)  :: somme
+
+      ii = 1
+      DO i = 1, Nsize
+         ip = indx(i)
+         IF (ip <= 0 .OR. ip > SIZE(col)) THEN
+            PRINT *, mod_name, " ERREUR : indx(", i, ") invalide", ip
+            STOP
+         END IF
+         somme = col(ip)
+         col(ip) = col(i)
+         IF (ii > 0) THEN
+            IF (optim3_plain) THEN
+               somme = somme - SUM(mat2(i, ii: i -1) * col(ii: i -1))
+            ELSE
+               DO j = ii, i -1
+                  somme = somme - mat2(i, j) * col(j)
+               END DO
+            END IF
+         ELSE
+            IF (somme == 0.0_dp) THEN
+               ii = i
+            END IF
+         END IF
+
+         col(i) = somme
+      END DO
+
+      DO i = Nsize, 1, -1
+         IF (optim4_plain) THEN
+            col(i) = (col(i) - SUM(mat2(i, i + 1: Nsize) * col(i + 1: Nsize))) / mat2(i, i)
+         ELSE
+            somme = col(i)
+            DO j = i + 1, Nsize
+               somme = somme - mat2(i, j) * col(j)
+            END DO
+            col(i) = somme / mat2(i, i)
+         END IF
+      END DO
+    END SUBROUTINE  luksb
+
+  END FUNCTION Inverse
+END MODULE algebra

Src1D/aretes.f90

+!!!  HIGH ORDER IN SPACE AND TIME DEFERRED CORRECTION (EXPLICIT) 
+!!!     RESIDUAL DISTRIBUTION METHOD 
+!!!  DESIGNED FOR THE SYSTEM GIVEN BY THE EULER EQUATIONS in 1D and 2D
+!!!
+!!!  Authors:
+!!!  Remi Abgrall (University of Zurich),
+!!!  Paola Bacigaluppi (University of Zurich),
+!!!  Svetlana Tokareva (University of Zurich)
+!!!  Institute of Mathematics and Institute of Computational Sciences
+!!!  University of Zurich
+!!!  July 10, 2018
+!!!  Correspondance:	remi.abgrall@math.uzh.ch
+!!!  ------------------------------------------
+MODULE arete_class
+  USE PRECISION
+  IMPLICIT NONE
+  TYPE, PUBLIC:: arete
+     INTEGER:: nsommets, itype, nvertex! nbre de dofs, type element: 1-> P1
+     ! 2-> B2
+     !3->P2
+     ! nombre de sommet dans cet element arete
+     LOGICAL:: bord
+     INTEGER:: jt1=-1, jt2 =-1 ! les deux elements de par et d'autre de l'arete.
+     INTEGER, DIMENSION(:,:), POINTER :: nu =>Null() ! nu( indice des elements, indice des voisins): on cherche a connaitre le numero local des 
+     ! points communs (puisque le maillage est confome)  aux deux elements sur cette face dans chacun des elements jt1 et jt2
+     REAL(dp), DIMENSION(:,:), POINTER   :: coor=>Null() ! il s'agit des coordonnees physique des dofs (communs) sur la face (commune)
+     REAL(dp)                           :: volume=0  !
+     REAL(dp)                           :: jump_flag=1.0 !(between 0 and 1 which gives the weight of the edge term
+     ! for this one check if there is a discontinuity around and take the maximum value)
+     REAL(dp),  DIMENSION(:), POINTER :: n   =>Null()   ! normales exterieures
+     INTEGER                        :: log    ! logique
+!!!!   quadrature de surface !
+     REAL(dp),   DIMENSION(:,:),POINTER :: quad =>Null()  ! point de quadrature 
+     REAL(dp),     DIMENSION(:),POINTER :: weight=>Null() ! poids 
+     INTEGER                        :: nquad  ! nbre de points de quadrature 
+!!! quadrature bord (dimension -1)
+     REAL(dp),   DIMENSION(:,:),POINTER :: quad_1 =>Null()  ! point de quadrature 
+     REAL(dp),     DIMENSION(:),POINTER :: weight_1=>Null() ! poids 
+     INTEGER                        :: nquad_1  ! nbre de points de quadrature 
+!!!
+   CONTAINS
+     PROCEDURE, PUBLIC:: aire=>aire_arete
+     PROCEDURE, PUBLIC:: quadrature=>quadrature_arete
+     PROCEDURE, PUBLIC:: normale=>normale_arete
+     !FINAL:: clean_arete
+
+  END TYPE arete
+
+
+CONTAINS
+
+  REAL(dp) FUNCTION aire_arete(e)
+    CLASS(arete), INTENT(in):: e
+    REAL(dp), DIMENSION(2):: a
+    a= e%coor(:,2)-e%coor(:,1)
+
+    aire_arete=SQRT(a(1)**2+ a(2)**2 )
+  END FUNCTION aire_arete
+
+  FUNCTION normale_arete(e) RESULT(n)
+    CLASS(arete), INTENT(in)::e
+    REAL(dp), DIMENSION(2):: n
+    INTEGER:: l, k1, k2
+    INTEGER, DIMENSION(3), PARAMETER:: ip1=(/2,3,1/)
+
+    n(1)=e%coor(2,2)-e%coor(2,1)
+    n(2)=e%coor(1,1)-e%coor(1,2)
+
+  END FUNCTION normale_arete
+
+
+
+  SUBROUTINE quadrature_arete(e)
+    CLASS(arete), INTENT(inout):: e
+    REAL(dp):: w,zo,xo,s
+    INTEGER:: nquad
+
+    !Gaussia formula, exact for polynomials of degree 5
+    PRINT*, "quadrature_arete"
+    STOP
+
+    e%nquad=3
+    nquad=e%nquad
+    ALLOCATE(e%quad(2,e%nquad),e%weight(e%nquad))
+
+    s=SQRT(0.6d0)
+    e%quad(1,1)=0.5d0*(1.0d0 - s)
+    e%quad(2,1)=0.5d0*(1.0d0 + s)
+    e%weight(1) = 5.0d0/18.0d0
+    e%quad(1,2)=0.5d0*(1.0d0 + s)
+    e%quad(2,2)=0.5d0*(1.0d0 - s)
+    e%weight(2) = 5.0d0/18.0d0
+    e%quad(1,3)=0.5d0
+    e%quad(2,3)=0.5d0
+    e%weight(3)=8.0d0/18.0d0
+
+
+
+  END SUBROUTINE quadrature_arete
+
+  SUBROUTINE clean_arete(e)
+    TYPE(arete):: e
+    IF (ASSOCIATED(e%nu)) NULLIFY(e%nu)
+    IF (ASSOCIATED(e%coor)) NULLIFY(e%coor)
+    IF (ASSOCIATED(e%quad)) NULLIFY(e%quad)
+    IF (ASSOCIATED(e%weight)) NULLIFY(e%weight)
+  END SUBROUTINE clean_arete
+END MODULE arete_class

Src1D/elements_1D.f90

+!!!  HIGH ORDER IN SPACE AND TIME DEFERRED CORRECTION (EXPLICIT) 
+!!!     RESIDUAL DISTRIBUTION METHOD 
+!!!  DESIGNED FOR THE SYSTEM GIVEN BY THE EULER EQUATIONS in 1D and 2D
+!!!
+!!!  Authors:
+!!!  Remi Abgrall (University of Zurich),
+!!!  Paola Bacigaluppi (University of Zurich),
+!!!  Svetlana Tokareva (University of Zurich)
+!!!  Institute of Mathematics and Institute of Computational Sciences
+!!!  University of Zurich
+!!!  July 10, 2018
+!!!  Correspondance:	remi.abgrall@math.uzh.ch
+!!!  ------------------------------------------
+MODULE element_class
+  ! In this module are listed all the tools linked to the basis functions and quadrature formulas
+  ! Structure of the module:
+  ! ELEMENT CLASS
+  ! - aire: area of the element
+  ! - normale: normal of the element
+  ! - base: basis function -- here we do have Lagrangian and Bernstein polynomials
+  ! - eval_function: evaluation of the solution via basis functions SUM(base(:)*(u(:)%u(l))
+  ! -  eval_der: evaluation of the derivative SUM( (u(:)%u(l) *grad(1,:)) 
+  ! -  eval_der2: evaluation of the second derivative SUM( (u(:)%u(l) *grad2(1,:))
+  ! - der_sec: alternative to eval_der2
+  ! - gradient: corresponds to the gradient of the basis function
+  ! - gradient2: corresponds to the second order gradient of the basis function
+  ! - quadrature:  collects all the points and weights for each typology of element
+  ! - base_ref:  computes the values of the basis functions at the physical dofs
+  ! - bary: defines the barycentric coordinates of the points
+  ! - clean: cleans the memory of the features defined within the ELEMENT
+
+  USE algebra
+  USE variable_def
+  USE overloading
+  USE PRECISION
+  IMPLICIT NONE
+
+
+  TYPE, PUBLIC:: element
+     INTEGER:: type_flux=5
+     INTEGER:: diag=-1
+     INTEGER:: nsommets, itype, nvertex ! nbre de dofs, type element:
+     !1->P1
+     !2->B2
+     !3->P2
+     !4->P3
+     !5->B2
+     REAL(dp),  DIMENSION(:), POINTER :: coor =>NULL() ! 
+     !   *-----*-----* for quadratic, *---*---*---* for cubic elements
+     !   1     3     2                1   3   4   2
+     REAL(dp), DIMENSION(:,:), POINTER:: x=> NULL()
+     INTEGER, DIMENSION(:), POINTER   :: nu =>NULL() ! local connectivity table, see above for location
+     ! For Bezier, this corresponds to the Greville points
+     REAL(dp)                           :: volume =0.0_dp  ! volume
+     REAL(dp),  DIMENSION(:), POINTER   :: n =>NULL()     ! external normals 
+     INTEGER                        :: log    ! logic element  : this for boundary conditions, if needed
+!!!!   quadrature de surface
+     REAL(dp),   DIMENSION(:,:),POINTER :: quad =>NULL()  ! point de quadrature 
+     REAL(dp),   DIMENSION(:)  ,POINTER :: weight=>NULL() ! poids 
+     INTEGER                            :: nquad=0  ! nbre de points de quadrature
+     REAL(dp),DIMENSION(:,:),POINTER:: base0=>NULL(),base1=>NULL()
+     INTEGER, DIMENSION(:), POINTER :: dof2ind
+!!! int nabla phi_sigma * phi_sigma'
+     REAL(dp), DIMENSION(:,:), POINTER:: coeff,coeff_b=> NULL(), mass=> NULL()
+   CONTAINS
+     PRIVATE
+     PROCEDURE, PUBLIC:: aire=>aire_element
+     PROCEDURE, PUBLIC:: gradient=>gradient_element
+     PROCEDURE, PUBLIC:: gradient2=>gradient2_element
+     PROCEDURE, PUBLIC:: average => average_element
+     PROCEDURE, PUBLIC:: base=>base_element
+     PROCEDURE, PUBLIC:: normale=>normale_element
+     PROCEDURE, PUBLIC:: quadrature=>quadrature_element
+     PROCEDURE, PUBLIC:: base_ref=>base_ref_element
+     PROCEDURE, PUBLIC:: eval_func=>eval_func_element
+     PROCEDURE, PUBLIC:: eval_der=>eval_der_element
+     PROCEDURE, PUBLIC:: eval_der2=>eval_der2_element
+     PROCEDURE, PUBLIC:: der_sec=>der_sec_element
+     PROCEDURE, PUBLIC:: eval_coeff=>eval_coeff_element
+     FINAL:: clean
+  END TYPE element
+
+CONTAINS
+
+  REAL(dp) FUNCTION aire_element(e) ! area of a element
+    CLASS(element), INTENT(in):: e
+    REAL(dp), DIMENSION(1):: a,b
+    a= e%coor(2)-e%coor(1)
+    aire_element=ABS ( a(1))
+  END FUNCTION aire_element
+
+  FUNCTION normale_element(e) RESULT(n) ! inward normals
+    CLASS(element), INTENT(in)::e
+    REAL(dp), DIMENSION(2):: n
+    INTEGER:: l
+    DO l=1,2
+       n(1)= 1.0_dp ! at e%nu(1), i.e. x_{i+1}
+       n(2)=-1.0_dp ! at e%nu(2), i.e. x_{i}
+    ENDDO
+  END FUNCTION normale_element
+
+  REAL(dp) FUNCTION base_element(e,k,x) ! basis functions
+    CHARACTER(Len=*), PARAMETER :: mod_name="base_element"
+    CLASS(element), INTENT(in):: e
+    INTEGER, INTENT(in):: k ! index of basis function
+    REAL(dp), DIMENSION(2), INTENT(in):: x ! barycentric coordinate
+
+    SELECT CASE(e%itype)
+    CASE(1) ! P1
+       SELECT CASE(k)
+       CASE(2)
+          base_element=x(1)
+       CASE(1)
+          base_element=x(2)
+       CASE default
+          PRINT*,  mod_name
+
+          PRINT*, "P1, numero base ", k
+          STOP
+       END SELECT
+
+    CASE(2) ! B2 Bezier
+       SELECT CASE(k)
+       CASE(1)
+          base_element=x(2)*x(2)!(1.0-x(1))*(1.0-x(1))
+       CASE(2)
+          base_element=x(1)*x(1)
+       CASE(3)
+          base_element=2.0_dp*x(1)*x(2)!(1.0-x(1))
+       CASE default
+          PRINT*,  mod_name
+
+          PRINT*, "B2, numero base ", k
+          STOP
+       END SELECT
+
+    CASE(3)! P2
+       SELECT CASE(k)
+       CASE(1)
+          base_element=x(2)*(x(2)-x(1) )!(1.0-x(1))*(1.00-2.0*x(1))
+       CASE(2)
+          base_element=x(1)*(x(1)-x(2) )!x(1)*(2.0*x(1)-1.00)
+       CASE(3)
+          base_element=4._dp*x(1)*x(2) !4.00*x(1)*(1.00-x(1))
+       CASE default
+          PRINT*,  mod_name
+
+          PRINT*, "P2, numero base ", k
+          STOP
+       END SELECT
+
+    CASE(4) ! P3
+       SELECT CASE(k)
+       CASE(1)
+          base_element = -0.5_dp*(3._dp*x(1)-1._dp)*(3.*x(1)-2._dp)*(x(1)-1._dp)
+       CASE(2)
+          base_element = 0.5_dp*x(1)*(3._dp*x(1)-1._dp)*(3._dp*x(1)-2._dp)
+       CASE(3)
+          base_element = 1.5_dp*x(1)*(3._dp*x(1)-2._dp)*(3._dp*x(1)-3._dp)
+       CASE(4)
+          base_element = -1.5_dp*x(1)*(3._dp*x(1)-1._dp)*(3._dp*x(1)-3._dp)
+       CASE default
+          PRINT*,  mod_name
+
+          PRINT*, "P3, numero base ", k
+          STOP
+       END SELECT
+
+    CASE(5) ! B3
+       SELECT CASE(k)
+       CASE(1)
+          base_element =x(2)**3
+       CASE(2)
+          base_element = x(1)**3
+       CASE(3)
+          base_element = 3.0_dp*x(1)*x(2)*x(2)
+       CASE(4)
+          base_element = 3.0_dp*x(1)*x(1)*x(2)
+       CASE default
+          PRINT*,  mod_name
+
+          PRINT*, "P3, numero base ", k
+          STOP
+       END SELECT
+    CASE(6) ! B4
+       SELECT CASE(k)
+       CASE(1)
+          base_element = x(2)**