Neumann Boundary condition in MeshWorker for Automatic Differentiation

[Jost Arndt]

Hi there,

so far I had a great learning curve, but now i seem to be stuck.

I try to solve my own little PDE, which is a system of multiple nonlinear time dependent equations and read through plenty of the Tutorials.

Using Dirichlet-Boundary conditions everything worked out finally, now I wanted to switch to Neuman boundary conditions and I get some issue I do not understand.

I followed mostly step-72 to understand how to use automatic differentiation. For Neumann boundary conditions I have now to integrate over the boundary of my domain, therefore I jumped right into the cell_worker lambda function (step-72) right after the for-loop over all quadrature points.

There I tried to stick to the example given in here: https://dealii.org/developer/doxygen/deal.II/classMeshWorker_1_1ScratchData.html and created a for loop:

for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
                        {
                                if (cell->face(f)->at_boundary())//(face->at_boundary())
                                {
                                        std::cout << "reached line 1032" << std::endl;
                                        const auto &fe_face_values = scratch_data.reinit(cell, f);
                                        std::cout << "reached line 1034" << std::endl;

                                         /* here integration using fe_face_value should happen*/

but the last cout never gets executed, but the code crashes. I get a long but cryptic (to me?) stacktrace but also the following error message:

An error occurred in line <3044> of file <./source/fe/fe_values.cc> in function                                                                                                                                          
dealii::FEValuesBase<dim, spacedim>::FEValuesBase(unsigned int, unsigned int, dealii::UpdateFlags, const dealii::Mapping<dim, spacedim>&, const dealii::FiniteElement<dim, spacedim>&) [with int dim = 2; int spacedim
 = 2]                                                                                                                                                                                                                    
The violated condition was:                                                                                                                                                                                              
    n_q_points > 0                                                                                                                                                                                                        
Additional information:                                                                                                                                                                                                  
    There is nothing useful you can do with an FEValues object when using                                                                                                                                                
    a quadrature formula with zero quadrature points!

From here I am quite unsure what I did wrong or how else to create the FEFaceValues Object? I would be really thankful about any hint!

[Wolfgang Bangerth]

The error likely means that MeshWorker hasn't initialized the FEFaceValues
object. Nobody knows any more how MeshWorker really works, and so we try to
discourage use of that framework. But I do know that MeshWorker allows you to
not only loop over the cells, but also the faces of a triangulation. That is
the recommended way to deal with face integrals (rather than dealing with
these face integrals as part of cell integration -- because you want to make
sure that you visit each face only once, rather than twice in your approach)
and if you tell MeshWorker that you also want to use a face worker, then it
should (?) initialize the FEFaceValues object as well.

Best
W.

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

[Jost Arndt]

Thanks for the quick response!

Would you suggest me to read more tutorials on how the MeshWorker works to fix this issue then or try to set up my program without a MeshWorker? I thought I had set it up correctly and tried already for a while to figure out what I did wrong. Also I Initially only started using the MeshWorker because it seemed to be recommended in Step-72 and I thought following the recommendations is probably the best and easiest way. Now I am not so sure anymore about it.

Best,

Jost

[Abbas Ballout]

Jost, 

For looping on the boundaries, you'll have to write an extra lambda amd pass that to the MeshWorker. 

You can look at step 12 or 74 for this.

Also following this cause I want to know what are the MeshWorker alternatives since I am relying heavily on it in my code

Best, 

Abbas

[Wolfgang Bangerth]

> Also following this cause I want to know what are the MeshWorker alternatives
> since I am relying heavily on it in my code

The problem with MeshWorker is that it was never well documented or tested,
and that the people who wrote it are no longer affiliated with the project. We
simply don't have the expertise to support this framework, fix bugs, write
better documentation, etc.

We agreed that for the next release, step-12b, step-16b, and step-39 will
likely be removed. The only thing I think that will remain from that namespace
is MeshWorker::loop(), but the DoFInfo, InfoBox, ... classes will hopefully
all disappear.

Best
Wolfgang