FE-Project
Loading...
Searching...
No Matches
Functions/Subroutines
scale_polynomial Module Reference

Module common / Polynomial. More...

Functions/Subroutines

real(rp) function, dimension(size(x), nord+1), public polynomial_genlagrangepoly (nord, x_lgl, x)
 A function to obtain the Lagrange basis functions related to the Gauss-Legendre-Lobatto (GLL) points.
real(rp) real(rp) function, dimension(n, k), public polynomial_gendlagrangepoly_lglpt (nord, x_lgl)
 Differential values of Lagrange basis functions at the GLL points.
real(rp) real(rp) function, dimension(m, k), public polynomial_gendlagrangepoly (nord, x_lgl, x_eval)
 Differential values of Lagrange basis functions defined at GLL nodes, evaluated at arbitrary points.
subroutine, public polynomial_genlegendrepoly_sub (nord, x, p)
 A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points.
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(size(x), nord+1), public polynomial_gendlegendrepoly (nord, x, p)
 A function to obtain differential values of Legendre polynomials which are evaluated at arbitrary points.
real(rp) function, dimension(nord+1), public polynomial_gengausslobattopt (nord)
 A function to calculate the Legendre-Gauss-Lobatto (LGL) points.
real(rp) function, dimension(nord+1), public polynomial_gengausslobattoptintweight (nord)
 A function to calculate the Gauss-Lobbato weights.
real(rp) function, dimension(nord), public polynomial_gengausslegendrept (nord)
 A function to calculate the Gauss-Legendre (GL) points.
real(rp) function, dimension(nord), public polynomial_gengausslegendreptintweight (nord)
 A function to calculate the Gauss-Legendre (GL) weights.

Detailed Description

Module common / Polynomial.

Description
A module to provide utilities for polynomials
Reference
Author
Yuta Kawai, Team SCALE

Function/Subroutine Documentation

◆ polynomial_genlagrangepoly()

real(rp) function, dimension(size(x), nord+1), public scale_polynomial::polynomial_genlagrangepoly ( integer, intent(in) nord,
real(rp), dimension(nord+1), intent(in) x_lgl,
real(rp), dimension(:), intent(in) x )

A function to obtain the Lagrange basis functions related to the Gauss-Legendre-Lobatto (GLL) points.

Parameters
NordOrder of Lagrange polynomial
x_lglPositions of GLL points
xPositions where the Lagrange basis functions are evaluated

Definition at line 68 of file scale_polynomial.F90.

69 implicit none
70
71 integer, intent(in) :: Nord
72 real(RP), intent(in) :: x_lgl(Nord+1)
73 real(RP), intent(in) :: x(:)
74 real(RP) :: l(size(x), Nord+1)
75
76 integer :: n, i
77 real(RP) :: P_lgl(Nord+1,Nord+1)
78 real(RP) :: P(size(x),Nord+1)
79 real(RP) :: Pr(size(x),Nord+1)
80
81 ! real(RP) :: w(Nord+1)
82 ! real(RP) :: int_w(Nord+1)
83 !---------------------------------------------------------------------------
84
85 p_lgl(:,:) = polynomial_genlegendrepoly(nord, x_lgl)
86 p(:,:) = polynomial_genlegendrepoly(nord, x)
87 pr(:,:) = polynomial_gendlegendrepoly(nord, x, p)
88
89 do n=1, nord+1
90 do i=1, size(x)
91 if ( abs(x(i)-x_lgl(n)) < 1e-15_rp ) then
92 l(i,n) = 1.0_rp
93 else
94 l(i,n) = &
95 (x(i) - 1.0_rp)*(x(i) + 1.0_rp)*pr(i,nord+1) &
96 / (dble(nord*(nord+1))*p_lgl(n,nord+1)*(x(i) - x_lgl(n)))
97 end if
98 end do
99 end do
100
101 !- Calculate interpolation coefficient based on barycentric form
102 ! Eq. (3.46) in Wang, Huybrechs & Vandewalle (2012): Explicit barycentric weights for polynomial interpolation in the roots or extrema of classical orthogonal polynomials
103 ! int_w(:) = Polynomial_GenGaussLobattoPtIntWeight( Nord )
104 ! do n=1, Nord+1
105 ! w(n) = (-1)**mod(n-1,2) * sqrt(int_w(n))
106 ! end do
107 ! do n=1, Nord+1
108 ! do i=1, size(x)
109 ! if ( abs(x(i)-x_lgl(n)) < 1E-16_RP ) then
110 ! l(i,n) = 1.0_RP
111 ! else
112 ! l(i,n) = &
113 ! ( w(n) / ( x(i) - x_lgl(n) ) ) &
114 ! / sum( w(:) / ( x(i) - x_lgl(:) ) )
115 ! end if
116 ! end do
117 ! end do
118
119 return

References polynomial_gendlegendrepoly(), and polynomial_genlegendrepoly().

Referenced by scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), scale_meshfield_filter_operation_base::meshfieldfilteroperationbase_apply_reconst1d_y_2(), scale_mesh_hierarchy_base::meshhierarchy_construct_pmg_mat1d(), scale_element_quadrilateral::quadrilateralelement_init(), and scale_element_siacfilter::siac_filter_init().

◆ polynomial_gendlagrangepoly_lglpt()

real(rp) real(rp) function, dimension(n,k), public scale_polynomial::polynomial_gendlagrangepoly_lglpt ( integer, intent(in) nord,
real(rp), dimension(nord+1), intent(in) x_lgl )

Differential values of Lagrange basis functions at the GLL points.

Parameters
NordOrder of Lagrange polynomial
x_lglPositions of GLL points

Definition at line 127 of file scale_polynomial.F90.

128 implicit none
129
130 integer, intent(in) :: Nord
131 real(RP), intent(in) :: x_lgl(Nord+1)
132 real(RP) :: lr(Nord+1, Nord+1) !> lr(n,k) = derivative of k-th Lagrange basis at x_lgl(n)
133
134 integer :: n, k
135 real(RP) :: P(Nord+1,Nord+1)
136 real(RP) :: lr_nn
137 !---------------------------------------------------------------------------
138
139 p(:,:) = polynomial_genlegendrepoly(nord, x_lgl)
140
141 do n=1, nord+1
142 lr_nn = 0.0_rp
143 do k=1, nord+1
144 if (k==1 .and. n==1) then
145 lr(k,n) = - 0.25_rp*dble(nord*(nord+1))
146 else if (k==nord+1 .and. n==nord+1) then
147 lr(k,n) = + 0.25_rp*dble(nord*(nord+1))
148 else if (k==n) then
149 lr(k,n) = 0.0_rp
150 else
151 lr(k,n) = p(n,nord+1)/(p(k,nord+1)*(x_lgl(n) - x_lgl(k)))
152 end if
153
154 if ( k /= n ) then
155 lr_nn = lr_nn + lr(k,n)
156 end if
157 end do
158 lr(n,n) = - lr_nn
159 end do
160
161 return

References polynomial_genlegendrepoly().

Referenced by scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), and scale_element_quadrilateral::quadrilateralelement_init().

◆ polynomial_gendlagrangepoly()

real(rp) real(rp) function, dimension(m,k), public scale_polynomial::polynomial_gendlagrangepoly ( integer, intent(in) nord,
real(rp), dimension(nord+1), intent(in) x_lgl,
real(rp), dimension(:), intent(in) x_eval )

Differential values of Lagrange basis functions defined at GLL nodes, evaluated at arbitrary points.

Parameters
NordPolynomial order
x_lglGLL interpolation nodes
NevalNumber of evaluation points
x_evalEvaluation points

Definition at line 171 of file scale_polynomial.F90.

172 implicit none
173 integer, intent(in) :: Nord
174 real(RP), intent(in) :: x_lgl(Nord+1)
175 real(RP), intent(in) :: x_eval(:)
176 real(RP) :: lr(size(x_eval),Nord+1) !> lr(m,k) = derivative of k-th Lagrange basis at x_eval(m)
177
178 real(RP) :: P(Nord+1,Nord+1)
179 real(RP) :: bw(Nord+1) ! Barycentric weights. Common normalization factor is unnecessary.
180
181 real(RP) :: a(Nord+1)
182 real(RP) :: S1, S2
183 real(RP) :: dx
184 real(RP) :: tol
185
186 integer :: k, j, m
187 integer :: knode
188 !-------------------------------------------------------------
189
190 p(:,:) = polynomial_genlegendrepoly(nord, x_lgl)
191 bw(:) = 1.0_rp / p(:,nord+1)
192
193 tol = 100.0_rp * epsilon(1.0_rp)
194
195 do m = 1, size(x_eval)
196
197 !- Check whether evaluation point coincides with a GLL node
198 knode = 0
199 do k=1, nord+1
200 if ( abs(x_eval(m)-x_lgl(k)) <= tol ) then
201 knode = k
202 exit
203 end if
204 end do
205
206 if ( knode /= 0 ) then
207 !* Evaluation exactly at a GLL node. Use the existing nodal formula.
208 do k = 1, nord+1
209 if ( k /= knode ) then
210 lr(m,k) = p(knode,nord+1) / ( p(k,nord+1) * ( x_lgl(knode)-x_lgl(k) ) )
211 else
212 lr(m,k) = 0.0_rp
213 end if
214 end do
215
216 ! Sum of derivatives of all basis functions is zero: sum_k l'_k(x) = 0
217 lr(m,knode) = -sum(lr(m,:))
218
219 else
220 !* General evaluation point
221 ! l'_k(x) = l_k(x) * [ ( sum_j w_j/(x - x_j,lgl)**2 ) / [sum_j a_j(x)] - 1 / ( x - x_k,lgl ) ]
222 ! where
223 ! l_k(x) = a_k(x) / sum_j a_j(x),
224 ! a_k(x) = w_k / ( x - x_k,lgl )
225
226 s1 = 0.0_rp; s2 = 0.0_rp
227 do j = 1, nord+1
228 dx = x_eval(m) - x_lgl(j)
229
230 a(j) = bw(j) / dx
231 s1 = s1 + a(j)
232 s2 = s2 + bw(j) / ( dx * dx )
233 end do
234
235 !　l'_k(x)
236 do k = 1, nord+1
237 dx = x_eval(m) - x_lgl(k)
238 lr(m,k) = ( a(k) / s1 ) * ( s2 / s1 - 1.0_rp / dx )
239 end do
240
241 end if
242
243 end do
244 return

References polynomial_genlegendrepoly().

◆ polynomial_genlegendrepoly_sub()

subroutine, public scale_polynomial::polynomial_genlegendrepoly_sub ( integer, intent(in) nord,
real(rp), dimension(:), intent(in) x,
real(rp), dimension(size(x), nord+1), intent(out) p )

A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points.

Parameters
NordOrder of Lagrange polynomial
xPositions where the Legendre polynomials are evaluated
PValues of the Legendre polynomials at x

Definition at line 253 of file scale_polynomial.F90.

254 implicit none
255
256 integer, intent(in) :: Nord
257 real(RP), intent(in) :: x(:)
258 real(RP), intent(out) :: P(size(x), Nord+1)
259
260 integer :: n
261 !---------------------------------------------------------------------------
262
263 if (nord < 0) then
264 log_error("Polynomial_GenLegendrePoly",*) "Nord must be larger than 0."
265 end if
266
267 p(:,1) = 1.0_rp
268 if (nord==0) return
269
270 p(:,2) = x(:)
271 do n=2, nord
272 p(:,n+1) = ( dble(2*n-1)*x(:)*p(:,n) - dble(n-1)*p(:,n-1) )/dble(n)
273 end do
274
275 return

Referenced by scale_file_common_meshfield::file_common_meshfield_put_field1d_cartesbuf(), scale_file_common_meshfield::file_common_meshfield_put_field2d_cartesbuf(), scale_file_common_meshfield::file_common_meshfield_put_field3d_cartesbuf(), and polynomial_genlegendrepoly().

◆ polynomial_genlegendrepoly()

real(rp) function, dimension(size(x), nord+1), public scale_polynomial::polynomial_genlegendrepoly ( integer, intent(in) nord,
real(rp), dimension(:), intent(in) x )

A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points.

Parameters
NordOrder of Lagrange polynomial
xPositions where the Legendre polynomials are evaluated
PValues of the Legendre polynomials at x

Definition at line 284 of file scale_polynomial.F90.

285 implicit none
286
287 integer, intent(in) :: Nord
288 real(RP), intent(in) :: x(:)
289 real(RP) :: P(size(x), Nord+1)
290 !---------------------------------------------------------------------------
291
292 call polynomial_genlegendrepoly_sub( nord, x(:), & ! (in)
293 p(:,:) ) ! (out)
294 return

References polynomial_genlegendrepoly_sub().

Referenced by scale_element_operation_tensorprod3d_interp::elementoptrtensorprod3dinterp_init(), scale_file_common_meshfield::file_common_meshfield_put_field2d_cubedsphere_cartesbuf(), scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), polynomial_gendlagrangepoly(), polynomial_gendlagrangepoly_lglpt(), polynomial_gengausslegendreptintweight(), polynomial_gengausslobattoptintweight(), polynomial_genlagrangepoly(), scale_element_quadrilateral::quadrilateralelement_init(), and scale_element_siacfilter::siac_filter_init().

◆ polynomial_gendlegendrepoly()

real(rp) function, dimension(size(x), nord+1), public scale_polynomial::polynomial_gendlegendrepoly ( integer, intent(in) nord,
real(rp), dimension(:), intent(in) x,
real(rp), dimension(:,:), intent(in) p )

A function to obtain differential values of Legendre polynomials which are evaluated at arbitrary points.

Parameters
NordOrder of Lagrange polynomial
xPositions where the Legendre polynomials are evaluated
PValues of the Legendre polynomials at x

Definition at line 303 of file scale_polynomial.F90.

304 implicit none
305
306 integer, intent(in) :: Nord
307 real(RP), intent(in) :: x(:)
308 real(RP), intent(in) :: P(:,:)
309 real(RP) :: GradP(size(x), Nord+1)
310
311 integer :: n
312 !---------------------------------------------------------------------------
313
314 if (nord < 0) then
315 log_error("Polynomial_GenDLegendrePoly",*) "Nord must be larger than 0."
316 end if
317
318 gradp(:,1) = 0.0_rp
319 if (nord == 0) return
320
321 gradp(:,2) = 1.0_rp
322 do n=2, nord
323 gradp(:,n+1) = 2.0_rp*x(:)*gradp(:,n) - gradp(:,n-1) + p(:,n)
324 end do
325
326 return

Referenced by scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), polynomial_gengausslegendreptintweight(), polynomial_genlagrangepoly(), and scale_element_quadrilateral::quadrilateralelement_init().

◆ polynomial_gengausslobattopt()

real(rp) function, dimension(nord+1), public scale_polynomial::polynomial_gengausslobattopt ( integer, intent(in) nord)

A function to calculate the Legendre-Gauss-Lobatto (LGL) points.

Parameters
NordOrder of Lagrange polynomial
ptsPosition of the LGL points

Definition at line 334 of file scale_polynomial.F90.

335 implicit none
336
337 integer, intent(in) :: Nord
338 real(RP) :: pts(Nord+1)
339
340 integer :: N1
341 !---------------------------------------------------------------------------
342
343 pts(1) = -1.0_rp; pts(nord+1) = 1.0_rp
344 if (nord==1) return
345
346 call gen_jacobigaussquadraturepts( 1, 1, nord-2, pts(2:nord) )
347 return

Referenced by scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), polynomial_gengausslobattoptintweight(), and scale_element_quadrilateral::quadrilateralelement_init().

◆ polynomial_gengausslobattoptintweight()

real(rp) function, dimension(nord+1), public scale_polynomial::polynomial_gengausslobattoptintweight ( integer, intent(in) nord)

A function to calculate the Gauss-Lobbato weights.

Parameters
NordOrder of Lagrange polynomial
int_weight_lglGauss-Lobbato weights

Definition at line 355 of file scale_polynomial.F90.

356 implicit none
357
358 integer, intent(in) :: Nord
359 real(RP) :: int_weight_lgl(Nord+1)
360
361 real(RP) :: lglPts1D(Nord+1)
362 real(RP) :: P1D_ori(Nord+1, Nord+1)
363 !---------------------------------------------------------------------------
364
365 lglpts1d(:) = polynomial_gengausslobattopt( nord )
366 p1d_ori(:,:) = polynomial_genlegendrepoly( nord, lglpts1d )
367
368 int_weight_lgl(:) = 2.0_rp/(dble(nord*(nord+1))*p1d_ori(:,nord+1)**2)
369
370 return

References polynomial_gengausslobattopt(), and polynomial_genlegendrepoly().

Referenced by scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_init(), scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_gen_intweight(), scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), and scale_element_quadrilateral::quadrilateralelement_init().

◆ polynomial_gengausslegendrept()

real(rp) function, dimension(nord), public scale_polynomial::polynomial_gengausslegendrept ( integer, intent(in) nord)

A function to calculate the Gauss-Legendre (GL) points.

Parameters
NordOrder of the Legendre polynomial
ptsPosition of the GL points

Definition at line 378 of file scale_polynomial.F90.

379 implicit none
380
381 integer, intent(in) :: Nord
382 real(RP) :: pts(Nord)
383 !---------------------------------------------------------------------------
384
385 call gen_jacobigaussquadraturepts( 0, 0, nord-1, pts(:) )
386 return

Referenced by scale_element_operation_tensorprod3d_interp::elementoptrtensorprod3dinterp_init(), scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), scale_meshfield_filter_operation_base::meshfieldfilteroperationbase_apply_reconst1d_y_2(), scale_mesh_hierarchy_base::meshhierarchy_construct_pmg_mat1d(), polynomial_gengausslegendreptintweight(), and scale_element_quadrilateral::quadrilateralelement_init().

◆ polynomial_gengausslegendreptintweight()

real(rp) function, dimension(nord), public scale_polynomial::polynomial_gengausslegendreptintweight ( integer, intent(in) nord)

A function to calculate the Gauss-Legendre (GL) weights.

Parameters
NordOrder of the Legendre polynomial
int_weight_glGauss-Legendre weights

Definition at line 394 of file scale_polynomial.F90.

395 implicit none
396
397 integer, intent(in) :: Nord
398 real(RP) :: int_weight_gl(Nord)
399
400 real(RP) :: glPts1D(Nord)
401 real(RP) :: P1D_ori(Nord, Nord+1)
402 real(RP) :: dP1D_ori(Nord, Nord+1)
403 !---------------------------------------------------------------------------
404
405 glpts1d(:) = polynomial_gengausslegendrept( nord )
406 p1d_ori(:,:) = polynomial_genlegendrepoly( nord, glpts1d )
407 dp1d_ori(:,:) = polynomial_gendlegendrepoly( nord, glpts1d, p1d_ori )
408
409 int_weight_gl(:) = 2.0_rp / ( (1.0_rp - glpts1d(:)**2) * dp1d_ori(:,nord+1)**2 )
410
411 return

References polynomial_gendlegendrepoly(), polynomial_gengausslegendrept(), and polynomial_genlegendrepoly().

Referenced by scale_element_operation_tensorprod3d_interp::elementoptrtensorprod3dinterp_init(), scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), scale_meshfield_filter_operation_base::meshfieldfilteroperationbase_apply_reconst1d_y_2(), scale_mesh_hierarchy_base::meshhierarchy_construct_pmg_mat1d(), and scale_element_quadrilateral::quadrilateralelement_init().

Generated by doxygen 1.14.0