Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add ReferenceCell::volume() and ::barycenter(). #13678

Merged
merged 5 commits into from May 6, 2022

Conversation

bangerth
Copy link
Member

@bangerth bangerth commented May 5, 2022

I've needed something like a midpoint quadrature in ASPECT, and that turned out to be surprisingly difficult to generate for general cell shapes. Knowing the volume and barycenter is a good first step.

/rebuild

@bangerth bangerth added this to the Release 9.4 milestone May 5, 2022
tests/grid/reference_cell_type_03.output Outdated Show resolved Hide resolved
@bangerth
Copy link
Member Author

bangerth commented May 6, 2022

So updated.

Comment on lines 59 to 76
std::unique_ptr<FiniteElement<dim>> fe;
if (reference_cell == ReferenceCells::Line)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Triangle)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Quadrilateral)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Tetrahedron)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Wedge)
fe = std::make_unique<FE_WedgeP<dim>>(1);
else if (reference_cell == ReferenceCells::Pyramid)
fe = std::make_unique<FE_PyramidP<dim>>(1);
else if (reference_cell == ReferenceCells::Hexahedron)
fe = std::make_unique<FE_Q<dim>>(1);

// Set up the objects to compute an integral on the reference cell
FEValues<dim> fe_values(*fe, *q, update_JxW_values);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for nit-picking, but shouldn't this work:

Suggested change
std::unique_ptr<FiniteElement<dim>> fe;
if (reference_cell == ReferenceCells::Line)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Triangle)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Quadrilateral)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Tetrahedron)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Wedge)
fe = std::make_unique<FE_WedgeP<dim>>(1);
else if (reference_cell == ReferenceCells::Pyramid)
fe = std::make_unique<FE_PyramidP<dim>>(1);
else if (reference_cell == ReferenceCells::Hexahedron)
fe = std::make_unique<FE_Q<dim>>(1);
// Set up the objects to compute an integral on the reference cell
FEValues<dim> fe_values(*fe, *q, update_JxW_values);
FE_Nothing<dim> fe;
// Set up the objects to compute an integral on the reference cell
FEValues<dim> fe_values(fe, *q, update_JxW_values);

Since you never evaluate any of the shape functions?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would hope this does not work ;) But there is a FE_Nothing version that takes a reference cell!

Comment on lines 59 to 73
std::unique_ptr<FiniteElement<dim>> fe;
if (reference_cell == ReferenceCells::Line)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Triangle)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Quadrilateral)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Tetrahedron)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Wedge)
fe = std::make_unique<FE_WedgeP<dim>>(1);
else if (reference_cell == ReferenceCells::Pyramid)
fe = std::make_unique<FE_PyramidP<dim>>(1);
else if (reference_cell == ReferenceCells::Hexahedron)
fe = std::make_unique<FE_Q<dim>>(1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here

Suggested change
std::unique_ptr<FiniteElement<dim>> fe;
if (reference_cell == ReferenceCells::Line)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Triangle)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Quadrilateral)
fe = std::make_unique<FE_Q<dim>>(1);
else if (reference_cell == ReferenceCells::Tetrahedron)
fe = std::make_unique<FE_SimplexP<dim>>(1);
else if (reference_cell == ReferenceCells::Wedge)
fe = std::make_unique<FE_WedgeP<dim>>(1);
else if (reference_cell == ReferenceCells::Pyramid)
fe = std::make_unique<FE_PyramidP<dim>>(1);
else if (reference_cell == ReferenceCells::Hexahedron)
fe = std::make_unique<FE_Q<dim>>(1);
FE_Nothing<dim> fe;

Copy link
Member

@kronbichler kronbichler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK apart from the comments regarding the test.

@bangerth
Copy link
Member Author

bangerth commented May 6, 2022

I had assumed that FE_Nothing doesn't work for anything other than hypercubes, and didn't know about the added ability of that class. Thanks for teaching me about it :-)

So fixed now!

else if (reference_cell == ReferenceCells::Pyramid)
q = std::make_unique<QGaussPyramid<dim>>(2);
else if (reference_cell == ReferenceCells::Hexahedron)
q = std::make_unique<QGauss<dim>>(2);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also picky, but we could use ReferenceCell::get_gauss_type_quadrature(2)) here for improved genericity.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would you allow me to do this as a follow-up for cloned tests? I think there is value in testing the equivalence of the output since what constitutes a "Gauss" formula is not unique for anything other than tensor-product formulas (and the function you propose may or may not return the formulas I use above).

else if (reference_cell == ReferenceCells::Pyramid)
q = std::make_unique<QGaussPyramid<dim>>(2);
else if (reference_cell == ReferenceCells::Hexahedron)
q = std::make_unique<QGauss<dim>>(2);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and here

@drwells drwells merged commit 754de58 into dealii:master May 6, 2022
@bangerth bangerth deleted the ref-cell branch May 6, 2022 19:59
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants