FE_DGQLegendre/FE_DGP: accept other polynomials #18177
Conversation
There was a problem hiding this comment.
The FE_Q class (I believe) does this by accepting a set of interpolation points. Do you want to use the new constructor to create modal elements, or why do you take a vector of polynomials?
Separately, don't you need to address the various functions that return nodal points? I think this should be tested too :-)
include/deal.II/fe/fe_dgp.h
Outdated
| /** | ||
| * Constructor for tensor product polynomials of given polynomials @p poly. | ||
| */ | ||
| FE_DGP(const std::vector<Polynomials::Polynomial<double>> &poly); | ||
|
|
There was a problem hiding this comment.
The comment was only copied from the other constructor, but I don't think it is clear at this point how the space is actually constructed. The given argument describes a one dimensional polynomial space. From this, we fill the multi-dimensional space by a truncated tensor product that picks the polynomials from the given space starting index zero in this list. So if one wants a complete polynomial degree as is suggested by the general class description, one needs to make sure that the polynomial space is sorted appropriately (to whatever space one wants to get after taking the tensor product).
In other words, I think this is a dangerous constructor because we do not enforce anything on a polynomial space, and if one hands in strange polynomials one will get a strange element. In particular, I can hand in a vector of degree much higher than the size of the polynomials, if I intend to implement some funny bubble space, but then FiniteElement::degree and similar are wrong and might produce confusing results. We could say it is up to the user to make this decision by making sure that the polynomial space is later matched with the right downstream code such as suitable quadrature formulas, but it should be clearly stated what the effects of the choices are at this point.
include/deal.II/fe/fe_dgq.h
Outdated
| /** | ||
| * Constructor for tensor product polynomials based on the | ||
| * given polynomials @p poly. | ||
| */ | ||
| FE_DGQLegendre(const std::vector<Polynomials::Polynomial<double>> &poly); |
There was a problem hiding this comment.
I don't like this constructor because it will yield a polynomial space that is not a Legendre space if the given general polynomial space is not a 1D Legendre space. I think we should rename this class into FE_DGQModal or similar (or FE_DGQArbitraryPolynomialSpace, because we are not tied to use a modal space, either, if we directly hand in the polynomials), implement all functionality there including the constructor you use here, and set up FE_DGQLegendre as a light-weight class on top of FE_DGQModal.
Also note the comments I gave in the other class: It is the user's responsibility to provide useful a useful polynomial space.
| FiniteElementData<dim>(get_dpo_vector(poly.size() - 1), | ||
| 1, | ||
| degree, | ||
| poly.size() - 1, | ||
| FiniteElementData<dim>::L2), |
There was a problem hiding this comment.
This is one of the points I tried to make in the comment above: We derive the degree from the size of the polynomials. Which clearly makes sense for a sane polynomial space. But by giving the space into the hand of the user, our internal variable names get very blurry.
There was a problem hiding this comment.
I think that could probably be addressed by adding an assertion to the body of constructor, right? Right now, you already assert that for each given polynomial, p.degree() <= degree. You will also need to assert that the max of the p.degree() is == degree.
|
@bangerth @kronbichler Thank you for the feedback! In the last commits, I made the following changes:
The class hierarchy of |
Indeed, there is such a constructor. But that is protected. Maybe someone knows the historical reason?
Yes, I would like to be able model elements but not necessarily based on Legendre polynomials. |
source/fe/fe_dgq.cc
Outdated
| @@ -1005,11 +1005,45 @@ FE_DGQArbitraryNodes<dim, spacedim>::clone() const | |||
|
|
|||
|
|
|||
|
|
|||
| // ---------------------------------- FE_DGQLegendre ------------------------- | |||
There was a problem hiding this comment.
| // ---------------------------------- FE_DGQLegendre ------------------------- | |
| // ---------------------------------- FE_DGQArbitraryPolynomialSpace ------------------------- |
I think you refer to the class Regarding the history: We historically had |
|
@kronbichler Thank you for the history lesson :) |
|
Regarding the doxygen warning. It complains about |
|
Hm, the error is strange, I see Could it be that the std::vector can be constructed from a lenght parameter |
| FiniteElementData<dim>(get_dpo_vector(poly.size() - 1), | ||
| 1, | ||
| degree, | ||
| poly.size() - 1, | ||
| FiniteElementData<dim>::L2), |
There was a problem hiding this comment.
I think that could probably be addressed by adding an assertion to the body of constructor, right? Right now, you already assert that for each given polynomial, p.degree() <= degree. You will also need to assert that the max of the p.degree() is == degree.
| class FE_DGQLegendre : public FE_DGQ<dim, spacedim> | ||
| class FE_DGQArbitraryPolynomialSpace : public FE_DGQ<dim, spacedim> |
There was a problem hiding this comment.
I've got to admit that this looks a bit backward. FE_DGQ should be the derived class, and FE_DGQArbitraryPolynomialSpace the base class. Do you know why it is this way around?
(not just Legendre basis)