scale_meshfield_filter_operation_base.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Filter operation base
!!
!! @par Description
!!           This module provides a base class to apply filter operations to MeshField data
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshfield_filter_operation_base
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_element_modalfilter, only: modalfilter
  use scale_localmesh_base, only: localmeshbase
  use scale_localmeshfield_base, only: localmeshfieldbase
  use scale_meshfieldcomm_base, only: &
    meshfieldcontainer
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Base type to represent filter operation
  type, public :: meshfieldfilteroperationbase
    integer :: operator_type
    integer :: hhalosize
 
    integer :: nnode_h1d_reconst
 
    !- Reconstruction matrices    
    real(rp), allocatable :: filtermat_h1d(:,:)
 
    !- Reconstruction matrices
    real(rp), allocatable :: minv_ml_tr(:,:)
    real(rp), allocatable :: minv_mc_tr(:,:)
    real(rp), allocatable :: minv_mr_tr(:,:)
    real(rp), allocatable :: intrpmat(:,:)
 
    !- Reconstruction matrices (2)
    real(rp), allocatable :: ml_tr(:,:)
    real(rp), allocatable :: mc_tr(:,:)
    real(rp), allocatable :: mr_tr(:,:)
 
    !- Modal Filter
    type(modalfilter) :: m_filter
 
    !- Interface Correction
    integer :: if_r = 8
    real(rp), allocatable :: if_gl(:)
    real(rp), allocatable :: if_gr(:)
 
  contains
    procedure, private :: prepair_filtermat => meshfieldfilteroperationbase_prepair_filter_matrix
    procedure, private :: prepair_reconstmat => meshfieldfilteroperationbase_prepair_reconstruct_matrix
    procedure, private :: prepair_reconstmat2 => meshfieldfilteroperationbase_prepair_reconstruct2_matrix
    procedure, private :: prepair_reconstmat2_gl => meshfieldfilteroperationbase_prepair_reconstruct2_gl_matrix
    procedure, private :: prepair_interfacecorrection => meshfieldfilteroperationbase_prepair_interface_correction
  end type meshfieldfilteroperationbase
 
  public :: meshfieldfilteroperationbase_init
  public :: meshfieldfilteroperationbase_final
 
  public :: meshfieldfilteroperationbase_apply_filter1d_x
  public :: meshfieldfilteroperationbase_apply_reconst1d_x
  public :: meshfieldfilteroperationbase_apply_reconst1d_x_2
  public :: meshfieldfilteroperationbase_apply_interface_correction1d_x
  public :: meshfieldfilteroperationbase_apply_filter1d_y
  public :: meshfieldfilteroperationbase_apply_reconst1d_y
  public :: meshfieldfilteroperationbase_apply_reconst1d_y_2
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  !
  private :: calc_filter_kenrnel
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
  integer, parameter, public :: filter_optrtype_convfilter           = 1
  integer, parameter, public :: filter_optrtype_reconstruct          = 2
  integer, parameter, public :: filter_optrtype_reconstruct2         = 3
  integer, parameter, public :: filter_optrtype_reconstruct2_gl      = 4
  integer, public, parameter :: filter_optrtype_interface_correction = 5
  integer, parameter, public :: filter_optrtype_modalfilter          = 6
 
contains
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_init( this, &
    Nnode_h1D,                                                      &
    FilterOptrType, FilterShape, FilterWidthFac, Nnode_h1D_reconst, &
    Nnode_h1D_GL, IF_r )
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    integer, intent(in) :: nnode_h1d
    character(*), intent(in) :: filteroptrtype
    character(*), intent(in) :: filtershape
    real(rp), intent(in) :: filterwidthfac
    integer, intent(in) :: nnode_h1d_reconst
    integer, intent(in), optional :: nnode_h1d_gl
    integer, intent(in), optional :: if_r
    !------------------------------
 
    select case(filteroptrtype)
    case ('ConvolFilter')
      call this%Prepair_FilterMat( nnode_h1d, filtershape, filterwidthfac )
      this%operator_type = filter_optrtype_convfilter
    case ('Reconstruction')
      call this%Prepair_ReconstMat( nnode_h1d, nnode_h1d_reconst )
      this%operator_type = filter_optrtype_reconstruct
      this%Nnode_h1D_reconst = nnode_h1d_reconst
    case ('Reconstruction2')
      call this%Prepair_ReconstMat2( nnode_h1d, nnode_h1d_reconst )
      this%operator_type = filter_optrtype_reconstruct2
      this%Nnode_h1D_reconst = nnode_h1d_reconst
    case ('Reconstruction2_GL')
      call this%Prepair_ReconstMat2_GL( nnode_h1d, nnode_h1d_reconst, nnode_h1d_gl )
      this%operator_type = filter_optrtype_reconstruct2_gl
      this%Nnode_h1D_reconst = nnode_h1d_reconst
    case ('InterfaceCorrection')
      call this%Prepair_InterfaceCorrection( nnode_h1d, if_r )
      this%operator_type = filter_optrtype_interface_correction
    case default
      log_info('MeshFieldFilterOperationBase_Init',*) "Unsupported filter operation is specified. Check!", trim(filteroptrtype)
      call prc_abort
    end select
    return
  end subroutine meshfieldfilteroperationbase_init
 
  subroutine meshfieldfilteroperationbase_final( this )
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    !------------------------------
 
    if ( allocated(this%FilterMat_h1D) ) then
      deallocate( this%FilterMat_h1D )
    end if
    if ( allocated(this%Minv_Ml_tr) ) then
      deallocate( this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat )
    end if
    if ( allocated(this%Ml_tr) ) then
      deallocate( this%Ml_tr, this%Mc_tr, this%Mr_tr )
    end if
    if ( allocated(this%IF_gL) ) then
      deallocate( this%IF_gL, this%IF_gR )
    end if
 
    return
  end subroutine meshfieldfilteroperationbase_final
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_filter1d_x( q, q0, Filter1D, Npx, Npy, Npz, Ne, NeA, Nnode_h1D )
    implicit none
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: nea
    integer, intent(in) :: nnode_h1d
    real(rp), intent(out) :: q(npx,npy,npz,nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: filter1d(nnode_h1d,-npx+1:2*npx)
 
    integer :: ke, pz, py, px, p
    real(rp) :: tmp
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,p,tmp) collapse(2)
    do ke=1, ne
    do pz=1, npz
      do py=1, npy
      do px=1, npx
        tmp = 0.0_rp
        do p=-npx+1, npx+npx
          tmp = tmp + filter1d(px,p) * q0(p,py,pz,ke)
        end do
        q(px,py,pz,ke) = tmp
      end do
      end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_filter1d_x
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_reconst1d_x( q, q0, Minv_Ml_tr, Minv_Mc_tr, Minv_Mr_tr, IntrpMat, &
    Npx, Npy, Npz, Ne, Npx_reconst, lmesh )
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: npx_reconst
    real(rp), intent(out) :: q(npx,npy,npz,lmesh%nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: minv_ml_tr(npx,npx_reconst)
    real(rp), intent(in) :: minv_mc_tr(npx,npx_reconst)
    real(rp), intent(in) :: minv_mr_tr(npx,npx_reconst)
    real(rp), intent(in) :: intrpmat(npx,npx_reconst)
 
    integer :: ke, pz, py, px
 
    integer :: i,j
 
    real(rp) :: tmp(npx_reconst)
    real(rp) :: s
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,i,j, tmp,s) collapse(2)
    do ke=1, ne
    do pz=1, npz
    do py=1, npy
      do i=1, npx_reconst
        s = 0.0_rp
        do j=1, npx
          s = s &
            + minv_ml_tr(j,i) * q0(j-npx,py,pz,ke) &
            + minv_mc_tr(j,i) * q0(    j,py,pz,ke) &
            + minv_mr_tr(j,i) * q0(j+npx,py,pz,ke)
        end do
        tmp(i) = s
      end do
      q(:,py,pz,ke) = matmul(intrpmat, tmp)
    end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_reconst1d_x
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_reconst1d_x_2( q, q0, Minv_Ml_tr, Minv_Mc_tr, Minv_Mr_tr, &
    Npx, Npy, Npz, Ne, Npx_reconst, lmesh )
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: npx_reconst
    real(rp), intent(out) :: q(npx,npy,npz,lmesh%nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: minv_ml_tr(npx,npx)
    real(rp), intent(in) :: minv_mc_tr(npx,npx)
    real(rp), intent(in) :: minv_mr_tr(npx,npx)
 
    integer :: ke, pz, py, px
 
    integer :: i,j
 
    real(rp) :: tmp(npx)
    real(rp) :: s
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,i,j, tmp,s) collapse(2)
    do ke=1, ne
    do pz=1, npz
    do py=1, npy
      do i=1, npx
        s = 0.0_rp
        do j=1, npx
          s = s &
            + minv_ml_tr(j,i) * q0(j-npx,py,pz,ke) &
            + minv_mc_tr(j,i) * q0(    j,py,pz,ke) &
            + minv_mr_tr(j,i) * q0(j+npx,py,pz,ke)
        end do
        tmp(i) = s
      end do
      q(:,py,pz,ke) = tmp(:)
    end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_reconst1d_x_2
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_interface_correction1d_x( q, q0, gL, gR, &
    Npx, Npy, Npz, Ne, lmesh )    
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: npx, npy, npz, ne    
    real(rp), intent(out) :: q(npx,npy,npz,lmesh%nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in)  :: gl(npx)
    real(rp), intent(in)  :: gr(npx)
 
    integer :: ke, px, py, pz
 
    real(rp) :: ql, qr
    real(rp) :: qpl, qpr
    real(rp) :: qstarl, qstarr
    real(rp) :: dql, dqr
    !-----------------------------------------------------
 
    !$omp parallel do private(ke, px,py,pz,qL,qR,qPL,qPR,qstarL,qstarR,dqL,dqR) collapse(3)
    do ke=1, ne
      do pz=1, npz
      do py=1, npy
        ! Current-element traces
        ql = q0(1,   py,pz,ke)
        qr = q0(npx, py,pz,ke)
        ! Neighbor traces
        qpl = q0(0,     py,pz,ke)
        qpr = q0(npx+1, py,pz,ke)
 
        ! Common interface values
        qstarl = 0.5_rp * ( ql + qpl )
        qstarr = 0.5_rp * ( qr + qpr )
 
        dql = qstarl - ql
        dqr = qstarr - qr
 
        do px=1, npx
          q(px,py,pz,ke) = q0(px,py,pz,ke) &
            + dql * gl(px)   &
            + dqr * gr(px)
        end do
 
      end do
      end do
    end do
 
    return
  end subroutine meshfieldfilteroperationbase_apply_interface_correction1d_x  
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_filter1d_y( q, q0, Filter1D, Npx, Npy, Npz, Ne, NeA, Nnode_h1D )
    implicit none
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: nea
    integer, intent(in) :: nnode_h1d
    real(rp), intent(out) :: q(npx,npy,npz,nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: filter1d(nnode_h1d,-npy+1:2*npy)
 
    integer :: ke, px, py, pz, p
    real(rp) :: tmp
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,p,tmp) collapse(2)
    do ke=1, ne
    do pz=1, npz
      do py=1, npy
      do px=1, npx
        tmp = 0.0_rp
        do p=-npy+1, npy+npy
          tmp = tmp + filter1d(py,p) * q0(px,p,pz,ke)
        end do
        q(px,py,pz,ke) = tmp
      end do
      end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_filter1d_y
 
  
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_reconst1d_y( q, q0, Minv_Ml_tr, Minv_Mc_tr, Minv_Mr_tr, IntrpMat, &
    Npx, Npy, Npz, Ne, Npy_reconst, lmesh )
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: npy_reconst
    real(rp), intent(out) :: q(npx,npy,npz,lmesh%nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: minv_ml_tr(npy,npy_reconst)
    real(rp), intent(in) :: minv_mc_tr(npy,npy_reconst)
    real(rp), intent(in) :: minv_mr_tr(npy,npy_reconst)
    real(rp), intent(in) :: intrpmat(npy,npy_reconst)
 
    integer :: ke, pz, py, px
 
    integer :: i,j
    real(rp) :: tmp(npy_reconst)
    real(rp) :: s
    real(rp) :: q0_l(npy), q0_c(npy), q0_r(npy)
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,i,j, tmp,s, &
    !$omp q0_L,q0_C,q0_R) collapse(2)
    do ke=1, ne
    do pz=1, npz
    do px=1, npx
      do j=1, npy
        q0_l(j) = q0(px,j-npy,pz,ke)
        q0_c(j) = q0(px,j    ,pz,ke)
        q0_r(j) = q0(px,j+npy,pz,ke)
      end do
      do i=1, npy_reconst
        s = 0.0_rp
        do j=1, npy
          s = s &
            + minv_ml_tr(j,i) * q0_l(j) &
            + minv_mc_tr(j,i) * q0_c(j) &
            + minv_mr_tr(j,i) * q0_r(j)
        end do
        tmp(i) = s
      end do
      q(px,:,pz,ke) = matmul(intrpmat, tmp)
    end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_reconst1d_y
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_apply_reconst1d_y_2( q, q0, Minv_Ml_tr, Minv_Mc_tr, Minv_Mr_tr, &
    Npx, Npy, Npz, Ne, Npy_reconst, lmesh )
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: npx, npy, npz, ne
    integer, intent(in) :: npy_reconst
    real(rp), intent(out) :: q(npx,npy,npz,lmesh%nea)
    real(rp), intent(in) :: q0(-npx+1:2*npx,-npy+1:2*npy,npz,ne)
    real(rp), intent(in) :: minv_ml_tr(npy,npy)
    real(rp), intent(in) :: minv_mc_tr(npy,npy)
    real(rp), intent(in) :: minv_mr_tr(npy,npy)
 
    integer :: ke, pz, py, px
 
    integer :: i,j
    real(rp) :: tmp(npy)
    real(rp) :: s
    real(rp) :: q0_l(npy), q0_c(npy), q0_r(npy)
    !-------------------------------------
 
    !$omp parallel do private(ke,px,py,pz,i,j, tmp,s, &
    !$omp q0_L,q0_C,q0_R) collapse(2)
    do ke=1, ne
    do pz=1, npz
    do px=1, npx
      do j=1, npy
        q0_l(j) = q0(px,j-npy,pz,ke)
        q0_c(j) = q0(px,j    ,pz,ke)
        q0_r(j) = q0(px,j+npy,pz,ke)
      end do
      do i=1, npy
        s = 0.0_rp
        do j=1, npy
          s = s &
            + minv_ml_tr(j,i) * q0_l(j) &
            + minv_mc_tr(j,i) * q0_c(j) &
            + minv_mr_tr(j,i) * q0_r(j)
        end do
        tmp(i) = s
      end do
      q(px,:,pz,ke) = tmp(:)
    end do
    end do
    end do
    return
  end subroutine meshfieldfilteroperationbase_apply_reconst1d_y_2
 
!-- Private routines -----------------
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_prepair_filter_matrix( this, Nnode_h1D, FilterShape, FilterWidthFac )
    use scale_polynomial, only: &
      polynomial_genlagrangepoly, &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    use scale_element_line, only: lineelement
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    integer, intent(in) :: nnode_h1d
    character(*), intent(in) :: filtershape
    real(rp), intent(in) :: filterwidthfac
 
    integer :: p1, p2
 
    real(rp) :: filterw
    real(rp) :: filterw2
    type(lineelement) :: elem1d
    type(lineelement) :: elem1d_intrp
 
    real(rp), allocatable :: lag(:,:)
    real(rp), allocatable :: filter_func(:)
 
    integer :: nintnode
    real(rp), allocatable :: r_int1d(:)
    real(rp), allocatable :: w_int1d(:)
 
    integer :: polyorder_h
    !--------------------------
 
    polyorder_h = nnode_h1d - 1
    call elem1d%Init( polyorder_h, .false. )
    call elem1d_intrp%Init( min(2*polyorder_h, 11), .false. )
    
    ! LOG_INFO("prepair_filter_matrix Pos:",*) elem1D%x1(:)
    ! LOG_INFO("prepair_filter_matrix IntWeight:",*) elem1D%IntWeight_lgl(:)
 
    allocate( this%FilterMat_h1D(elem1d%Np,-elem1d%Np+1:elem1d%Np+elem1d%Np) )
 
    nintnode = min(2*polyorder_h, 11)
    ! NGLnodes = elem1D_intrp%Np
    
    allocate( r_int1d(nintnode), w_int1d(nintnode) )
    r_int1d(:) = polynomial_gengausslegendrept( nintnode )
    w_int1d(:) = polynomial_gengausslegendreptintweight( nintnode )
    ! r_int1D(:) = elem1D_intrp%x1(:)
    ! w_int1D(:) = elem1D_intrp%IntWeight_lgl(:)
 
    allocate( filter_func(nintnode) )
 
    allocate( lag(nintnode,elem1d%PolyOrder+1) )
    lag(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, r_int1d(:) )
 
    this%FilterMat_h1D(:,:) = 0.0_rp
 
    filterw = filterwidthfac * 2.0_rp / real(elem1d%Np,kind=rp)
    this%hHaloSize = elem1d%Np
 
    ! LOG_INFO("prepair_filter_matrix FilterW:",*) filterW
 
    do p1=1, elem1d%Np
      !--
      call calc_filter_kenrnel( filter_func, &
        filtershape, r_int1d(:) - elem1d%x1(p1), filterw, nintnode )      
      do p2=1, elem1d%Np
        this%FilterMat_h1D(p1,p2) = sum( w_int1d(:) * lag(:,p2) * filter_func(:) )
      end do
      filterw2 = sum( w_int1d(:) * filter_func(:) )
 
      !--
      call calc_filter_kenrnel( filter_func, &
        filtershape, r_int1d(:) - 2.0_rp - elem1d%x1(p1), filterw, nintnode )
      
      ! this%FilterMat_h1D(p1,0) = sum( elem1D_intrp%IntWeight_lgl(:) * filter_func(:) )
      do p2=1, elem1d%Np
        this%FilterMat_h1D(p1,-elem1d%Np+p2) = sum( w_int1d(:) * lag(:,p2) * filter_func(:) )
      end do
      filterw2 = filterw2 + sum( w_int1d(:) * filter_func(:) )
 
      !--
      call calc_filter_kenrnel( filter_func, &
        filtershape, r_int1d(:) + 2.0_rp - elem1d%x1(p1), filterw, nintnode )
      
      ! this%FilterMat_h1D(p1,elem1D%Np+1) = sum( elem1D_intrp%IntWeight_lgl(:) * filter_func(:) )
      do p2=1, elem1d%Np
        this%FilterMat_h1D(p1,elem1d%Np+p2) = sum( w_int1d(:) * lag(:,p2) * filter_func(:) )
      end do
      filterw2 = filterw2 + sum( w_int1d(:) * filter_func(:) )
      
      ! LOG_INFO("prepair_filter_matrix",*) p1, ":", filterW2, ":", this%FilterMat_h1D(p1,:)
      this%FilterMat_h1D(p1,:) = this%FilterMat_h1D(p1,:) / filterw2 ! normalization
    end do
 
    call elem1d%Final()
    call elem1d_intrp%Final()
    return
  end subroutine meshfieldfilteroperationbase_prepair_filter_matrix
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_prepair_reconstruct_matrix( this, Nnode_h1D, Nnode_h1D_reconst )
    use scale_polynomial, only: &
      polynomial_genlagrangepoly, &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    use scale_element_line, only: lineelement
    use scale_element_modalfilter, only: modalfilter
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    integer, intent(in) :: nnode_h1d
    integer, intent(in) :: nnode_h1d_reconst
 
    integer :: p1, p2
 
    type(lineelement) :: elem1d
    type(lineelement) :: elem1d_reconst
 
    real(rp), allocatable :: lagr_l(:,:)
    real(rp), allocatable :: lagr_c(:,:)
    real(rp), allocatable :: lagr_r(:,:)
 
    real(rp), allocatable :: rec_lagr_l(:,:)
    real(rp), allocatable :: rec_lagr_c(:,:)
    real(rp), allocatable :: rec_lagr_r(:,:)
 
    integer :: nintnode
    real(rp), allocatable :: r_int1d(:)
    real(rp), allocatable :: w_int1d(:)
 
    real(rp) :: m_h1d_l(nnode_h1d_reconst,nnode_h1d)
    real(rp) :: m_h1d_c(nnode_h1d_reconst,nnode_h1d)
    real(rp) :: m_h1d_r(nnode_h1d_reconst,nnode_h1d)
 
    real(rp) :: minv(nnode_h1d_reconst,nnode_h1d_reconst)
    real(rp) :: tmpm(nnode_h1d_reconst,nnode_h1d)
 
    real(rp), allocatable :: x_int(:)
    real(rp) :: x_c(nnode_h1d)
 
    integer :: polyorder
    integer :: polyorder_reconst
    type(modalfilter) :: modalfilter1d
    !--------------------------
 
    polyorder = nnode_h1d - 1
    polyorder_reconst = nnode_h1d_reconst - 1
 
    call elem1d%Init( polyorder, .false. )    
    call elem1d_reconst%Init( polyorder_reconst, .false. )    
    call modalfilter1d%Init(elem1d_reconst, 0.0_rp, 2e1_rp, 16)
 
    this%hHaloSize = elem1d%Np
 
    ! NIntNode = elem1D%Np 
    nintnode = ceiling( 0.5_rp * ( polyorder + polyorder_reconst ) ) + 1
 
    allocate( r_int1d(nintnode), w_int1d(nintnode) )
    allocate( x_int(nintnode) )
    r_int1d(:) = polynomial_gengausslegendrept( nintnode )
    w_int1d(:) = polynomial_gengausslegendreptintweight( nintnode )
 
    !-
    allocate( lagr_l(nintnode,elem1d%PolyOrder+1), rec_lagr_l(nintnode,polyorder_reconst+1) )
    allocate( lagr_c(nintnode,elem1d%PolyOrder+1), rec_lagr_c(nintnode,polyorder_reconst+1) )
    allocate( lagr_r(nintnode,elem1d%PolyOrder+1), rec_lagr_r(nintnode,polyorder_reconst+1) )
    
    !-
!    x_int(:) = -1.0_RP + 0.25_RP * ( 1.0_RP + r_int1D(:) )
    x_int(:) = -1.0_rp + (1.0_rp / 3.0_rp ) * ( 1.0_rp + r_int1d(:) )
    rec_lagr_l(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
!    x_int(:) = 0.0_RP + 0.5_RP * ( 1.0_RP + r_int1D(:) )
    x_int(:) = r_int1d(:)
    lagr_l(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
    !-
!    x_int(:) = -0.5_RP + 0.5_RP * ( 1.0_RP + r_int1D(:) )
    x_int(:) = -1.0_rp/3.0_rp + (1.0_rp / 3.0_rp ) * ( 1.0_rp + r_int1d(:) )
    rec_lagr_c(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
    x_int(:) = r_int1d(:)
    lagr_c(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
    !-
!    x_int(:) = 0.5_RP + 0.25_RP * ( 1.0_RP + r_int1D(:) )
    x_int(:) = +1.0_rp/3.0_rp + (1.0_rp / 3.0_rp ) * ( 1.0_rp + r_int1d(:) )
    rec_lagr_r(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
    ! x_int(:) = -1.0_RP + 0.5_RP * ( 1.0_RP + r_int1D(:) )
    x_int(:) = r_int1d(:)
    lagr_r(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
    !$omp parallel do collapse(2)
    do p2=1, elem1d%Np
      do p1=1, elem1d_reconst%Np
        ! M_h1D_l(p1,p2) = 0.25_RP * sum( w_int1D(:) * lagr_l(:,p2) * rec_lagr_l(:,p1) )
        ! M_h1D_c(p1,p2) = 0.5_RP  * sum( w_int1D(:) * lagr_c(:,p2) * rec_lagr_c(:,p1) )
        ! M_h1D_r(p1,p2) = 0.25_RP * sum( w_int1D(:) * lagr_r(:,p2) * rec_lagr_r(:,p1) )
        m_h1d_l(p1,p2) = 1.0_rp/3.0_rp * sum( w_int1d(:) * lagr_l(:,p2) * rec_lagr_l(:,p1) )
        m_h1d_c(p1,p2) = 1.0_rp/3.0_rp * sum( w_int1d(:) * lagr_c(:,p2) * rec_lagr_c(:,p1) )
        m_h1d_r(p1,p2) = 1.0_rp/3.0_rp * sum( w_int1d(:) * lagr_r(:,p2) * rec_lagr_r(:,p1) )
      end do
    end do
    
    allocate( this%Minv_Ml_tr(elem1d%Np,nnode_h1d_reconst) )
    allocate( this%Minv_Mc_tr(elem1d%Np,nnode_h1d_reconst) )
    allocate( this%Minv_Mr_tr(elem1d%Np,nnode_h1d_reconst) )
    allocate( this%IntrpMat(elem1d%Np,nnode_h1d_reconst) )
 
    minv(:,:) = elem1d_reconst%invM(:,:)
    ! Minv(:,:) = matmul( modalFilter1D%FilterMat, Minv )
    
    tmpm(:,:) = matmul(minv, m_h1d_l)
    this%Minv_Ml_tr(:,:) = transpose(tmpm)
 
    tmpm(:,:) = matmul(minv, m_h1d_c)
    this%Minv_Mc_tr(:,:) = transpose(tmpm)
 
    tmpm(:,:) = matmul(minv, m_h1d_r)
    this%Minv_Mr_tr(:,:) = transpose(tmpm)
 
!    x_c(:) = -0.5_RP + 0.5_RP * ( 1.0_RP + elem1D%x1(:) )
    x_c(:) = - 1.0_rp/3.0_rp + (1.0_rp/3.0_rp) * ( 1.0_rp + elem1d%x1(:) )
    this%IntrpMat(:,:) = polynomial_genlagrangepoly( elem1d_reconst%PolyOrder, elem1d_reconst%x1, x_c(:) )
 
    !-
    call elem1d%Final()
    call elem1d_reconst%Final()
    call modalfilter1d%Final()
    return
  end subroutine meshfieldfilteroperationbase_prepair_reconstruct_matrix
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_prepair_reconstruct2_matrix( this, Nnode_h1D, Nnode_h1D_reconst )
    use scale_const, only: &
      eps => const_eps
    use scale_polynomial, only: &
      polynomial_genlagrangepoly, &
      polynomial_gengausslegendrept, &
      polynomial_gengausslegendreptintweight
    use scale_element_line, only: lineelement
    use scale_element_modalfilter, only: modalfilter
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    integer, intent(in) :: nnode_h1d
    integer, intent(in) :: nnode_h1d_reconst
 
    integer :: p0, p1, p2
    real(rp) :: x0
    real(rp) :: xr_rec, xl_rec
    real(rp) :: xr, xl
    real(rp) :: coef_l, coef_c, coef_r
 
    type(lineelement) :: elem1d
    type(lineelement) :: elem1d_reconst
 
    real(rp), allocatable :: lagr_l(:,:)
    real(rp), allocatable :: lagr_c(:,:)
    real(rp), allocatable :: lagr_r(:,:)
 
    real(rp), allocatable :: rec_lagr_l(:,:)
    real(rp), allocatable :: rec_lagr_c(:,:)
    real(rp), allocatable :: rec_lagr_r(:,:)
 
    integer :: nintnode
    real(rp), allocatable :: r_int1d(:)
    real(rp), allocatable :: w_int1d(:)
 
    real(rp) :: m_h1d_l(nnode_h1d_reconst,nnode_h1d)
    real(rp) :: m_h1d_c(nnode_h1d_reconst,nnode_h1d)
    real(rp) :: m_h1d_r(nnode_h1d_reconst,nnode_h1d)
 
    real(rp) :: minv(nnode_h1d_reconst,nnode_h1d_reconst)
 
    real(rp) :: minv_ml(nnode_h1d_reconst,nnode_h1d,nnode_h1d)
    real(rp) :: minv_mc(nnode_h1d_reconst,nnode_h1d,nnode_h1d)
    real(rp) :: minv_mr(nnode_h1d_reconst,nnode_h1d,nnode_h1d)
 
    real(rp) :: intrpmat(1,nnode_h1d_reconst)
    real(rp) :: tmpmat(1,nnode_h1d)
 
    real(rp), allocatable :: x_int(:)
    real(rp) :: x_c(1)
 
    integer :: polyorder
    integer :: polyorder_reconst
    type(modalfilter) :: modalfilter1d
    !--------------------------
 
    polyorder = nnode_h1d - 1
    polyorder_reconst = nnode_h1d_reconst - 1
 
    call elem1d%Init( polyorder, .false. )    
    call elem1d_reconst%Init( polyorder_reconst, .false. )    
    call modalfilter1d%Init(elem1d_reconst, 0.0_rp, 1d3, 16)
 
    this%hHaloSize = elem1d%Np
 
    ! NIntNode = elem1D%Np 
    nintnode = ceiling( 0.5_rp * ( polyorder + polyorder_reconst ) ) + 1
 
    allocate( r_int1d(nintnode), w_int1d(nintnode) )
    allocate( x_int(nintnode) )
    r_int1d(:) = polynomial_gengausslegendrept( nintnode )
    w_int1d(:) = polynomial_gengausslegendreptintweight( nintnode )
 
    !-
    allocate( lagr_l(nintnode,elem1d%PolyOrder+1), rec_lagr_l(nintnode,polyorder_reconst+1) )
    allocate( lagr_c(nintnode,elem1d%PolyOrder+1), rec_lagr_c(nintnode,polyorder_reconst+1) )
    allocate( lagr_r(nintnode,elem1d%PolyOrder+1), rec_lagr_r(nintnode,polyorder_reconst+1) )
    
    !-
    minv(:,:) = elem1d_reconst%invM(:,:)
    ! Minv(:,:) = matmul( modalFilter1D%FilterMat, Minv )
 
    do p0=1, elem1d%Np
      x0 = elem1d%x1(p0)
 
      xl_rec = - 1.0_rp
      xr_rec = xl_rec + 0.5_rp * ( 1.0_rp - x0 )
      coef_l = 0.5_rp * ( xr_rec - xl_rec )
      x_int(:) = xl_rec + coef_l * ( 1.0_rp + r_int1d(:) )
      rec_lagr_l(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
      xl = x0
      xr = 1.0_rp
      x_int(:) = xl + 0.5_rp * ( xr - xl ) * ( 1.0_rp + r_int1d(:) )
      lagr_l(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
      !-
      xl_rec = xr_rec
      xr_rec = xl_rec + 0.5_rp * 2.0_rp
      coef_c = 0.5_rp * ( xr_rec - xl_rec )
      x_int(:) = xl_rec + coef_c * ( 1.0_rp + r_int1d(:) )
      rec_lagr_c(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
      x_int(:) = r_int1d(:)
      lagr_c(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
      !-
      xl_rec = xr_rec
      xr_rec = xl_rec + 0.5_rp * ( x0 + 1.0_rp )
      coef_r = 0.5_rp * ( xr_rec - xl_rec )
      x_int(:) = xl_rec + coef_r * ( 1.0_rp + r_int1d(:) )
      rec_lagr_r(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int(:) )
 
      xl = -1.0_rp
      xr = x0
      x_int(:) = xl + 0.5_rp * ( xr - xl ) * ( 1.0_rp + r_int1d(:) )
      lagr_r(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int(:) )
 
      !$omp parallel do collapse(2)
      do p2=1, elem1d%Np
        do p1=1, elem1d_reconst%Np
          m_h1d_l(p1,p2) = coef_l * sum( w_int1d(:) * lagr_l(:,p2) * rec_lagr_l(:,p1) )
          m_h1d_c(p1,p2) = coef_c * sum( w_int1d(:) * lagr_c(:,p2) * rec_lagr_c(:,p1) )
          m_h1d_r(p1,p2) = coef_r * sum( w_int1d(:) * lagr_r(:,p2) * rec_lagr_r(:,p1) )
        end do
      end do
 
      minv_ml(:,:,p0) = matmul( minv, m_h1d_l )
      minv_mc(:,:,p0) = matmul( minv, m_h1d_c )
      minv_mr(:,:,p0) = matmul( minv, m_h1d_r )
    end do
 
    x_c(:) = 0.0_rp
    intrpmat(:,:) = polynomial_genlagrangepoly( elem1d_reconst%PolyOrder, elem1d_reconst%x1, x_c(:) )
 
    allocate( this%Ml_tr(elem1d%Np,elem1d%Np) )
    allocate( this%Mc_tr(elem1d%Np,elem1d%Np) )
    allocate( this%Mr_tr(elem1d%Np,elem1d%Np) )
 
    do p0=1, nnode_h1d
      tmpmat(:,:) = matmul(intrpmat, minv_ml(:,:,p0))
      this%Ml_tr(:,p0) = tmpmat(1,:)
 
      tmpmat(:,:) = matmul(intrpmat, minv_mc(:,:,p0))
      this%Mc_tr(:,p0) = tmpmat(1,:)
 
      tmpmat(:,:) = matmul(intrpmat, minv_mr(:,:,p0))
      this%Mr_tr(:,p0) = tmpmat(1,:)      
    end do
 
    !-
    call elem1d%Final()
    call elem1d_reconst%Final()
    call modalfilter1d%Final()
    return
  end subroutine meshfieldfilteroperationbase_prepair_reconstruct2_matrix
 
!OCL SERIAL
subroutine meshfieldfilteroperationbase_prepair_reconstruct2_gl_matrix( this, &
    Nnode_h1D, Nnode_h1D_reconst, Nnode_h1D_GL )
 
  use scale_polynomial, only: &
    polynomial_genlagrangepoly, &
    polynomial_gengausslegendrept, &
    polynomial_gengausslegendreptintweight
  use scale_element_line, only: lineelement
 
  implicit none
 
  class(meshfieldfilteroperationbase), intent(inout) :: this
  integer, intent(in) :: nnode_h1d
  integer, intent(in) :: nnode_h1d_reconst
  integer, intent(in), optional :: nnode_h1d_gl
 
  integer :: pg, p1, p2
 
  real(rp) :: x0
  real(rp) :: xr_rec, xl_rec
  real(rp) :: xr, xl
  real(rp) :: coef_l, coef_c, coef_r
 
  type(lineelement) :: elem1d
  type(lineelement) :: elem1d_reconst
 
  real(rp), allocatable :: lagr_l(:,:)
  real(rp), allocatable :: lagr_c(:,:)
  real(rp), allocatable :: lagr_r(:,:)
 
  real(rp), allocatable :: rec_lagr_l(:,:)
  real(rp), allocatable :: rec_lagr_c(:,:)
  real(rp), allocatable :: rec_lagr_r(:,:)
 
  integer :: nintnode
  real(rp), allocatable :: r_int1d(:)
  real(rp), allocatable :: w_int1d(:)
  real(rp), allocatable :: x_int(:)
 
  real(rp) :: m_h1d_l(nnode_h1d_reconst,nnode_h1d)
  real(rp) :: m_h1d_c(nnode_h1d_reconst,nnode_h1d)
  real(rp) :: m_h1d_r(nnode_h1d_reconst,nnode_h1d)
 
  real(rp) :: minv(nnode_h1d_reconst,nnode_h1d_reconst)
 
  real(rp) :: minv_ml(nnode_h1d_reconst,nnode_h1d)
  real(rp) :: minv_mc(nnode_h1d_reconst,nnode_h1d)
  real(rp) :: minv_mr(nnode_h1d_reconst,nnode_h1d)
 
  real(rp) :: intrpmat(1,nnode_h1d_reconst)
  real(rp) :: tmpmat(1,nnode_h1d)
 
  real(rp) :: x_c(1)
 
  integer :: polyorder
  integer :: polyorder_reconst
 
  integer :: nnode_h1d_gl_
  real(rp), allocatable :: x_gl(:)
  !------------------------------------------------------------
 
  polyorder         = nnode_h1d         - 1
  polyorder_reconst = nnode_h1d_reconst - 1
 
  if ( present(nnode_h1d_gl) ) then
    nnode_h1d_gl_ = nnode_h1d_gl
  else
    nnode_h1d_gl_ = nnode_h1d
  end if
 
  ! Original DG element:
  !     DOFs are located at LGL nodes.
  call elem1d%Init( polyorder, .false. )
 
  this%hHaloSize = elem1d%Np
 
 
  ! Polynomial space used for patch reconstruction.
  call elem1d_reconst%Init( polyorder_reconst, .false. )
 
  ! GL points at which reconstructed values are required.
  allocate( x_gl(nnode_h1d_gl_) )
  x_gl(:) = polynomial_gengausslegendrept( nnode_h1d_gl_ )
 
  ! Quadrature used to construct projection matrices.
  ! Keep the same choice as Reconstruction2.
  !
  nintnode = ceiling( &
       0.5_rp * real(polyorder + polyorder_reconst,kind=rp) ) + 1
 
  allocate( r_int1d(nintnode) )
  allocate( w_int1d(nintnode) )
  allocate( x_int(nintnode) )
 
  r_int1d(:) = polynomial_gengausslegendrept( nintnode )
  w_int1d(:) = polynomial_gengausslegendreptintweight( nintnode )
 
  allocate( lagr_l(nintnode,nnode_h1d) )
  allocate( lagr_c(nintnode,nnode_h1d) )
  allocate( lagr_r(nintnode,nnode_h1d) )
 
  allocate( rec_lagr_l(nintnode,nnode_h1d_reconst) )
  allocate( rec_lagr_c(nintnode,nnode_h1d_reconst) )
  allocate( rec_lagr_r(nintnode,nnode_h1d_reconst) )
 
  ! Inverse mass matrix in reconstruction space.
  minv(:,:) = elem1d_reconst%invM(:,:)
 
  ! The target point is always mapped onto the center of the reconstructed patch.
  x_c(1) = 0.0_rp
  intrpmat(:,:) = polynomial_genlagrangepoly( elem1d_reconst%PolyOrder, elem1d_reconst%x1, x_c )
 
 
  allocate( this%Ml_tr(nnode_h1d,nnode_h1d_gl_) )
  allocate( this%Mc_tr(nnode_h1d,nnode_h1d_gl_) )
  allocate( this%Mr_tr(nnode_h1d,nnode_h1d_gl_) )
  
  this%Ml_tr(:,:) = 0.0_rp
  this%Mc_tr(:,:) = 0.0_rp
  this%Mr_tr(:,:) = 0.0_rp
 
  !- Loop over arbitrary-order GL nodes.
  !
 
  do pg = 1, nnode_h1d_gl_
 
    x0 = x_gl(pg)
 
    !==========================================================
    ! Left element contribution
    !
    ! Original-coordinate interval: [x0, 1]
    ! is mapped to the left part of the reconstructed patch.
    !==========================================================
 
    xl_rec = -1.0_rp
    xr_rec = xl_rec + 0.5_rp * (1.0_rp - x0)
 
    coef_l = 0.5_rp * (xr_rec - xl_rec)
 
    x_int(:) = xl_rec &
             + coef_l * (1.0_rp + r_int1d(:))
 
    rec_lagr_l(:,:) = polynomial_genlagrangepoly( &
         polyorder_reconst, &
         elem1d_reconst%x1, &
         x_int )
 
    xl = x0
    xr = 1.0_rp
 
    x_int(:) = xl &
             + 0.5_rp * (xr-xl) * (1.0_rp+r_int1d(:))
 
    lagr_l(:,:) = polynomial_genlagrangepoly( &
         elem1d%PolyOrder, &
         elem1d%x1, &
         x_int )
 
    !==========================================================
    ! Center element contribution
    !==========================================================
 
    xl_rec = xr_rec
    xr_rec = xl_rec + 1.0_rp
 
    coef_c = 0.5_rp * (xr_rec-xl_rec)
 
    x_int(:) = xl_rec &
             + coef_c * (1.0_rp+r_int1d(:))
 
    rec_lagr_c(:,:) = polynomial_genlagrangepoly( &
         polyorder_reconst, &
         elem1d_reconst%x1, &
         x_int )
 
    x_int(:) = r_int1d(:)
 
    lagr_c(:,:) = polynomial_genlagrangepoly( &
         elem1d%PolyOrder, &
         elem1d%x1, &
         x_int )
 
    !==========================================================
    ! Right element contribution
    !
    ! Original-coordinate interval: [-1, x0]
    ! is mapped to the right part of the reconstructed patch.
    !==========================================================
 
    xl_rec = xr_rec
    xr_rec = xl_rec + 0.5_rp * (x0 + 1.0_rp)
    coef_r = 0.5_rp * (xr_rec-xl_rec)
    x_int(:) = xl_rec + coef_r * (1.0_rp+r_int1d(:))
    rec_lagr_r(:,:) = polynomial_genlagrangepoly( polyorder_reconst, elem1d_reconst%x1, x_int )
 
    xl = -1.0_rp
    xr = x0
    x_int(:) = xl + 0.5_rp * (xr-xl) * (1.0_rp+r_int1d(:))
    lagr_r(:,:) = polynomial_genlagrangepoly( elem1d%PolyOrder, elem1d%x1, x_int )
 
    !==========================================================
    ! Cross mass matrices
    !==========================================================
 
    do p2 = 1, nnode_h1d
      do p1 = 1, nnode_h1d_reconst
        m_h1d_l(p1,p2) = coef_l * sum( w_int1d(:) * lagr_l(:,p2) * rec_lagr_l(:,p1) )
        m_h1d_c(p1,p2) = coef_c * sum( w_int1d(:) * lagr_c(:,p2) * rec_lagr_c(:,p1) )
        m_h1d_r(p1,p2) = coef_r * sum( w_int1d(:) * lagr_r(:,p2) * rec_lagr_r(:,p1) )
      end do
    end do
 
    !==========================================================
    ! Reconstruction coefficients
    !==========================================================
 
    minv_ml(:,:) = matmul( minv, m_h1d_l )
    minv_mc(:,:) = matmul( minv, m_h1d_c )
    minv_mr(:,:) = matmul( minv, m_h1d_r )
 
    !==========================================================
    ! Evaluate reconstructed polynomial at patch center.
    !
    ! The center corresponds to x0 in the original
    ! central element.
    !==========================================================
 
    tmpmat(:,:) = matmul( intrpmat, minv_ml )
    this%Ml_tr(:,pg) = tmpmat(1,:)
 
    tmpmat(:,:) = matmul( intrpmat, minv_mc )
    this%Mc_tr(:,pg) = tmpmat(1,:)
 
    tmpmat(:,:) = matmul( intrpmat, minv_mr )
    this%Mr_tr(:,pg) = tmpmat(1,:)
  end do
 
  !------------------------------------------------------------
 
  deallocate( lagr_l, lagr_c, lagr_r )
  deallocate( rec_lagr_l, rec_lagr_c, rec_lagr_r )
  deallocate( r_int1d, w_int1d, x_int )
 
  call elem1d%Final()
  call elem1d_reconst%Final()
  return
end subroutine meshfieldfilteroperationbase_prepair_reconstruct2_gl_matrix  
 
!OCL SERIAL
  subroutine meshfieldfilteroperationbase_prepair_interface_correction( this, Nnode_h1D, IF_r )
    use scale_element_line, only: lineelement
    implicit none
    class(meshfieldfilteroperationbase), intent(inout) :: this
    integer, intent(in) :: nnode_h1d
    integer, intent(in), optional :: if_r
 
    type(lineelement) :: elem1d
    !---------------------------------------------
 
    if ( present(if_r) ) then
      this%IF_r = if_r
    else
      this%IF_r = nnode_h1d
    end if
    call elem1d%Init( nnode_h1d-1, .false. )
 
    this%hHaloSize = elem1d%Np
 
    allocate( this%IF_gL(elem1d%Np) )
    allocate( this%IF_gR(elem1d%Np) )
    this%IF_gL(:) = ( 0.5_rp * ( 1.0_rp - elem1d%x1(:) ) )**this%IF_r
    this%IF_gR(:) = ( 0.5_rp * ( 1.0_rp + elem1d%x1(:) ) )**this%IF_r
 
    call elem1d%Final()
    return
  end subroutine meshfieldfilteroperationbase_prepair_interface_correction
 
!OCL SERIAL
  subroutine calc_filter_kenrnel( filter_kernel, &
      FilterShape, x, filter_width, Np )
    implicit none
    integer, intent(in) :: np
    real(rp), intent(out) :: filter_kernel(np)
    character(*), intent(in) :: filtershape
    real(rp), intent(in) :: x(np)
    real(rp), intent(in) :: filter_width
    !------------------------------------------------
 
    select case(filtershape)
    case( "GAUSSIAN")
      filter_kernel(:) = exp( - (x(:)/filter_width)**2 )
    case( "TOPHAT")
      filter_kernel(:) = 0.5_rp * ( sign(1.0_rp, x(:) + 0.5_rp * filter_width) - sign(1.0_rp, x(:) - 0.5_rp * filter_width) )
    case default
      log_error("MeshFieldFilterOperation3D_calc_filter_kernel",*) "The specified FilterShape is not supported. Check!", filtershape
      call prc_abort
    end select
    return
  end subroutine calc_filter_kenrnel
 
end module scale_meshfield_filter_operation_base