Adaptivity on lower dim elems in higher dim space

[permcody]

libMesh has a bug that prevents us from using adaptivity with lower dimensional elements in a mesh with a higher dimension number.

Currently we us a workaround where we set the mesh dim == elem dim even if the elements aren't limited to the lower dimensional space. The issue stems from inspecting the wrong dim value in the adaptivity routines.

Cross-referencing: https://github.com/idaholab/moose/issues/4530

[pbauman]

@permcody - could you elaborate a bit? @cahaynes and myself will care about this very much soon. In particular, were asserts tripped or was it different failure modes? I scanned MeshRefinement, but nothing stuck out immediately (I didn't try very hard though) - any pointers? @cahaynes and I can chip in on this.

[friedmud]

Asserts get tripped. It's in the EquationSystems::reinit() path... On Wed, Feb 4, 2015 at 7:18 PM Paul T. Bauman notifications@github.com wrote:

@permcody https://github.com/permcody - could you elaborate a bit? @cahaynes https://github.com/cahaynes and myself will care about this very much soon. In particular, were asserts tripped or was it different failure modes? I scanned MeshRefinement, but nothing stuck out immediately (I didn't try very hard though) - any pointers? @cahaynes https://github.com/cahaynes and I can chip in on this.

— Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/433#issuecomment-72969859.

[pbauman]

Thanks for the pointer @friedmud. Do you guys have a running test that catches this? I'm curious if #442 maybe helped - there was two spots where a _dim got changed to elem->dim() because it was the right thing.

[permcody]

@andrsd was the one that ran into this. I'll work with him to get a test case set up.

[pbauman]

No rush. If we hit it before y'all have a chance, I'll let you know.

[roystgnr]

So here's a question for everyone: when you're mixing lower and higher dimensional elements, is it ever possible for the centroids of two different elements to overlap? If the answer is anything other than an unqualified "no" then there's some LocationMap-based internals that I'm going to need to rewrite...

[jwpeterson]

On Mon, Mar 2, 2015 at 2:30 PM, roystgnr notifications@github.com wrote:

So here's a question for everyone: when you're mixing lower and higher dimensional elements, is it ever possible for the centroids of two different elements to overlap? If the answer is anything other than an unqualified "no" then there's some LocationMap-based internals that I'm going to need to rewrite...

I can imagine a layer of Quads cutting perfectly through the middle of a layer of Hexes... so not an unqualified "no", but I wonder if this a realistic use case.

John

[pbauman]

For me: no. The lower dimensional element will always overlap on the boundary of the higher dimensional element for any case of interest (to me) that I can imagine at the moment. So possibly the centroid of the lower dimensional element with the centroid of the side of the higher dimensional one, but not the centroid of the interior.

[roystgnr]

@pbauman, the recent elasticity tests you've added to GRINS cover the case where a lower-dimensional element shares multiple nodes with a higher-dimensional element? Could we add a variant of one of those with, say, a single uniform refinement? I want to make sure I'm not regressing the multiple-dimensions case when I add the next set of fixes for #436.

[pbauman]

the recent elasticity tests you've added to GRINS cover the case where a lower-dimensional element shares multiple nodes with a higher-dimensional element?

Correct, but only 1D connected to 2D.

Could we add a variant of one of those with, say, a single uniform refinement? I want to make sure I'm not regressing the multiple-dimensions case when I add the next set of fixes for #436.

Will do. In theory, the libMesh unit test in tests/mesh should pick up regressions on the connectivity - I have a uniform refinement test in there.

[roystgnr]

Thanks!

[friedmud]

Yes: we will definitely have the case of 2D elements sharing centroid a of 3D elements and even 1D sharing centroids with 2D and 3D elements... On Mon, Mar 2, 2015 at 4:55 PM roystgnr notifications@github.com wrote:

Thanks!

— Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/433#issuecomment-76834459.

[roystgnr]

Thanks, @friedmud. Out of preemptive fear, I came up with a line of attack last night, and I think "LocationMap" is about to get relegated to "HistoricalCuriosityMap".

Are there any test cases in MOOSE yet that would hit the overlapping-centroids case? If not, I think I can trick one out of GRINS by starting with one of Paul's "membranes plus stiffeners" problems and moving the stiffener to the diagonal of the membrane...

[friedmud]

We don't currently have any in the test suite that I'm aware of. It would be pretty easy to make one though. Let me know if you need one...

Derek

On Tue, Mar 3, 2015 at 9:39 AM roystgnr notifications@github.com wrote:

Thanks, @friedmud https://github.com/friedmud. Out of preemptive fear, I came up with a line of attack last night, and I think "LocationMap" is about to get relegated to "HistoricalCuriosityMap".

Are there any test cases in MOOSE yet that would hit the overlapping-centroids case? If not, I think I can trick one out of GRINS by starting with one of Paul's "membranes plus stiffeners" problems and moving the stiffener to the diagonal of the membrane...

— Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/433#issuecomment-76956730.

[roystgnr]

Could someone point me to a failing case? I seem to be getting good AMR results on 2D-in-3D manifolds myself, and I don't know whether that means Paul or I fixed it or whether that means I'm just not hitting the same code paths.

[roystgnr]

@andrsd sent me a stack trace that clarifies the remaining problem:

Stack frames: 12
0: 0   libmesh_opt.0.dylib                 0x0000000109c0ee55 libMesh::print_trace(std::__1::basic_ostream<char, std::__1::char_traits<char> >&)
+ 1237
1: 1   libmesh_opt.0.dylib                 0x0000000109c0d269 libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*) +
329
2: 2   libmesh_opt.0.dylib                 0x0000000109f6bbbc libMesh::QGauss::init_2D(libMesh::ElemType, unsigned int) + 796
3: 3   libmesh_opt.0.dylib                 0x0000000109ce7a76 libMesh::FE<3u, (libMesh::FEFamily)0>::reinit(libMesh::Elem const*, unsigned int,
double, std::__1::vector<libMesh::Point, std::__1::allocator<libMesh::Point> > const*, std::__1::vector<double, std::__1::allocator<double> >
const*) + 262
4: 4   libmesh_opt.0.dylib                 0x0000000109c2ad92 libMesh::JumpErrorEstimator::reinit_sides() + 66
5: 5   libmesh_opt.0.dylib                 0x0000000109c2a42f libMesh::JumpErrorEstimator::estimate_error(libMesh::System const&,
libMesh::ErrorVector&, libMesh::NumericVector<double> const*, bool) + 3951
6: 6   libmoose-opt.0.dylib                0x00000001092c37e4 Adaptivity::adaptMesh() + 164
7: 7   libmoose-opt.0.dylib                0x0000000109330582 FEProblem::adaptMesh() + 146
8: 8   libmoose-opt.0.dylib                0x00000001094ab5de Transient::incrementStepOrReject() + 46
9: 9   libmoose-opt.0.dylib                0x00000001094ab552 Transient::execute() + 146
10: 10  relap-7-opt                         0x00000001087b955e main + 110
11: 11  libdyld.dylib                       0x00007fff8d6f05c9 start + 1
[0] ../src/quadrature/quadrature_gauss_2D.C, line 1292, compiled nodate at notime
Element type not supported!:27

It looks as if I was succeeding because the mesh_dimension matched the highest elem dim. JumpErrorEstimator is assuming that all elements in the mesh share the same dimension. This is going to be a very intrusive fix (ripping out half the JumpErrorEstimator guts and replacing them with two new @pbauman -improved FEMContext objects), but I don't think there'll be any obstacles; it's analogous to what I've already done in System::project_vector().

[roystgnr]

@andrsd, would you mind trying out that same test with #587?

[andrsd]

I'll try to get to it this afternoon and report back the result...

[andrsd]

With this I get:

Assertion `!_side_fe[dim].empty()' failed.

Stack frames: 13
0: 0   libmesh_devel.0.dylib               0x000000010a03675a libMesh::print_trace(std::__1::basic_ostream<char, std::__1::char_traits<char> >&) + 234
1: 1   libmesh_devel.0.dylib               0x000000010a03251c libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*) + 396
2: 2   libmesh_devel.0.dylib               0x000000010a9351bb libMesh::FEMContext::side_fe_reinit() + 427
3: 3   libmesh_devel.0.dylib               0x000000010a07bf20 libMesh::JumpErrorEstimator::estimate_error(libMesh::System const&, libMesh::ErrorVector&, libMesh::NumericVector<double> const*, bool) + 5504
4: 4   libmoose-devel.0.dylib              0x00000001086460ee Adaptivity::adaptMesh() + 1822
5: 5   libmoose-devel.0.dylib              0x000000010871c61b FEProblem::adaptMesh() + 139
6: 6   libmoose-devel.0.dylib              0x000000010896a257 Transient::incrementStepOrReject() + 199
7: 7   libmoose-devel.0.dylib              0x0000000108969f2a Transient::execute() + 330
8: 8   libmoose-devel.0.dylib              0x000000010885c33c MooseApp::executeExecutioner() + 524
9: 9   libmoose-devel.0.dylib              0x000000010885d0b1 MooseApp::run() + 129
10: 10  relap-7-devel                       0x0000000105e55953 main + 323
11: 11  libdyld.dylib                       0x00007fff8d6f05c9 start + 1
12: 12  ???                                 0x0000000000000003 0x0 + 3
[0] ../src/systems/fem_context.C, line 1373, compiled nodate at notime

[roystgnr]

Looks like there's a missing entry in MeshBase::elem_dimensions? How is your mesh getting built?

[andrsd]

Sorry, I did not have time to look into this lately. In RELAP-7, we generate our 1D and 2D elements from a parametrical descriptions. We have "pipes" that are 1D and rectangles (2D). So, we build all those ourselves and call add_elem(). Do we need to fill in MeshBase::elem_dimensions ourselves or do we need to call something?

[roystgnr]

No worries; I'm not going to get to this before my vacation either, and there's definitely work to be done on the libMesh side still. Even after we find and fix whatever is causing RELAP-7 to fail where fem_system_ex4 succeeds, fem_system_ex4 fails shortly thereafter when the AMR solution projection can't handle the mixed dimensions, and that's going to be a bit tricky for me to fix.

After the add_elem() calls, do you call prepare_for_use()? That should be all that's necessary to get the dimensions cached correctly.

[andrsd]

After the add_elem() calls, do you call prepare_for_use()?

Yes, prepare_for_use() gets called after our mesh is built.

[roystgnr]

@andrsd - now that mixed dim refinement is working in fem_system_ex4, I'd like to fix whatever problem remains in the interaction with relap. Which test should I be modifying/running?

[andrsd]

@roystgnr It is adapt_test.i sitting in tests/misc/adapt (devel branch)

[jwpeterson]

Ping @roystgnr, Rich just asked me about this multi-dimensional adaptive refinement issue. I can confirm that the test in question is still failing in RELAP-7. Do you think you can add this somewhere at the top of your TODO list? I will see if I can replicate it in a pure MOOSE or libMesh code, but no promises that that will be possible.

[roystgnr]

Damn; will do. libMesh master head + moose devel head + relap-7 devel head?

[jwpeterson]

That will show the failure yes. The test is currently skipped, so you just have to run it by hand to see the segfault.

[roystgnr]

Safe to close this now that we've traced the RELAP-7 problem down to app level? We've got increasing test coverage of these cases in GRINS.

[jwpeterson]

:+1: Hurray

[pbauman]

We've got increasing test coverage of these cases in GRINS.

Should still be increasing too with some AMR cases in the near future (CC @cahaynes).