!> anisotorpic mesh adaptation
module anisotropic
  use paramets
  use main_data  ! contains "type(mesh) ::  grid"   for computation
  use mesh_oper
  use problem_oper
  use ama_L2interpol
  use set_solution
  use mesh_adaptation
!  use AMA_interpol
  !use AMA_estims
  use AMAdata
  !use angener77
  use angener90

  implicit none

  public:: AdaptMesh_Angener
  public:: ReadingOfParamets
  public:: ANI_Set_Params
  public:: SetAngenerTriang
  public:: SetGrid
  public:: GenerateAMAmetric

contains
  !> Anisotropic mesh adaptation useing the web variant of ANGENER 3.2 in Fortran 77
  subroutine AdaptMesh_Angener(metric, identical_AMA_grids)
    logical, intent(in) :: metric ! =1 metric generated by F77 ANGENER, =0 in AnisotErrorEstimates
    logical, intent(inout) :: identical_AMA_grids ! end of mesh adaptation
    integer :: ifig, ifig1, i, icon
    integer :: id1, id2, nconstr
    integer :: notUsed_i, notUsed_ii
    character*1  ch
    character(len=50) :: not_used
    character(len=50) :: subMeshType


    !!allocate(AMA)  - allocated in mainAD.f90

    AMA%adapt_level = state%space%adapt%adapt_level

    print*
    print*,'# Calling ANGENER', ndim

    !call system('free')

    ! logical variant fr a test
    AMA%test = .false.

    !vizialization aftear each set of local operations: mesh0, mesh1, ...
    AMA%ifig = -1  !     ...  no vizualization
    !AMA%ifig = 0   !  ...     vizualization
    !if(AMA%adapt_level == 4) AMA%ifig = 0   !     ...     vizualization

    AMA%ifig1 = 9

    open(AMA%ifig1, file='AMA_opers', status='UNKNOWN', position='append')

    !call ReadingOfParametsOLD(ifv, icrack, iwa, iwall )
    call ReadingOfParamets( metric )

    !write(*,*) ityp, ifv,pos, pos1, numel, epsilon1, p, Re, xte1, yte1, xte2, yte2
    !write(*,*) AMA%ityp, AMA%ifv,AMA%pos, AMA%pos1, AMA%numel, AMA%epsilon1, AMA%p, AMA%Re, &
    !     AMA%xte(:,:)


    call ANI_Set_Params( )

    ! estimate of the size of the new mesh

    !if(AMA%ityp .ne. 0) then
       !     standard
       AMA%melem = max(20*AMA%numel,10*AMA%nelem)
       AMA%mpoin = max(10*AMA%numel,10*AMA%npoin)
       AMA%mbelm = max(40*int(sqrt(1.*AMA%numel)),10*AMA%nbelm)
    !else
    !   !     fine isotropic
    !   AMA%melem = max(2*AMA%numel,7*AMA%nelem)
    !   AMA%mpoin = max(1*AMA%numel,7*AMA%npoin)
    !   AMA%mbelm = max(4*int(sqrt(1.*AMA%numel)),7*AMA%nbelm)
    !endif


    !reading of node position constrains
    if(len(lines_file) > 3) then
       icon = 19
       open(icon,file=lines_file, status='old')
       read(icon, *) AMA%nconstr

       if( AMA%nconstr > 0) then
          read(icon, *) subMeshType

          if ( subMeshType=='NONE' .or. subMeshType=='convex' .or. subMeshType=='boundary') then
              ! DO NOTHING
              read(icon, *) notUsed_i
          else if ( subMeshType == 'convex2') then
               !print*, 'lines_file of type convex2 contains an additional parameter!'
               read(icon, *) notUsed_i, notUsed_ii
          else
            print*, 'unknown type of lines_file in AdaptMesh_Angener'
            stop 'Please, control number of lines in your file!'
          endif

          nconstr = AMA%nconstr
          ! we need two constrains in order to keep the nodes, the second one is identical
          if(AMA%nconstr == 1)  nconstr = 2


          allocate(AMA%constr(1:nconstr, 1:2, 1:2))
          allocate(AMA%iconstr(1:nconstr))

          do i=1,  AMA%nconstr
             read(icon, *) AMA%iconstr(i), AMA%constr(i, 1, 1:2), AMA%constr(i, 2, 1:2)
             !write(*, '(a10, i5, 4es12.4)') &
             !     'Constrs:', AMA%iconstr(i), AMA%constr(i, 1, 1:2), AMA%constr(i, 2, 1:2)
          enddo

          if(AMA%nconstr == 1) then
             AMA%nconstr = 2
             AMA%iconstr(2) = 2
             AMA%constr(2, 1, 1:2) = AMA%constr(1, 2, 1:2)
             AMA%constr(2, 2, 1:2) = AMA%constr(1, 1, 1:2)
             print*,'Attention HERE, it is really necessary ?? anisot.f90'
          endif

          !print*,'###################################################################'
          !do i=1,  AMA%nconstr
          !   print*, AMA%iconstr(i), AMA%constr(i, 1, 1:2), AMA%constr(i, 2, 1:2)
          !   write(38, *) AMA%constr(i, 1, 1:2)
          !   write(38, *) AMA%constr(i, 2, 1:2)
          !enddo
          !print*,'###################################################################'


       endif

       close(icon)

       if(state%space%adapt%adapt_level == 0) then
          open(icon,file="ama-constrains", status='unknown')
          do i=1,  AMA%nconstr
             write(icon, *) AMA%constr(i, 1, 1:2)
             write(icon, *) AMA%constr(i, 2, 1:2)
             write(icon,'(x)')
          enddo
          close(icon)
       endif


    else
       AMA%nconstr = 0
    endif



    !! when comented, then some troubles in ANGENER caused by the gfortran optimization??
    !! replace commons by something else ???
    !write(*, '(a6,3i8, 30i6)') 'MAX:',AMA%melem,AMA%mpoin,AMA%mbelm,AMA%maxdeg,ndim,AMA%ipoint, AMA%nbp

    allocate( AMA%x(AMA%mpoin),AMA%y(AMA%mpoin), AMA%lnd(1:AMA%melem, 1:3), &
         AMA%lbn(1:AMA%mbelm,1:2), AMA%ibc(1:AMA%mbelm), AMA%itc(1:AMA%mbelm), &
         AMA%icyc(1:AMA%mpoin, 1:AMA%maxdeg) )
    
    allocate( AMA%xold(1:AMA%mpoin), AMA%yold(1:AMA%mpoin) )
    
    allocate( AMA%lndold(1:AMA%melem, 1:3),&
         AMA%iae(1:AMA%melem, 1:3),&
         AMA%wp(1:AMA%mpoin,1: ndim+1), &
         AMA%w(1:AMA%melem,1:ndim+1),&
         AMA%wpold(1:AMA%mpoin,1:ndim+1), &
         AMA%lnd1(AMA%melem,3), AMA%iae1(AMA%melem,3),&
         AMA%rga(1:AMA%mpoin), AMA%rgb(1:AMA%mpoin), AMA%rgc(1:AMA%mpoin), &
         AMA%nserr(1:AMA%melem*3,1:2),&
         AMA%ibb(1:AMA%mpoin, 1:3+AMA%nconstr),&
         AMA%iba(1:AMA%melem, 1:3), AMA%xb(AMA%ipoint),AMA%yb(AMA%ipoint), &
         AMA%ibpoin(0:AMA%nbp),&
         AMA%ibp(1:AMA%mpoin,1:2), AMA%itrans(1:AMA%melem),&
         AMA%iaegr(1:AMA%melem,1:4) )


    ! local arrays
    allocate(AMA%loc_M(1:3), AMA%loc_A(1:3, 1:4) )

    ! setting of input arrays for ANGENER
    call SetAngenerTriang(  )

    !print*,'########## ANGENER SKIPPED in anisot.f90'
    !!do i=1, grid%nelem
    !!   grid%elem(i)%ama_p = 0.
    !!enddo


    if(.not. metric) then
       call GenerateAMAmetric( )

    else
       ! ANGENER metric, p order is fixed
       AMA%wpold(:, ndim+1) = state%space%deg
    endif



    ! anisotropic mesh adaptation
    call ANGENER_90(metric, ndim)

    !print*,'# ANGENER finished'
    !stop

    !print*,'SW3x1'
    !pause

    !print*, 'Allocation of gridN in anisot.f90 122'



    ! allocate(gridN, source = grid )
    allocate ( MeshAMA_t :: gridN )

    select type (gridN)
      type is ( MeshAMA_t )

      class default
         stop 'other type of mesh'
    end select

    call SetGrid(gridN )

    call EvalAMAparams(gridN)
    ! identical or almost identical grids
    i = state%space%adapt%adapt_level

    id1 = 0
    id2 = -1
    !id1 = max(1, int( 0.005 * state%space%adapt%AMAhistory(i, 2) ))
    !id2 = max(1, int( 0.005 * state%space%adapt%AMAhistory(i, 3) ))
    identical_AMA_grids = .false.

    !if(i >= 2) then
    !   print*,'### Variation of the Shp space_______________________________'
    !   write(*,'(2(a6,i7,a4, i4))') &
    !        'V1:',abs(state%space%adapt%AMAhistory(i, 2) - state%space%adapt%AMAhistory(i-1, 2)),'<?',id1, &
    !        ', V2:',abs(state%space%adapt%AMAhistory(i, 3) - state%space%adapt%AMAhistory(i-1, 3)),'<?',id2
    !   write(*,'(2(a6,i7,a4, i4))') &
    !        'V1:',abs(state%space%adapt%AMAhistory(i, 2) - state%space%adapt%AMAhistory(i-2, 2)),'<?',id1, &
    !        ', V2:',abs(state%space%adapt%AMAhistory(i, 3) - state%space%adapt%AMAhistory(i-2, 3)),'<?',id2
    !   print*,'_____________________________________________________________'
    !endif

    if(i >= 1) then
       if(abs( state%space%adapt%AMAhistory(i, 2) - state%space%adapt%AMAhistory(i-1, 2)) <= id1 .and. &
            abs(state%space%adapt%AMAhistory(i, 3) - state%space%adapt%AMAhistory(i-1, 3)) <= id2) then
          identical_AMA_grids = .true.
       endif
    endif
    if(i >= 2) then
       if( abs(state%space%adapt%AMAhistory(i, 2) - state%space%adapt%AMAhistory(i-2, 2)) <= id1 .and. &
            abs(state%space%adapt%AMAhistory(i, 3) - state%space%adapt%AMAhistory(i-2, 3)) <=id2 ) then
          identical_AMA_grids = .true.
       endif
    endif

    if(identical_AMA_grids) then
       write(*,*)'# The new hp-grid is (almost) identical with the previous ones'
       write(*,*)'# End of computation, state%space%adapt%AMAhistory (level, npoin, nelem, dof):'
       write(*,'(a8,i5,3i8)') 'AMAlevel',i, state%space%adapt%AMAhistory(i, 1:3)
       write(*,'(a8,i5,3i8)') 'AMAlevel',i-1, state%space%adapt%AMAhistory(i-1, 1:3)
       write(*,'(a8,i5,3i8)') 'AMAlevel',i-2, state%space%adapt%AMAhistory(i-2, 1:3)
       write(*,'(a8,i5,3i8)') 'AMAlevel',i-3, state%space%adapt%AMAhistory(i-3, 1:3)
       return
    endif

    print*,'# Passing of data from grid to gridN, deallocation of grid'
!    !! 1st neighbours, 2nd curved
    !call SeekNeighbours(gridN)
    if(state%modelName == 'pedes' ) call DeallocatePedestrianEikonal( )

    call gridN%seekNeighbours( )
    !call seekNeighboursNew(gridN )

    gridN%curved_deg = grid%curved_deg
    call SeekCurvedBoundary(gridN)


    !print*,'SW3x1.a'
    !pause
!
!    print*,'before Reprepare: (in AdaptMesh_angener)'
    !call system('free')

    call ReprepareProblem(gridN, grid)

    !already in reprepare problem
    !call gridN%setSubmesh( lines_file ) ! setConvexSubmesh( lines_file )
!    call grid%plotSubmesh( 1, '../DWR/insideSubmesh.gnu')

    !print*,'AFTER Reprepare:'
    !call system('free')

    !print*,'# PW CONSTANT Interpolation on the new mesh'

    !print*,'# PW polynomial Simple  Interpolation on the new mesh'
    !call SimpleInterpolDGsolution(gridN, grid)


    !print*,'SW3x2'
    !   pause

    !print*,'# PW polynomial Advanced (adaptive) Interpolation on the new mesh'
    !FR already done in ReprepareProblem()
    !call AdvancedInterpolDGsolution(gridN, grid)

   !    print*,'SW3x3'
   !    pause

    !print*,'# PW polynomial Interpolation on the new mesh'
    !call InterpolDGsolution(gridN, grid)

    !call system('free')


    call DeallocateGrid(grid)
    deallocate (grid)

    !call system('free')


    grid => gridN


    ! dealocation of ANGENER arrays
    deallocate(AMA%lndold, AMA%iaegr,  AMA%xold, AMA%yold, &
         AMA%x, AMA%y, AMA%lnd, AMA%lbn, AMA%ibc, AMA%itc, AMA%icyc, &
         AMA%wpold, AMA%w, AMA%wp, AMA%iae, AMA%rga, AMA%rgb, AMA%rgc, &
         AMA%nserr, AMA%xb, AMA%yb, AMA%ibpoin, AMA%ibp, &
         AMA%iba, AMA%ibb, AMA%lnd1, AMA%iae1, AMA%itrans, AMA%iwall)
    deallocate(AMA%loc_M, AMA%loc_A )

    if(AMA%nconstr > 0) deallocate( AMA%constr, AMA%iconstr)

    !print*,'          end of ANGENER'
    !call system('free')

    !deallocate(AMA)

    !! plotting of recomputed solution
    !state%space%adapt%adapt_level = state%space%adapt%adapt_level + 1
    !call WriteProgressOutput( 'ST' )
    !state%space%adapt%adapt_level = state%space%adapt%adapt_level - 1

    !do i=1,grid%nelem
    !   write(31,*) grid%elem(i)%xc(:), grid%elem(i)%deg
    !enddo



    close(AMA%ifig1)
    !print*,' # end subroutine AdaptMesh_Angener'
    !stop

  end subroutine AdaptMesh_Angener


  !> reading of data from file paramet
  !> parameter metric:
  ! =1 metric generated by F77 ANGENER, =0 in AnisotErrorEstimates
  subroutine ReadingOfParamets(metric )
    logical, intent(in) :: metric
    integer :: isid, j

    if(metric) then
       isid = 58

       open(isid,file='paramet',status='old')

       read(isid,*) AMA%ityp, AMA%icrack
       read(isid,*) AMA%ndim
       read(isid,*) AMA%ifv
       read(isid,*) AMA%pos, AMA%posW, AMA%pos1
       read(isid,*) AMA%numel
       read(isid,*) AMA%epsilon1
       read(isid,*) AMA%p

    else

       AMA%ityp = 0
       AMA%icrack = 0
       AMA%ndim = ndim
       AMA%ifv = 1

       AMA%maximal_angle_optimal = 165. / 180 * pi
       !AMA%maximal_angle_optimal = 170. / 180 * pi

       !AMA%maximal_angle = 150. / 180 * pi
       AMA%maximal_angle = 165. / 180 * pi
       !AMA%maximal_angle = 170. / 180 * pi
       !AMA%maximal_angle = 200. / 180 * pi ! NO constrain


       !print*,'AMA%maximal_angle_optimal = ', AMA%maximal_angle_optimal

       AMA%pos1 = 0.0

       if(state%modelName == 'scalar' .or.state%modelName == '2eqs' .or.state%modelName == 'porous') then
          ! used for paper ESCO14,  AMAtdp for moving Front
          !AMA%pos  = 0.0005
          !AMA%posW = 0.0005

          AMA%pos  = 0.001  ! used for ESCO18
          AMA%posW = 0.001

          !AMA%pos  = 0.002
          !AMA%posW = 0.002

          !AMA%pos  = 0.005
          !AMA%posW = 0.005

          !AMA%pos  = 0.01
          !AMA%posW = 0.01

          !AMA%pos  = 0.02
          !AMA%posW = 0.02

          !AMA%pos  = 0.025
          !AMA%posW = 0.025

          !AMA%pos  = 0.05
          !AMA%posW = 0.05

          !AMA%pos  = 0.1
          !AMA%posW = 0.1

          if( state%modelName == 'porous' ) then
             AMA%pos = 0.02
             AMA%posW = 0.02
          endif


          if(state%space%adapt%adapt_type == 'Ihp') then
             AMA%pos = 0.12   ! 0.12
             AMA%posW = 0.15  ! 0.15
             !AMA%posW = 0.15
          endif
       else
          ! anis_hp for NSe
          AMA%pos = 0.0001       ! criterion 2|K| /sum l^2 > pos
          AMA%posW = 0.0001

          !AMA%pos = 0.025       ! criterion 2|K| /sum l^2 > pos
          !AMA%posW = 0.025

          !AMA%pos = 0.01       ! criterion 2|K| /sum l^2 > pos
          !AMA%posW = 0.01

          AMA%pos1 =  AMA%posW *1.   ! criterion (l1 + l2 )/l3 > (1 + pos1)

          ! used for paper ESCO14,  AMAtdp for shock-vortex interaction
          !AMA%pos = 0.005
          !AMA%posW = 0.01

          if(state%space%adapt%adapt_type == 'Ihp') then
             !AMA%pos = 0.12
             !AMA%posW = 0.15

             AMA%pos = 0.02
             AMA%posW = 0.02
          elseif( state%space%adapt%adapt_method == 'ANI')  then
             AMA%pos = 0.01
             AMA%posW = 0.01
          endif
       endif


       AMA%numel = 200
       AMA%epsilon1 = 1E+20
       AMA%p = 1E+05

    endif

    write(*,'(a18, 4(a10, es9.2))') 'ANGENER params:', &
         ', A%pos =', AMA%pos, &
         ', A%posW =', AMA%posW, &
         ', A%pos1 =', AMA%pos1, &
         ', A%max_a =', AMA%maximal_angle_optimal

    if(AMA%ndim /= ndim) then
       print*,' nonconsistency with ndim in *.ini and paramet'
       print*,ndim, AMA%ndim
       stop
    endif


    if( AMA%ifv  /= 1) then
       print *,'Sorry, incorect number ifv in paramet, only ifv == 1 is possible'
       print *,'Ifv in paramet =',AMA%ifv
       stop
    endif


    if(AMA%ityp .ge. 3) then
       read(isid,*) AMA%Re
       read(isid,*) AMA%xte(1,1)
       read(isid,*) AMA%xte(1,2)
       read(isid,*) AMA%xte(2,1)
       read(isid,*) AMA%xte(2,2)
       read(isid,*) AMA%iwa

       allocate(AMA%iwall(1:AMA%iwa) )

       do j=1,AMA%iwa
          read(isid,*) AMA%iwall(j)
       enddo
    else
       AMA%Re = 0.
       AMA%xte(1,1) = 0.
       AMA%xte(1,2) = 0.
       AMA%xte(2,1) = 1.
       AMA%xte(2,2) = 0.
       AMA%iwa = 0
       allocate(AMA%iwall(1:1) )  ! only for a unique deallocation
    endif

    if(metric) close(isid)

  end subroutine ReadingOfParamets



  subroutine ANI_Set_Params( )
    integer :: i

    AMA%npoin = grid%npoin
    AMA%nelem = grid%nelem
    AMA%nbelm = grid%nbelm

    AMA%xper(1:2, 1:2) =  grid%xper(1:2,1:2)
    AMA%iper(1:2, 1:2) =  grid%iper(1:2,1:2)

    AMA%maxdeg = 35

    AMA%nbp = curvbound%nbp
    AMA%ipoint = 0
    if(curvbound%nbp > 0) then
       do i= 1, curvbound%nbp
          AMA%ipoint = AMA%ipoint + size(curvbound%CBP(i)%pnt, 1)
       enddo
    endif


  end subroutine ANI_Set_Params


  !> filling of arrays for ANGENER
  subroutine SetAngenerTriang(  )
    integer :: i, k, k1, iflag, j, j1, ipoin
    integer, dimension(:,:), pointer :: lnd, lbn
    integer, dimension(:), pointer :: ibc, itc

    lbn => AMA%lbn(1:AMA%mbelm,1:2)
    ibc => AMA%ibc(1:AMA%mbelm)
    itc => AMA%itc(1:AMA%mbelm)
    lnd => AMA%lnd(1:AMA%melem, 1:3)

    !print*,'AMA: nodes coordinates'
    do k=1,AMA%npoin
       AMA%x(k) = grid%x(k,1)
       AMA%y(k) = grid%x(k,2)
       !read (ivlt,*) x(k),y(k)
    enddo

    !print*,'AMA: triangles'
    do k=1,AMA%nelem
       lnd(k,1:3) = grid%elem(k)%face(idx, 1:3)
       !read (ivlt,*) lnd(k,1),lnd(k,2),lnd(k,3)
    enddo

    AMA%npoinold = AMA%npoin
    AMA%nelemold = AMA%nelem

    AMA%ibp(1:AMA%npoin,1:2) = 0
    AMA%ibb(1:AMA%npoin, 3 ) = 0   ! interior constraint
    AMA%ibb(1:AMA%npoin, 4: ) =0   ! interior lines for corner nodes, AMA%ibb(1:AMA%npoin, 3 ) = - 1

    AMA%xold(1:AMA%npoin) = AMA%x(1:AMA%npoin)
    AMA%yold(1:AMA%npoin) = AMA%y(1:AMA%npoin)

    AMA%lndold(1:AMA%nelem, 1:3) = lnd(1:AMA%nelem, 1:3)

!      iflag = 2
    iflag = 0

    AMA%nbc = 0   !     maximal index of components

    !print*,'AMA: boundary edges'
    do k=1,AMA%nbelm
       !read (ivlt,*) lbn(k,1),lbn(k,2),ibc(k)
       lbn(k,1:2) = grid%b_edge(k)%lbn(1:2)
       ibc(k) = grid%b_edge(k)%ibc
       AMA%nbc = max(AMA%nbc, ibc(k) )
    enddo

    do k=1,AMA%nbelm
       if(iflag .eq. 2) then
          print *,'I am here, TRY it!!!!'
          do i=1,AMA%iwa
             if(ibc(k) .eq. AMA%iwall(i) ) then
                AMA%ibp(lbn(k,1),1) = -1
                AMA%ibp(lbn(k,2),1) = -1
                goto 10
             endif
11           continue
          enddo
       endif
       if(k .eq. 1) then
          k1 = AMA%nbelm
       else
          k1 = k -1
       endif
       if(ibc(k) .ne. ibc(k1)) AMA%ibb(lbn(k,1),3) = -1
       !if(AMA%ibb(lbn(k,1),3) == -1) &
       !     write(*,'(a8, 2i5,a2, 5i5, 2es12.4)') &
       !     'IBC C', lbn(k, 1:2), '|',k, k1, ibc(k), ibc(k1), AMA%ibb(lbn(k,1),3), &
       !     AMA%x(lbn(k,1)), AMA%y(lbn(k,1))
10     continue
    enddo
    !close(ivlt)

    !do i=1,AMA%npoin
    !   if(AMA%ibb(i,3) .eq. -1) write(95,*) AMA%x(i),AMA%y(i),i,'&^%'
    !enddo

    !print*,'AMA: seeking of itc array'
    do  k=1,AMA%nbelm
       do  i=1,AMA%nelem
          do  j=1,3
             j1 = mod(j,3) +1
             if(lnd(i,j) .eq. lbn(k,1) .and. lnd(i,j1) .eq. lbn(k,2) ) then
                itc(k) = i
                goto 681
             endif
          enddo
       enddo
681    continue
    enddo


    !print*,'AMA: element averaging', AMA%ityp
    if(AMA%ityp .ne. 0 .and. AMA%ityp .ne. 3 .and. AMA%ityp .ne. -1 .and. AMA%ityp .ne. -10) then
       do i=1,AMA%nelem
          call  Eval_aver_w_Elem(grid%elem(i), AMA%w(i,1:ndim) )
       enddo
    else
       ! uniform metric
       AMA%w(:,:) = 1.
    endif

    !print*,'AMA: curved boundary'
    AMA%ibpoin(0) = 0
    !read(iprof,*) nbppp
    do i=1, AMA%nbp
       !read(iprof,*) ipoin
       ipoin = size(curvbound%CBP(i)%pnt, 1)
       AMA%ibpoin(i) = AMA%ibpoin(i-1) + ipoin
       do j=1,ipoin
          !read(iprof,*) xb(AMA%ibpoin(i-1)+j),yb(AMA%ibpoin(i-1)+j)
          AMA%xb(AMA%ibpoin(i-1)+j) = curvbound%CBP(i)%pnt(j,1)
          AMA%yb(AMA%ibpoin(i-1)+j) = curvbound%CBP(i)%pnt(j,2)
       enddo
    enddo
    !close(iprof)

    !print *,'The reading of datas done'

  end subroutine SetAngenerTriang

  ! setting the data from a ANGENER generated grid to new grid
  subroutine SetGrid(gridN )
    class( mesh), intent(inout) :: gridN
    real, dimension(:,:), pointer :: w
    integer, dimension(:), pointer :: itrans
    class(element), pointer :: elem
    integer :: i, k, Tdeg, dof

    itrans => AMA%itrans(1:AMA%melem)
    w => AMA%w( 1:AMA%melem, 1:ndim+1)

    gridN%npoin = AMA%npoin
    gridN%nelem = AMA%nelem
    gridN%nbelm = AMA%nbelm
    gridN%nbc = AMA%nbc

    !print*,'>>>>>> A1', gridN%npoin, gridN%nelem, gridN%nbelm
    gridN%periodic = grid%periodic
    gridN%xper(1:2, 1:2)  = grid%xper(1:2, 1:2)
    gridN%iper(1:2, 1:2)  = grid%iper(1:2, 1:2)

    allocate(gridN%x(gridN%npoin,2) )
    allocate(gridN%xcur(1:gridN%npoin, 1:nbDim))

    !print*,'>>>>>> A2'

    do i=1,gridN%npoin
       gridN%x(i,1) = AMA%x(i)
       gridN%x(i,2) = AMA%y(i)
    enddo

    allocate(gridN%elem(1:gridN%nelem) )

    !print*,'>>>>>> A5'

    do i=1, gridN%nelem
      elem => gridN%elem(i)
      elem%flen = 3
      elem%type = 3
      elem%HGnode = .false.
      elem%RGtype  = 'N'
      elem%RGhistory = 0
      elem%RGlevel = 0
      elem%RGindex = 0

      !print*,'>>>>>> B1',i,max_face,elem%flen, w(itrans(i), :)

      allocate(elem%face(1:max_face, 1:elem%flen))
      elem%face(idx,1:3) = AMA%lnd1(i, 1:3)
      elem%i = i

      !elem%deg = state%space%deg
      elem%deg =  int( w(itrans(i), ndim + 1) + 0.5)
      elem%Tdeg = state%time%deg

      call elem%initElementDof()

      elem%hsplit = 0
      elem%psplit = 0

      elem%to_recompute = .true.


      Tdeg = state%time%deg+1
      if(state%time%disc_time == 'STDG') Tdeg =1

      !print*,'>>>>>> B2',i,max_face,elem%flen, Tdeg, elem%deg, elem%dof, ndim

      allocate(elem%w(0:Tdeg, 1:ndim* elem%dof) )

      !print*,'>>>>>> B3',i,max_face,elem%flen

      call SetOneElementConstantIC(elem, w(itrans(i), 1:ndim) )

      !print*,'>>>>>> B4',i,max_face,elem%flen

      if (state%time%disc_time == 'STDG') then
         allocate(elem%wST(1:ndim, 1:elem%dof, 1:elem%Tdof ) )
         allocate(elem%wActual(1:ndim, 1:elem%dof_plus), source = 0.0 )
         allocate(elem%rhsST(1:ndim, 1:elem%dof_plus, 1:elem%Tdof_plus ))
         allocate( elem%wSTfin(1:ndim, 1:elem%dof) ) !in problem.f90
     endif

      ! taking values of the solution
      !allocate(gridN%elem(i)%wS(1:1,1:ndim))

      ! passing of the solution from the Amesh to mesh
      !do j=1,ndim
      !   gridN%elem(i)%wS(1, j) = state%BC(1)%ww(j) ! doscasne, jen aby to probehlo
      !enddo
   enddo

   if ( state%space%estim_space == 'DWR' ) then

      call allocateZST( gridN, .false.)
      if (state%dwr_p) then
        call allocateWSTplus( gridN)
        call allocateZSTplus( gridN)
      end if
   end if

   !print*,'>>>>>> A6 at anisot.f90'

   ! boundary segments

   allocate(gridN%b_edge(1:gridN%nbelm) )
   do i=1, gridN%nbelm
      gridN%b_edge(i)%lbn(1:2) = AMA%lbn(i, 1:2)
      gridN%b_edge(i)%ibc = AMA%ibc(i)
   enddo

   gridN%ElemSupports = .false.
   
   !print*,'>>>>>> A7'

 end subroutine SetGrid


 !> generate metric for AMA
 !> recomputation from elem%rgabc to ANGENER's arrray rga, rgb, rgc
 !> evaluation of wp in mesh vertexes !!!
 subroutine GenerateAMAmetric( )
   real, dimension(:), allocatable :: supp, wi
   class(element), pointer :: elem
   integer :: nelem, npoin, i,j,k, imt, is
   real, dimension(1:3) :: rgabc
   real :: rnorm, epsilon2, eps1
   character(len=15) :: file1, file2
   character(len=5) :: ch5
   real, dimension(:), pointer :: rga, rgb, rgc
   real, dimension(:,:), pointer :: wp, wpold

   wp    => AMA%wp(   1:AMA%mpoin,1:ndim+1)
   wpold => AMA%wpold(1:AMA%mpoin,1:ndim+1)

   rga => AMA%rga( 1:AMA%mpoin )
   rgb => AMA%rgb( 1:AMA%mpoin )
   rgc => AMA%rgc( 1:AMA%mpoin )

   npoin = grid%npoin
   !nelem = grid%nelem

   allocate (supp(1:npoin),  wi(1:ndim ) )

   ! weighted average in nodes
   wp(:,:) = 0.
   wpold(:,:) = 0.
   supp(:) = 0.
   rga(:) = 0.
   rgb(:) = 0.
   rgc(:) = 0.

   do i=1, grid%nelem
      elem => grid%elem(i)

      call  Eval_aver_w_Elem(elem, wi(1:ndim) ) ! w average on elements

      do j=1, elem%flen
         k = elem%face(idx, j)
         ! solution in points
         wp(k, 2:ndim+1) = wp(k, 2:ndim+1) + wi(1:ndim) * elem%area

         wpold(k, ndim+1) = wpold(k, ndim+1) + (elem%deg+elem%ama_p) * elem%area
         !wpold(k, ndim+1) = wpold(k, ndim+1) + (elem%deg+elem%psplit) * elem%area

                  ! elem%rg* already generated(e.g., Hessian matrix)
         rga(k) = rga(k) + elem%rgabc(1) * elem%area
         rgb(k) = rgb(k) + elem%rgabc(2) * elem%area
         rgc(k) = rgc(k) + elem%rgabc(3) * elem%area

         ! uniform metric
         !rga(k) = 0.
         !rgb(k) = 0.
         !rgc(k) = 0.

         ! weights
         supp(k) = supp(k) + elem%area

      enddo
   enddo

   ! normalization by weights
   do k=1,npoin
      wp(k, 2:ndim+1) = wp(k, 2:ndim+1) / supp(k)
      wp(k, 1) = wp(k, 2)
      wpold(k, ndim+1) = wpold(k, ndim+1)  / supp(k)
      rga(k) = rga(k) / supp(k)
      rgb(k) = rgb(k) / supp(k)
      rgc(k) = rgc(k) / supp(k)

      !!print*,'.ded39d3',k,  wpold(k, ndim+1)
   enddo

   ! file output



   !do i=1, grid%nelem
   !   elem => grid%elem(i)
   !   write(200+state%space%adapt%adapt_level, *)  elem%xc(1:2), elem%deg+elem%psplit

   !   do j=1, elem%flen + 1
   !      k = elem%face(idx, mod(j, elem%flen)+ 1)
   !      write(200+state%space%adapt%adapt_level, *)  &
   !           grid%x(k, 1:2), wp(k, 1:ndim + 1), wpold(k, ndim + 1), i,j,k
   !   enddo
   !   write(200+state%space%adapt%adapt_level, '(x)' )
   !   write(200+state%space%adapt%adapt_level, '(x)' )
   !   write(200+state%space%adapt%adapt_level, '(x)' )
   !enddo


   ! "lower" and "upper" regularization of metrices
   !rnorm = 1.29903810567 *AMA%numel/state%space%domain_volume
   !print*,'!!!!', AMA%numel, state%space%domain_volume, rnorm
   !epsilon2 = AMA%epsilon1/ AMA%p
   !epsilon2 = max (1.,epsilon2)

   !do i=1, npoin
   !   eps1 = AMA%epsilon1/(epsilon2+max(rga(i),rgc(i)))
   !
   !   !rgb(i) = 0.
   !   !rga(i) = (1. + eps1*abs(rga(i)))*rnorm
   !   !rgc(i) = (1. + eps1*abs(rgc(i)))*rnorm
   !
   !   !write(24, '(a6, i5, 5es14.6)') 'metr:',i, rga(i), rgb(i), rgc(i), rga(i)*rgc(i) - rgb(i)**2
   !enddo

   ! drawing of ellipses
   !file1 = 'metrixV00000'
   !is = 0
   !if(state%space%adapt%adapt_level > 0) is = int(log(1. * state%space%adapt%adapt_level)/log(10.))
   !
   !write( ch5, '(i5)' ) state%space%adapt%adapt_level  ! change the format if num_size /= 5 !!!
   !file1(12-is: 12)  = ch5(5-is:5)
   !imt = 118
   !
   !open(imt,file=file1,status='unknown')
   !do  i=1,npoin
   !   rgabc(1)  = rga(i)
   !   rgabc(2)  = rgb(i)
   !   rgabc(3)  = rgc(i)
   !   call DrawEllips(imt, rgabc(1:3), grid%x(i,1:2) )
   !enddo
   !close(imt)


   ! FILE OUTPUT
   ! file1 = 'metrixF00000'
   ! is = 0
   ! if(state%space%adapt%adapt_level > 0) is = int(log(1. * state%space%adapt%adapt_level)/log(10.))

   ! write( ch5, '(i5)' ) state%space%adapt%adapt_level  ! change the format if num_size /= 5 !!!
   ! file1(12-is: 12)  = ch5(5-is:5)
   ! imt = 118

   ! open(imt,file=file1,status='unknown')
   ! do  i=1,npoin
   !    write(imt,*), rga(i), rgb(i), rgc(i)
   ! enddo
   ! close(imt)

   deallocate(supp, wi)

 end subroutine GenerateAMAmetric




end module anisotropic
