extrude triangulation with n_slices = 1

[Greg Wang]

Hi there,

We want to construct a 3D triangulation by extruding a 2D triangulation (one that potentially contains hanging nodes) and we only want one slice/layer of mesh on the extrusion direction.

Looking around in the GridGenerator namespace led me to the extrude_triangulation function. It’s doing everything we desire except for that (a) the slices/layers on the extrusion direction has to be at least two and (b) the 2D mesh must be a coarse mesh. I’m wondering if there are tips on getting around these two restrictions.

Originally I was using GridGenerator::hyper_rectangle teamed with anisotropic refinement cut_xy. But this ceases to work after the code was re-implemented with parallel::distributed::Triangulation because “this class does not support anisotropic refinement, because it relies on the p4est library that does not support this” [1].

Thanks a lot in advance for any insights and suggestions!

Best,
Greg

---

[1]
https://www.dealii.org/current/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html

​

[Wolfgang Bangerth]

Greg:

> We want to construct a 3D triangulation by extruding a 2D triangulation
> (one that potentially contains hanging nodes) and we only want one
> slice/layer of mesh on the extrusion direction.
>
> Looking around in the GridGenerator namespace led me to the
> extrude_triangulation function. It’s doing everything we desire except
> for that (a) the slices/layers on the extrusion direction has to be at
> least two

I think this is poorly described in the documentation. The number of
slices = the number of cell layers plus one. Two slices => one layer of
cells. I've fixed this here:
https://github.com/dealii/dealii/pull/15028

> and (b) the 2D mesh must be a coarse mesh. I’m wondering if
> there are tips on getting around these two restrictions.

This you can't get around.

> Originally I was using GridGenerator::hyper_rectangle teamed with
> anisotropic refinement cut_xy. But this ceases to work after the code
> was re-implemented with parallel::distributed::Triangulation because
> “this class does not support anisotropic refinement, because it relies
> on the p4est library that does not support this” [1].

Yes, and this is true regardless of how you generate the mesh. You can't
create an anisotropically refined mesh for p::d::T. This also means that
you cannot extrude a refined 2d mesh into 3d that has only one layer of
cells. You can extrude a coarse mesh, though, and then refine the
resulting mesh -- it will then have more than one layer in z-direction
in some locations, however.

Best
W.

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

[Greg Wang]

Hi Wolfgang,

Thanks a lot for clarifying!

I decided to modify the code to adopt a p::shared::T model with METIS and
realized anisotropic refinement doesn’t seem to work for this case either. The
error comes from AffineConstraints’ unacceptance of any refinement cases that
are not isotropic. The same error can be reproduced for my serial code as well.
If I change the RefinementCase, when it’s, say, cut_yz, to cut_xyz, hence going
back to isotropic refinement, both codes run without a problem.

My code is implementing a space-time IP-HDG method. It would be nice to have
anisotropic refinements in order to separate local time stepping and spatial
refinement within a space-time slab. Here is the part of the code where
AffineConstraints gets involved and it’s pretty much copied verbatim from step-51:

constraints.clear(); DoFTools::make_hanging_node_constraints(dof_handler, constraints); std::map<types::boundary_id, const Function<dim> *> boundary_functions; Solution<dim> solution_function(nu); boundary_functions[0] = &solution_function; VectorTools::interpolate_boundary_values(dof_handler, boundary_functions, constraints); constraints.close();

And here is the runtime error message:

An error occurred in line <1102> of file <path/to/deal.II/9.4.1-petsc-hypre/source/source/dofs/dof_tools_constraints.cc> in function
void dealii::DoFTools::internal::make_hp_hanging_node_constraints(const dealii::DoFHandler&, dealii::AffineConstraints&) [with int dim = 3; int spacedim = 3; number = double]
The violated condition was:
cell->face(face)->refinement_case() == RefinementCase::isotropic_refinement
Additional information:
You are trying to use functionality in deal.II that is currently not
implemented. In many cases, this indicates that there simply didn’t
appear much of a need for it, or that the author of the original code
did not have the time to implement a particular case. If you hit this
exception, it is therefore worth the time to look into the code to
find out whether you may be able to implement the missing
functionality. If you do, please consider providing a patch to the
deal.II development sources (see the deal.II website on how to
contribute).

I’m wondering if it would be ill-advised to simply remove the assertion and
re-compile the library. If so, I’m thinking about going a bit deeper and see if
I can come up with a patch. In that case, I would really appreciate some
insights in terms of AffineConstraints’ incompatibilities, if any, with hanging
nodes created by anisotropic refinement.

Thanks,
Greg

​

[Wolfgang Bangerth]

Greg:

> I decided to modify the code to adopt a p::shared::T model with METIS and
> realized anisotropic refinement doesn’t seem to work for this case either. The
> error comes from AffineConstraints’ unacceptance of any refinement cases that
> are not isotropic. The same error can be reproduced for my serial code as well.
> If I change the RefinementCase, when it’s, say, cut_yz, to cut_xyz, hence going
> back to isotropic refinement, both codes run without a problem.

That's interesting. What's the element you are using?

Anisotropic refinement is "supposed" to work, but it is a poorly tested part
of the library and I am not greatly surprised (though I find it disappointing)
that there are cases that don't work as expected.

> And here is the runtime error message:
>
> An error occurred in line <1102> of file
> <path/to/deal.II/9.4.1-petsc-hypre/source/source/dofs/dof_tools_constraints.cc> in function
> void dealii::DoFTools::internal::make_hp_hanging_node_constraints(const
> dealii::DoFHandler&, dealii::AffineConstraints&) [with int dim = 3; int
> spacedim = 3; number = double]
> The violated condition was:
> cell->face(face)->refinement_case() == RefinementCase::isotropic_refinement
> Additional information:
> You are trying to use functionality in deal.II that is currently not
> implemented. In many cases, this indicates that there simply didn’t
> appear much of a need for it, or that the author of the original code
> did not have the time to implement a particular case. If you hit this
> exception, it is therefore worth the time to look into the code to
> find out whether you may be able to implement the missing
> functionality. If you do, please consider providing a patch to the
> deal.II development sources (see the deal.II website on how to
> contribute).
>
> I’m wondering if it would be ill-advised to simply remove the assertion and
> re-compile the library.

Yes :-) Someone placed the assertion there to make sure things that are not
implemented are not actually executed. The way forward is to understand what
is necessary to actually implement the missing functionality. This would
probably start with a minimal example that illustrates the error and nothing else.

The case you seem to end up in is concerned with different elements (the hp
case) *and* hanging nodes. It is conceivable that the case where you have
anisotropic refinement and the same element on both sides is relatively easy
to implement, and so maybe that's the special case to consider first.

[Greg Wang]

Hi Wolfgang,

Thanks a lot for your suggestions! My code is actually not using hp but simply
FE_DGQ<3>(1) and FE_FaceQ<3>(1). The related source code shows that the
boolean of dof_handler.get_fe_collection().hp_constraints_are_implemented()
decides whether internal::make_hp_hanging_node_constraints(dof_handler, constraints) should be called and for the FE_FaceQ class,
make_hp_hanging_node_constraints
always
returns true. This makes sense since the 2D implementation of
make_oldstyle_hanging_node_constraints seems to be looping over cells and
their 1D faces to locate hanging nodes, which won’t apply to FE_FaceQ<3>(1).

As a reference, I wrote a stripped-down version applying IPHDG to step-4’s
Poisson problem and reproduced the same error message. Demo-wise, four cases
were ran with attached pngs numbered accordingly:

1. isotropic refinement + no hanging nodes;
2. isotropic refinement + hanging nodes;
3. anisotropic refinement + no hanging nodes;
4. anisotropic refinement + hanging nodes (error).

Would you consider as a way out of this rewriting the 2D version of
make_oldstyle_hanging_node_constraints by essentially replacing the cell/face
pairs with face/line pairs?

Best regards,
Greg

​

[Wolfgang Bangerth]

Hi Greg:

> Thanks a lot for your suggestions! My code is actually not using hp but
> simply
> FE_DGQ<3>(1) and FE_FaceQ<3>(1). The related source code shows that the
> boolean of dof_handler.get_fe_collection().hp_constraints_are_implemented()
> decides whether internal::make_hp_hanging_node_constraints(dof_handler,
> constraints) should be called and for the FE_FaceQ class,
> make_hp_hanging_node_constraints
> always

> <https://dealii.org/current/doxygen/deal.II/fe__face_8cc_source.html#l00279>
> returns true.

Ah yes, I remember now.

> This makes sense since the 2D implementation of
> make_oldstyle_hanging_node_constraints seems to be looping over cells and
> their 1D faces to locate hanging nodes, which won’t apply to FE_FaceQ<3>(1).

We should just nuke the old-style stuff. It's been ~20 years since we
built finite elements that way, and the hp-function is the way to go.

> As a reference, I wrote a stripped-down version applying IPHDG to step-4’s
> Poisson problem and reproduced the same error message. Demo-wise, four cases
> were ran with attached pngs numbered accordingly:
>

> 1. isotropic refinement + no hanging nodes;
> 2. isotropic refinement + hanging nodes;
> 3. anisotropic refinement + no hanging nodes;
> 4. anisotropic refinement + hanging nodes (error).

>
> Would you consider as a way out of this rewriting the 2D version of
> make_oldstyle_hanging_node_constraints by essentially replacing the
> cell/face
> pairs with face/line pairs?

No, the old-style way just needs to go. What is missing is that the hp
function needs to learn how to do anisotropic refinement.

Since you're already looking into all of the right places: I have
forgotten how anisotropic constraints work. The place where we currently
generate the error you are seeing suggests that the function in question
simply can't do anisotropic refinement at all. But we have programs that
do anisotropic refinement. How do they end up in the old-style function?
In other words, how does this function
https://github.com/dealii/dealii/blob/master/source/dofs/dof_tools_constraints.cc#L1789-L1808
not always go to the hp case given that for almost all finite elements,
hp_constraints_are_implemented() is true?

Or is it true that these days anisotropic refinement only works for DG
elements where constraints are not an issue? I should really know the
answers to these questions, but not a lot of people use anisotropic
refinement and so a lot of the knowledge behind it got lost...

[Greg Wang]

Hi Wolfgang,

Or is it true that these days anisotropic refinement only works
for DG elements where constraints are not an issue?

I grepped anisotropic refinement related funtion names from
tutorial programs and only found step-30. Indeed as you observed,
it’s a DG finite element with L2 conformity, hence sparing
AffineConstraint. It appears get_dpo_vector is implemented
relatively
trivially
for these FEs such that n_dofs_per_face would return 0. The
make_hp_hanging_node_constraints function would have therefore
ignored
them anyways.

No, the old-style way just needs to go. What is missing is that
the hp function needs to learn how to do anisotropic
refinement.

After going through the hp-paper and the
make_hp_hanging_node_constraints function, it seems that, at
least for my current purposes, a lot of code won’t be triggered
since my dofhandlers are not hp capable and hence all affine
constraints fall under the simple case in the hp paper. Though
I’m not sure how susceptible the simple case code would be to the
complications anisotropic refinement brings in step-30. I’ll
spend more time looking into it.

Best regards,
Greg

​

[Wolfgang Bangerth]

On 4/11/23 08:26, Greg Wang wrote:
> I grepped anisotropic refinement related funtion names from
> tutorial programs and only found step-30. Indeed as you observed,
> it’s a DG finite element with L2 conformity, hence sparing
> AffineConstraint. It appears get_dpo_vector is implemented
> relatively
> trivially

> <https://dealii.org/current/doxygen/deal.II/fe__dgq_8cc_source.html#l00183>

> for these FEs such that n_dofs_per_face would return 0. The
> make_hp_hanging_node_constraints function would have therefore
> ignored

> <https://dealii.org/current/doxygen/deal.II/dof__tools__constraints_8cc_source.html#l01096>

> them anyways.
>
> No, the old-style way just needs to go. What is missing is that
> the hp function needs to learn how to do anisotropic
> refinement.
>
> After going through the hp-paper and the
> make_hp_hanging_node_constraints function, it seems that, at
> least for my current purposes, a lot of code won’t be triggered
> since my dofhandlers are not hp capable and hence all affine

> constraints fall under the /simple case/ in the hp paper. Though

> I’m not sure how susceptible the simple case code would be to the
> complications anisotropic refinement brings in step-30. I’ll
> spend more time looking into it.

Greg:
thanks for poking around and reminding me what works and what doesn't. I
believe that what you found means that the "new-style" (hp-style) function we
like to use just doesn't support anisotropic refinement at all. We should
probably put a better assertion error message into the relevant places.

I suspect that making it work would require quite a lot of work because the
finite elements simply do not provide the kind of information this function
would require if one wanted to implement anisotropic refinement. I don't know
how much time you want to spend on it -- if you were interested, I could try
and point out what finite elements would need to provide, and you could see
how one would implement that. But I don't want to give the impression that
that would be an easy fix.

What you could try, if you wanted, is to see what it would take to make your
case work with the old function. It is conceivable that that, too, would
require a lot of work -- I just don't know that function well enough any more.

I'm sorry I don't have better news :-(

Best
Wolfgang