Commit 4f7fe35f by rabgra

### correction of bugs, add scheme 8

parent 9ae45a58
 MODULE quickhull !-------------------------------------------------- ! compute the convex hull of a family of points !Ref: https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#Pseudo-code ! 2D, complexity O(nlog(n)) ! algo: ! Input: a list P of points in the plane. ! !Precondition: There must be at least 3 points. ! !Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate). ! !Initialize U and L as empty lists. !The lists will hold the vertices of upper and lower hulls respectively. ! !for i = 1, 2, ..., n: ! while L contains at least two points and the sequence of last two points ! of L and the point P[i] does not make a counter-clockwise turn: ! remove the last point from L ! append P[i] to L ! !for i = n, n-1, ..., 1: ! while U contains at least two points and the sequence of last two points ! of U and the point P[i] does not make a counter-clockwise turn: ! remove the last point from U ! append P[i] to U ! !Remove the last point of each list (it's the same as the first point of the other list). !Concatenate L and U to obtain the convex hull of P. !Points in the result will be listed in counter-clockwise order. !__________________________________________________________ USE PRECISION IMPLICIT NONE PRIVATE PUBLIC:: test_quick REAL(DP), PARAMETER:: tol=-SQRT(EPSILON(0.0_dp)) REAL(dp), PARAMETER:: tol_hull=1.e-4!0._dp!-epsilon(0.0_dp) CONTAINS SUBROUTINE tri(x,y) IMPLICIT NONE REAL(DP), DIMENSION(:,:), INTENT(in) :: x REAL(DP), DIMENSION(:,:), INTENT(out) :: y REAL(DP), DIMENSION(2, SIZE(x, dim=2)):: buff LOGICAL, DIMENSION(SIZE(x,dim=2)) :: ind INTEGER :: m, j, k, i m=SIZE(x,dim=2) ind=.TRUE. k=0 DO i=1, m j=MINLOC(x(1,:),1, ind) k=k+1 buff(:, k)=x(:,j) ind(j)=.FALSE. ENDDO y=buff END SUBROUTINE tri REAL(dp) FUNCTION Cross(v1,v2,v3) ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------- ! INPUT VARIABLES REAL(DP),INTENT(IN) :: v1(2) !< input vector 1 REAL(DP),INTENT(IN) :: v2(2) !< input vector 2 REAL(DP),INTENT(IN) :: v3(2) !< input vector 3 !----------------------------------------------- ! OUTPUT VARIABLES REAL(DP) :: w1(2), w2(2) !----------------------------------------------- ! LOCAL VARIABLES !=============================================== w1=(v2-v1)/( SQRT( dot_PRODUCT( (v2-v1) , (v2-v1) ) ) -tol ) w2=(v3-v1)/( SQRT( dot_PRODUCT( (v3-v1) , (v3-v1) ) ) -tol ) Cross= w1(1) * w2(2) - w1(2) * w2(1) END FUNCTION Cross SUBROUTINE ConvHull(nPoints,Points_o,nHull,Hull)!, ilower, Lower, iupper, upper) ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE CHARACTER(LEN = *), PARAMETER :: mod_name ="ConvHull" !------------------------------------------------ ! INPUT VARIABLES INTEGER,INTENT(IN) :: nPoints REAL(DP),INTENT(IN) :: Points_o(2,0:nPoints-1) !------------------------------------------------ ! OUTPUT VARIABLES INTEGER,INTENT(OUT) :: nHull ! NOTE: allocate Hull always one point greater than Points, because we save the first value twice REAL(DP),INTENT(OUT) :: Hull(2,0:nPoints) !------------------------------------------------ REAL(DP) :: Points(2,0:nPoints-1) ! LOCAL VARIABLES REAL(DP) :: Lower(2,0:nPoints-1) REAL(DP) :: Upper(2,0:nPoints-1) INTEGER :: iLower,iUpper INTEGER :: i !================================================ IF(nPoints.LE.1)THEN Hull = Points nHull = nPoints ELSE CALL tri(Points_o,Points) iLower = 0 Lower = -HUGE(1._dp) DO i=0,nPoints-1 DO WHILE(iLower.GE.2.AND.Cross(Lower(:,iLower-2),Lower(:,iLower-1),Points(:,i)).LE.0._dp) Lower(:,iLower) = -HUGE(1.) iLower = iLower - 1 END DO Lower(:,iLower) = Points(:,i) iLower = iLower + 1 END DO iUpper = 0 Upper = HUGE(1._dp) DO i=nPoints-1,0,-1 DO WHILE(iUpper.GE.2.AND.Cross(Upper(:,iUpper-2),Upper(:,iUpper-1),Points(:,i)).LE.0._dp) Upper(:,iUpper) = HUGE(1._dp) iUpper = iUpper - 1 END DO Upper(:,iUpper) = Points(:,i) iUpper = iUpper + 1 END DO iLower = iLower-1 iUpper = iUpper-1 nHull = iLower+iUpper+1 ! NOTE: Initialize Hull with zeros Hull = 0._dp ! NOTE: save values in Hull Hull(:,0 :iLower -1) = Lower(:,0:iLower-1) Hull(:,iLower:iLower+iUpper-1) = Upper(:,0:iUpper-1) ! NOTE: save first value twice Hull(:, iLower+iUpper ) = Hull (:,0 ) END IF END SUBROUTINE ConvHull LOGICAL FUNCTION is_in (Nhull, Hull, x, lk, jt) IMPLICIT NONE ! The parameter tol must be negative to be less strict that <0. CHARACTER(LEN = *), PARAMETER :: mod_name ="is_in" INTEGER, INTENT(in):: Nhull, lk, jt REAL(DP), DIMENSION(2,0:Nhull-1), INTENT(in):: Hull REAL(DP), DIMENSION(2), INTENT(in):: x REAL(DP), DIMENSION(Nhull-1):: prod INTEGER:: l DO l=1, Nhull-1 prod(l)=cross(x, hull(:,l-1), hull(:,l)) ENDDO IF (MINVAL(prod)*MAXVAL(prod).LE.-tol_hull) THEN! this is essential is_in=.FALSE. ELSE is_in=.TRUE. ENDIF END FUNCTION is_in LOGICAL FUNCTION test_quick(nPoints,Points,x,n,jt, lk, flag) IMPLICIT NONE CHARACTER(LEN = *), PARAMETER :: mod_name ="test_quick" ! INPUT VARIABLES INTEGER,INTENT(IN) :: nPoints,n, jt, lk REAL(DP),INTENT(IN) :: Points(2,0:nPoints-1) REAL(DP), DIMENSION(2,n), INTENT(in):: x LOGICAL, INTENT(in), OPTIONAL:: flag INTEGER :: nHull ! NOTE: allocate Hull always one point greater than Points, because we save the first value twice REAL(DP) :: Hull(2,0:nPoints) INTEGER:: i LOGICAL:: test CALL ConvHull(npoints, Points, nHull,Hull) IF (PRESENT(flag)) THEN PRINT*, nhull DO i=0, nhull-1 WRITE(3,*) hull(:,i) ENDDO ENDIF DO i=1, n test=is_in(Nhull,Hull,x(:,i),lk, jt) IF (.NOT.test) THEN test_quick=.FALSE. RETURN ENDIF ENDDO test_quick=.TRUE. END FUNCTION test_quick END MODULE quickhull
 ... @@ -16,7 +16,7 @@ MODULE geom ... @@ -16,7 +16,7 @@ MODULE geom USE param2d USE param2d USE Boundary USE Boundary IMPLICIT NONE IMPLICIT NONE PRIVATE private ! GMSH_ELE(i, j): i = #Ndofs, j = # Verts ! GMSH_ELE(i, j): i = #Ndofs, j = # Verts INTEGER, DIMENSION(31, 2), PARAMETER :: GMSH_ELE = & INTEGER, DIMENSION(31, 2), PARAMETER :: GMSH_ELE = & ... @@ -28,10 +28,16 @@ MODULE geom ... @@ -28,10 +28,16 @@ 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 !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/) (/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 CONTAINS ... @@ -101,7 +107,7 @@ CONTAINS ... @@ -101,7 +107,7 @@ CONTAINS ALLOCATE (Physname(Nphys)) ALLOCATE (Physname(Nphys)) DO i=1,Nphys DO i=1,Nphys READ(UNIT,*) dummy,idummy,PhysName(idummy) READ(UNIT,*) dummy,idummy,PhysName(idummy) PRINT*, TRIM(ADJUSTL(physName(idummy))) PRINT*, trim(adjustl(physName(idummy))) ENDDO ENDDO READ(UNIT,*) !!$EndPhysicalNames READ(UNIT,*) !!$EndPhysicalNames ... @@ -130,6 +136,10 @@ CONTAINS ... @@ -130,6 +136,10 @@ CONTAINS ALLOCATE( ele(i)%NU(GMSH_ELE(ele_type, 1)) ) ALLOCATE( ele(i)%NU(GMSH_ELE(ele_type, 1)) ) ALLOCATE( ele(i)%VV(GMSH_ELE(ele_type, 2)) ) 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) ) ALLOCATE( v_dummy(n_tag-1) ) ... @@ -142,7 +152,7 @@ CONTAINS ... @@ -142,7 +152,7 @@ CONTAINS ele(i)%VV = ele(i)%NU( 1:GMSH_ELE(ele_type, 2) ) ele(i)%VV = ele(i)%NU( 1:GMSH_ELE(ele_type, 2) ) ! ele(i)%bc_tag = v_dummy(1) DEALLOCATE(v_dummy) DEALLOCATE(v_dummy) ... @@ -202,7 +212,7 @@ CONTAINS ... @@ -202,7 +212,7 @@ CONTAINS PRINT*, 'N_b = ', N_b PRINT*, 'N_b = ', N_b Mesh%Nt = N_ele - N_b Mesh%Nt = N_ele - N_b ! Nullify(Mesh%e) ALLOCATE (Mesh%e(Mesh%nt)) ALLOCATE (Mesh%e(Mesh%nt)) ... @@ -237,7 +247,7 @@ CONTAINS ... @@ -237,7 +247,7 @@ CONTAINS ALLOCATE(Mesh%e(ii)%base_at_dofs( SIZE(ele(i)%nu), SIZE(ele(i)%nu))) 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)%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)%base_ref() CALL Mesh%e(ii)%grad_ref() CALL Mesh%e(ii)%grad_ref() ... @@ -296,9 +306,9 @@ CONTAINS ... @@ -296,9 +306,9 @@ CONTAINS ENDDO ENDDO ! If third order we look for ! dans le cas ordre 3 on cherche IF (Mesh%period) THEN IF (Mesh%period) THEN IF (SIZE(ele(1)%nu).GT.2) THEN !assume square IF (SIZE(ele(1)%nu).GT.2) THEN !assume carre ALLOCATE(voisin(Mesh%ns,2)) ALLOCATE(voisin(Mesh%ns,2)) voisin=-1 voisin=-1 DO i=1, N_b DO i=1, N_b ... @@ -322,7 +332,7 @@ CONTAINS ... @@ -322,7 +332,7 @@ CONTAINS IF ((i1==j1.AND.i2==j2).OR.(i1==j2.AND.i2==j1)) THEN IF ((i1==j1.AND.i2==j2).OR.(i1==j2.AND.i2==j1)) THEN Mesh%per(ele(i)%nu(3))=ele(i)%nu(3) Mesh%per(ele(i)%nu(3))=ele(i)%nu(3) EXIT ! we are on a side that maps onto itself EXIT ! on est sur un cote qui se mape lui meme ENDIF ENDIF l11=voisin(j1,1) l11=voisin(j1,1) ... @@ -363,12 +373,55 @@ CONTAINS ... @@ -363,12 +373,55 @@ CONTAINS ! MassMatrLumped is the global lumped mass matrix (PETSc) ! MassMatrLumped is the global lumped mass matrix (PETSc) ! MatLumped is the global lumped mass matrix (Fortran) ! MatLumped is the global lumped mass matrix (Fortran) ! MatLumpedInv is its inverse (using Inverse function from algebra.f90) ! 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 SUBROUTINE calculateMassMatrLumped(DATA, Mesh) IMPLICIT NONE ALLOCATE(Mesh%aires(1:Mesh%ns)) CHARACTER(LEN = *), PARAMETER :: mod_name ="calculateMassMatrLumped" 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" TYPE(donnees), INTENT(in) :: DATA TYPE(donnees), INTENT(in) :: DATA TYPE(maillage), INTENT(inout) :: Mesh TYPE(maillage), INTENT(inout) :: Mesh TYPE(element) :: e TYPE(element) :: e ... @@ -440,11 +493,90 @@ CONTAINS ... @@ -440,11 +493,90 @@ CONTAINS END SUBROUTINE volume END SUBROUTINE volume END SUBROUTINE calculateMassMatrLumped 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 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(Mesh) SUBROUTINE coord_kinetic_new(Mesh) ! valid for quad & triangles ! valid for quad & triangles IMPLICIT NONE IMPLICIT NONE CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic_new" CHARACTER(LEN = *), PARAMETER :: mod_name ="coord_kinetic_new" ... @@ -458,29 +590,29 @@ CONTAINS ... @@ -458,29 +590,29 @@ CONTAINS INTEGER :: ndim=2 INTEGER :: ndim=2 vals=0._dp vals=0._dp ! e%nsommet=3,triangle ! e%nsommet=3,triangle VALS(1 ,1,3)=1._qp/6._qp VALS(1 ,1,3)=1._qp/6._qp VALS(2:3,1,3)=1._qp/12._qp VALS(2:3,1,3)=1._qp/12._qp VALS(2 ,2,3)=VALS(1,1,3) VALS(2 ,2,3)=VALS(1,1,3) VALS(1 ,2,3)=VALS(2,1,3) VALS(1 ,2,3)=VALS(2,1,3) VALS(3 ,2,3)=VALS(2,1,3) VALS(3 ,2,3)=VALS(2,1,3)