Commit 52a6e305 authored by lmical's avatar lmical
Browse files

Revert "I've read, commented, cleaned and (to a certain extent) corrected the...

Revert "I've read, commented, cleaned and (to a certain extent) corrected the 1D code. All the changes are in the README present in Src1D and we can discuss if integrating them or not in the main branch. Most of the changes (epsecially the critical ones) have been discussed with Davide. In particular I added ._DP to all the assignments involving real(DP) variables, I removed the DP at the exponents when we elevate sth to the power of some integer, this is the only case where one should not put DP (for example x**2, here it is better not to put DP). mood was totally in SP (it wasn't simply missing ._DP, the variables were al declaread as REAL(SP)). I applied the BC before and after the udating in fin() (otherwise we make a mistake when we want to impose strong Dirichlet, in the README it is explained better). I think we should import the corrections in the main branch, maybe not all the comments even if I think that this version may be very useful to the ones who are starting to work on the code for the first time. I did it also for the 1d code but in a less cleaned way so it is better not to upload it, or at least not now. Best :)"

This reverts commit f3d858e2.

I uploaded my changes in the main XD
parent f3d858e2
......@@ -11,12 +11,6 @@
!!! July 10, 2018
!!! Correspondance: remi.abgrall@math.uzh.ch
!!! ------------------------------------------
!*----------------------------------------------------------------------------------------
!*Module to pass from coefficients to values and viceversa
!*contr=coefficient
!*cons=values (cons stays for conserved variables)
!*----------------------------------------------------------------------------------------
MODULE Model
! This module allows to pass from control to physical variables and viceversa
......@@ -26,33 +20,21 @@ MODULE Model
IMPLICIT NONE
CONTAINS
!*--------------------------------------------------------------------------------------
!*coeff->values
!*--------------------------------------------------------------------------------------
FUNCTION Control_to_Cons(u,e) RESULT (u_cons)
! crom control point, compute the values a physical dofs
TYPE(Pvar), DIMENSION(:), INTENT(in):: u !*coeff in input
TYPE(Pvar), DIMENSION(:), INTENT(in):: u
TYPE(element), INTENT(in):: e
TYPE(Pvar),DIMENSION(SIZE(u,dim=1)):: u_cons !*values as output
TYPE(Pvar),DIMENSION(SIZE(u,dim=1)):: u_cons
INTEGER:: k,l
u_cons=u !*Set u_cons=u
u_cons=u
SELECT CASE(e%itype)
CASE(1,3,4,7,11,12,13,14) ! Lagrange
!*In this case there is nothing to do because values=coeff
CASE(2,5,6) ! Bezier
!*In this case we need to do a projection
DO l=1, e%nsommets !*Loop on the DoFs
u_cons(l) = e%eval_func(u(:),e%x(:,l)) !*He's storing in u_cons the pointwise evaluations of the solution in the DoFs
!*eval_func_element(u,y) takes in input the coefficients of the solution (or a general function) in the DoFs and the barycentric coordinates y and evaluates the solution in y
!*u->coefficients in the different nodes
!*y->barycentric coordinates of the point
!*e%x are the barycentric coordinates of the DoFs
!*2 indices x(:,:)
!*Second one referred to the DoF of which we consider the barycentric coordinates 1:nsommets
!*First one referred to the 2 barycentric coordinates 1:2
DO l=1, e%nsommets
u_cons(l) = e%eval_func(u(:),e%x(:,l))
ENDDO
CASE default
......@@ -61,49 +43,33 @@ CONTAINS
END SELECT
END FUNCTION Control_to_cons
!*--------------------------------------------------------------------------------------
!*values->coeff
!*--------------------------------------------------------------------------------------
FUNCTION Cons_to_Control(u_cons,e) RESULT (u)
! from values at dofs, compute control points
TYPE(Pvar), DIMENSION(:), INTENT(in):: u_cons !*values in input
TYPE(Pvar), DIMENSION(:), INTENT(in):: u_cons
TYPE(element), INTENT(in):: e
TYPE(Pvar),DIMENSION(SIZE(u_cons,dim=1)):: u !*coefficients as output
TYPE(Pvar),DIMENSION(SIZE(u_cons,dim=1)):: u
INTEGER:: l, k
u=u_cons !*He's copying the values in the output variable
u=u_cons
SELECT CASE(e%itype)
CASE(1,3,4,7,11,12,13,14) ! Lagrange
!*In this case there is nothing to do because coeff=values
CASE(2,5,6)! cubic bezier
!*Projection needed
DO k=1, n_vars !*Loop on the components of the vector u
DO k=1, n_vars
u(:)%u(k) = MATMUL(e%base1,u_cons%u(k))
!*NB: u_cons%u(k) has the component k in all the DoFs, we mean u_cons(:)%u(k), it's the same
!*RMK, from the fields of element
!*base0 matrix of the values of the basis functions at the physical DoFs
!*base1 inverse of base0
!*base0:
!*First index DoF
!*Second index basis function
!*e%base0=[ phi_1(x_1) phi_2(x_1) ... phi_N(x_1)
!* phi_1(x_2) phi_2(x_2) ... phi_N(x_2)
!* . .
!* . .
!* . .
!* phi_1(x_N) phi_2(x_N) ... phi_N(x_n) ]
!*
!*NB: e%base0*vectorofthecoefficient=vectorofthevalues
!*vectorofthecoefficients=inv(e%base0)*vectorofthevalues=e%base1*vectorofthevalues
ENDDO
CASE default
PRINT*, "erreur dans Model/Control_to_Cons"
STOP
END SELECT
!enddo
END FUNCTION Cons_to_Control
END MODULE Model
The final time in DATA is not used, it is in init_bc
also, not all the mentioned schemes are present
MENTIONED SCHEMES
scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
ACTUAL SCHEMES
-1 lxf
4 galerkin+jump
5 psi+jump
->main_dec.f90
!*TYPE(Pvar), DIMENSION(:,:), ALLOCATABLE:: resJ !*NOT USED
!*REAL(dp), DIMENSION(3):: x=0._dp !*NOT USED
!*INTEGER:: nb_args, n, lp, liter !*NOT USED
INTEGER:: Impre !*AS IF IT WAS NOT USED
!*TYPE(Pvar),DIMENSION(1):: FL,FR !*NOT USED
!*INTEGER, DIMENSION(:), ALLOCATABLE :: ShIndArray!*NOT USED
Var%Ncells=Mesh%ns !*PROBABLY NOT NEEDED !*Var%Ncells=Mesh%ns is always 4096, probably an old structure not used anymore
DELETED Var%Ncells & Mesh%ns and there are no problems :)
!*TYPE(Pvar), DIMENSION(:),ALLOCATABLE:: u_b, u_c !*NOT NEEDED
!*ALERT A subroutine test is called (from timestepping) to perform the mood stabilization. It's basically all in SP, the real are declared as REAL and not REAL(DP)
!*Also in the 1D case the boundary residuals of the different subtimesteps should be combined through the thetas.
!*In this case it is not a big problem because we have
!*STRONG IMPOSITION: nothing to combine (no boundary update so no boundary residuals)
!*PERIODIC: Var%un(1) and Var%un(Mesh%nt) (node residuals of the boundary nodes (first and last) and so combined) are given to the first and the last node.
!*In 2d it is more critical: the weak boundary conditions generate boundary residuals that must be combined through the thetas and this must be corrected in the 2d code
!*ALERT
I report a problem on the strong imposition of the boundary conditions in 1d.Basically BC is called before updating the solution with in fin() 
DO is=1, Mesh%ndofs
       var%ua(is,k-1)=Var%up(is,k-1)-Var%un(is)*Mesh%aires(is) !*Updating
ENDDO
Whenever we want to impose strong Dirichlet (see for example case(4) in init_bc_euler), BC modifies Var%ua through
Var%ua(k-1,i1)%u=whatever we want
But then when we update we lose this information (because we set  var%ua(is,k-1)=Var%up(is,k-1)-Var%un(is)*Mesh%aires(is))
I've talked a bit with Davide about it:One possibility is to call BC before and after the updating. This works but I think we lose a bit in "linearity" and "clarity" of the code. So I suggest to think a bit about it and to find a solution.
NB:->Setting Var%un=0 doesn't work because later we would still havevar%ua(is,k-1)=Var%up(is,k-1)-Var%un(is)*Mesh%aires(is)=Var%up(is,k-1) 
->so you may think that we could modify Var%up and set Var%un=0This doesn't work either because at the next subtimestep you would have a different value for Var%up
For the moment I suggest to apply BCs before and after the updating but as I wrote, this is not a proper solution if we want to fix the code forever.
!*n in Pas_de_temps to be declared with ._DP
n_theta=0 !*<- It is a vector of integers, no need to declare it as DP
->elements_1D
TONS OF ._DP MISSING (for example in gradient_element)
e%log !*NOT NEEDED, I commented and nothing changes
e%n !*INTERNAL normal
e%coeff_b=> NULL() !*NOT USED, not initialized, they can give segmentation fault in fact
aire_element(e)
b is useless
a can be defined as a scalar
normale_element(e)
l useless and also the loop
eval_func_element
x, alpha, beta, a ,b, c, aa, cc, flag
PROCEDURE, PUBLIC:: der_sec=>der_sec_element !*OLD, IT CAN BE REMOVED !*APPARENTLY an alternative to eval_der2
!*In practice much less complete, it lacks P2, P3
!*Moreover for B3 it gives as output sth of dimension 2 which is not so well clear what is
!*FURTHER it takes as input a vector of real and not of PVar. Definitely an old feature to remove
PROCEDURE, PUBLIC:: der2_poly=>der2_poly_element !*OLD, IT CAN BE REMOVED !*APPARENTLY another alternative to eval_der2
!*In practice much less complete, it lacks P2, P3
!*Moreover for B3 it gives as output sth of dimension 4 which is not so well clear what is
!*FURTHER it takes as input a vector of real and not of PVar. Definitely an old feature to remove
FUNCTION gradient_element(e,k,x)
fy,fz useless
average_element(e,u)
i useless
base_ref_element(e) !*ALERT
Depending on the compiler not specifying ._DP can be dangerous
eval_coeff_element_old(e) !*It can be removed !*IT IS NOT USED HERE AND NOT EVEN A PUBLIC PROCEDURE
!*ANYWAY IT IS NOT WRONG
eval_coeff_element(e)
eps !*NOT NEEDED
!*NOT NEEDED
REAL(dp), DIMENSION(2,2):: xp
xp(1,1)=0.0_dp; xp(2,1)=1.0_dp; xp(1,2)=1.0_dp; xp(2,2)=0.0_dp
->aretes !*Arete is the class of the boundaries of the elements which are faces in 3d, segments in 2d, points in 1d,
!*In this case (1d) it is almost useless
!*We have a lot of useless fields which just make sense for higher dimension and are not be used/initialized here
THREE FIELDS ONLY ARE IMPORTANT, THE OTHERS CAN BE DELETED
!*I)bord !*GOOD !*It tells whether the arete is on the boundary of the domain or not
!*NB:
!*Mesh%edge(1)%bord=T
!*Mesh%edge(Mesh%nsegmt)%bord=Mesh%edge(Mesh%nt+1)%bord=T
!*The others are F
!*II) jt1=-1, jt2 =-1 ! les deux elements de par et d'autre de l'arete. !*GOOD
!*The two (in 1d) elements sharing the arete
!*NB: The aretes at the boundary of the domain will have just one element !*IN THIS CASE IT IS SET jt1=jt2=that element !*Clearly this happens for the first and the last one
!*Mesh%edge(1)%jt1=Mesh%edge(1)%jt2=1
!*Mesh%edge(nsegmt)%jt1=Mesh%edge(nsegmt)%jt2=Mesh%nt
!*For the rest we have
!*Mesh%edge(indi)%jt1=indi-1
!*Mesh%edge(indi)%jt2=indi
!*NB:TO REMOVE THE FOLLOWING STRUCTURES YOU HAVE TO CLEAN AND RECOMPILE
INTEGER:: nsommets, itype, nvertex! nbre de dofs, type element: 1-> P1 !*USELESS
!*nsommets !*-10, clearly not used, by logic it should be the number of DoFs in the "edge" which would be 1 in the 1d case
!*itype -1, clearly not used, by logic it should be the itype
!*-1, clearly not used, by logic it should be the number of vertices in the "edge" which would be 1 in the 1d case
INTEGER, DIMENSION(:,:), POINTER :: nu !*USELESS
REAL(dp), DIMENSION(:,:), POINTER :: coor !*USELESS !*NOT EVEN INITIALIZED -> If you try to print you get segmentation fault
volume and jump_flag !*USELESS, moreover they were declared without._DP
REAL(dp), DIMENSION(:), POINTER :: n =>Null() !*USELESS !*NOT EVEN INITIALIZED -> If you try to print you get segmentation fault
!*INTEGER :: log ! logique !*Always 0 !*USELESS
ALL THE QUADRATURE STRUCTURES
->variable_def_euler
ALERT: TONS OF ._DP MISSING SPECIALLY IN THE PARAMETERS
For example pi or gamma
ALSO REMOVE THE DP FROM THE EXPONENT WHEN IT IS A NATURAL **2 AND NOT ._DP
FUNCTION roe_eul(e,u,n) RESULT (J) NOT NEEDED
Inherited from 2d
In flux: FUNCTION flux_eul(Var,x) RESULT(f)
REAL(dp), DIMENSION(n_dim), INTENT(in) :: x !*NOT NEEDED OR AT LEAST OPTIONAL
FUNCTION source_euler(e,x,test) RESULT(s)
Attempt to introduce the source.
It is not bad but probably not used in fact it is identically 0
I leave it even if I think that the source must be in utils
FUNCTION AbsJacobian_eul(Var,x) RESULT(J)
!*REAL(dp) :: Vi2, u, H, eps, p !*NOT NEEDED
lvectors
!*REAL(dp) :: rho, u, p, a, E, H, gmm1, eps
!*NB lvectors are got through inversion, it is better to directly define them to avoid the inverision of the matrix every time
PROCEDURE, PUBLIC:: Nmat => Nmat_eul !*NOT NEEDED AND NOT EVEN PROPERLY DEFINED AT THE END IT IS SET 0._DP
!*These two last functions are drafts of the positive and the negative parts of the Jacobian but they are not public and not used
-min_mat_eul
-max_mat_eul
->timestepping
Adjust the indentation of edge_main_update and add a comment where you say that
even if in fin() edge_main_update is called only for DATA%ischema==5 .OR. DATA%ischema==4 it is called also in mood so the IFs make sense
!*Anyway note how in the "standard" case where we do not have mood we enter in these IFs
!*ALERT
!*THIS IS NOT EVEN COMPILED, I GUESS IT BELONGS TO THE PREVIOUS VERSION OF THE DeC, WHERE THE OPERATOR L^1 WAS MORE COMPLICATED.
!*ANYWAY IT IS NOT USED
!*I SMELL THAT IT COULD ALSO BE WRONG -> Before the updating ua and up are equal so
!*up1(l,lp)-ua1(l,lp)
!*I DO NOT DELETE
!*DAVIDE TOLD ME THAT IT IS IMPORTANT FOR SO I LEAVE IT
!*ANYWAY TO BE USED YOU MUST DO STH ON ua OR up OTHERWISE THEY ARE THE SAME
#if (1==0)
!Gauss-Seidel like
DO l=1,e1%nsommets
DO lp=1, k-1
u1(l)=u1(l)-GAMMA(lp,DATA%iordret-1,DATA%iordret)* (up1(l,lp)-ua1(l,lp))
ENDDO
ENDDO
DO l=1,e2%nsommets
DO lp=1,k-1
u2(l)=u2(l)-GAMMA(lp,DATA%iordret-1,DATA%iordret)* (up2(l,lp)-ua2(l,lp))
ENDDO
ENDDO
#endif
SUBROUTINE test(k_iter,Debug,Var, mesh, DATA, flux_mood) !*ALERT IT WAS MOSTLY IN SP, I'M CHANGING IT INTO DP
->scheme.f90
FUNCTION galerkin
is not used anymore. (Now we use galerkin_new.)
I tried to comment it and it runs without problems,
I tried to rename it and it runs without problems.
so actually galerkin is not used anymore
I leave it because it computes the galerkin residuals with the integration by parts (in galerkin_new it performs the integration integrating directly phi div f without passing the derivative on the test function) and it may be useful to have it already implmented even if it is not used.
NB: IT MUST BE CORRECTED
I REMOVED
u_loc = e%eval_func(u ,e%quad(:,iq))
AND THE FIRST ALTERNATIVE OF THE IF DELETING THE IF AND TAKING JUST THE #ELSE
#if (1==1)
! not quadrature free
fluxloc=u_loc%flux(xx) ! warning: mismatch between bary coord and e%coord
BECAUSE the flux is nonlinear and you cannot expect to compute it on the combination of the u in time thorugh the thetas and get the combination of the fluxes in the different subtimesteps through the thetas
NB: This problem exists also in 2d
lxf
!*NOT USED !*u_diff(i)=(u(i)-vbar(i))
The loop on i is useless
jump
Some ._DP missing
Add the DP
ubar1%u(l)= SUM(uu1%u(l))/REAL(SIZE(uu1))
ubar2%u(l)= SUM(uu2%u(l))/REAL(SIZE(uu2))
jump_old NOT USED, not even public. FUrther I made it public and I tested it, it doesn't work so I deleted it.
NB:I moved a bit some functions and subroutines to make a logic flow: for example, now I have
galerkin
galerkin_new
lxf
psi_galerkin
jump
limit
limit_char
limit_var
->geometry
dx=DATA%Length/REAL(nt,DP) <-DP
._DP missing in the definition of the stuff for internal DoFs
k=0 and several k=k+1 not used
Very misleading comments in the construction of the edges
Mesh%aires(e%nu)=Mesh%aires(e%nu)+e%volume/REAL(e%nsommets,DP)<- DP
->init_bc_euler
Wrong unit of measures in the comments
They should be
[kg/m^3]
[kg/m^3*m/s]
[J/m^3]
I corrected them
Let's add a comment saying that we give the IC in term of conserved variables (and as values that must be later projected into the space of the coefficients)
usual problem with dp
->utils
NOT USED for Euler
usual problem with dp
->algebra
some missing ._DP
I was tired and I just made sure that there were no missing DP or variables in SP
->postprocessing
sol(l)IC(e%coor(l)-DATA%temps,DATA) only good for linear advection with speed=1
I substituted it with
sol(l)=0._DP !*In case it is known, this is the exact solution which one can use to have feedbacks or to make convergence analyses. For example in the linear advection with speed a=1 we have sol(l)IC(e%coor(l)-DATA%temps,DATA)
There is a redundant
u(l)=Var%ua(e%nu(l),0)
NOT USED
DO k=1, e%nsommets
base(k)=e%base(k,e%x(:,l))
grad(:,k)=e%gradient(k,e%x(:,l))
ENDDO
base and grad actually useless, I deleted them
In some files we write the same information. This is due to the fact that some of them are overwritten every time (the "last" files), while some are written for the first time.
DESPITE THIS, some file could be redundant, check.
......@@ -11,11 +11,6 @@
!!! July 10, 2018
!!! Correspondance: remi.abgrall@math.uzh.ch
!!! ------------------------------------------
!*----------------------------------------------------------------------------------------
!*This module (only) provides the function inverse to invert a square matrix
!*----------------------------------------------------------------------------------------
MODULE algebra
USE PRECISION
IMPLICIT NONE
......@@ -63,8 +58,8 @@ CONTAINS
CALL ludcmp(mat2, Nsize, indx, d) !-- indx OUT
DO j = 1, Nsize
col = 0.0_DP
col(j) = 1.0_DP
col = 0.0
col(j) = 1.0
CALL luksb(mat2, Nsize, indx, col) !-- col IN OUT; indx IN
Mat1(:, j) = col(: )
END DO
......@@ -100,7 +95,7 @@ CONTAINS
WRITE( *, *) mod_name, " ERREUR : Matrice singuliere, i==", i
STOP
END IF
vv(i) = 1.0_DP / Big
vv(i) = 1.0 / Big
END DO
DO j = 1, Nsize
......
......@@ -11,100 +11,37 @@
!!! July 10, 2018
!!! Correspondance: remi.abgrall@math.uzh.ch
!!! ------------------------------------------
!*---------------------------------------------------------------------------------
!*Arete is the class of the boundaries of the elements which are faces in 3d, segments in 2d, points in 1d,
!*In this case (1d) it is almost useless
!*We will have a lot of useless fields which just make sense for higher dimension and will not be used/initialized here !*THEY WILL BE MARKED BY USELESS AND IT WILL BE SPECIFIED WHETHER THEY CAN BE SAFELY CANCELED OR NOT !*(NB: If they are used I write GOOD)
!*---------------------------------------------------------------------------------
!*SPOILER: THREE FIELDS ONLY ARE IMPORTANT, THE OTHERS CAN BE DELETED
!*I)bord !*GOOD !*It tells whether the arete is on the boundary of the domain or not
!*NB:
!*Mesh%edge(1)%bord=T
!*Mesh%edge(Mesh%nsegmt)%bord=Mesh%edge(Mesh%nt+1)%bord=T
!*The others are F
!*II) jt1=-1, jt2 =-1 ! les deux elements de par et d'autre de l'arete. !*GOOD
!*The two (in 1d) elements sharing the arete
!*NB: The aretes at the boundary of the domain will have just one element !*IN THIS CASE IT IS SET jt1=jt2=that element !*Clearly this happens for the first and the last one
!*Mesh%edge(1)%jt1=Mesh%edge(1)%jt2=1
!*Mesh%edge(nsegmt)%jt1=Mesh%edge(nsegmt)%jt2=Mesh%nt
!*For the rest we have
!*Mesh%edge(indi)%jt1=indi-1
!*Mesh%edge(indi)%jt2=indi
!*---------------------------------------------------------------------------------
MODULE arete_class
USE PRECISION
IMPLICIT NONE
TYPE, PUBLIC:: arete
!*---SUCCEEDED IN DELETING---
!*INTEGER:: nsommets, itype, nvertex! nbre de dofs, type element: 1-> P1 !*USELESS ---SUCCEEDED IN DELETING---
INTEGER:: nsommets, itype, nvertex! nbre de dofs, type element: 1-> P1
! 2-> B2
!3->P2
! nombre de sommet dans cet element arete
!*nsommets !*-10, clearly not used, by logic it should be the number of DoFs in the "edge" which would be 1 in the 1d case
!*itype -1, clearly not used, by logic it should be the itype
!*-1, clearly not used, by logic it should be the number of vertices in the "edge" which would be 1 in the 1d case
LOGICAL:: bord !*GOOD !*It tells whether the arete is on the boundary of the domain or not
!*NB:
!*Mesh%edge(1)%bord=T
!*Mesh%edge(Mesh%nsegmt)%bord=Mesh%edge(Mesh%nt+1)%bord=T
!*The others are F
INTEGER:: jt1=-1, jt2 =-1 ! les deux elements de par et d'autre de l'arete. !*GOOD
!*The two elements sharing the arete
!*NB: The aretes at the boundary of the domain will have just one element !*IN THIS CASE IT IS SET jt1=jt2=that element !*Clearly this happens for the first and the last one
!*Mesh%edge(1)%jt1=Mesh%edge(1)%jt2=1
!*Mesh%edge(nsegmt)%jt1=Mesh%edge(nsegmt)%jt2=Mesh%nt
!*For the rest we have
!*Mesh%edge(indi)%jt1=indi-1
!*Mesh%edge(indi)%jt2=indi
!*---SUCCEEDED IN DELETING---
!*INTEGER, DIMENSION(:,:), POINTER :: nu =>Null() ! nu( indice des elements, indice des voisins): on cherche a connaitre le numero local des !*USELESS ---SUCCEEDED IN DELETING--- !*0, clearly not used, by logic it should be the vector of the global indices of the DoFs in the "edge" which would be just the global index of the node corresponding to that "edge" in the 1d case
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
!*---SUCCEEDED IN DELETING---
!*REAL(dp), DIMENSION(:,:), POINTER :: coor=>Null() ! il s'agit des coordonnees physique des dofs (communs) sur la face (commune) !*USELESS ---SUCCEEDED IN DELETING--- !*NOT EVEN INITIALIZED, when you try to print it you may get a segmentation fault
!*By logic it should be the coordinates of the DoFs contained by the edge (coordinates of a single point in 1d)
!*volume and jump_flag USELESS ---SUCCEEDED IN DELETING---
!*NB: volume and jump_flag were declared without ._DP
!*I guess they are useless but in any case...
!*Anyway
!*volume and jump_flag are never changed, they stay 0 and 1
!*I guess they are not used
!*volume should be the measure of the arete (one point in 1d, so trivially 0)
!*jump_flag??????????
!*---SUCCEEDED---!*REAL(dp) :: volume=0._DP !
!*---SUCCEEDED---!*REAL(dp) :: jump_flag=1.0_DP !(between 0 and 1 which gives the weight of the edge term
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)
!*USELESS ---SUCCEEDED IN DELETING---
!*REAL(dp), DIMENSION(:), POINTER :: n =>Null() ! normales exterieures !*USELESS
!*Actually in 2d it makes sense to define the normal to an arete because it is a segment, in 1d we have a point so it makes no sense !*SEGMENTATION FAULT, NOT EVEN DEFINED
!*USELESS ---SUCCEEDED IN DELETING---
!*INTEGER :: log ! logique !*Always 0
!*USELESS ---SUCCEEDED IN DELETING---
!*The quadrature doesn't make sense for an arete in 1d (point)
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
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
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 !*USELESS !*---SUCCEEDED IN DELETING---
!*PROCEDURE, PUBLIC:: quadrature=>quadrature_arete
!*PROCEDURE, PUBLIC:: normale=>normale_arete !*USELESS !*---SUCCEEDED IN DELETING---
CONTAINS
PROCEDURE, PUBLIC:: aire=>aire_arete
PROCEDURE, PUBLIC:: quadrature=>quadrature_arete
PROCEDURE, PUBLIC:: normale=>normale_arete
!FINAL:: clean_arete
END TYPE arete
......@@ -112,70 +49,60 @@ MODULE arete_class
CONTAINS
!*aire_arete and normale_arete are clearly bidimensional oriented
!*They can be removed. Moreover they involve coor which is not even initialized (if you try to print it you get a segmentation fault)
!*!*-------------------
!*!*In 2d it would be the length of the boundary element !*USELESS ---SUCCEEDED IN DELETING---
!*!*-------------------
!*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
!*!*-------------------
!*!*In 2d it would be the normal to the boundary element ---SUCCEEDED IN DELETING---
!*!*-------------------
!*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) ---SUCCEEDED IN DELETING---
!*IF (ASSOCIATED(e%coor)) NULLIFY(e%coor) ---SUCCEEDED IN DELETING---
!*IF (ASSOCIATED(e%quad)) NULLIFY(e%quad)
!*IF (ASSOCIATED(e%weight)) NULLIFY(e%weight)