Ordering of polynomials in FE_DGQLegendre<3>

[vachanpo...@gmail.com]

Dear All,

I wasn't able to figure out the ordering of polynomials stored in FE_DGQLegendre<3>. If I am correct, this class is constructed by FE_Poly's constructor using TensorProductPolynomials of Polynomials::Legendre. If this is indeed so, then I wasn't able to find the ordering mentioned in TensorProductPolynomials as well.

So, like always, is the ordering lexicographic (since these are tensor product basis functions)? Or is it based on the order of basis function (say, lowest to highest)?

Thanking in anticipation,

Vachan

[Wolfgang Bangerth]

Vachan,
I don't know, but here is the way to find out assuming that you are interested
in the local enumeration (i.e., the order of shape functions within a single
cell): You create such a finite element and then output the values of the i'th
shape function via
FiniteElement::shape_value(i,p)
on a regular grid of points on the reference cell [0,1]^d.

Alternatively: For DG elements, the local enumeration equals the global
enumeration. You can generate a mesh with one cell, associated a DoFHandler
with it, and then use DataOut::build_patches() with a relatively large number
of subdivisions to output the finite element field associated with a solution
vector that has only one entry (i.e., a unit vector). This is how the pictures
of the shape functions that are shown in several of the finite element
documentation pages are generated.

Of course, if you learn how the shape functions are numbered, it would be most
excellent if you could write a short patch that adds this piece of information
to the documentation!

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[vachanpo...@gmail.com]

Wolfgang,

Thanks a lot for your reply. What I actually need is a change of basis from Lagrange polynomials (nodal) to Legendre polynomials (modal). I then want to know the coefficients of certain modes.

So, if there is any straightforward way to do this in dealii, I would proceed with that. I, could not find any such functions and hence planned on doing it manually, using the shape functions. However, even using any inbuilt functions would only address my issue partly if the ordering is not clear.

One workaround I think is to use 1d basis functions. Since 3d basis functions are tensor products of these, and the ordering of shape functions is fixed for FE_DGQ, I can use 1d Lagrange polynomials to get the required modes. The 1d Lagrange polynomials have a definite order, as can be seen from Legendre::generate_complete_basis().

I will also try your suggestion to figure out the ordering for a 3d element once I get some time.

Thanks again!

[vachanpo...@gmail.com]

I am extremely sorry, the ordering is actually mentioned in the documentation of TensorProductPolynomials. It is lexicographic. I didn't look at this carefully. I suppose the ordering would then be the same for FE_DGQLegendre.

Apologies for my oversight.

[Wolfgang Bangerth]

On 5/24/21 11:52 PM, vachanpo...@gmail.com wrote:
> I am extremely sorry, the ordering is actually mentioned in the documentation
> of TensorProductPolynomials

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassTensorProductPolynomials.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C7305b7f6099345ea6abd08d91f415786%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637575188746128694%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=d4t5afMecx2mhH%2FwxLFUGEEAHltSOgvy2mFOtcZoLO4%3D&reserved=0>.

> It is lexicographic. I didn't look at this carefully. I suppose the ordering
> would then be the same for FE_DGQLegendre.

It doesn't have to be. For example, for FE_Q, we also build on
TensorProductPolynomials but the ordering is not lexicographic. So it would
still be of interest to document the order of shape functions if you end up
finding out what it is!

Best
W.

[Wolfgang Bangerth]

Vachan

> Thanks a lot for your reply. What I actually need is a change of basis from
> Lagrange polynomials (nodal) to Legendre polynomials (modal). I then want to
> know the coefficients of certain modes.
>
> So, if there is any straightforward way to do this in dealii, I would proceed
> with that. I, could not find any such functions and hence planned on doing it
> manually, using the shape functions. However, even using any inbuilt functions
> would only address my issue partly if the ordering is not clear.

I don't recall whether there is an easy way to achieve what you are looking
for, but some elements definitely do something like what you are trying to
achieve. For example, the difference between FE_Q and FE_DGQ is, in essence,
just a change of basis where the basis functions are permuted. Similarly, the
difference between FE_Q and FE_QHierarchical is similar to the nodal -> modal
change you are interested in. Finally, there is also the case of the FE_Q
constructor that receives a Quadrature object as argument and that then
computes a basis change.

You might want to look into how all of these are implemented. Most of this
kind of functionality exists in some kind of helper function that might be
useful to you.

Best
W.

[vachan potluri]

Wolfgang,

 

It doesn't have to be. For example, for FE_Q, we also build on
TensorProductPolynomials but the ordering is not lexicographic. So it would
still be of interest to document the order of shape functions if you end up
finding out what it is!

Noted. So I have verified this with the following code.

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/base/point.h>

#include <deal.II/fe/fe_dgq.h>

#include <iostream>

#include <cmath>

#include <vector>

using namespace dealii;

int main(int argc, char **argv)

{

const int degree = 4;

const int n_poly1 = degree+1;

const FE_DGQLegendre<1> fe1(degree);

const FE_DGQLegendre<3> fe3(degree);

const QGauss<3> quad(degree+1);

const int n_qp = quad.size();

const std::vector<Point<3>>& points = quad.get_points();

const std::vector<double>& weights = quad.get_weights();

for(int k=0; k<n_poly1; k++){

for(int j=0; j<n_poly1; j++){

for(int i=0; i<n_poly1; i++){

int index_3d = i + j*n_poly1 + k*n_poly1*n_poly1;

// compute error between 'index_3d'-th polynomial of 'fe3' and the product of

// i-th, j-th and k-th polynomials of 'fe1'

double error = 0;

for(int q=0; q<n_qp; q++){

error += weights[q]*std::pow(

fe3.shape_value(index_3d, points[q]) - (

fe1.shape_value(i, Point<1>(points[q][0]))*

fe1.shape_value(j, Point<1>(points[q][1]))*

fe1.shape_value(k, Point<1>(points[q][2]))

),

2

);

}

std::cout << "Tensor indices: " << i << " " << j << " " << k << ", "

<< "3d index: " << index_3d << ", error: " << error << "\n";

}

}

}

}

And the output shows all errors to be 0! If this is ok, I will create a pr with a patch to the documentation shortly.

I don't recall whether there is an easy way to achieve what you are looking
for, but some elements definitely do something like what you are trying to
achieve. For example, the difference between FE_Q and FE_DGQ is, in essence,
just a change of basis where the basis functions are permuted. Similarly, the
difference between FE_Q and FE_QHierarchical is similar to the nodal -> modal
change you are interested in. Finally, there is also the case of the FE_Q
constructor that receives a Quadrature object as argument and that then
computes a basis change.
You might want to look into how all of these are implemented. Most of this
kind of functionality exists in some kind of helper function that might be
useful to you.

Thanks a lot for providing these details. For now, I have done this by myself, but I will keep this in mind.

Thanks very much!

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4d767eb0-cf66-b11a-8791-55deab3666de%40colostate.edu.

[Wolfgang Bangerth]

On 5/25/21 11:26 PM, vachan potluri wrote:
> And the output shows all errors to be 0! If this is ok, I will create a pr
> with a patch to the documentation shortly.

Yes, this little program makes sense to me. So that means that indeed shape
functions are lexicographically ordered.

We're looking forward to the patch! :-)

[vachanpo...@gmail.com]

Please have a look at pr12331.

[Wolfgang Bangerth]

On 5/28/21 1:33 AM, vachanpo...@gmail.com wrote:
> Please have a look at pr12331

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F12331&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C753c3f8501374cdabf3108d921aadf6b%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637577840820012778%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=NoTMRHX2wlHZrzn69U5T4LBuru7v9Ofpi1QJPAmyD5o%3D&reserved=0>.

Excellent -- like all patches, we really appreciate yours!
Best
Wolfgang