diff --git a/README.md b/README.md index c8d3088..0aacaee 100644 --- a/README.md +++ b/README.md @@ -152,6 +152,39 @@ K_2D = finite_element_matrix(l_i,l_j,coordinate_x1, ``` where we have two coordinates, `x1` (with polynomials `l_i` and `l_j`) and `x2` (with polynomials `l_k` and `l_n`). Note that the `@enum` variables used to select the Lagrange polynomials are the same between the two coordinates. The matrix created above is indexed in the order `K_2D[i,k,j,n]`, indexing over the first polynomial `l_i` of `x1`, then the first polynomial `l_k` of `x2`, then the second polynomial `l_j` of `x1` with the final index representing the polynomial `l_n` of `x2`. The same adaptive quadrature keyword arguments are used as in the 1D case. + +# Integration limits + +For some function inputs $`\rho(x)`$, the domain over which $`\rho(x)`$ is nonzero may not coincide with the element limits $`[x_l, x_u]`$. To handle these cases efficiently with Gauss-Legendre quadrature points, we provide a 1D quadrature option `GLSpecifiedLimits(x_min, x_max)` which can limit the range of integration where necessary without introducing the need to define $`\rho(x)`$ where it is zero, e.g., with a Heaviside function. Note that the values `x_min` and `x_max` should be supplied in the physical coordinate normalisation for the coordinate $`x`$=`coordinate`, rather than in the reference normalisation on the reference domain $`[-1,1]`$ defined for each element. For example, for a $`\rho(x)`$ which is nonzero between $`[x_{-},x_{+}]`$ = `[x_min, x_max]`, we compute the integral +```math +A_{ij} = \int^{\textrm{min}(x_u,x_{+})}_{\textrm{max}(x_l,x_-)} P_i(x) Q_j(x) \rho(x) d x. +``` +For the integral +```math +M_{ij} = \int^{\textrm{min}(x_u,x_{+})}_{\textrm{max}(x_l,x_-)} \Phi_i(x) \Phi_j(x) \rho(x) d x, +``` +the corresponding syntax for the matrix element $`M_{ij}`$ is as follows. +``` +finite_element_matrix(lagrange_x,lagrange_x,coordinate, + kernel_function=(x -> rho(x)), + quadrature_option=GLSpecifiedLimits(x_min,x_max)) +``` +For matrices of integrals in two dimensions, we use a similar syntax, as follows. +``` +finite_element_matrix(lagrange_x,lagrange_x,x1_coordinate, + lagrange_x,lagrange_x,x2_coordinate, + kernel_function=((x1,x2)->rho(x1,x2)), + quadrature_option_x1=GLSpecifiedLimits(x1_min,x1_max), + quadrature_option_x2=GLSpecifiedLimits(x2_min,x2_max)) +``` +The default quadrature options for integrals in the range $`[x_{l},x_{u}]`$ are as follows. +``` +quadrature_option=Default1DQuadrature() +quadrature_option_x1=Default1DQuadrature() +quadrature_option_x2=Default1DQuadrature() +``` + + # Examples There are several examples of taking first and second derivatives diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl index 9d4a9fc..7c23f4a 100644 --- a/src/FiniteElementMatrices.jl +++ b/src/FiniteElementMatrices.jl @@ -11,13 +11,25 @@ using FastGaussQuadrature: gausslegendre export lagrange_x, d_lagrange_dx, finite_element_matrix, - ElementCoordinates + ElementCoordinates, + GLSpecifiedLimits @enum LagrangeFunctionType begin lagrange_x d_lagrange_dx end +abstract type Abstract1DQuadrature end +struct Default1DQuadrature <: Abstract1DQuadrature end +struct GLSpecifiedLimits <: Abstract1DQuadrature + # minimum of range of integration + # in physical grid v = s x + c + # with x in [-1, 1] + v_min::Float64 + # maximum of range of integration + v_max::Float64 +end + struct ElementCoordinates # precomputed data for calculating # the Lagrange polynomials, including @@ -132,8 +144,9 @@ function _finite_element_matrix( kernel_function::TFunction=((v -> 1.0)), additional_quadrature_points::Int64=0, # argument to permit single implementation of - adaptive_quadrature_points::Int64=0 - ) where TFunction + adaptive_quadrature_points::Int64=0, + quadrature_option::TQuad=Default1DQuadrature() + ) where {TFunction, TQuad <: Abstract1DQuadrature} lpoly_data = coordinate.lpoly_data ngrid = length(coordinate.lpoly_data.x_nodes) scale = coordinate.scale @@ -146,7 +159,7 @@ function _finite_element_matrix( # nquad chosen for exact results for default inputs # with kernel = 1.0 and zero additional quadrature points nquad = ngrid + additional_quadrature_points + adaptive_quadrature_points - zz, wz = gausslegendre(nquad) + zz, wz = quadrature1D(quadrature_option,nquad,coordinate) # compute integral # int P_i(z) Q_j(z) poly(z) s d z # with poly(z) = (s z + c)^power @@ -193,8 +206,9 @@ function _finite_element_matrix( # rather than of the reference coordinate z on [-1,1] kernel_function::TFunction=((v -> 1.0)), additional_quadrature_points::Int64=0, - adaptive_quadrature_points::Int64=0 - ) where TFunction + adaptive_quadrature_points::Int64=0, + quadrature_option::TQuad=Default1DQuadrature(), + ) where {TFunction, TQuad <: Abstract1DQuadrature} lpoly_data = coordinate.lpoly_data ngrid = length(coordinate.lpoly_data.x_nodes) scale = coordinate.scale @@ -208,7 +222,7 @@ function _finite_element_matrix( # nquad chosen for exact results for default inputs # with kernel = 1.0 and zero additional quadrature points nquad = 2*ngrid + additional_quadrature_points + adaptive_quadrature_points - zz, wz = gausslegendre(nquad) + zz, wz = quadrature1D(quadrature_option,nquad,coordinate) # compute integral # int P_i(z) Q_j(z) S_k(z) poly(z) d z # with poly(z) = (s z + c)^power @@ -270,8 +284,12 @@ function _finite_element_matrix( kernel_function::TFunction=((v1,v2) -> 1.0), additional_quadrature_points_x1::Int64=0, additional_quadrature_points_x2::Int64=0, - adaptive_quadrature_points::Int64=0 - ) where TFunction + adaptive_quadrature_points::Int64=0, + quadrature_option_x1::TQuad1=Default1DQuadrature(), + quadrature_option_x2::TQuad2=Default1DQuadrature(), + ) where {TFunction, + TQuad1 <: Abstract1DQuadrature, + TQuad2 <: Abstract1DQuadrature} # coordinate x1 data lpoly_data_x1 = coordinate_x1.lpoly_data ngrid_x1 = length(coordinate_x1.lpoly_data.x_nodes) @@ -293,9 +311,9 @@ function _finite_element_matrix( # nquad chosen for exact results for default inputs # with kernel = 1.0 and zero additional quadrature points nquad_x1 = ngrid_x1 + additional_quadrature_points_x1 + adaptive_quadrature_points - zz_x1, wz_x1 = gausslegendre(nquad_x1) + zz_x1, wz_x1 = quadrature1D(quadrature_option_x1,nquad_x1,coordinate_x1) nquad_x2 = ngrid_x2 + additional_quadrature_points_x2 + adaptive_quadrature_points - zz_x2, wz_x2 = gausslegendre(nquad_x2) + zz_x2, wz_x2 = quadrature1D(quadrature_option_x2,nquad_x2,coordinate_x2) # compute integral # \int \int \left(P1_i(z_1) Q1_j(z_1) P2_i(z_1) Q2_j(z_2) # kernel(s_1 z_1 + c_1, s_2 z_2 + c_2) \right) s_1 s_2 d z_1 d z_2 @@ -333,5 +351,44 @@ function _finite_element_matrix( return matrix end +function quadrature1D(::Default1DQuadrature, nquad::Int64, coordinate::ElementCoordinates) + # default quadrature running from [-1,1] + zz, wz = gausslegendre(nquad) + return zz, wz +end +function quadrature1D(qopt::GLSpecifiedLimits, nquad::Int64, coordinate::ElementCoordinates) + # limits specified for integration + v_min = qopt.v_min + v_max = qopt.v_max + # limits based on the range supported by the Lagrange polynomials + v_lower_limit = coordinate.shift - coordinate.scale + v_upper_limit = coordinate.shift + coordinate.scale + # some checks on the inputs + if !(v_min < v_upper_limit) + error("invalid integration range: v_min >= v_upper_limit") + elseif !(v_max > v_lower_limit) + error("invalid integration range: v_max =< v_lower_limit") + elseif !(v_min < v_max) + error("invalid integration range: v_min >= v_max") + end + # obtain valid limits of integration + # in the physical range [v_lower_limit,v_upper_limit] + v_min = max(v_min, v_lower_limit) + v_max = min(v_max, v_upper_limit) + # in the reference range [-1,1] + x_min = (v_min - coordinate.shift)/coordinate.scale + x_max = (v_max - coordinate.shift)/coordinate.scale + # default quadrature running from [-1,1] + zz, wz = gausslegendre(nquad) + # quadrature running from [x_min, x_max] + zx = zeros(Float64, nquad) + wx = zeros(Float64, nquad) + x_scale = 0.5*(x_max - x_min) + x_shift = 0.5*(x_max + x_min) + @. zx = x_scale*zz + x_shift + @. wx = x_scale*wz + return zx, wx +end + end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 9299d0a..3f3aec3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,8 @@ using Test: @testset, @test using FiniteElementMatrices: lagrange_x, d_lagrange_dx, finite_element_matrix, - ElementCoordinates + ElementCoordinates, + GLSpecifiedLimits using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre using LinearAlgebra: lu, ldiv!, mul! @@ -52,9 +53,9 @@ function test_first_derivative(;nodes::node_type=GLL, # check S_ij = P_ji to cover tests where # derivative acts on test function @test isapprox(S,transpose(P), atol=2.0e-15) - result = Array{Float64,1}(undef,ngrid) - dummy = Array{Float64,1}(undef,ngrid) - err = Array{Float64,1}(undef,ngrid) + result = Array{Float64,1}(undef,ngrid) + dummy = Array{Float64,1}(undef,ngrid) + err = Array{Float64,1}(undef,ngrid) # function to test f = Array{Float64,1}(undef,ngrid) df = Array{Float64,1}(undef,ngrid) @@ -96,9 +97,9 @@ function test_first_derivative_nonlinear_operators( # check Y100_ijk = Y001_kji @test isapprox(permutedims(Y100, [3,2,1]),Y001,atol=2.0e-15) # check that a first derivative can be carried out with Y100 - result = Array{Float64,1}(undef,ngrid) + result = Array{Float64,1}(undef,ngrid) dummy = Array{Float64,1}(undef,ngrid) - err = Array{Float64,1}(undef,ngrid) + err = Array{Float64,1}(undef,ngrid) # function to test f = Array{Float64,1}(undef,ngrid) df = Array{Float64,1}(undef,ngrid) @@ -135,14 +136,14 @@ function test_second_derivative(;nodes::node_type=GLL, M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate) S = -finite_element_matrix(d_lagrange_dx,lagrange_x,0,coordinate) K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,coordinate) - result = Array{Float64,1}(undef,ngrid) - dummy = Array{Float64,1}(undef,ngrid) - err = Array{Float64,1}(undef,ngrid) + result = Array{Float64,1}(undef,ngrid) + dummy = Array{Float64,1}(undef,ngrid) + err = Array{Float64,1}(undef,ngrid) # function to test f = Array{Float64,1}(undef,ngrid) df = Array{Float64,1}(undef,ngrid) d2f = Array{Float64,1}(undef,ngrid) - # to test the weak first derivative, + # to test the weak first derivative, # choose a fn that vanishes at x = +-1 # to avoid including boundary terms for i in 1:ngrid @@ -150,19 +151,19 @@ function test_second_derivative(;nodes::node_type=GLL, f[i] = sin(pi*(y-shift)/scale) df[i] = (pi/scale)*cos(pi*(y-shift)/scale) end - + # test the performance of a first derivative - # taken using the "weak" methods via + # taken using the "weak" methods via # integration by parts luM = lu(M) mul!(dummy,S,f) ldiv!(result,luM,dummy) - + @. err = result - df maxerr = maximum(abs.(err)) @test maxerr < atol - - # to test the weak second derivative, + + # to test the weak second derivative, # choose a fn that where the derivative vanishes at x = +-1 # to avoid including boundary terms for i in 1:ngrid @@ -171,14 +172,14 @@ function test_second_derivative(;nodes::node_type=GLL, d2f[i] = -((pi/scale)^2)*cos(pi*(y-shift)/scale) end # test the performance of a second derivative - # taken using the "weak" methods via + # taken using the "weak" methods via # integration by parts mul!(dummy,K,f) ldiv!(result,luM,dummy) @. err = result - d2f maxerr = maximum(abs.(err)) @test maxerr < atol - + end function test_nonlinear_operators(;nodes::node_type=GLL, @@ -195,9 +196,9 @@ function test_nonlinear_operators(;nodes::node_type=GLL, Y100 = finite_element_matrix(d_lagrange_dx,lagrange_x,lagrange_x,0,coordinate) Y010 = finite_element_matrix(lagrange_x,d_lagrange_dx,lagrange_x,0,coordinate) Y001 = finite_element_matrix(lagrange_x,lagrange_x,d_lagrange_dx,0,coordinate) - result = Array{Float64,1}(undef,ngrid) + result = Array{Float64,1}(undef,ngrid) dummy = Array{Float64,1}(undef,ngrid) - err = Array{Float64,1}(undef,ngrid) + err = Array{Float64,1}(undef,ngrid) # function to test f = Array{Float64,1}(undef,ngrid) u = Array{Float64,1}(undef,ngrid) @@ -330,6 +331,62 @@ function test_nonpolynomial_kernel(;nodes::node_type=GLL, #println(maxerr) end +function test_limited_integration_range(;nodes::node_type=GLL, + ngrid::Int64=5, + y0::Float64=-1.1, + y1::Float64=1.3, + atol::Float64=2.0e-13) + x = reference_nodes(nodes,ngrid) + scale = 0.5*(y1-y0) + shift = 0.5*(y0+y1) + coordinate = ElementCoordinates(x,scale,shift) + # midpoint splitting the integral + ymid = y0 + (1.0/3.0)*(y1 - y0) + # test of operators with two polynomials + M_total = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate) + M_lower_range = finite_element_matrix(lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(y0,ymid)) + M_upper_range = finite_element_matrix(lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(ymid,y1)) + M_err = deepcopy(M_total) + # test the total matches the exact result + @. M_err = M_total - M_lower_range - M_upper_range + maxerr = maximum(abs.(M_err)) + #println(maxerr) + @test maxerr < atol + + # test of operators with three polynomials + M_nl_total = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,0,coordinate) + M_nl_lower_range = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(y0,ymid)) + M_nl_upper_range = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(ymid,y1)) + M_nl_err = deepcopy(M_nl_total) + # test the total matches the exact result + @. M_nl_err = M_nl_total - M_nl_lower_range - M_nl_upper_range + maxerr_nl = maximum(abs.(M_nl_err)) + #println(maxerr_nl) + @test maxerr_nl < atol + + # test of 2D operators with two polynomials + M2D_total = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate, + lagrange_x,lagrange_x,0,coordinate) + M2D_range1 = finite_element_matrix(lagrange_x,lagrange_x,coordinate, + lagrange_x,lagrange_x,coordinate, + quadrature_option_x1=GLSpecifiedLimits(y0,ymid), quadrature_option_x2=GLSpecifiedLimits(y0,ymid)) + M2D_range2 = finite_element_matrix(lagrange_x,lagrange_x,coordinate, + lagrange_x,lagrange_x,coordinate, + quadrature_option_x1=GLSpecifiedLimits(ymid,y1), quadrature_option_x2=GLSpecifiedLimits(y0,ymid)) + M2D_range3 = finite_element_matrix(lagrange_x,lagrange_x,coordinate, + lagrange_x,lagrange_x,coordinate, + quadrature_option_x1=GLSpecifiedLimits(y0,ymid), quadrature_option_x2=GLSpecifiedLimits(ymid,y1)) + M2D_range4 = finite_element_matrix(lagrange_x,lagrange_x,coordinate, + lagrange_x,lagrange_x,coordinate, + quadrature_option_x1=GLSpecifiedLimits(ymid,y1), quadrature_option_x2=GLSpecifiedLimits(ymid,y1)) + M2D_err = deepcopy(M2D_total) + # test the total matches the exact result + @. M2D_err = M2D_total - M2D_range1 - M2D_range2 - M2D_range3 - M2D_range4 + maxerr_2D = maximum(abs.(M2D_err)) + #println(maxerr_2D) + @test maxerr_2D < atol +end + function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL, ngrid::Int64=20, func1::Function=(x -> sin(x)), @@ -514,7 +571,7 @@ function runtests() y0 = -0.5, y1=1.0, ngrid=n) end - + println("test nonlinear operators:") for nodes in nodes_list n = 25 @@ -541,6 +598,12 @@ function runtests() test_2D_linear_operators(;nodes=nodes, ngrid_x1=4, ngrid_x2=5, atol=1.0e-13) end + + println("test limited integration range:") + for nodes in nodes_list + println(" -- test: $(nodes) ngrid=5 poly") + test_limited_integration_range(;nodes = nodes, ngrid=5) + end end end