Skip to content

GitLab

Commit 80c342c5 authored by rabgra's avatar rabgra
Browse files

remove some comments and un-necessary things

parent 216aa3ea
Hide whitespace changes
Inline Side-by-side
Src2D2/geom.f90
View file @ 80c342c5
......@@ -16,7 +16,7 @@ MODULE geom
USE param2d
USE Boundary
IMPLICIT NONE
private
PRIVATE
! GMSH_ELE(i, j): i = #Ndofs, j = # Verts
INTEGER, DIMENSION(31, 2), PARAMETER :: GMSH_ELE = &
......@@ -28,16 +28,10 @@ MODULE geom
!1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
(/2,3,4,4,8,6,5,2,3,4,4,8,6,5,1,4,8,6,5,3,3,3,3,3,3,2,2,2,4,4,4/)
INTERFACE coord_kinetic
MODULE PROCEDURE coord_kinetic_new
END INTERFACE coord_kinetic
INTERFACE calculateMassMatrLumped
MODULE PROCEDURE calculateMassMatrLumped_new
END INTERFACE calculateMassMatrLumped
public:: ReadMeshGMSH2, coord_kinetic, calculateMassMatrLumped
PUBLIC:: ReadMeshGMSH2, coord_kinetic, calculateMassMatrLumped
CONTAINS
......@@ -107,7 +101,7 @@ CONTAINS
ALLOCATE (Physname(Nphys))
DO i=1,Nphys
READ(UNIT,*) dummy,idummy,PhysName(idummy)
PRINT*, trim(adjustl(physName(idummy)))
PRINT*, TRIM(ADJUSTL(physName(idummy)))
ENDDO
READ(UNIT,*) !!$EndPhysicalNames
......@@ -136,10 +130,6 @@ CONTAINS
ALLOCATE( ele(i)%NU(GMSH_ELE(ele_type, 1)) )
ALLOCATE( ele(i)%VV(GMSH_ELE(ele_type, 2)) )
! if (ele_type == 26) then
! print*, 'ele_type = ', ele_type
! print*, 'GMSH_ELE(ele_type, 1) = ', GMSH_ELE(ele_type,1), ', GMSH_ELE(ele_type, 2) = ', GMSH_ELE(ele_type, 2)
! end if
ALLOCATE( v_dummy(n_tag-1) )
......@@ -152,7 +142,7 @@ CONTAINS
ele(i)%VV = ele(i)%NU( 1:GMSH_ELE(ele_type, 2) )
! ele(i)%bc_tag = v_dummy(1)
DEALLOCATE(v_dummy)
......@@ -212,7 +202,7 @@ CONTAINS
PRINT*, 'N_b = ', N_b
Mesh%Nt = N_ele - N_b
! Nullify(Mesh%e)
ALLOCATE (Mesh%e(Mesh%nt))
......@@ -247,7 +237,7 @@ CONTAINS
ALLOCATE(Mesh%e(ii)%base_at_dofs( SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
ALLOCATE(Mesh%e(ii)%inv_base_at_dofs( SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
allocate(Mesh%e(ii)%grad_at_dofs(n_dim,SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
ALLOCATE(Mesh%e(ii)%grad_at_dofs(n_dim,SIZE(ele(i)%nu), SIZE(ele(i)%nu)))
CALL Mesh%e(ii)%base_ref()
CALL Mesh%e(ii)%grad_ref()
......@@ -306,9 +296,9 @@ CONTAINS
ENDDO
! dans le cas ordre 3 on cherche
! If third order we look for
IF (Mesh%period) THEN
IF (SIZE(ele(1)%nu).GT.2) THEN !assume carre
IF (SIZE(ele(1)%nu).GT.2) THEN !assume square
ALLOCATE(voisin(Mesh%ns,2))
voisin=-1
DO i=1, N_b
......@@ -332,7 +322,7 @@ CONTAINS
IF ((i1==j1.AND.i2==j2).OR.(i1==j2.AND.i2==j1)) THEN
Mesh%per(ele(i)%nu(3))=ele(i)%nu(3)
EXIT ! on est sur un cote qui se mape lui meme
EXIT ! we are on a side that maps onto itself
ENDIF
l11=voisin(j1,1)
......@@ -373,55 +363,12 @@ CONTAINS
! MassMatrLumped is the global lumped mass matrix (PETSc)
! MatLumped is the global lumped mass matrix (Fortran)
! MatLumpedInv is its inverse (using Inverse function from algebra.f90)
SUBROUTINE calculateMassMatrLumped_old(DATA, Mesh)
implicit none
CHARACTER(LEN = *), PARAMETER :: mod_name = "calculateMassMatrLumped_old"
TYPE(donnees), INTENT(in) :: DATA
TYPE(maillage), INTENT(inout) :: Mesh
TYPE(element) :: e
INTEGER, ALLOCATABLE :: iss(:)
INTEGER:: jt, i
ALLOCATE(Mesh%aires(1:Mesh%ns))
Mesh%aires = 0.0_dp
! Loop through mesh elements
DO jt=1, Mesh%nt
! get element
e=Mesh%e(jt)
ALLOCATE(iss(e%nsommets))
! mapping between the local ang global DOFs
iss(1:e%nsommets)=e%nu(1:e%nsommets)
! Loop through the DOFs of the element e
DO i = 1,e%nsommets
! update the entry of the lumped mass matrix (Fortran)
Mesh%aires(iss(i)) = Mesh%aires(iss(i)) + e%volume/REAL(e%nsommets,dp)
END DO ! i
DEALLOCATE(iss)
END DO ! jt
#ifndef parallel
Mesh%aires=1.0_dp/Mesh%aires
#endif
END SUBROUTINE calculateMassMatrLumped_old
SUBROUTINE calculateMassMatrLumped_new(DATA, Mesh)
implicit none
CHARACTER(LEN = *), PARAMETER :: mod_name ="calculateMassMatrLumped_new"
SUBROUTINE calculateMassMatrLumped(DATA, Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="calculateMassMatrLumped"
TYPE(donnees), INTENT(in) :: DATA
TYPE(maillage), INTENT(inout) :: Mesh
TYPE(element) :: e
......@@ -493,90 +440,11 @@ CONTAINS
END SUBROUTINE volume
END SUBROUTINE calculateMassMatrLumped_new
SUBROUTINE coord_kinetic_old(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic_old"
TYPE(maillage), INTENT(inout):: Mesh
REAL(dp), DIMENSION(:,:), ALLOCATABLE:: vals
REAL(dp), DIMENSION(:,:), ALLOCATABLE:: coords
REAL(dp):: airs
TYPE(element):: e
TYPE(frontiere):: efr
INTEGER:: jt, l
INTEGER:: ndim=2
SELECT CASE (Mesh%e(1)%itype)
CASE(1)
ALLOCATE( vals(3,3) )
VALS(1 ,1)=1._qp/6._qp
VALS(2:3,1)=1._qp/12._qp
VALS(2 ,2)=VALS(1,1)
VALS(1 ,2)=VALS(2,1)
VALS(3 ,2)=VALS(2,1)
VALS(1:2,3)=VALS(2,1)
VALS(3 ,3)=VALS(1,1)
CASE(2)
ALLOCATE(vals(3,6) )
VALS(1,1)=1._qp/10._qp
VALS(2:3,1)=1._qp/30._qp
VALS(2,2)=VALS(1,1)
VALS(1,2)=VALS(2,1)
VALS(3,2)=VALS(2,1)
VALS(1:2,3)=VALS(2,1)
VALS(3,3)=VALS(1,1)
VALS(3,4)=1._qp/30._qp
VALS(2,5)=VALS(3,4)
VALS(1,6)=VALS(3,4)
VALS(1:2,4)=1._qp/15._qp
VALS(2:3,6)=1._qp/15._qp
VALS(1,5)=1._qp/15._qp
VALS(3,5)=1._qp/15._qp
CASE default
PRINT*, mod_name
PRINT*, "Kinetic: this case does not exist"
PRINT*, "type=", Mesh%e(1)%itype
END SELECT
ALLOCATE(coords(ndim, Mesh%ndofs))
coords=0._dp
DO jt=1, Mesh%nt
ALLOCATE(Mesh%e(jt)%yy(ndim,e%nsommets))
e=Mesh%e(jt)
DO l=1,e%nsommets
e%yy(1,l)=SUM(VALS(1:3,l)*e%coor(1,1:3) )
e%yy(2,l)=SUM(VALS(1:3,l)*e%coor(2,1:3) )
coords(:,e%nu(l))=coords(:,e%nu(l))+ e%yy(:,l) *e%volume
ENDDO
ENDDO
coords(1,:)=coords(1,:) * mesh%aires
coords(2,:)=coords(2,:) * mesh%aires
DO jt=1, Mesh%nt
ALLOCATE(Mesh%e(jt)%y(ndim,e%nsommets))
e=Mesh%e(jt)
DO l=1, e%nsommets
e%y(:,l)=coords(:,e%nu(l))
ENDDO
ENDDO
! for the boundary conditions
END SUBROUTINE calculateMassMatrLumped
DO jt=1, Mesh%Nsegfr
ALLOCATE(Mesh%fr(jt)%y(ndim, Mesh%fr(jt)%nsommets))
efr=Mesh%fr(jt)
DO l=1, efr%nsommets
efr%y(:,l)=coords(:,efr%nu(l))
ENDDO
ENDDO
DEALLOCATE(coords)
END SUBROUTINE coord_kinetic_old
SUBROUTINE coord_kinetic_new(Mesh)
SUBROUTINE coord_kinetic(Mesh)
! valid for quad & triangles
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic_new"
......@@ -590,29 +458,29 @@ CONTAINS
INTEGER :: ndim=2
vals=0._dp
! e%nsommet=3,triangle
VALS(1 ,1,3)=1._qp/6._qp
VALS(2:3,1,3)=1._qp/12._qp
VALS(2 ,2,3)=VALS(1,1,3)
VALS(1 ,2,3)=VALS(2,1,3)
VALS(3 ,2,3)=VALS(2,1,3)
VALS(1:2,3,3)=VALS(2,1,3)
VALS(3 ,3,3)=VALS(1,1,3)
! e%nsommet=6,triangle
VALS(1,1,6)=1._qp/10._qp
VALS(2:3,1,6)=1._qp/30._qp
VALS(2,2,6)=VALS(1,1,6)
VALS(1,2,6)=VALS(2,1,6)
VALS(3,2,6)=VALS(2,1,6)
VALS(1:2,3,6)=VALS(2,1,6)
VALS(3,3,6)=VALS(1,1,6)
VALS(3,4,6)=1._qp/30._qp
VALS(2,5,6)=VALS(3,4,6)
VALS(1,6,6)=VALS(3,4,6)
VALS(1:2,4,6)=1._qp/15._qp
VALS(2:3,6,6)=1._qp/15._qp
VALS(1,5,6)=1._qp/15._qp
VALS(3,5,6)=1._qp/15._qp
! e%nsommet=3,triangle
VALS(1 ,1,3)=1._qp/6._qp
VALS(2:3,1,3)=1._qp/12._qp
VALS(2 ,2,3)=VALS(1,1,3)
VALS(1 ,2,3)=VALS(2,1,3)
VALS(3 ,2,3)=VALS(2,1,3)
VALS(1:2,3,3)=VALS(2,1,3)
VALS(3 ,3,3)=VALS(1,1,3)
! e%nsommet=6,triangle
VALS(1,1,6)=1._qp/10._qp
VALS(2:3,1,6)=1._qp/30._qp
VALS(2,2,6)=VALS(1,1,6)
VALS(1,2,6)=VALS(2,1,6)
VALS(3,2,6)=VALS(2,1,6)
VALS(1:2,3,6)=VALS(2,1,6)
VALS(3,3,6)=VALS(1,1,6)
VALS(3,4,6)=1._qp/30._qp
VALS(2,5,6)=VALS(3,4,6)
VALS(1,6,6)=VALS(3,4,6)
VALS(1:2,4,6)=1._qp/15._qp
VALS(2:3,6,6)=1._qp/15._qp
VALS(1,5,6)=1._qp/15._qp
VALS(3,5,6)=1._qp/15._qp
!!!
ALLOCATE(coords(ndim, Mesh%ndofs))
......@@ -631,7 +499,7 @@ CONTAINS
coords(:,e%nu(l))=coords(:,e%nu(l))+ e%yy(:,l) *e%volume
ENDDO
CASE(4) ! quadrangle
DO l=1, e%nsommets
DO l=1, e%nsommets
CALL machin(l,e,z)
e%yy(:,l)=z/e%volume
coords(:,e%nu(l))=coords(:,e%nu(l))+ z
......@@ -673,22 +541,14 @@ CONTAINS
REAL(dp), DIMENSION(2) :: z
REAL(dp), DIMENSION(2,2) :: jac
REAL(dp), DIMENSION(2) :: a,b,c
y(:)=0._dp; vol=0._dp
DO iq=1, e%nquad
x=e%quad(1:2,iq)
z=e%iso(x)
! a=e%coor(:,2)-e%coor(:,1)
! b=e%coor(:,4)-e%coor(:,1)
! c=e%coor(:,3)-e%coor(:,2)+e%coor(:,1)-e%coor(:,4)
! z=e%coor(:,1)+ ( a+c*x(2) )*x(1) + b*x(2) ! location of the point thanks to iso parametric transformation
! jac(1,1)=a(1)+x(2)*c(1)
! Jac(1,2)=b(1)+x(1)*c(1)
! Jac(2,1)=a(2)+x(2)*c(2)
! Jac(2,2)=b(2)+x(1)*c(2)
Jac=e%d_iso(x)
z=e%iso(x)
Jac=e%d_iso(x)
det=ABS( Jac(1,1)*Jac(2,2) - Jac(1,2)*Jac(2,1) )
phi=e%base(l,x(1:2))
......@@ -697,13 +557,13 @@ Jac=e%d_iso(x)
ENDDO
END SUBROUTINE machin
END SUBROUTINE coord_kinetic_new
END SUBROUTINE coord_kinetic
SUBROUTINE coord_kinetic_ent(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic entropy"
......@@ -714,12 +574,12 @@ Jac=e%d_iso(x)
SUBROUTINE voisinage(Mesh)
IMPLICIT NONE
CHARACTER(LEN = *), PARAMETER :: mod_name = "voisinage"
CHARACTER(LEN = *), PARAMETER :: mod_name = "voisinage"
TYPE(Maillage), INTENT(inout):: mesh
type(element):: e
TYPE(element):: e
INTEGER:: jt, k, is
integer, dimension(:), allocatable:: is2jt, compte
INTEGER, DIMENSION(:), ALLOCATABLE:: is2jt, compte
ALLOCATE(is2jt(Mesh%ns), compte(Mesh%ns) )
is2jt=0
compte=0
......@@ -729,7 +589,7 @@ Jac=e%d_iso(x)
is2jt(e%nu(k))=is2jt(e%nu(k))+1
ENDDO
ENDDO
allocate (Mesh%vois(Mesh%ns))!dofs) )
ALLOCATE (Mesh%vois(Mesh%ns))!dofs) )
DO is=1, Mesh%ns!dofs
ALLOCATE(Mesh%vois(is)%nvois(is2jt(is)), Mesh%vois(is)%loc(is2jt(is)) )
Mesh%vois(is)%nvois=0
......@@ -745,12 +605,12 @@ Jac=e%d_iso(x)
mesh%vois( e%nu(k) )%loc( compte(e%nu(k)) )=k
ENDDO
ENDDO
do is=1, Mesh%ns!dofs
mesh%vois(is)%nbre=size(Mesh%vois(is)%nvois)
enddo
DO is=1, Mesh%ns!dofs
mesh%vois(is)%nbre=SIZE(Mesh%vois(is)%nvois)
ENDDO
deallocate(compte, is2jt)
DEALLOCATE(compte, is2jt)
END SUBROUTINE voisinage
END MODULE geom
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment