scale_element_siacfilter Module Reference

module FElib / Element/ SIAC filter More...

Data Types
type	siac_filter
	Derived type to provide a SIAC filter. More...

Functions/Subroutines
subroutine	siac_filter_init (this, r, l, x_pts_per_elem, elem1d)
	Initialize a object for SIAC filter.

Detailed Description

module FElib / Element/ SIAC filter

Description

A module to provide a Smoothness-Increasing Accuracy-Increasing (SIAC) filter

Reference

• Cockburn et al. 2003: Enhanced accuracy by post-processing for finite element methods for hyperbolic equations. Mathematics of Computation, 72(242), 577–506.

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

siac_filter_init()

subroutine scale_element_siacfilter::siac_filter_init	(	class(siac_filter), intent(inout)	this,
		integer, intent(in)	r,
		integer, intent(in)	l,
		real(rp), dimension(:), intent(in)	x_pts_per_elem,
		class(lineelement), intent(in)	elem1d )

Initialize a object for SIAC filter.

Parameters

r	the number of B-splines - 1
l	the order of B-spline

Definition at line 82 of file scale_element_SIACfilter.F90.

    use scale_polynomial, only: &
      polynomial_genlagrangepoly, &
      polynomial_genlegendrepoly
    implicit none
    class(SIAC_filter), intent(inout) :: this
    integer, intent(in) :: r
    integer, intent(in) :: l
    real(RP), intent(in) :: x_pts_per_elem(:)
    class(LineElement), intent(in) :: elem1D
 
    integer :: i
    integer :: p
    integer :: m
 
    real(RP), allocatable :: IntrpMat_dummy(:,:)
    real(RP) :: x0, x1, x2
    real(RP), allocatable :: int_x_tmp(:)
 
    real(RP), allocatable :: P1D_ori(:,:)
    !--------------------------------------
 
    this%spline_r = r
    this%spline_num = r + 1
    this%spline_ord = l
    
    this%KernelHalfW = ceiling(0.5_rp * real(r+l,kind=rp))
 
    this%Npts_per_elem = size(x_pts_per_elem)
 
    !---
    allocate( this%x_pts_per_elem(this%Npts_per_elem) )
    this%x_pts_per_elem(:) = x_pts_per_elem(:)
 
    this%NintGLPt = ceiling( ( r + l ) / 2.0_rp )
    allocate( intrpmat_dummy(this%NintGLPt,elem1d%Np) )
    allocate( this%int_x(this%NintGLPt), this%int_w(this%NintGLPt) )
 
    intrpmat_dummy(:,:) = elem1d%GenIntGaussLegendreIntrpMat( this%NintGLPt, this%int_w, this%int_x )
 
    allocate( this%IntrpMat(this%NintGLPt,elem1d%Np,2,this%Npts_per_elem) )
    allocate( int_x_tmp(this%NintGLPt) )
    allocate( p1d_ori(this%NintGLPt,elem1d%Np) )
    do i=1, this%Npts_per_elem
      x0 = - 1.0_rp; x2 = 1.0_rp
      x1 = x_pts_per_elem(i)
 
      int_x_tmp(:) = x0 + 0.5_rp * (x1 - x0) * ( 1.0_rp + this%int_x(:) )
      p1d_ori(:,:) = polynomial_genlegendrepoly( elem1d%PolyOrder, int_x_tmp(:) )
      do p=1, elem1d%Np
        p1d_ori(:,p) = p1d_ori(:,p) * sqrt(real(p-1,kind=rp) + 0.5_rp)
      end do
      this%IntrpMat(:,:,1,i) = matmul( p1d_ori, elem1d%invV )
      ! this%IntrpMat(:,:,1,i) = Polynomial_GenLagrangePoly( elem1D%PolyOrder, elem1D%x1, int_x_tmp(:) )
      
      int_x_tmp(:) = x1 + 0.5_rp * (x2 - x1) * ( 1.0_rp + this%int_x(:) )
      p1d_ori(:,:) = polynomial_genlegendrepoly( elem1d%PolyOrder, int_x_tmp(:) )
      do p=1, elem1d%Np
        p1d_ori(:,p) = p1d_ori(:,p) * sqrt(real(p-1,kind=rp) + 0.5_rp)
      end do
      this%IntrpMat(:,:,2,i) = matmul( p1d_ori, elem1d%invV )
      ! this%IntrpMat(:,:,2,i) = Polynomial_GenLagrangePoly( elem1D%PolyOrder, elem1D%x1, int_x_tmp(:) )
    end do
 
    !---
    allocate( this%kernel_func_coef(0:r) )
    call calculate_kernel_func_coef( this%kernel_func_coef, &
        r, l, this%int_x, this%int_w, this%NintGLPt )
    
    allocate( this%kernel_func(this%NintGLPt,2,-this%KernelHalfW:this%KernelHalfW,this%Npts_per_elem) )
 
    do i=1, this%Npts_per_elem
      call construct_kernel_func( this%kernel_func(:,:,:,i), &
        r, l, this%kernel_func_coef, this%KernelHalfW, &
        this%int_x, this%NintGLPt, x_pts_per_elem(i) )
    end do
 
    log_info("SIAC_filter_Init",*) "r, l=", r, l        
    log_info("SIAC_filter_Init",*) "KernelFunc coef:", this%kernel_func_coef
    do i=1, this%Npts_per_elem
      log_info("SIAC_filter_Init",*) "--- KernelFunc xi=", x_pts_per_elem(i)
      do m=-this%KernelHalfW,this%KernelHalfW
        log_info("SIAC_filter_Init",*) this%kernel_func(:,1,m,i), ":", this%kernel_func(:,2,m,i)
      end do
    end do
 
    do i=1, this%Npts_per_elem
      log_info("SIAC_filter_Init",*) "--- InterpMat xi=", x_pts_per_elem(i)
      log_info("SIAC_filter_Init",*) "L", this%IntrpMat(1,:,1,i)
      log_info("SIAC_filter_Init",*) "L", this%IntrpMat(2,:,1,i)
      log_info("SIAC_filter_Init",*) "R", this%IntrpMat(1,:,2,i)
      log_info("SIAC_filter_Init",*) "R", this%IntrpMat(2,:,2,i)
    end do
 
    return

References scale_polynomial::polynomial_genlagrangepoly(), and scale_polynomial::polynomial_genlegendrepoly().