module target_functional_mod
  use blocks_integ
  use eval_sol
  use eval_jumps
  use lapack_oper
  use mesh_mod
  use mesh_oper
  use model_oper
  use paramets
  use st_interpol
  use stdgm_mod
  use weight_fun_mod
  use set_solution
  !  use sort_mod

   implicit none

   type, public,abstract :: Abstr_Target_functional_t

   contains
      procedure(initTarFunc), deferred :: init
      procedure(findSupp), deferred :: findSupp
      !      procedure(evalWeightFunction), deferred :: evalWeightFunction
!      final :: delete_Target_functional
   end type Abstr_Target_functional_t

   abstract interface
   subroutine initTarFunc( this, xy_coord, grid, weight_i )
      import :: Abstr_Target_functional_t
      import :: nbDim
      import :: mesh
      class(  Abstr_Target_functional_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i
   end subroutine initTarFunc

   subroutine findSupp( this, grid )
      import ::  Abstr_Target_functional_t
      import ::  mesh
      class(  Abstr_Target_functional_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
   end subroutine findSupp

   end interface

   type, extends( Abstr_Target_functional_t ) :: Target_functional_t
      integer :: id ! id of the functional
      character(len=20) :: name ! name of the target quantity
      real :: Ju ! target quantity of the computed solution
      real :: Ju_BiCG ! target quantity resulting from BiCG method
      real :: Ju_exact
!      real, dimension(:,:), allocatable :: dJ_phi ! array - for linear J equals J(phi), for nonlinear J'[u](phi)
      logical :: linear ! target functional is linear -> no need of J'(u)(.)
      logical :: boundary ! target functional involves boundary integration
      ! not boundary -> target functional involves volume integration
      logical :: integralAver ! the functional is divided by |\om_J|
      logical :: time_dependent ! target functional involves time integration
      integer, allocatable, dimension(:) :: supp ! indices of elements in support of the target functional J
      integer, allocatable, dimension(:) :: suppFace ! 1:nelemSupport
      !  - local index of the face in the element with index J%supp(:) in supp(J) ! FOR BOUNDARY J ONLY
      integer :: isupp ! # of elements in support, i.e. size of the supp array
      real :: vol_supp ! volume of the support of the functional
      real, allocatable, dimension(:) :: xy_coord ! coordinates of the point value in tarFunc for id=, or x_ref for Momentum coeff
      class( mesh ), allocatable :: grid ! only for point val
      real :: eps                              ! diameter for pointvalue of tarFunc for id=
      integer :: iWeight ! weighting function \psi to be used as J(v) = int v*\psi, i.e. j_\om or j_gamma resp.

   contains
      procedure :: clean => cleanTargetFunctional ! nulify all parts which change after adaptation
      procedure :: init => init_tarFunc
      procedure :: initGrid => initGrid_tarFunc ! initializes the dual mesh - support of the functional
      procedure :: findSupp => findSupp_tarFunc
      procedure :: evalWeightFunction !=> evalWeightFunction_tarFunc
      procedure :: evalWeightFunctionVector
      procedure :: computeJu_exact
      procedure :: computeJu
      procedure :: delete => delete_Target_functional

   end type Target_functional_t

   !> integration of the solution over a subdomain \f$ \Omega_s \f$ of \f$ \Omega\f$
   !> dwr%id = 4
   type, extends( Target_functional_t ) :: U_over_subdomain_t

   contains
      procedure :: computeJu => computeJu_u_over_subdomain
      procedure :: computeJu_exact => computeJu_exact_u_over_subdomain
!      procedure :: evalWeightFunction => evalWeightFunction_u_over_subdomain
      procedure :: findSupp => findSupp_u_over_subdomain
      procedure :: init => init_u_over_subdomain

   end type U_over_subdomain_t

   !> integration of the solution over a subdomain \f$ \Omega_s\f$  of \f$ \Omega\f$
   !> dwr%id = 5
   type, extends( Target_functional_t ) :: dudx_t
      integer :: dx ! 1:nbDim - derivative with respect to x (1) or y (2)

   contains
      procedure :: computeJu => computeJu_dudx
      procedure :: computeJu_exact => computeJu_exact_dudx
!      procedure :: evalWeightFunction => evalWeightFunction_dudx
      procedure :: findSupp => findSupp_dudx
      procedure :: init => init_dudx

!      procedure :: setRHS => setRHS_dudx

   end type dudx_t

   !> value of the solution over part of the boundary (Neumann)
   !> dwr%id = 1
   !> will not work if one of the elements contains more than ONE edge in supp(J)
   type, extends( Target_functional_t ) :: BoundaryValue_t

   contains
      procedure :: computeJu => computeJu_BoundaryValue
      procedure :: computeJu_exact => computeJuExact_BoundaryValue
!      procedure :: evalWeightFunction ! => evalWeightFunction_BoundaryValue
      procedure :: findSupp => findSupp_BoundaryValue
      procedure :: init => init_BoundaryValue

   end type BoundaryValue_t

   !> flux of the solution over part of the boundary (Dirichlet)
   !> J(v) = \int j_D * A(v,grad(v)*n
   !> IT HAS TO BE MODIFIED IN ORDER TO GET ADJOINT CONSISTENCY
   !> JJ(v) = J(v) - \sum_{\partial K \cap Gamma_D) ( \int_{-}( \sigma(v-u_d)*j_D ) + \int_{+}((\sigma + b*n)(v-u_D)j_D ) )
   !> dwr%id = 2
   !> will not work if one of the elements contains more than ONE edge in supp(J)
   !> then we solve DP: a(v,z) = JJ'_u(v) , since JJ'_u(u-u_h) = JJ(u) - JJ(u_h)
   type, extends( Target_functional_t ) :: BoundaryFlux_t
   contains
      procedure :: computeJu => computeJu_BoundaryFlux
      procedure :: computeJu_exact => computeJuExact_BoundaryFlux
      procedure :: findSupp => findSupp_BoundaryFlux
      procedure :: init => init_BoundaryFlux


   end type BoundaryFlux_t

   !> Drag and Lift Computation for the Euler equations
   !> dwr%id = 10 - Drag, 11 - Lift, 12 - Momentum
   type, extends( Target_functional_t ) :: DragAndLift_t
      !integer :: dlm ! 1 - drag, 2 - lift, 3 - momentum
      character(len=20) :: dlm ! 'drag', 'lift', 'momentum'
      real, dimension(:), allocatable :: theta ! (0,theta1, theta2,0)

   contains
      procedure :: computeJu => computeJu_DragAndLift
      procedure :: computeJu_exact => computeJu_exact_DragAndLift
      procedure :: findSupp => findSupp_DragAndLift
      procedure :: init => init_DragAndLift
      procedure :: setTheta
      procedure :: evalWeightFunctionVector => evalWeightFunctionVector_Euler
      procedure :: testBC_Euler_W

   end type DragAndLift_t


contains

   subroutine init_tarFunc( this, xy_coord, grid, weight_i )
      class( Target_functional_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      stop 'init_tarFunc should not be called - abstract'

   end subroutine init_tarFunc

   subroutine initGrid_tarFunc( this, gridfile )
      class( Target_functional_t ), intent(inout) :: this
      character(len=50), intent(in) :: gridfile     ! file with dual mesh

      stop 'initGrid_tarFunc should not be called - abstract'

   end subroutine initGrid_tarFunc

   subroutine findSupp_tarFunc( this , grid)
      class( Target_functional_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid

      stop 'findSupp_tarFunc should not be called - abstract'

   end subroutine findSupp_tarFunc

   subroutine cleanTargetFunctional( this )
      class( Target_functional_t ) :: this

      this%Ju = 0.0 ! target quantity of the computed solution
      this%Ju_exact = 0.0 ! ??
      if (allocated(this%supp)) &
         deallocate(this%supp)

      this%isupp = 0
      !this%vol_supp = ??
      !grid ! only for point val

   end subroutine cleanTargetFunctional

   subroutine delete_Target_functional( this )
    class( Target_functional_t ) :: this

    if (allocated( this%supp ) ) &
      deallocate( this%supp )
    if (allocated( this%xy_coord ) ) &
      deallocate( this%xy_coord )
    if (allocated( this%grid ))then
      deallocate( this%grid )
      print*, 'Control deallocation of DWR%grid.'
    endif

   end subroutine delete_Target_functional

   subroutine computeJu_exact( this, grid )
      class( Target_functional_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid

      stop 'computeJu_exact should not be called - abstract'

   end subroutine computeJu_exact

   subroutine computeJu( this, grid )
      class( Target_functional_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid

      stop 'computeJu should not be called - abstract'

   end subroutine computeJu



   !> evaluate the value of J_\om, or J_D or J_N
   !> same routine for all types of functionals
   !> time independent - for stationary problems only !!!
   function evalWeightFunction( this, x, t) result(f)
    class( Target_functional_t ),intent(in) :: this
    real, dimension(1:nbDim), intent(in) :: x
    real, intent(in) :: t
    real, dimension(1:ndim) :: f
    real :: aver

    f(1:ndim) = 0.0

        ! divide the functional value by |\om_J|
    if (this%integralAver) then
       aver = this%vol_supp
       !print*,'AVER CHAHGED'
       !aver = aver / 12 * 0.04  ! scaling for Ajay functional
       !print*,'aver = ', aver


    else
      aver = 1.0 ! do not compute average
    endif

    !print*, 'HERE evalWeightFunction_u_over_subdomain', aver, this%vol_supp, this%iWeight

    !if (t == state%time%FinTime) then
    f( 1 ) = weightFun( x, this%iWeight, state%model%Re1 ) / aver
!    else
!      stop 'evalWeightFunction: CONTROL nonzero only in final time'
!    endif

   end function evalWeightFunction

   !> compute weight function in quadrature nodes
   !> for Euler case
   !> NOT TESTED: for scalar case: ie > 0 index of an edge, ie = 0 volume
   !>
   subroutine evalWeightFunctionVector( this, elem, ie, Qdof, fi )
        class( Target_functional_t ), intent(inout) :: this
        class( element ), intent(inout) :: elem
        integer, intent(in) :: ie ! index of the edge
        integer, intent(in) :: Qdof ! quad nodes on the edge
        real, dimension(1:Qdof, 1:ndim), intent(inout) :: fi
        integer :: j
        real :: t

        t = state%time%finTime

        if (ie == 0) then !volume
          if (Qdof /= elem%Qdof) &
            stop "wrong Qdof in evalWeightFunctionVector"
        else !boundary
          if (Qdof /= elem%face(fGdof, ie )) &
            stop "wrong Qdof in evalWeightFunctionVector"
        end if
        do j = 1, Qdof
          fi(j, 1:ndim) = evalWeightFunction( this, elem%xi( ie,j,1:nbDim), t)
        end do

        print*, "evalWeightFunctionVector for SCALAR case was not tested!"
   end subroutine evalWeightFunctionVector





!!!!!!!!!!! BOUNDARY VALUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


   subroutine init_BoundaryValue( this, xy_coord, grid, weight_i)
      class( BoundaryValue_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      this%linear = .true.
      this%boundary = .true.
      this%time_dependent = .false.
      this%name = 'BoundaryValue'

      this%iWeight = weight_i !
!      1 ==constant 1
!      2 ! weight fun psi = y
!      3 ! weight fun psi = y*sin
!      print*, 'We use weight function 1/y in target functional over dOm!'
      this%integralAver = .false. !do not divide by |\om_J|
   end subroutine init_BoundaryValue

   !> marks the elements of the support of J, counts them and computes also volume of this support (edges)
   !> support = element which have iSubmesh == 1
   subroutine findSupp_BoundaryValue( this, grid )
      class( BoundaryValue_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
      class( element ), pointer :: elem
      integer :: nelem, i, j
      integer, dimension(:), allocatable :: temp_supp, temp_suppFace
      real :: area

      nelem = grid%nelem
      j = 0
      area = 0.0

      allocate( temp_supp(1:nelem), source = -1 )
      allocate( temp_suppFace(1:nelem), source = -1 )

      do i = 1, nelem
         elem => grid%elem(i)
         if (elem%iSubMesh == 1) then
            j = j+1
            ! add the index of the element with support
            temp_supp(j) = i
            ! add the local index of the edge in the support
            temp_suppFace(j) = elem%iSubmeshFace
            ! add the length of the edge
            area = area + elem%dn( elem%iSubmeshFace )
         end if
      end do !i

      this%isupp = j

      if (allocated(this%supp) ) &
         deallocate(this%supp)

      if (allocated( this%suppFace) ) &
         deallocate(this%suppFace)

      allocate( this%supp(1:this%isupp), source = 0 )
      allocate( this%suppFace(1:this%isupp), source = 0)
      this%supp(1:this%isupp) = temp_supp(1:this%isupp)
      this%suppFace(1:this%isupp) = temp_suppFace(1:this%isupp)

      deallocate( temp_supp , temp_suppFace )

      this%vol_supp = area !area ! NOT TESTED, NOT USED
      !print*, 'Volume of the support of the boundary functional:', area


   end subroutine findSupp_BoundaryValue

   subroutine computeJuExact_BoundaryValue( this, grid )
      class( BoundaryValue_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time
      integer :: iFace

      ! for some problems the precomputed solution is more precise even with known solution
      if ( state%model%known_sol ) then

         !      print*, 'F@R: computeJu_exact - set time ( temporarily finTime )'
         time = state%time%finTime

         this%Ju_exact = 0.0

         do i = 1, this%isupp
            elem => grid%elem( this%supp(i) )
            iFace = this%suppFace(i) ! local index of the edge we need to integrate over

            Qdof = elem%face(fGdof, iFace )

            allocate( wi(1:Qdof, 1:ndim), source = 0.0 )
            allocate( f(1:Qdof, 1:ndim), source = 0.0 )

            ! eval w, f in integ nodes, 0 - endpoint
            do j = 1, Qdof
               ! integ nodes, time, w, ityp = 1 (exact solution)
               call Set_Model_Data( elem%xi( iFace ,j,1:nbDim) , time , wi(j,1:ndim), 1 )
               ! xi on the face !
               f(j,1:ndim) = this%evalWeightFunction(elem%xi( iFace ,j,1:nbDim), time)
            enddo

            !write(*,'(a8,2i5, 400es12.4)') 'HJu 1:', i, Qdof, f(1:Qdof, 1:ndim)
            !write(*,'(a8,2i5, 400es12.4)') 'HJu 2:', i, Qdof, wi(1:Qdof, 1:ndim)
            !print*

            f(1:Qdof,1:ndim) = f(1:Qdof, 1:ndim) * wi(1:Qdof, 1:ndim)

            ! integrate over the edge iFace - normalizes by the length of the edge
            call IntegrateFunctionEdge( elem, iFace, f(1:Qdof, 1), local_Ju )
            this%Ju_exact = this%Ju_exact + local_Ju

            deallocate(f, wi)

         end do !i
      else if (state%model%icase == 63 ) then   ! battery, approximate value
         this%Ju_exact = 2.731436273280E+03 ! +-5E-00, AMAhp, DOF = 311429

      else if (state%model%icase == 74 ) then
         ! idiff =15, iconv = 13, iexact = 38, delta = 0.05, gamma = 0.05, val1 = 1.E-9
         ! subgrid 0.625 - 0.75
         this%Ju_exact = 0.048544 ! AMAhp, 0.04865 !uniform P3 with 65000 elements

      else if (state%model%icase == 76 ) then
         ! gamma = 0.01, val1 = 0.0, delta = 0.05
!         this%Ju_exact = 1.17997900 ! P4, uniform 16 000
!         this%Ju_exact = 1.17997895 ! P4, adaptive 4 000
         !idiff = 15 iconv = 17 iexact = 71 this%ireac = 1 (param1=3),
         this%Ju_exact = 0.32402676943309316 ! AMAh P4 , 5202 elements  ! Harriman
         this%Ju_exact = 0.32402676958 ! +- 5E-11  AMAhp ,43000 DoF  ! Harriman

      else if ( state%model%icase == 80  ) then ! Carpio 2013
         ! J_3(u)
         !this%Ju_exact = 0.074082   ! +- 1E-6  HG tests,      Carpio  0.0739210
         this%Ju_exact = 0.07408122  ! +- 1E-8     AMA MULTI-tests    Carpio  0.0739210 SISC
         this%Ju_exact = 0.07408120962  ! +- 1E-10     AMA MULTI-tests    Carpio  0.0739210

      else if ( state%model%icase == 81  ) then ! Carpio 2013
         ! J_2(u)
         !this%Ju_exact = 3.970304   ! +- 1E-6     HG tests,      Carpio  3.96992104
         !this%Ju_exact = 3.970295   ! +- 1E-6     AMA tests,     Carpio  3.96992104
         !this%Ju_exact = 3.970300   ! +- 1E-6     AMA tests,     Carpio  3.96992104
         !this%Ju_exact = 3.970298   ! +- 1E-6     AMA tests,     Carpio  3.96992104
         !this%Ju_exact = 3.9703045   ! +- 1E-6     AMA tests,    Carpio  3.96992104
         this%Ju_exact  = 3.9703044  ! +- 1E-7  AMA MULTI-tests,  Carpio  3.96992104  SISC
         this%Ju_exact  = 3.97030506  ! +- 2E-8  AMA MULTI-tests, Carpio  3.96992104
         this%Ju_exact  = 3.97030506625  ! +- 1E-10  AMA MULTI-tests, Carpio  3.96992104

      else
         this%Ju_exact = 0.0
         print*, '# DWR: The exact solution is not a priori known for icase' , &
            state%model%icase,'. Ju_exact cannot be computed!'
      endif

   end subroutine computeJuExact_BoundaryValue

   subroutine computeJu_BoundaryValue( this, grid )
      class( BoundaryValue_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof, iFace
      real :: local_Ju
      real :: time
      real, dimension(:), allocatable :: Jloc

      time = state%time%finTime
      this%Ju = 0.0

      allocate(Jloc (1:this%isupp), source = 0.0 )

      do i = 1, this%isupp
         elem => grid%elem( this%supp(i) )
         iFace = this%suppFace(i) ! local index of the edge we need to integrate over
         Qdof = elem%face(fGdof, iFace )


         allocate(wi(1:Qdof,1:ndim), source = 0.0)
         allocate( f(1:Qdof, 1:ndim) , source = 0.0 )

         ! eval w in !FACE! integ nodes, 0 - endpoint
         call EvalwSTEdge(elem, iFace, elem%TQnum, 0, wi(1:Qdof,1:ndim), .false. )

         do j = 1,Qdof
            ! xi on the face !
            f(j,1:ndim) = this%evalWeightFunction(elem%xi( iFace ,j,1:nbDim), time)
         end do !j


         f(1:Qdof,1:ndim) = f(1:Qdof, 1:ndim) * wi(1:Qdof, 1:ndim)

         ! integrate over the edge iFace
         call IntegrateFunctionEdge( elem, iFace, f(1:Qdof, 1), local_Ju )
         this%Ju = this%Ju + local_Ju

         Jloc(i) = local_Ju

         !write(*,'(a15, i5,2es22.14, 30es12.4)') '####  J(u_h)',i, this%Ju, local_Ju
         deallocate(f, wi)

      end do

      print*
      !write(*,'(a15, i5,2es22.14, 30es12.4)') '####  J(u_h)',-99, &
      !     this%Ju, sum(Jloc( 1:this%isupp) ), this%Ju - sum(Jloc(1:this%isupp) )

      deallocate(Jloc)

   end subroutine computeJu_BoundaryValue


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! BOUNDARY FLUX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   subroutine init_BoundaryFlux( this, xy_coord, grid, weight_i)
      class( BoundaryFlux_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      this%linear = .false.
      this%boundary = .true.
      this%time_dependent = .false.
      this%name = 'BoundaryFlux'

      print*, ''
      print*, 'Boundary Flux target functional have not been tested yet!'
      print*, ''

      this%iWeight = weight_i
      !1 =  constant 1
!      this%iWeight = 2 ! weight fun psi = y
!      this%iWeight = 3 ! weight fun psi = y*sin
!      print*, 'We use weight function 1/y in target functional over dOm!'
      this%integralAver = .false. !do not divide by |\om_J|

   end subroutine init_BoundaryFlux

   !> marks the elements of the support of J, counts them and computes also volume of this support (edges)
   !> support = element which have iSubmesh == 1
   !> the same as findSupp_BoundaryFlux - these two should have a parent boundary object
   subroutine findSupp_BoundaryFlux( this, grid )
      class( BoundaryFlux_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
      class( element ), pointer :: elem
      integer :: nelem, i, j
      integer, dimension(:), allocatable :: temp_supp, temp_suppFace
      real :: area

      nelem = grid%nelem
      j = 0
      area = 0.0

      allocate( temp_supp(1:nelem), source = -1 )
      allocate( temp_suppFace(1:nelem), source = -1 )

      do i = 1, nelem
         elem => grid%elem(i)
         if (elem%iSubMesh == 1) then
            j = j+1
            ! add the index of the element with support
            temp_supp(j) = i
            ! add the local index of the edge in the support
            temp_suppFace(j) = elem%iSubmeshFace
            ! add the length of the edge
            area = area + elem%dn( elem%iSubmeshFace )
         end if
      end do !i

      this%isupp = j

      if (allocated(this%supp) ) &
         deallocate(this%supp)

      if (allocated( this%suppFace) ) &
         deallocate(this%suppFace)

      allocate( this%supp(1:this%isupp), source = 0 )
      allocate( this%suppFace(1:this%isupp), source = 0)
      this%supp(1:this%isupp) = temp_supp(1:this%isupp)
      this%suppFace(1:this%isupp) = temp_suppFace(1:this%isupp)

      deallocate( temp_supp , temp_suppFace )

      this%vol_supp = area !area ! NOT TESTED, NOT USED

   end subroutine findSupp_BoundaryFlux

   !> compute J(u) = \int \mA(u,grad(u))*n * j_D dS
   !> J(uExact) = JJ(uExact) ie the modified functional equals J
   subroutine computeJuExact_BoundaryFlux( this, grid )
      class( BoundaryFlux_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: f, Rflux
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time
      integer :: iFace

      ! for some problems the precomputed solution is more precise even with known solution
      ! works only for SCALAR problems
      if ( state%model%known_sol ) then

         time = state%time%finTime
         this%Ju_exact = 0.0

         do i = 1, this%isupp
            elem => grid%elem( this%supp(i) )
            iFace = this%suppFace(i) ! local index of the edge we need to integrate over
            Qdof = elem%face(fGdof, iFace )

            ! compute diffusive flux
            allocate( Rflux(1:Qdof,1:ndim), source = 0.0 )
            allocate( f(1:Qdof, 1:ndim), source = 0.0 )

            ! compute A(u,grad(u))*n
            call EvalSolutionEdge_NeumannBoundaryFlux(elem, iFace, Set_R_s_scalar, Rflux(1:Qdof,1:ndim))

            ! the added part in modified J is 0 for exact solution :-)

            ! eval w, f in integ nodes, 0 - endpoint
            do j = 1, Qdof
               ! xi on the face !
               f(j,1:ndim) = this%evalWeightFunction( elem%xi( iFace ,j,1:nbDim), time)
            enddo

            f(1:Qdof,1:ndim) = f(1:Qdof, 1:ndim) * Rflux(1:Qdof, 1:ndim)

            ! integrate over the edge iFace - normalizes by the length of the edge
            call IntegrateFunctionEdge( elem, iFace, f(1:Qdof, 1), local_Ju )
            this%Ju_exact = this%Ju_exact + local_Ju

            deallocate(f, Rflux)

         end do !i
      else
         this%Ju_exact = 0.0
         print*, '# DWR: The exact solution is not a priori known for icase' , &
            state%model%icase,'. Ju_exact cannot be computed!'
      endif

   end subroutine computeJuExact_BoundaryFlux

   !> compute the functional for the discrete solution
   !> we work with the modified functional JJ:
   !> J(v) = \int j_D * A(v,grad(v)*n
   !> JJ(v) = J(v) - \sum_{\partial K \cap Gamma_D) ( \int_{-}( \sigma(v-u_d)*j_D ) + \int_{+}((\sigma + b*n)(v-u_D)j_D ) )
   subroutine computeJu_BoundaryFlux( this, grid )
      class( BoundaryFlux_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi, Re_1, Rflux
      real, allocatable, dimension(:,:) :: jD, bTimesN, jump
      real, allocatable, dimension(:,:,:) :: R_s_inner, Dwi
      integer :: i, j, Qdof, dof, iFace, l
      !logical :: partialPlus ! \partial K^+
      real :: local_Ju
      real :: time, sigma

      time = state%time%finTime
      this%Ju = 0.0

      do i = 1, this%isupp
         ! this edge should be part of Dirichlet BC
         elem => grid%elem( this%supp(i) )
         iFace = this%suppFace(i) ! local index of the edge we need to intagrate over
         Qdof = elem%face(fGdof, iFace )

         ! compute weighting function
         allocate(jD(1:Qdof,1:ndim) , source = 0.0)
         do j = 1, Qdof
               ! xi on the face !
               jD(j,1:ndim) = this%evalWeightFunction( elem%xi( iFace ,j,1:nbDim), time)
         enddo

         !!!! A(v,grad(v)*n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! amount  of diffusivity \epsilon = 1/Re
         allocate( Re_1(1:iRe,1:Qdof), source = state%model%Re1 )
         ! compute flux, w and Dw from inside
         allocate( wi(1:Qdof, 1:ndim ), Dwi(1:Qdof, 1:ndim, 1:nbDim) )
         allocate( R_s_inner(1:Qdof, 1:nbDim, 1:ndim) )
         allocate( Rflux(1:Qdof,1:ndim), source = 0.0 )

         call Eval_w_Edge(elem, iFace, wi(1:Qdof,1:ndim), .false.)
         call Eval_Dw_Edge(elem, iFace, Dwi(1:Qdof, 1:ndim,1:nbDim), .false.)
         call Set_R_s_scalar( ndim, nbDim, iRe, Qdof, wi(1:Qdof,1:ndim), &
                              Dwi(1:Qdof, 1:ndim, 1:nbDim), Re_1(1:iRe,1:Qdof), &
                              R_s_inner(1:Qdof, 1:nbDim, 1:ndim) , &
                              elem%xi(iFace ,1:Qdof, 1:nbDim) )

         ! compute the flux
         do l=1,ndim
           Rflux(1:Qdof, l) = matmul( R_s_inner(1:Qdof, 1:nbDim, l),elem%n(iFace,1:nbDim) ) &
                            / elem%dn(iFace)
         end do !l
         deallocate( wi, Dwi, R_s_inner, Re_1 )

         ! compute the added part, differs if \gamma \in \partial K+ or \partial K-
         ! compute b*n and partial PLUS/MINUS
         allocate( bTimesN(1:Qdof,1:ndim), source = 0.0 )
         ! b*n is in there only for \partial K^PLUS ! - last argument -> max(bn,0)
         call computeBtimesN( elem, iFace, Qdof, bTimesN(1:Qdof,1:ndim), .true.)

         ! add sigma
         sigma = elem%d_gamma * state%model%Re1 / elem%dn(iFace)
         ! b*n is only in \partial K^+
         ! partialPlus not needed - computeBtimesN(...,.true.) gives 0 for b*n<0
         !if (partialPlus) then
         bTimesN(1:Qdof,1:ndim) = bTimesN(1:Qdof,1:ndim) + sigma
         !else
         !   bTimesN(1:Qdof,1:ndim) = sigma
         !endif

         ! compute jump (w_h - u_D)
         allocate( jump(1:Qdof,1:ndim), source = 0.0 )
         call ElementEdgeJump(elem,  iFace, jump(1:Qdof, 1:ndim) )

         ! put together
         jD(1:Qdof,1:ndim) = jD(1:Qdof,1:ndim)* &
            ( Rflux(1:Qdof,1:ndim) - bTimesN(1:Qdof,1:ndim)*jump(1:Qdof,1:ndim) )

         ! integrate over the edge iFace
         call IntegrateFunctionEdge( elem, iFace, jD(1:Qdof, 1), local_Ju )
         this%Ju = this%Ju + local_Ju

         deallocate(jD, jump, bTimesN, Rflux)
      end do

   end subroutine computeJu_BoundaryFlux






!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! NOT WORKING !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! POINT VALUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   subroutine init_point_value( this, xy_coord, grid, weight_i)
!      class( Point_value_t ), intent(inout) :: this
!      real, dimension(1:nbDim), intent(in) :: xy_coord
!      class( mesh ), intent(in) :: grid
!      character(len=50) :: gridfile     ! file with dual mesh
!      integer, intent(in) :: weight_i ! index for weighting function j_om
!
!      this%linear = .true.
!      this%boundary = .false.
!      this%time_dependent = .false.
!      this%name = 'PointValue'
!      this%eps = 0.1!0.01 * grid%diam
!
!!      this%eps = 0.02 !max( 0.1 , 0.5 * h)
!      print*, 'Epsilon =', this%eps
!!      print*, 'h,eps:', h, this%eps
!      allocate( this%xy_coord(1,1:nbDim) , source = 0.0 )
!      this%xy_coord(1,1:2) = xy_coord(1:2)
!
!
!      gridfile = '../Grids/gridPointVal.grid'
!
!      call this%initGrid(gridfile)
!!      stop 'AFTER initGrid'
!
!   end subroutine init_point_value
!
!
!   subroutine initGrid_point_value( this, gridfile )
!      class( Point_value_t ), intent(inout) :: this
!      character(len=50), intent(in) :: gridfile     ! file with dual mesh
!      integer :: nelem, npoin,nbelm
!      real :: eps
!      real, dimension(1:5,1:nbDim) :: x
!      real, dimension(1:2,1:nbDim) :: xper
!      integer, dimension(1:2,1:nbDim) :: iper
!      integer, dimension(:), allocatable :: lenloc
!      integer, dimension(:,:), allocatable :: iloc
!      integer :: i
!
!      allocate( mesh :: this%grid )
!
!      print*, 'TODO : deinit target Grid in the end of computation!'
!      ! read the template square
!      call this%grid%read( gridfile )
!
!!      print*, 'xy_coord and eps:', this%xy_coord(1,1:nbDim), this%eps
!      call this%grid%shiftMesh( 1, this%xy_coord(1,1:nbDim), this%eps )
!      call this%grid%seekNeighbours()
!
!      ! CURVED ???
!      this%grid%b_edge(1:this%grid%nbelm)%icurv = 0
!      this%grid%elem(1:this%grid%nelem)%ibcur = -1
!      this%grid%elem(1:this%grid%nelem)%jcur = 0
!      this%grid%elem(1:this%grid%nelem)%deg_cur = 1
!
!      call this%grid%computeGeometry()
!
!      call this%grid%plot('../DWR/meshTargetJ.gnu')
!
!      ! controls !
!
!!      deallocate( lenloc, iloc )
!
!   end subroutine initGrid_point_value
!
!   !> marks the elements of the support of J, counts them and computes also volume of this support
!   !> we also compute the value of this%eps from the parameters of the mesh
!   subroutine findSupp_point_value( this, grid )
!      class( Point_value_t ), intent(inout) :: this
!      class( mesh ), intent(in) :: grid
!      class( element ), pointer :: elem
!      integer :: nelem, i, j, dof, Tdof, k
!      integer, dimension(:,:), allocatable :: temp_supp
!      real :: dist
!      logical :: inSupp
!      integer :: nsq, max_tri
!      real, dimension(1:2,1:2) :: rmx
!      integer, dimension(:, :, :), allocatable:: itri
!      integer, dimension(:, :), allocatable :: ntri
!      real, dimension(1:2) :: x_bary ! not used - only for the subroutine call
!
!!      integer, dimension(1:3) :: A = (/ 1, 5, 2 /)
!!      integer, dimension(4) :: B = (/ 1, -50, 2, 0 /)
!!      integer, dimension(:), allocatable :: C
!!
!!
!!      call heapsort(A)
!!      write(*,*)'Sorted array :',A
!!      C = MergeArrays(A,B)
!!      print*, C
!
!
!
!      type(intersect)			:: inter
!      integer, dimension(:), allocatable	:: triangles! numbers of triangles in gridN which have
!      integer :: NumTri
!
!   end subroutine findSupp_point_value
!
!
!      !> marks the elements of the support of J, counts them and computes also volume of this support
!   subroutine findSupp_point_value_old( this, grid )
!      class( Point_value_t ), intent(inout) :: this
!      class( mesh ), intent(in) :: grid
!      class( element ), pointer :: elem
!      integer :: nsq, max_tri
!      real, dimension(1:2,1:2) :: rmx
!      integer, dimension(:, :, :), allocatable:: itri
!      integer, dimension(:, :), allocatable :: ntri
!      real, dimension(1:2) :: x_bary ! not used - only for the subroutine call
!
!      print*, ' HERE findSupp_point_value_old!'
!
!      ! grid, nsq = size of the square, frames of the domain, nelem, # of elems in the i,j square, indices of the elements in squares
!
!      call FindTriangleIndexs(grid, nsq, rmx, max_tri, ntri, itri )
!
!      call FindTriangleCoords( grid, elem, this%xy_coord(1,1:nbDim), x_bary(1:2), &
!            nsq, max_tri, ntri, itri, rmx(1:2,1:2) )
!
!      this%vol_supp = elem%area
!      this%isupp = 1
!
!      if (allocated(this%supp) ) &
!         deallocate(this%supp)
!
!      allocate( this%supp(1:this%isupp), source = 0 )
!      this%supp(1:this%isupp) = elem%i
!
!      ! => the point xy_coord is in the element elem !
!      this%eps = elem%diam
!
!   end subroutine findSupp_point_value_old
!
!!   !> computes the dual rhs in a point
!!   !> used when \f$  1/|B_{\epsilon}| \int_{B_{\epsilon}} u dx  \f$
!!   !> together with findSupp_point_value_old
!!   function evalWeightFunction_pointVal( this, x, t) result(f)
!!    class( Point_value_t ),intent(in) :: this
!!    real, dimension(1:nbDim), intent(in) :: x
!!    real, intent(in) :: t
!!    real, dimension(1:ndim) :: f
!!
!!    f(1:ndim) = 0.0
!!
!!    if (t == state%time%FinTime) then
!!      f(1) = 1.0 / (this%vol_supp)
!!    else
!!      stop 'evalWeightFunction_pointVal: CONTROL nonzero only in final time'
!!    endif
!!
!!   end function evalWeightFunction_pointVal
!
!   subroutine computeJu_exact_pointVal_old( this )
!      class( Point_value_t ), intent(inout) :: this
!!      class(mesh), intent(in) :: grid
!      real, dimension(1:ndim) :: w
!
!      print*, 'computeJu_exact_pointVal trule the value u(x,y) not integral!'
!
!      if ( state%time_dependent ) then
!         stop 'computeJu_exact_pointVal not implemented for time-dependent problems!'
!      endif
!
!      call Set_Model_Data( this%xy_coord(1,1:nbDim), state%time%FinTime , w, 1 )
!
!      this%Ju_exact = w( 1 )
!
!   end subroutine computeJu_exact_pointVal_old
!
!   !> computes the value J(u)
!   subroutine computeJu_pointVal_old( this, grid )
!      class( Point_value_t ), intent(inout) :: this
!      class(mesh), intent(in) :: grid
!      class( element ), pointer :: elem
!      real, allocatable, dimension(:) :: wi
!      real, allocatable, dimension(:,:) :: f
!      integer :: i, j, Qdof, dof
!      real :: local_Ju
!      real :: time
!
!      print*, 'computeJu_pointVal: Ju ~ on the choice of the pointVal approximation!!!'
!      print*, 'F@R: computeJu_pointVal - set time ( temporarily finTime )'
!      time = state%time%finTime
!
!      this%Ju = 0.0
!
!      do i = 1, this%isupp
!         elem => grid%elem( this%supp(i) )
!         Qdof = elem%Qdof
!
!         allocate(wi(1:Qdof), source = 0.0)
!         ! eval w in integ nodes, 0 - endpoint
!         wi(1:Qdof) = Eval_whST_iCoord_Elem( elem, elem%TQnum, 0, 1)
!         ! we use only wi(:,1)
!         !          call IntegrateVectorFunction2( elem, wi(1:Qdof,1:ndim), local_Ju)
!
!         allocate( f(1:Qdof, 1:ndim) , source = 0.0 )
!         do j = 1,Qdof
!            f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
!         end do !j
!
!         local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof,1), wi(1:Qdof) )
!
!         this%Ju = this%Ju + local_Ju
!
!         deallocate(f, wi)
!
!      end do
!
!   end subroutine computeJu_pointVal_old
!
!   ! computes the integral of the exact solution (one of its dimension from 1:ndim) over the elements of the support of J
!   subroutine computeJu_exact_pointVal( this, grid )
!      class( Point_value_t ), intent(inout) :: this
!      class(mesh), intent(in) :: grid
!      class( element ), pointer :: elem
!      real :: w
!
!      stop 'not done'
!      print*, 'computeJu_exact_pointVal trule the value u(x,y) not integral!'
!!
!!      if ( state%time_dependent ) then
!!         stop 'computeJu_exact_pointVal not implemented for time-dependent problems!'
!!      endif
!!
!!      call Set_Model_Data( this%xy_coord, state%time%FinTime , w, 1 )
!!
!!      do i = 1, this%isupp
!!         elem => grid%elem( this%supp(i) )
!!
!!         do j = 1, elem%Qdof
!!            call Set_Model_Data( this%xy_coord, state%time%FinTime , w, 1 )
!!             elem%xi(0, j, 1:nbDim)
!!
!!         end do !j
!!
!!      enddo
!!      this%Ju_exact = w( 1 )
!
!   end subroutine computeJu_exact_pointVal
!
!   !> computes the value J(u)
!   subroutine computeJu_pointVal( this, grid )
!      class( Point_value_t ), intent(inout) :: this
!      class(mesh), intent(in) :: grid
!      class( element ), pointer :: elem
!      real, allocatable, dimension(:) :: wi
!      real, allocatable, dimension(:,:) :: f
!      integer :: i, j, Qdof, dof
!      real :: local_Ju
!      real :: time
!
!      print*, 'computeJu_pointVal: Ju ~ on the choice of the pointVal approximation!!!'
!
!      print*, 'F@R: computeJu_pointVal - set time ( temporarily finTime )'
!      time = state%time%finTime
!
!      this%Ju = 0.0
!
!      do i = 1, this%isupp
!         elem => grid%elem( this%supp(i) )
!         Qdof = elem%Qdof
!
!         allocate(wi(1:Qdof), source = 0.0)
!         ! eval w in integ nodes, 0 - endpoint
!         wi(1:Qdof) = Eval_whST_iCoord_Elem( elem, elem%TQnum, 0, 1)
!         ! we use only wi(:,1)
!         !          call IntegrateVectorFunction2( elem, wi(1:Qdof,1:ndim), local_Ju)
!
!         allocate( f(1:Qdof, 1:ndim) , source = 0.0 )
!         do j = 1,Qdof
!            f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
!         end do !j
!
!         local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof, 1), wi(1:Qdof) )
!
!         this%Ju = this%Ju + local_Ju
!
!         deallocate(f, wi)
!
!      end do
!
!   end subroutine computeJu_pointVal












!!! U OVER SUBDOMAIN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


   subroutine init_u_over_subdomain( this, xy_coord, grid, weight_i)
      class( U_over_subdomain_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      this%linear = .true.
      this%boundary = .false.
      this%time_dependent = .false.
      this%name = 'UoverSubdomain'

      this%iWeight = weight_i
      ! 1 = constant 1
      this%integralAver = .true. !divide by |\om_J|
      !!this%integralAver = .false. !divide by |\om_J|

   end subroutine init_u_over_subdomain

   !> marks the elements of the support of J, counts them and computes also volume of this support
   !> support = element which have iSubmesh == 1
   subroutine findSupp_u_over_subdomain( this, grid )
      class( U_over_subdomain_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
      class( element ), pointer :: elem
      integer :: nelem, i, j
      integer, dimension(:), allocatable :: temp_supp
      real :: area


      if ( grid%elem(1)%iSubMesh == 0 ) then
         stop 'The subdomain partition was not done yet in findSupp_u_over_subdomain!'
      endif

      nelem = grid%nelem
      j = 0
      area = 0.0

      allocate( temp_supp(1:nelem), source = -1 )

      do i = 1, nelem
         elem => grid%elem(i)
         if (elem%iSubMesh == 1) then
            j = j+1
            temp_supp(j) = i
            area = area + elem%area
        end if
      end do !i

      this%isupp = j

      if (allocated(this%supp) ) &
         deallocate(this%supp)

      allocate( this%supp(1:this%isupp), source = 0 )
      this%supp(1:this%isupp) = temp_supp(1:this%isupp)

      deallocate( temp_supp )

      this%eps = elem%r_ins
      this%vol_supp = area

   end subroutine findSupp_u_over_subdomain

   subroutine computeJu_exact_u_over_subdomain( this, grid )
      class( U_over_subdomain_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time

      ! for some problems the precomputed solution is more precise even with known solution
      ! convection problem isca 29
      if ( state%model%icase == 29 ) then
         !this%Ju_exact = 0.114
         this%Ju_exact = 0.1175   ! obtained using DWR_AMA estimates with HGh adapt, P_3, DOF = 183580
         this%Ju_exact = 0.11745997905 ! obtained using DWR_AMA estimates with AMAhp adapt, P_3, DOF = 74938
      ! peaks
      else if ( state%model%icase == 64 ) then
         ! connected with twoSquaresIn250
         this%Ju_exact = 4.6769935E-3
         ! connected with twoSquaresIn164
!         this%Ju_exact = 1.379437E-2
      ! unknown solution
      ! LL shaped peak
      else if ( state%model%icase == 65 ) then
         this%Ju_exact = 2.138 * 0.001

      ! CROSS domain problem, Ainsworth, Rankin 2012
      ! exact = 0.01630471454734821 *25 (25: area of the subdomain = 0.04)
      else if ( state%model%icase == 70 ) then ! domain cross
         this%Ju_exact = 0.40761786368370525 ! sent from prof Rankin
         !this%Ju_exact = 0.4076211063069 ! Ajay modification
         
      else if ( state%model%icase == 71 ) then ! Forchheimer flow
         !imaterial = 2
         this%Ju_exact = 0.49114690441416270 ! P3 5-adapt cycles
         ! 0.496669 P3
         ! 0.496625
         ! 0.49139655025293516
!      else if ( state%model%icase == 74 .and. state%model%Re == 100. ) then
!         ! exact computed with P=8 on mesh with 5000elems, alpha = 200
!         this%Ju_exact = 0.157058818600
      else if ( state%model%icase == 80 ) then ! Carpio
        if ( state%model%Re /= 1000 ) print*, 'wrong size of diffusion coef Re!'
        !this%Ju_exact = 0.20305828 ! from the Carpio's article 2013
        !this%Ju_exact = 0.2031415  !+-1E-6 from HG strong refinement, DWR_AMA_CV_HG
        this%Ju_exact = 0.20314158  !+-1E-8 from AMA multi-refinement, SISC
        this%Ju_exact = 0.2031415665  !+-2E-9 from AMA multi-refinement
        this%Ju_exact = 0.2031415667117  !+-1E-12 from AMA multi-refinement


        !     ! we know the exact solution
      else if ( state%model%known_sol ) then

         !      print*, 'F@R: computeJu_exact - set time ( temporarily finTime )'
         time = state%time%finTime

         this%Ju_exact = 0.0

         do i = 1, this%isupp
            elem => grid%elem( this%supp(i) )
            Qdof = elem%Qdof

            allocate(wi(1:Qdof,1:ndim), source = 0.0)
            allocate( f(1:Qdof, 1:ndim) , source = 0.0 )
            ! eval w in integ nodes, 0 - endpoint
            do j = 1, Qdof
               ! integ nodes, time, w, ityp = 1 (exact solution)
               ! we use only wi(:,1)
               call Set_Model_Data( elem%xi(0,j,1:nbDim) , time , wi(j,1:ndim), 1 )

               f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
            enddo

            local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof, 1), wi(1:Qdof, 1) )
            this%Ju_exact = this%Ju_exact + local_Ju

            deallocate(f, wi)

         end do !i
      else
         this%Ju_exact = 0.0
         print*, '# DWR: The exact solution is not a priori known for icase' , &
            state%model%icase,'. Ju_exact cannot be computed!'
      endif

   end subroutine computeJu_exact_u_over_subdomain

   subroutine computeJu_u_over_subdomain( this, grid )
      class( U_over_subdomain_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:) :: local
      integer, allocatable, dimension(:) :: i_local
      real, allocatable, dimension(:) :: wi
      real, allocatable, dimension(:,:) :: f
      real, dimension(:,:), allocatable :: Fxi
      type(volume_rule), pointer :: V_rule

      integer :: i, j, Qdof, dof, ifile, ifile1, ifile2
      real :: local_Ju, val
      real :: time

      ifile  = 600 + state%space%adapt%adapt_level
      ifile1 = 700 + state%space%adapt%adapt_level
      ifile2 = 800 + state%space%adapt%adapt_level

!      write(debug,*) 'F@R: computeJu_pointVal - set time ( temporarily finTime )'
      time = state%time%finTime

      this%Ju = 0.0

      allocate(local(1: this%isupp ), i_local(1: this%isupp ) )

      do i = 1, this%isupp
         elem => grid%elem( this%supp(i) )
         Qdof = elem%Qdof

         allocate(wi(1:Qdof), source = 0.0)
         allocate( f(1:Qdof, 1:ndim) , source = 0.0 )
         ! eval w in integ nodes, 0 - endpoint, 1 - first component (density)
         wi(1:Qdof) = Eval_whST_iCoord_Elem( elem, elem%TQnum, 0, 1)
         ! we use only wi(:,1)
         !   call IntegrateVectorFunction2( elem, wi(1:Qdof,1:ndim), local_Ju)

         ! only for drawing test
         !V_rule => state%space%V_rule(elem%Qnum)
         !allocate( Fxi(1:Qdof, 1:2) )
         !call ComputeF(elem, Qdof, V_rule%lambda(1:Qdof, 1:2), Fxi(1:Qdof,1:nbDim) )

         do j = 1,Qdof
            f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
            !write(ifile1, *) Fxi(j, 1:2), f(j,1:ndim), wi(j)
         end do !j
!         print*, "eval weight fun =" , f(1,:)

         !write(ifile2, *) elem%xc(:)

         local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof, 1), wi(1:Qdof) )
         this%Ju = this%Ju + local_Ju

         !local(i) = local_Ju
         !i_local(i) = i

         !call PlotElemFunction3D(ifile, elem, elem%dof, elem%wST( 1:, 1:elem%dof, 1))
         !deallocate(Fxi)

         deallocate(f, wi)

      end do

      ! if(state%space%adapt%adapt_level == 0) then
      !    open(38, file = "vals_Ju", action ="write", status="replace")
      ! else
      !    open(38, file = "vals_Ju", action ="write", status="UNKNOWN", position = "append")
      ! endif

      ! write(38, * ) ' original this%Ju = ', this%Ju
      ! write(38, * ) ' sum      this%Ju = ', sum(local(:) )

      ! call order_estims( this%isupp, i_local(1: this%isupp), local(1: this%isupp) )

      ! write(38, * ) ' ordering this%Ju = ', sum(local(:) )

      ! val = 0.
      ! do i=1, this%isupp
      !    val = val + local(i)
      ! enddo
      ! write(38, * ) ' increase this%Ju = ', val

      ! val = 0.
      ! do i= this%isupp, 1, -1
      !    val = val + local(i)
      ! enddo
      ! write(38, * ) ' decrease this%Ju = ', val

      ! write(38,*)' __________________________________________', state%space%adapt%adapt_level
      ! close(38)
      ! deallocate( local, i_local)

   end subroutine computeJu_u_over_subdomain








!!!! DUDX_T !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


   subroutine init_dudx( this, xy_coord, grid, weight_i)
      class( dudx_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      this%linear = .true.
      this%boundary = .false.
      this%time_dependent = .false.
      this%name = 'dudx'
      this%iWeight = weight_i
      this%integralAver = .false. !.true. !divide by |\om_J|

      ! du/dx -- 1
      this%dx = 1

   end subroutine init_dudx

   !> marks the elements of the support of J, counts them and computes also volume of this support
   !> support = element which have iSubmesh == 1
   !> copied fro u_over_subdomain_t
   subroutine findSupp_dudx( this, grid )
      class( dudx_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
      class( element ), pointer :: elem
      integer :: nelem, i, j
      integer, dimension(:), allocatable :: temp_supp
      real :: area

      if ( grid%elem(1)%iSubMesh == 0 ) then
         stop 'The subdomain partition was not done yet in findSupp_u_over_subdomain!'
      endif

      nelem = grid%nelem
      j = 0
      area = 0.0

      allocate( temp_supp(1:nelem), source = -1 )

      do i = 1, nelem
         elem => grid%elem(i)
         if (elem%iSubMesh == 1) then
            j = j+1
            temp_supp(j) = i
            area = area + elem%area
        end if
      end do !i

      this%isupp = j

      if (allocated(this%supp) ) &
         deallocate(this%supp)

      allocate( this%supp(1:this%isupp), source = 0 )
      this%supp(1:this%isupp) = temp_supp(1:this%isupp)

      deallocate( temp_supp )

      this%eps = elem%r_ins
      this%vol_supp = area

   end subroutine findSupp_dudx

   subroutine computeJu_exact_dudx( this, grid )
      class( dudx_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time

      ! do we know the exact solution
      if ( state%model%known_sol ) then

         time = state%time%finTime

         this%Ju_exact = 0.0

         do i = 1, this%isupp
            elem => grid%elem( this%supp(i) )
            Qdof = elem%Qdof

            allocate(wi(1:Qdof,1:ndim), source = 0.0)
            allocate( f(1:Qdof, 1:ndim) , source = 0.0 )
            ! eval w in integ nodes, 0 - endpoint
            do j = 1, Qdof
               ! integ nodes, time, w, ityp = 1 (exact solution)
               ! we use only wi(:,1) , 8 - derivative du/dx
               call Set_Model_Data( elem%xi(0,j,1:nbDim) , time , wi(j,1:ndim), 8 )

               f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
            enddo

            local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof, 1), wi(1:Qdof, 1) )
            this%Ju_exact = this%Ju_exact + local_Ju

            deallocate(f, wi)

         end do !i
      ! convection problem isca 68
      else if ( state%model%icase == 68 ) then
         ! connected with submesh square8TwoTriangles
!         this%Ju_exact = 1.585090814E-3
         this%Ju_exact = 1.58509081387E-3 ! p7, AMAhp
!                        1.5850908130823156E-3 ! p7, HGh
!                        1.5850908137 ! p7 256 elements
!          0.00158509081389 !
!           .12345678901234 ~ E-14
      else
         this%Ju_exact = 0.0
         print*, '# DWR: The exact solution is not a priori known for icase' , state%model%icase,'. Ju_exact cannot be computed!'
      endif

   end subroutine computeJu_exact_dudx

   subroutine computeJu_dudx( this, grid )
      class( dudx_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time

      if (ndim > 1) &
         stop 'Problem in computeJu_dudx wi( 1 ,1:Qdof) with ndim>1!'

!      if ( state%dual ) &
!         stop 'the primal solution is not saved in elem%wST'

      write(debug,*) 'F@R: computeJu_pointVal - set time ( temporarily finTime )'
      time = state%time%finTime

      this%Ju = 0.0

      do i = 1, this%isupp
         elem => grid%elem( this%supp(i) )
         Qdof = elem%Qdof

         allocate(wi(1:ndim, 1:Qdof), source = 0.0)
         allocate( f(1:Qdof, 1:ndim) , source = 0.0 )

         ! eval dw/dx in integ nodes, 0 - endpoint
         wi( 1:ndim ,1:Qdof) = evalSTfunInIntTime_spaceDer( elem, ndim, elem%dof, elem%Tdof, &
            elem%wST(1:ndim,1:elem%dof,1:elem%Tdof), 0, elem%tQnum, this%dx )

         ! we use only wi(:,1)
         !          call IntegrateVectorFunction2( elem, wi(1:Qdof,1:ndim), local_Ju)
         do j = 1,Qdof
            f(j,1:ndim) = this%evalWeightFunction(elem%xi(0,j,1:nbDim), time)
         end do !j

         local_Ju = EvalL2ScalarProduct( elem, f(1:Qdof, 1), &
            wi(1, 1:Qdof) )

         this%Ju = this%Ju + local_Ju

         deallocate(f, wi)

      end do

   end subroutine computeJu_dudx

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!! DRAG AND LIFT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   subroutine init_DragAndLift( this, xy_coord, grid, weight_i)
      class( DragAndLift_t ), intent(inout) :: this
      real, dimension(1:nbDim), intent(in) :: xy_coord
      class( mesh ), intent(in) :: grid
      integer, intent(in) :: weight_i ! index for weighting function j_om

      this%linear = .false.
      this%boundary = .true.
      this%time_dependent = .false.
      this%name = 'dragAndLift'

      if (weight_i == 10) then
         this%dlm = "drag"
         print*, '# Drag target functional'
      else if (weight_i == 11 ) then
         this%dlm = "lift"
         print*, '# Lift target functional'
      else if (weight_i == 12 ) then
         this%dlm = "momentum"
         print*, '# Momentum target functional'
      else
         stop "unknown kind of DragAndLift target functional!"
      endif

      ! used for momentum coef only
      allocate( this%xy_coord (1:nbDim), source = xy_coord(1:nbDim) )

      if (ndim /= 4) &
         stop  'The drag and lift is only for ndim=4!'

   end subroutine init_DragAndLift

   !> computes the weighting function in the directional derivative of Jh
   !> fi = \vartheta^T \mH, where the matrix \mH depends on the choice of the discretization of the slip BC
   !> which is set by parameter state%model%discretizationOfSlipBC_Euler == 1, 2, 3
   !> 1) direct use of \mP_W -> \mH = \mP_W(w,n)
   !> 2) Hartmans approach -> \mH = \mP_W(u_Gamma,n) \mU , u_Gamma,\mU come from the boundary value operator (see subroutine UpdateMirror)
   !> 3) use of the Vijayasundaram flux -> \mH = \mP^+ + \mP^-\mM
   !> vartheta depends on whether Drag, Lift, or Momentum is chosen as J
   subroutine evalWeightFunctionVector_Euler( this, elem, ie, Qdof, fi )
        class( DragAndLift_t ), intent(inout) :: this
        class( element ), intent(inout) :: elem
        integer, intent(in) :: ie ! index of the edge
        integer, intent(in) :: Qdof ! quad nodes on the edge
        real, dimension(1:Qdof, 1:ndim), intent(inout) :: fi
        real :: c_infty
        real, dimension(:,:), allocatable :: theta, Q, G, wi, nc, ncu
        real, dimension(:,:,:), allocatable :: H
        integer :: j, k

        if (Qdof /= elem%face(fGdof, ie)) &
          stop "problem with quadrature nodes in evalWeightFunctionEuler"

        ! compute THETA
        allocate( theta(1:Qdof, 1:ndim), source = 0.0 )
        call this%setTheta( Qdof, elem%xi( ie ,1:Qdof,1:nbDim), theta(1:Qdof, 1:ndim) )

        !!!! compute \mH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        allocate( H(1:Qdof,1:ndim,1:ndim) , source = 0.0 )
        call computeBoundaryFluxMatrix_Euler(elem, ie, Qdof, H(1:Qdof,1:ndim,1:ndim) )
        ! compute \theta * \mH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        do j = 1,Qdof
          fi(j,1:ndim) = matmul( theta(j,1:ndim) , H(j,1:ndim,1:ndim) )
        end do

        ! integrate (put to RHS) is done outside this function
        deallocate(theta, H )

   end subroutine evalWeightFunctionVector_Euler


   subroutine setTheta( this, Qdof, xi, theta)
      class( DragAndLift_t ), intent(inout) :: this
      integer, intent(in) :: Qdof
      real, dimension(1:Qdof,1:nbDim), intent(in) :: xi
      real, dimension(1:Qdof,1:ndim), intent(out) :: theta
      integer :: j
      real :: c_infty
      real, dimension(:,:), allocatable :: Q, G
      ! FR_Euler is L_ref == 1 ???
      c_infty = 0.5 * state%rho_infty * state%v_infty * state%v_infty
      !print*, "C infinity in setTheta = ", c_infty

      theta(1:Qdof,1:ndim) = 0.0

      if (this%dlm == "drag") then
        theta(1:Qdof,2) = cos(state%alpha_infty) / c_infty
        theta(1:Qdof,3) = sin(state%alpha_infty) / c_infty

      else if (this%dlm == "lift") then
        theta(1:Qdof,2) = -sin(state%alpha_infty) / c_infty
        theta(1:Qdof,3) = cos(state%alpha_infty) / c_infty
      else if (this%dlm == "momentum") then
        allocate( Q(1:nbDim,1:nbDim), G(1:nbDim,1:nbDim), source = 0.0)

        Q(1,1:2) = (/ cos(state%alpha_infty), sin(state%alpha_infty)/)
        Q(2,1:2) = (/-sin(state%alpha_infty), cos(state%alpha_infty)/)

        G(1,1:2) = (/ 0, -1 /)
        G(2,1:2) = (/ 1, 0 /)
        do j = 1,Qdof
        theta(j,2:3) = matmul( xi(j,1:nbDim) - this%xy_coord(1:nbDim)  , &
                        matmul(G(1:nbDim,1:nbDim), Q(1:nbDim,1:nbDim)) )  / c_infty
        end do
        deallocate(Q,G)
      else
        stop "unknown kind of DragAndLift target functional in setTheta!"
      end if

   end subroutine setTheta

   !> marks the elements of the support of J, counts them and computes also volume of this support (edges)
   !> support = element which have iSubmesh == 1
   !> the support is on the elements which have an edge on the slip boundary part \Gamma_W
   subroutine findSupp_DragAndLift( this, grid )
      class( DragAndLift_t ), intent(inout) :: this
      class( mesh ), intent(in) :: grid
      class( element ), pointer :: elem
      integer :: nelem, i, j, ib, ie
      integer, dimension(:), allocatable :: temp_supp, temp_suppFace
      real :: area

      nelem = grid%nelem
      j = 0
      area = 0.0

      allocate( temp_supp(1:nelem), source = -1 )
      allocate( temp_suppFace(1:nelem), source = -1 )

      do ib =1,grid%nbelm
       if(grid%b_edge(ib)%BC == 0 ) then  ! impermeable walls = profiles
          j = j+1
          elem => grid%elem(grid%b_edge(ib)%itc)
          ie = grid%b_edge(ib)%jtc

          temp_supp(j) = elem%i
          temp_suppFace(j) = ie

          area = area + elem%dn( ie )
       end if
      end do

      this%isupp = j

      if (allocated(this%supp) ) &
         deallocate(this%supp)
      if (allocated( this%suppFace) ) &
         deallocate(this%suppFace)

      allocate( this%supp(1:this%isupp), source = 0 )
      allocate( this%suppFace(1:this%isupp), source = 0)
      this%supp(1:this%isupp) = temp_supp(1:this%isupp)
      this%suppFace(1:this%isupp) = temp_suppFace(1:this%isupp)

      deallocate( temp_supp , temp_suppFace )
      this%vol_supp = area !area ! NOT TESTED, NOT USED
!      print*, 'findSupp_DragAndLift: Size of the profile (not tested):', area


   end subroutine findSupp_DragAndLift

   subroutine computeJu_exact_DragAndLift( this, grid )
      class( DragAndLift_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi
      real, allocatable, dimension(:,:) :: f
      integer :: i, j, Qdof, dof
      real :: local_Ju
      real :: time

      ! do we know the exact solution
      if ( state%model%known_sol ) then
         this%Ju_exact = 0.0
      !        convection problem isca 68
      !else if ( state%model%icase == 68 ) then
      else
         this%Ju_exact = 0.0
         print*, '# DWR: The exact solution is not a priori known for icase' , state%model%icase,'. Ju_exact cannot be computed!'
      endif

   end subroutine computeJu_exact_DragAndLift

   subroutine computeJu_DragAndLift( this, grid )
      class( DragAndLift_t ), intent(inout) :: this
      class(mesh), intent(in) :: grid
      class( element ), pointer :: elem
      real, allocatable, dimension(:,:) :: wi, theta, nc, wiAver, wiMir
      real, allocatable, dimension(:,:) :: f
      real, allocatable, dimension(:,:,:,:) :: Ppm
      integer :: i, j, Qdof, dof, ie, Qnum
      real :: local_Ju
      real :: time

      this%Ju = 0.0

!      print*, "Inconstistent adjoint discretization TEST in computeJu_DragAndLift"
!      print*, "---- put back for correct calculation ---"

      do i = 1, this%isupp
        elem => grid%elem( this%supp(i) )
        ie = this%suppFace(i) ! boundary face

        Qdof = elem%face(fGdof, ie)
        ! solution on the edge
        allocate(wi(1:Qdof,1:ndim))
!        call Eval_w_Edge(elem, ie, wi, .false.)
!        print*, "norm w = " , norm2(wi)

        ! FR_Euler - there were errors with what is saved in elem%w(0,:)
!        call Transfer_wST_to_w_Elem(elem, 0, elem%TQnum)
!        call Eval_w_Edge(elem, ie, wi, .false.)
!        print*, "2nd norm w = " , norm2(wi)

        ! second possibility
        call EvalwSTEdge(elem, ie, elem%TQnum, 0, wi, .false.)
!        print*, "Second way norm w = " , norm2(wi)

        ! setting of outer normals in integration nodes
        allocate(nc(1:Qdof, 1:nbDim) )
        if(elem%ibcur > 0) then
           nc(1:Qdof,1:nbDim) = elem%nc(1:Qdof,1:nbDim)
        else
           nc(1:Qdof,1) = elem%n(ie,1)
           nc(1:Qdof,2) = elem%n(ie,2)
        endif

        allocate( theta(1:Qdof,1:ndim), source = 0.0 )
        call this%setTheta(Qdof, elem%xi(ie,1:Qdof,1:nbDim), &
          theta(1:Qdof,1:ndim) )

        allocate( f(1:Qdof,1:ndim), source = 0.0 )

        if (state%model%discretizationOfSlipBC_Euler == 1) then

!          if (elem%i == 1707) then
!            print*, "This way Jh should equal directly computed drag! TODO"
!            print*, "Element : " , elem%i
!            print*, "Sol w (1) = ", elem%wST
!            print*, "Sol w (1) = ", wi(1:Qdof, 1:ndim)
!            !print*, "Pressure = " , pressure(ndim, wi(1,1:ndim))
!            do j = 1, Qdof
!                print*, "Pressure : j  ", pressure(ndim, wi(j,1:ndim))
!                print*, "F = ", pressure(ndim, wi(j,1:ndim)) * nc(j,1:nbDim)
!
!            end do
!            print*, "Theta:" , theta(1,1:ndim)
!          end if

          do j = 1, Qdof
            !print*, "Pressure : j  ", pressure(ndim, wi(j,1:ndim))
            f(j, 2:3) = pressure(ndim, wi(j,1:ndim)) * nc(j,1:nbDim)
          end do


        else if (state%model%discretizationOfSlipBC_Euler == 2) then
          ! compute u_Gamma(w) = w - (w*n)n
          ! length of n is not important here
          call UpdateMirror(ndim,Qdof,wi(1:Qdof,1:ndim),nc(1:Qdof, 1:nbDim))
          do j = 1, Qdof
            f(j, 2:3) = pressure(ndim, wi(j,1:ndim)) * nc(j,1:nbDim)
          end do
        else if (state%model%discretizationOfSlipBC_Euler == 3) then
          ! average wi = (wi + \Mir(w)) / 2
          allocate( wiAver(1:Qdof,1:ndim), source = wi(1:Qdof, 1:ndim) )
          allocate(Ppm(1:Qdof, 1:nbDim, 1:ndim, 1:ndim), source = 0.0 )

          call UpdateMirror(ndim,Qdof,wiAver(1:Qdof,1:ndim),nc(1:Qdof, 1:nbDim))

          call Set_Ppm_Euler(ndim, nbDim, Qdof, wiAver(1:Qdof,1:ndim), nc(1:Qdof,1:nbDim), &
         elem%xi(ie, 1:Qdof, 1:nbDim),  &
         Ppm(1:Qdof,1:2, 1:ndim, 1:ndim), 1./elem%area, elem )

         allocate( wiMir(1:Qdof,1:ndim), source = wi(1:Qdof, 1:ndim) )
         call Mirror_W(ndim,Qdof,wiMir(1:Qdof,1:ndim),nc(1:Qdof, 1:nbDim))

         do j =1, Qdof
           f(j, 1:ndim) = matmul( Ppm(j,1, 1:ndim, 1:ndim) , wi(j,1:ndim) ) &
           + matmul( Ppm(j,2, 1:ndim, 1:ndim) , wiMir(j,1:ndim) )
         end do

         deallocate(Ppm, wiAver, wiMir)
        else
          stop "Problem in computeJu_DragAndLift"
        end if

        ! for test inconsistency
!        ! FR_Euler temporarily - for test an adjoint inconsistent discretization
!        do j = 1, Qdof
!            !print*, "Pressure : j  ", pressure(ndim, wi(j,1:ndim))
!            f(j, 2:3) = pressure(ndim, wi(j,1:ndim)) * nc(j,1:nbDim)
!        end do

        !!! integrate f * theta
        do j =1,Qdof
          f(j, 1) = dot_product( f(j,1:ndim), theta(j,1:ndim) )
        end do
        local_Ju = 0.0

        ! TODO FR_Euler - integration of the function, normal n is not UNIT
        ! new way
        Qnum = elem%face(fGnum,ie)
        local_Ju = dot_product( f(1:Qdof,1), &
                state%space%G_rule(Qnum)%weights(1:Qdof))

        ! old way
        !call IntegrateFunctionEdge( elem, ie, f(1:Qdof,1), local_Ju)


        this%Ju = this%Ju + local_Ju
!        if (elem%i == 1707) then
!            print*, "Local j(u) = " , local_Ju
!        endif

        deallocate(theta, wi, f, nc)
      end do

      print*, "---- kappa1 (gamma1 ) = " , state%model%kappa1
      print*, "Value of Jh(wh) = ", this%Ju
      print*, "Value od directly computed drag = " , state%cDLM(state%time%iter,1) , &
        " lift = " , state%cDLM(state%time%iter,2)
      print*, "----"



   end subroutine computeJu_DragAndLift


   !> test how the boundary condition: n1*z2+n2*z3 = n*\vartheta is fulfilled
   subroutine testBC_Euler_W( this, elem, bcZerror )
    class( DragAndLift_t ), intent(inout) :: this
    class( element ), intent(inout) :: elem
    real, intent(out) :: bcZerror
    real, dimension(:,:), allocatable :: zi, nc, ncu, wi, theta
    integer :: j, k, Qdof, dof
    real, dimension(:), allocatable :: bcW, bcZ
    real :: bcWerror

    dof = elem%dof

    do j = 1, elem%flen
      k = elem%face(neigh,j)
      if( k <= 0 .and. elem%iBC(j) == 0 ) then

        Qdof = elem%face(fGdof,j)
        allocate(nc(1:Qdof, 1:nbDim) )

        if(elem%ibcur > 0) then  ! UNIT normal
            nc(1:Qdof,1) = elem%nc(1:Qdof,1) / elem%dnc(1:Qdof)
            nc(1:Qdof,2) = elem%nc(1:Qdof,2)  / elem%dnc(1:Qdof)
        else
            nc(1:Qdof,1) = elem%n(j,1) / elem%dn(j)
            nc(1:Qdof,2) = elem%n(j,2) / elem%dn(j)
        endif

        allocate( wi(1:Qdof,1:ndim), source = 0.0 )
        allocate( zi(1:Qdof,1:ndim), source = 0.0 )
        call Eval_w_Edge( elem, j, wi(1:Qdof,1:ndim), .false.)
        call Eval_z_Edge( elem, j, zi(1:Qdof,1:ndim), .false.)

        allocate( theta(1:Qdof,1:ndim), source = 0.0 )
        call this%setTheta(Qdof, elem%xi(j,1:Qdof,1:nbDim), &
          theta(1:Qdof,1:ndim) )

        allocate( bcZ(1:Qdof), source = 0.0 )
        allocate( bcW(1:Qdof), source = 0.0 )
        do k= 1,Qdof
            bcZ(k) = zi(k,2)*nc(k,1) + zi(k,3)*nc(k,2) &
                   - theta(k,2)*nc(k,1) - theta(k,3)*nc(k,2)
            bcW(k) = wi(k,2)*nc(k,1) + wi(k,3)*nc(k,2)
        end do

        bcZerror =  EvalScalarProdEdge(elem, j, bcZ(1:Qdof), bcZ(1:Qdof))
        bcZerror = sqrt(bcZerror)

        bcWerror =  EvalScalarProdEdge(elem, j, bcW(1:Qdof), bcW(1:Qdof))
        bcWerror = sqrt(bcWerror)

        if ( state%getP_mod() == 0) then
!             if (bcWerror > 0.05 .or. bcZerror > 0.05) then
!               print*, "testBC_Euler_W: The BC on Gamma_W is inconsistent on element = ", elem%i
!               print*, "wi = " ,  maxval(abs(bcW(:)))  , " zi - theta =" , maxval(abs(bcZ(:)))
!             end if

!            if ( maxval( abs( bcW(:) ) ) > 0.05 .or. maxval(abs(bcZ(:)) ) > 0.05) then
!               print*, "testBC_Euler_W: The BC on Gamma_W is inconsistent on element = ", elem%i
!               print*, "wi = " ,  maxval(abs(bcW(:)))  , " zi - theta =" , maxval(abs(bcZ(:)))
!               !k = maxloc(abs(bcZ(:)))
!!               do k = 1, Qdof
!!                  print*," zi - theta =" , abs(bcZ(k))
!!                  print*, "Z2,Z3, rhs:", zi(k,2)*nc(k,1) , zi(k,3)*nc(k,2) , &
!!                  theta(k,2)*nc(k,1) + theta(k,3)*nc(k,2)
!!                  print*, "n1,n2,Z2,Z3", nc(k,1) , nc(k,2) , zi(k,2), zi(k,3)
!!
!!               end do
!               elem%eta(dwr_mark,1) =  maxval(abs(bcZ(:)) )
!               ! for testing - COMMENT IN FUTURE
!!               elem%eta(29,1) =  elem%i
!            end if
        endif

      deallocate(zi, wi,  nc, bcZ, bcW, theta)
      end if
    end do

  end subroutine testBC_Euler_W



end module target_functional_mod

