project_solution fails with infinite elements

[BalticPinguin]

Dear all,
I have a bugreport and a related question about intended usage/performance:
In my program I have an input-function (i.e. I can compute its value at any point) and which needs to be projected into the FEM-space.
Currently, I am solving a FEM-system M*x=f where M is just the mass matrix and f is computed by integrating the input-quantity.

Now I wanted to change this, using the project_vector class, aiming at a speed-up.
The question-part now is: How does this function work? It looks quite complex, and seems to do something similar, but on an element-wise level.
Am I right? I see that this saves me from constructing the large matrix; but can anyone tell me something about the performance in terms of time? I am just curious...

Now to the bug: My mesh contains infinite elements,leading to the project_solution part to break with the traceout:

ERROR: Bad ElemType = INFPRISM12 for SECOND order approximation!
Stack frames: 18
0: libMesh::print_trace(std::ostream&)
1: libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*)
2: /home/tm162/bin/libmesh/dev_libmesh/lib/libmesh_dbg.so.0(+0x1b2c7a0) [0x7ffff60c87a0]
3: libMesh::FE<3u, (libMesh::FEFamily)0>::n_dofs(libMesh::ElemType, libMesh::Order)
4: libMesh::FE<3u, (libMesh::FEFamily)0>::n_shape_functions(libMesh::ElemType, libMesh::Order)
5: libMesh::FE<3u, (libMesh::FEFamily)0>::map(libMesh::Elem const*, libMesh::Point const&)
6: libMesh::FE<3u, (libMesh::FEFamily)0>::inverse_map(libMesh::Elem const*, libMesh::Point const&, double, bool)
7: libMesh::FE<3u, (libMesh::FEFamily)0>::inverse_map(libMesh::Elem const*, std::__debug::vector<libMesh::Point, std::allocator<libMesh::Point> > const&, std::__debug::vector<libMesh::Point, std::allocator<libMesh::Point> >&, double, bool)
8: libMesh::FE<3u, (libMesh::FEFamily)0>::edge_reinit(libMesh::Elem const*, unsigned int, double, std::__debug::vector<libMesh::Point, std::allocator<libMesh::Point> > const*, std::__debug::vector<double, std::allocator<double> > const*)
9: libMesh::FEMContext::edge_fe_reinit()
10: libMesh::GenericProjector<libMesh::FEMFunctionWrapper<std::complex<double> >, libMesh::FEMFunctionWrapper<libMesh::VectorValue<std::complex<double> > >, std::complex<double>, libMesh::VectorSetAction<std::complex<double> > >::operator()(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&) const
11: void libMesh::Threads::parallel_for<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, libMesh::GenericProjector<libMesh::FEMFunctionWrapper<std::complex<double> >, libMesh::FEMFunctionWrapper<libMesh::VectorValue<std::complex<double> > >, std::complex<double>, libMesh::VectorSetAction<std::complex<double> > > >(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, libMesh::GenericProjector<libMesh::FEMFunctionWrapper<std::complex<double> >, libMesh::FEMFunctionWrapper<libMesh::VectorValue<std::complex<double> > >, std::complex<double>, libMesh::VectorSetAction<std::complex<double> > > const&)
12: libMesh::System::project_vector(libMesh::NumericVector<std::complex<double> >&, libMesh::FEMFunctionBase<std::complex<double> >*, libMesh::FEMFunctionBase<libMesh::VectorValue<std::complex<double> > >*, int) const
13: libMesh::System::project_vector(libMesh::NumericVector<std::complex<double> >&, libMesh::FunctionBase<std::complex<double> >*, libMesh::FunctionBase<libMesh::VectorValue<std::complex<double> > >*, int) const
14: libMesh::System::project_solution(libMesh::FunctionBase<std::complex<double> >*, libMesh::FunctionBase<libMesh::VectorValue<std::complex<double> > >*) const
15: /scratch/tm162/bin/FreeWilly/FrWll-dbg(+0xc179a) [0x55555561579a]

I think, the issue here is around stack # 9: the FEMContext should look carefully, whether the fe_reinit or inf_reinit should be called.
Another way to look at it is that the _edge_fe is filled via FEAbstract::build(dim, fe_type), which is not sensitive to it.

I must admit that I am a bit stuck at this point, because I am not very familiar with the classes involved and I fear that it is not easy to generalize this part.
Probably someone can comment on this and give me some direction how to fix this?

Many thanks in advance,
Hubert

[roystgnr]

I must admit that I am a bit stuck at this point, because I am not very familiar with the classes involved and I fear that it is not easy to generalize this part.

That's a tragic understatement. There's currently a massive rewrite in process in #1938 because the current code fails to give the right solution for some corner cases with subdomain-restricted variables.

Probably someone can comment on this and give me some direction how to fix this?

For starters, can you set up a unit test for us?

At first glance it looks like you've identified the correct root problem, FEMContext isn't infinite-element-aware and needs to be made so. I'm not sure how simple this will be to do efficiently without slowing down the common case, but I suspect allocating InfFE objects on the fly might be reasonable.

(Actually, IMHO making InfFE an alternative class rather than a helper class was a silly mistake which breaks a natural abstraction, but I don't think there's a quick fix for that so we can just keep using the standard idiom for now)

[BalticPinguin]

Thank you for that answer; I hope I manage to make a test by tomorrow.

[BalticPinguin]

Ok; I have modified the miscellaneous_ex1 now to use the project_vector() function to create the rhs, giving me the same as posted above. I am not sure doing so is physically meaningfull, but I think it is somewhat similar to what was done before and was not too much of work.
You can find the code here
Is it convenient as is?

Just to avoid miscommunication: Unless told otherwise, I will not work on this issue and rather use my old scheme in the meantime, at least until #1938 is done. However, I am willing to contribute work and, having some direction, can implement some generalisations.
Best regards

[BalticPinguin]

great; issue is resolved now.
Many thanks!