Bug when calling EquationSystems::reinit with NodeElem and discontinuous variables

[dknez]

1. Make or read in a mesh, and add a NodeElem to it. Make an EquationSystems object, es, based on the mesh.

2. Add a system to es that only has discontinuous basis functions, e.g. CONSTANT MONOMIALs, or FIRST L2_LAGRANGEs.

3. Do es.init(); followed by es.reinit(); On the call to reinit, I get the stack trace below. I gather the issue is that reinit fails because the NodeElem has no variables on it associated with the system added in point 2 (though it may have variables in other systems in es).

Assertion `var < this->n_vars(s)' failed.
var = 0
this->n_vars(s) = 0


Stack frames: 17
0: libMesh::print_trace(std::ostream&)
1: libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*)
2: libMesh::DofObject::n_comp(unsigned int, unsigned int) const
3: libMesh::DofMap::_dof_indices(libMesh::Elem const&, int, std::__debug::vector<unsigned int, std::allocator<unsigned int> >&, unsigned int, libMesh::Node const* const*, unsigned int, unsigned long&) const
4: libMesh::DofMap::dof_indices(libMesh::Elem const*, std::__debug::vector<unsigned int, std::allocator<unsigned int> >&) const
5: libMesh::FEMContext::pre_fe_reinit(libMesh::System const&, libMesh::Elem const*)
6: libMesh::GenericProjector<libMesh::OldSolutionValue<double, &(void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const)>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &(void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double> >(unsigned int, libMesh::Point const&, double&, double) const)>, double, libMesh::VectorSetAction<double> >::operator()(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&) const
7: void libMesh::Threads::parallel_for<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, libMesh::GenericProjector<libMesh::OldSolutionValue<double, &(void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const)>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &(void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double> >(unsigned int, libMesh::Point const&, double&, double) const)>, double, libMesh::VectorSetAction<double> > >(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, libMesh::GenericProjector<libMesh::OldSolutionValue<double, &(void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const)>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &(void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double> >(unsigned int, libMesh::Point const&, double&, double) const)>, double, libMesh::VectorSetAction<double> > const&)
8: libMesh::System::project_vector(libMesh::NumericVector<double> const&, libMesh::NumericVector<double>&, int) const
9: libMesh::System::project_vector(libMesh::NumericVector<double>&, int) const
10: libMesh::System::restrict_vectors()
11: libMesh::System::prolong_vectors()
12: libMesh::EquationSystems::reinit_solutions()
13: libMesh::EquationSystems::reinit()
14: ./example-dbg() [0x491074]
15: __libc_start_main
16: ./example-dbg() [0x48fe09]
[0] /home/dknez/software/libmesh_install/dbg_real/libmesh/include/libmesh/dof_object.h, line 784, compiled Apr 19 2018 at 13:22:05

[roystgnr]

Could you set this up in a unit test I can cherry-pick?

[roystgnr]

We do currently call n_comp(sysnum, var) on every node of an element unless we know (due to subdomain restrictions, or because it's a SCALAR) that var doesn't live on that element at all. If var is piecewise constant then that should still be a valid call, it's just supposed to return 0 on nodes.

DofMap::reinit() isn't trying to do anything fancy here, either, just looping over all elements and all nodes and calling set_n_vars_per_group() on each...

[dknez]

I made a PR with the unit test. It replicates the issue for me. Let me know if you need anything further.

[roystgnr]

Let me know if anything like this crops up again. That one fix was sufficient for the one unit test, but clearly we don't have enough test coverage of NodeElem usage.

[dknez]

I'll keep a lookout, thanks.

FYI, I use NodeElem a lot as a way to couple unconnected nodes using GhostingFunctor.

[lindsayad]

In #4345 I'm trying to add a debugging check that makes sure discontinuous variables only live on one value of elem->dim(). Discontinuous variables are generally meant for discontinuous Galerkin methods and generally speaking, inserting dofs associated with different element dimensions is going to lead to lead to singular systems of equations. I'm wondering if it is ok if I associate the CONSTANT MONOMIAL variable in this test with a subdomain tied to the 1D elements?