    # Isogeometric analysis with piece-wise constant test functions

We focus on finite element method computations. We show that systems of linear equations resulting from Galerkin discretizations with C1 test functions are algebraically equivalent to systems of linear equations resulting from discretizations with piece-wise constant test functions. We show that test functions that fulfill the partition of unity can be eliminated from the system of linear equations. The rows of the system matrix can be combined, and the test functions can be sum up to 1 using the partition of unity property at the quadrature points. The test functions in higher continuity IGA can be set to piece-wise constants without losing any accuracy. It is equivalent to testing with piece-wise constant basis functions, with supports span over some parts of the domain. This observation has the following consequences. The numerical integration cost can be reduced because we do not need to evaluate the test functions since they are equal to 1. This reduction of test functions can be performed for an arbitrary linear differential operator resulting from the Galerkin method applied to a PDE where we use C^1 continuity basis functions. This observation is valid for any basis functions preserving the partition of unity property. It is independent of the problem dimension and geometry of the computational domain. It also can be used in time-dependent problems, e.g., in the explicit dynamics computations, where we can reduce the cost of generation of the right-hand side. In particular, this observation works for any computational problem that can be solved using the Galerkin method with isogeometric analysis with higher continuity test functions. The Galerkin approximation with C1 continuity is equivalent to linear combinations of the collocations at points and with weights resulting from applied quadrature, over the spans defined by supports of test functions.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The main result of this paper can be summarized as follows. We focus on finite element method discretization, with continuity basis functions, e.g. quadratic B-splines utilized in isogeometric analysis (IGA) [1; 2; 3]. Let us focus our attention on the one-dimensional Laplace equation for the simplicity of the presentation. In this case, the Galerkin method involves the integrals (assuming zero boundary condition also for simplicity). If we use higher continuity basis functions, e.g., continuity B-splines, the approximation is an element of a subspace of , and if we integrate exactly, with proper numerical quadrature, these integrals are the same. In other words, the matrix of the system of linear equations resulting from the Galerkin method

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−∫∇Bx1,p∇Bx1,pdx⋯−∫∇Bx1,p∇BxNx,pdx−∫∇Bx2,p∇Bx1,pdx⋯−∫∇Bx2,p∇BxNx,pdx⋮⋮⋮−∫∇BxNx,p∇Bx1,pdx⋯−∫∇BxNx,p∇BxNx,pdx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣∫RHS(x)Bx1(x)dx∫RHS(x)Bx2(x)dx⋮∫RHS(x)BxNx(x)dx⎤⎥ ⎥ ⎥ ⎥ ⎥⎦

with continuity basis functions have identical double precision values as the system not integrated by parts

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫Bx1,pΔBx1,pdx⋯∫Bx1,pΔBxNx,pdx∫Bx2,pΔBx1,pdx⋯∫Bx2,pΔBxNx,pdx⋮⋮⋮∫BxNx,pΔBx1,pdx⋯∫BxNx,pΔBxNx,pdx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣∫RHS(x)Bx1(x)dx∫RHS(x)Bx2(x)dx⋮∫RHS(x)BxNx(x)dx⎤⎥ ⎥ ⎥ ⎥ ⎥⎦

The matrices as well as the right-hand-sides of both system are equal. The fluxes between elements are zero when we employ discretization. It does not matter which method we use for the generation of the system on the computer, and the resulting floating-point values will be the same (up to double precision round-off errors).

The second observation is that the system where we test the Laplace equation with B-splines is algebraically equivalent to the one where we have set the test functions to 1 (we, however, integrate over the spans of the piece-wise constant test functions).

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pΔBx1,pdx⋯∫IBx1,pΔBxNx,pdx∫IBx2,pΔBx1,pdx⋯∫IBx2,pΔBxNx,pdx⋮⋮⋮∫IBxNx,pΔBx1,pdx⋯∫IBxNx,pΔBxNx,pdx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pRHS(x)dx∫IBx2,pRHS(x)dx⋮∫IBxNx,pRHS(x)dx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

where if and if . The detailed derivation of this algebraic equivalence is present later in the paper. It is based on the idea of combining the rows of the matrix. The rows are combined in such a way that test functions sum up to 1, using the partition of unity property.

This observation has the following important consequences. First, the numerical integration cost will be reduced, since we do not need to integrate the test functions (we do not need to evaluate the test B-splines at quadrature points). Second, this equivalence does not depend on the selected quadrature points. Third, we will show that this algebraic equivalence is true if we replace the Laplacian by any partial differential operator resulting from a PDE that can be solved with approximations with basis functions, preserving the partition of unity property. We selected B-splines for the simplicity of the presentation but there are several other options for discretization available [7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19]. The critical here is the partition of unity property. Fourth, this equivalence is also independent of the dimension of the problem, and it works in two or three-dimensions, or in the space-time formulations. Fifth, this equivalence is independent of the geometry of the computational domain and of the Jacobian of the transformation of the patch of elements into the master patch. Recently, it is also possible to extend the continuity between patches of elements , so the equivalence with peice-wise constants also co be extended there. Sixth, we end up with the integrals using the values of the trial functions at the quadrature points.

It is like combining the collocation points [4; 5] at quadrature points with quadrature weights, over the spans of piece-wise constant test functions. The quadrature and the spans of the piece-wise constant test functions define the locations of the collocation points. Several collocation points are combined into one equation by the integration operator.

This observation speeds up also the explicit simulations with IGA, since the integration of the right-hand side is cheaper.

Summing up, the test functions for IGA with basis functions are indeed piece-wise constant. The same logic applies to any basis functions that are globally and preserves the partition of unity property. The method can be easily verified by any finite element code using basis functions, preserving the partition of unity. We just implement in the element routine the discretization of the PDE, and we set the test functions to 1, in the matrix, and the right-hand-side. We, however, keep the spans of test functions in the integrals. We can also reduce the order of the quadrature since test functions are equal to 1. The numerical result will be identical.

## 2 One dimensional case

Let us focus on the general PDE in the following form

 Fu=RHS (1)

defined over interval. We partition the interval into finite elements. Let us use the Galerkin method with continuity of the discretization. We have the one dimensional B-spline basis functions

 {Bxi,p(x)}i=1,...,Nx (2)

where . We approximate the solution . We also test with B-splines.

If we have continuity of the basis functions and we use the exact quadratures during the integration, then the fact, if we integrate by parts or not, does not matter, the values in the matrix are the same, before or after the integration. So let us focus on B-splines and test our PDE with B-splines, and we do not integrate by parts.

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫Bx1,pF(Bx1,p)dx⋯∫Bx1,pF(BxNx,p)dx∫Bx2,pF(Bx1,p)dx⋯∫Bx2,pF(BxNx,p)dx⋮⋮⋮∫BxNx,pF(Bx1,p)dx⋯∫BxNx,pF(BxNx,p)dx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣∫RHS(x)Bx1dx∫RHS(x)Bx2dx⋮∫RHS(x)BxNxdx⎤⎥ ⎥ ⎥ ⎥ ⎥⎦

Let us select any quadrature with points and weights , resulting in the exact numerical integration. At a given quadrature point we have non-zero B-spline functions.

We take our system of linear equations, and we add the sum of rows to the first row. We also add the sum of rows to the second row. Similarly, we add the sum of rows to the row . Finally, we add the sum of rows to the row . We get the algebraically equivalent system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑m=1,p∫Bxm,pF(Bx1,p)dx⋯∑m=1,...,p∫Bxm,pF(BxNx,p)dx∑m=1,p+1∫Bxm,pF(Bx1,p)dx⋯∑m=1,p+1∫Bxm,pF(BxNx,p)dx⋮⋮⋮∑m=Nx−p,...,Nx∫Bxm,pF(Bx1,p)dx⋯∑m=Nx−p,...,Nx∫Bxm,pF(BxNx,p)dx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦ =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑m=1,...,p∫RHS(x)Bxm,pdx∑m=1,...,p+1∫RHS(x)Bxm,pdx⋮∑m=Nx−p,Nx∫RHS(x)Bxm,pdx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

Now, we integrate numerically using our exact quadrature points and weights, so we can move the summation inside

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owo[∑m=1,pBxm,p(xo)]F(Bx1,p(xo))Jac(xo)⋯∑owo[∑m=1,pBxm,p(xo)]F(BxNx,p(xo))Jac(xo)∑owo[∑m=1,p+1Bxm,p(xo)]F(Bx1,p(xo))Jac(xo)⋯∑owo[∑m=1,p+1Bxm,p(xo)]F(BxNx,p(xo))Jac(xo)⋮⋮⋮∑owo[∑m=Nx−p,...,NxBxm,p]F(Bx1,p(xo))Jac(xo)⋯∑owo[∑m=Nx−p,...,NxBxm,p(xo)]F(BxNx,p(xo))Jac(xo)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owoRHS(xo)[∑m=1,...,pBxm,p(xo)]Jac(xo)∑owoRHS(xo)[∑m=1,...,p+1Bxm,p(xo)]Jac(xo)⋮∑owoRHS(xo)[∑m=Nx−p,NxBxm,p(xo)]Jac(xo)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

At a given quadrature point, the B-splines sum up to 1, from the partition of unity property

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owoF(Bx1,p(xo))Jac(xo)⋯∑owoF(BxNx,p(xo))Jac(xo)∑owoF(Bx1,p(xo))Jac(xo)⋯∑owoF(BxNx,p(xo))Jac(xo)⋮⋮⋮∑owoF(Bx1,p(xo))Jac(xo)⋯∑owoF(BxNx,p(xo))Jac(xo)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣∑owoRHS(xo)Jac(xo)∑owoRHS(xo)Jac(xo)⋮∑owoRHS(xo)Jac(xo)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦

We can come back to the integrals, that are computed over the domains of the piece-wise constant test functions

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pF(Bx1,p)dx⋯∫IBx1,pF(BxNx,p)dx∫IBx2,pF(Bx1,p)dx⋯∫IBx2,pF(BxNx,p)dx⋮⋮⋮∫IBxNx,pF(Bx1,p)dx⋯∫IBxNx,pF(BxNx,p)dx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣u1u2⋮uNx⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pRHSdx∫IBx2,pRHSdx⋮∫IBxNx,pRHSdx⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

The domains of the original B-spline basis functions defines now the range of the collocation points computed at quadrate points and weights that span over the domains, see Figure 1. There, we show that

 ∫IBx5,2(x)Bx4,2(x)dx=∑iwiIBx5,2(xi)Bx4,2(xi)= w3Bx4,2(x3)+w4Bx4,2(x4)+w5Bx4,2(x5)+w6Bx4,2(x6) Figure 1: The domains of the original B-spline test functions define the range of the collocation points. The collocation points and weights are selected by the quadrate span over the domain.

## 3 Two dimensional case

We repeat our considerations in the two-dimensional case. We start from the general form of the PDE

 Fu=RHS (3)

where we discretize with continuity B-splines, and we do not integrate by parts (or we start from the Galerkin formulation to show the well-posedness, but then we discretize with higher continuity functions, and we integrate back by parts). We have the global system of linear equations

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫(Bx1,pBy1,p)F(Bx1,pBy1,p)∫(Bx1,pBy1,p)F(Bx2,pBy1,p)⋯∫(Bx1,pBy1,p)F(BxNx,pByNy,p)∫(Bx2,pBy1,p)F(Bx1,pBy1,p)∫(Bx2,pBy1,p)F(Bx2,pBy1,p)⋯∫(Bx2,pBy1,p)F(BxNx,pByNy,p)⋮⋮⋮⋮∫(BxNx,pByNy,p)F(Bx1,pBy1,p)∫(BxNx,pByNy,p)F(Bx2,pBy1,p)⋯∫(BxNx,pByNy,p)F(BxNx,pByNy,p)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢⎣u1,1u2,1⋮uNx,Ny⎤⎥ ⎥ ⎥ ⎥ ⎥⎦
 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫RHS(x,y)Bx1,pBy1,p∫RHS(x,y)Bx2,pBy1,p⋮∫RHS(x,y)BxNx,pByNy,p⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

We consider a quadrature with points and weights . At a given point from the selected quadrature, we have non-zero B-spline basis functions in one direction.

Since each row in the global matrix corresponds to one test function we can number them .

Now, we will replace the global matrix by a new one obtained by adding linear combinations of rows. In other words, we replace each row of the element matrix by a linear combination of other rows, in such a way that the test B-splines in one direction will sum up later to 1.

Namely, we add to the row the sum of rows number

 (min(1,i−p),j;k,l),...,(max(Nx,i+p),j;k,l) (4)

We get the equivalent global system

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑m=1,...,p+1∫(Bxm,pBy1,p)F(Bx1,pBy1,p)⋯∑m=1,...,p+1∫(Bxm,pBy1,p)F(BxNx,pByNy,p)⋮⋮⋮∑m=i−p,...,i+p∫(Bxm,pByj,p)F(Bx1,pBy1,p)⋯∑m=i−p,...,i+p∫(Bxm,pByj,p)F(BxNx,pByNy,p)⋮⋮⋮∑m=Nx−p,...,Nx∫(Bxm,pByNy,p)F(Bx1,pBy1,p)⋯∑m=Nx−p,...,Nx∫(Bxm,pByNy,p)F(BxNx,pByNy,p)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1,1⋮ui,j⋮uNx,Ny⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑m=1,...,p+1∫RHS(x,y)Bxm,p(x)By1,p(y)⋮∑m=i−p,...,i+p∫RHS(x,y)Bxm,p(x)Byj,p(y)⋮∑m=Nx−p,...,Nx∫RHS(x,y)Bxm,p(x)ByNy,p(y)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

Now, we compute the integrals exactly by using numerical quadrature rule

 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owoRHS(xo,yo)[∑m=1,...,p+1Bxm,p(xo)]By1,p(yo)Jac(xo,yo)⋮∑owoRHS(xo,yo)[∑m=i−p,...,i+pBxm,p(xo)]Byj,p(yo)Jac(xo,yo)⋮∑owoRHS(xo,yo)[∑m=Nx−p,...,NxBxm,p(xo)]ByNy,p(yo)Jac(xo,yo)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

At a given quadrature point, we sum up all the B-splines in one direction, and from the partition of unity property, they sum up to 1. So these summation terms they disappear

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owo(By1,p(yo))F(Bx1,p(xo)By1,p(yo))Jac(xo,yo)⋯⋮⋮∑owo(Byj,p(yo))F(Bx1,p(xo)By1,p(yo))Jac(xo,yo)⋯⋮⋮∑owo(ByNy,p(yo))F(Bx1,p(xo)By1,p(yo))Jac(xo,yo)⋯⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1,1⋮ui,j⋮uNx,Ny⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑owoRHS(xo,yo)By1,p(yo)Jac(xo,yo)⋮∑owoRHS(xo,yo)Byj,p(yo)Jac(xo,yo)⋮∑owoRHS(xo,yo)ByNy,p(yo)Jac(xo,yo)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

Now, we can come back to the integral, and we compute them over the domains of the piece-wise constant test functions

 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pRHS(x,y)By1,p(y)⋮∫IBxi,pRHS(x,y)Byj,p(y)⋮∫IBxNx,pRHS(x,y)ByNy,p(y)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

Now, we repeat the same logic with respect to the one-dimensional B-spline basis functions in the direction, to get

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pIBy1,pF(Bx1,pBy1,p)⋯∫IBx1,pIBy1,pF(BxNx,pByNy,p)⋮⋮⋮∫IBxi,pIByj,pF(Bx1,pBy1,p)⋯∫IBxi,pIByj,pF(BxNx,pByNy,p)⋮⋮⋮∫IBxNx,pIByNy,pF(Bx1,pBy1,p)⋯∫IBxNx,pIByNy,pF(BxNx,pByNy,p)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1,1⋮ui,j⋮uNx,Ny⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫IBx1,pIBy1,pRHS(x,y)⋮∫IBxi,pIByj,pRHS(x,y)⋮∫IBxNx,pIByNy,pRHS(x,y)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

## 4 Examples

### 4.1 Isogeometric L2 projections

First example is the L2 projection problem.

 u=RHS

which in the weak form is

 (v,u)=(v,RHS)

solved over .

We define the B-spline basis for trial and test and we discretize in the standard Galerkin way

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫(Bx1,pBy1,pBz1,p)(Bx1,pBy1,pBz1,p)⋯∫(Bx1,pBy1,pBz1,p)(BxNx,pByNy,pBzNz,p)∫(Bx2,pBy1,pBz1,p)(Bx1,pBy1,pBz1,p)⋯∫(Bx2,pBy1,pBz1,p)(BxNx,pByNy,pBzNz,p)⋮⋮⋮∫(BxNx,pByNy,pBzNz,p)(Bx1,pBy1,pBz1,p)⋯∫(BxNx,pByNy,pBzNz,p)(BxNx,pByNy,pBzNz,p)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢⎣u1,1,1u2,1,1⋮uNx,Ny,Nz⎤⎥ ⎥ ⎥ ⎥ ⎥⎦
 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∫RHS(x,y,z)Bx1,p(x)By1,p(y)Bz1,p(z)∫RHS(x,y,z)