Default quadrature order for simplex elements

[dladd]

We've been running into some odd behaviour using the default quadrature scheme for 2D simplex elements (triangles). I've created an edited the diffusion_example case to demonstrate the issue. This now has boundary conditions set to 0 in the lower left and 1 in the upper right. The figures below are at the end of the solution (this goes from t=0 to 1 with dt=0.001)

2D Linear Lagrange, default basis settings

Here we see a relatively smooth solution and the linear solver converges within relative tolerance (RTOL) in 7 iterations in the last several timesteps

Refining the mesh from a 10x10 to a 50x50 grid improves the solution and the linear solver now converges within 6 iterations.

2D Simplex, default basis settings: quadrature order = 3

The solution looks pretty rough compared to the Lagrange case and the linear solver now takes 29 iterations to converge.

Refining the mesh from 10x10 to 50x50 definitely improves the solution but the solver is now taking a whopping 172 iterations.

2D Simplex, manually specify quadrature order = 2

Solution at 10x10 looks similar to the Lagrange case and the solver now converges in 5 iterations

Refining to 50x50 again improves the solution and is converging in 12 iterations.

From this I would think the quadrature order should default to 2 rather than 3 for the 2D simplex case. Can anyone comment on how these were determined or when we should be tweaking them? I know quadrature is trickier to define for the Simplex case.

Thanks!

[dladd]

I've looked into this a little bit more and from Basis_GaussPointsCalculate, it appears that Gauss_Simplex is using the formulae from:

Liu, Yen and Vinokur, Marcel. "Exact Integrations of Polynomials and Symmetric Quadrature Formulas over Arbitrary Polyhedral Grids", Journal of Computational Physics, 140:122-147 (1998).

to calculate area coordinate expressions of Gauss point locations and weights. Then in Basis_AreaToXiCoordinates, these get converted to xi coordinates.

QO = 1

For triangles with a quadrature order (QO) = 1, there will be one Gauss point and Gauss_Simplex will place it at (in area coordinates):

a_c_1 = (1/3 , 1/3, 1/3)

converting with Basis_AreaToXiCoordinates appears to throw out the last area coordinate and define xi coordinates as (1 - area coordinates):

xi_c_1 = (2/3, 2/3)

QO = 2

Now we expect 3 Gauss points, which Gauss_Simplex places:

a_c_1 = (0, 1/2, 1/2) a_c_2 = (1/2, 0, 1/2) a_c_3 = (1/2, 1/2, 0)

and Basis_AreaToXiCoordinates converts to:

xi_c_1 = (1, 1/2) xi_c_2 = (1/2, 1) xi_c_3 = (1/2, 1/2)

QO = 3

Now we expect 4 Gauss points, which Gauss_Simplex places:

a_c_1 = (1/3, 1/3, 1/3) a_c_2 = (17, -8, -8) a_c_3 = (-8, 17, -8) a_c_4 = (-8, -8, 17)

and Basis_AreaToXiCoordinates converts to:

xi_c_1 = (2/3, 2/3) xi_c_2 = (-16, 9) xi_c_3 = (9, -16) xi_c_4 = (9, 9)

Which seems pretty odd to me... Maybe these aren't actually xi_coordinates and they get turned into something more reasonable by scaling with the weighting terms eventually? Will keep investigating.

[chrispbradley]

No, they do not seem correct. I recently went over those Gauss points and weights and though that I had corrected those bugs. However, I've just gone back to the paper for triangles with order 3 alpha1 is 25.0??? This can't be correct as all the coordinates should be within the triangle???

[chrispbradley]

It is a little hard to tell if they are talking about a Lobotto scheme or not but Equation (38) of the paper defines alpha as 2/(n+2) which would give 2/5 rather than 25? It also gives the correct w_alpha of (n+2)^2/4n(n+1) or 25/48. However it gives w_c as -n^2/n(n+1) which would be -9/12 and not -9/16 from table 3? Typos in the paper? Could try changing the alpha and w_c and see if it works?

[dladd]

Thanks @chrispbradley- that does sound worth a try. Does the area-to-xi coordinate conversion throwing out the last component seem right to you?

Also, is the xi orientation for simplices similar to this figure in the documentation? If so I would think even xi_c = (2/3, 2/3) or (1, 1/2) could be placed outside of a triangle.

[chrispbradley]

Hi @dladd, yes, area coordinates have a redundancy (they sum to a certain value). The Gauss points are stored as xi points (so that the number of xi's are consistent between quads and tets) and thus the 4 area coordinates are converted to 3 xi coordinates. The figure you link to is orientation that cmgui uses. We use the Zienkiewicz orientation e.g., https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/lineartriangle.svg https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/lineartetrahedron.svg https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/quadratictriangle.svg https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/quadratictetrahedron.svg https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/cubictriangle.svg https://github.com/OpenCMISS/documentation/blob/develop/notes/OpenCMISSNotes/svgs/BasisFunctions/cubictetrahedron.svg

[dladd]

Thanks @chrispbradley - that's fair enough re: the area-to-xi. On the xi orientation side what I am worried about is something like:

[dladd]

Sorry that shouldn't actually happen- thinking about simplicies is messing with my head...

[dladd]

Bugfix pull #158 adjusts alpha and w as @chrispbradley suggested above. With those changes, the "2D Simplex, default basis settings: QO = 3" case in my example now converges in 5 iterations (previously 29) at a 10x10 element resolution and 12 iterations (previously 172) at 50x50.

I have also now gone through and sanity checked the computed xi coordinate locations for simplicies (lines, triangles, and tets) for orders 2 through 5 (spreadsheet here for the curious). Looks like order=3 triangles was the only obvious issue- the rest are at least placed with xi coords between 0 and 1 with no repeats.