77 use scale_prc
78 implicit none
79 class(ElementOptrTensorProd3DInterp), intent(inout) :: this
80 type(ElementBase3D), intent(in), target :: elem
81 integer, intent(in) :: Polyorder_h_ip
82
83 real(RP) :: r_int1D_i(Polyorder_h_ip+1)
84 real(RP) :: r_int1Dw_i(Polyorder_h_ip+1)
85 real(RP) :: P_int1D_ori_h(Polyorder_h_ip+1,elem%PolyOrder_h+1)
86 real(RP) :: Vint1D(Polyorder_h_ip+1,elem%PolyOrder_h+1)
87
88 real(RP) :: P1D_gl_h(Polyorder_h_ip+1,Polyorder_h_ip+1)
89 real(RP) :: P1D_gll_h(elem%Nnode_h1D,Polyorder_h_ip+1)
90 real(RP) :: Vint1D_gl(elem%Nnode_h1D,Polyorder_h_ip+1)
91
92 integer :: p1, p1_
93 integer :: n_, l_
94
95 type(LineElement) :: elem1D
96
97
98 this%elem => elem
99 this%PolyOrder_h_intrp = polyorder_h_ip
100 this%NnodeH1D_GL = polyorder_h_ip + 1
101 this%Np3D_hGLvGLL = this%NnodeH1D_GL**2 * elem%Nnode_v
102 this%Np3D = elem%Np
103
104 allocate( this%IntrpMat1D_gll2gl(this%NnodeH1D_GL,elem%Nnode_h1D) )
105 allocate( this%IntrpMat1D_gll2gl_tr(elem%Nnode_h1D,this%NnodeH1D_GL) )
106
107
108 call elem1d%Init( elem%PolyOrder_h, .false. )
112
113 do p1_=1,this%NnodeH1D_GL
114 n_ = p1_
115 do p1=1, elem%Nnode_h1D
116 l_ = p1
117 vint1d(n_,l_) = p_int1d_ori_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
118 end do
119 end do
120 this%IntrpMat1D_gll2gl(:,:) = matmul( vint1d, elem1d%invV )
121 this%IntrpMat1D_gll2gl_tr(:,:) = transpose(this%IntrpMat1D_gll2gl)
122
123
124 allocate( this%V1D_GL(this%NnodeH1D_GL,polyorder_h_ip+1) )
125 allocate( this%inv_V1D_GL(polyorder_h_ip+1,this%NnodeH1D_GL) )
126 allocate( this%IntrpMat1D_gl2gll(elem%Nnode_h1D,this%NnodeH1D_GL) )
127 allocate( this%IntrpMat1D_gl2gll_tr(this%NnodeH1D_GL,elem%Nnode_h1D) )
128
131
132 do p1_=1, this%NnodeH1D_GL
133 n_ = p1_
134 do p1=1, polyorder_h_ip+1
135 l_ = p1
136 this%V1D_GL(n_,l_) = p1d_gl_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
137 end do
138 end do
140
141 vint1d_gl(:,:) = 0.0_rp
142 do p1_=1,this%elem%Nnode_h1D
143 n_ = p1_
144 do p1=1, this%NnodeH1D_GL
145 l_ = p1
146 vint1d_gl(n_,l_) = p1d_gll_h(p1_,p1) * sqrt( real(p1-1,kind=rp) + 0.5_rp )
147 end do
148 end do
149 this%IntrpMat1D_gl2gll(:,:) = matmul( vint1d_gl, this%inv_V1D_GL )
150 this%IntrpMat1D_gl2gll_tr(:,:) = transpose(this%IntrpMat1D_gl2gll(:,:))
151
152
153 call elem1d%Final()
154 return
module FElib / Element / line
Module common / Linear algebra.
real(rp) function, dimension(size(a, 1), size(a, 2)), public linalgebra_inv(a)
Calculate a inversion of matrix A.
Module common / Polynomial.
real(rp) function, dimension(nord), public polynomial_gengausslegendreptintweight(nord)
A function to calculate the Gauss-Legendre (GL) weights.
real(rp) function, dimension(size(x), nord+1), public polynomial_genlegendrepoly(nord, x)
A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points.
real(rp) function, dimension(nord), public polynomial_gengausslegendrept(nord)
A function to calculate the Gauss-Legendre (GL) points.
Derived type representing a line element.