hp::FECollection runtime bug in Mixed Mesh

[xjtukklaa]

Hello everyone

I am solving elasticity equations with mixed mesh via dealii. A bug was triggered very accidentally. The dofs seems not correct and fe is not correctly initializated. By looking at the source code of dealii, I believe I should have found the problem.

The order in which hp::FECollection is initialized can lead to bugs. This is the correct oder. if switching the oder of PyramidP and SimplexP, it will trigger the bug. if initializing the right oder as follow, it will works correct in release mode.

fe(FESystem<dim>(FE_PyramidP<dim>(degree)^3), 
    FESystem<dim>(FE_SimplexP<dim>(degree)^3), 
    FESystem<dim>(FE_Q<dim>(degree)^3))

I think the reason as follow. In dealii source code, source/fe/fe_pyramid_p.cc<196> and source/fe/fe_simplex_p.cc<804>. It seems that the pyramid cell will check tetrahedron cell and hexahedron cell. however, the tetrahedron cell can't check pyramid cell. then trigger the bug. my os is wsl Ubuntu.

[bangerth]

@xjtukklaa Do you think you could come up with a small but complete program that shows the bug? Including all #include files etc., something we can compile and run to see the error.

[xjtukklaa]

@bangerth hello and yes, here https://github.com/xjtukklaa/dealii-Mixed-Mesh-Test.git is source code and mesh file with some short description.

[bangerth]

@xjtukklaa Your test does not actually run on my system. After fixing compiling with a trivial patch, I run into the problem that I cannot seem to read the mesh file because it contains entities deal.II does not understand. Does this actually work on your system?

When I strip down the program to basically this...

void Mixed_Mesh::run()
{
  FESystem<dim> fe(FESystem<dim>(FE_PyramidP<dim>(degree)^3), 
                   FESystem<dim>(FE_SimplexP<dim>(degree)^3), 
                   FESystem<dim>(FE_Q<dim>(degree)^3));
}

int main()
{
  deallog.depth_console(2);
  Mixed_Mesh Mixed_Mesh_Problem;
  Mixed_Mesh_Problem.run();
  return 0;
}

...then I get the following error message:

--------------------------------------------------------
An error occurred in line <3282> of file </home/bangerth/p/deal.II/1/dealii/include/deal.II/fe/fe.h> in function
    const dealii::ComponentMask& dealii::FiniteElement<<anonymous>, <anonymous> >::get_nonzero_components(unsigned int) const [with int dim = 3; int spacedim = 3]
The violated condition was: 
    ::dealii::deal_II_exceptions::internals::compare_less_than(i, this->n_dofs_per_cell())
Additional information: 
    Index 12 is not in the half-open range [0,12).

Stacktrace:
-----------
#0  /home/bangerth/p/deal.II/1/install/lib/libdeal_II.g.so.9.6.0-pre: dealii::FiniteElement<3, 3>::get_nonzero_components(unsigned int) const
#1  /home/bangerth/p/deal.II/1/install/lib/libdeal_II.g.so.9.6.0-pre: std::vector<dealii::ComponentMask, std::allocator<dealii::ComponentMask> > dealii::FETools::Compositing::compute_nonzero_components<3, 3>(std::vector<dealii::FiniteElement<3, 3> const*, std::allocator<dealii::FiniteElement<3, 3> const*> > const&, std::vector<unsigned int, std::allocator<unsigned int> > const&, bool)
#2  ./step-1: std::vector<dealii::ComponentMask, std::allocator<dealii::ComponentMask> > dealii::FETools::Compositing::compute_nonzero_components<3, 3>(std::initializer_list<std::pair<std::unique_ptr<dealii::FiniteElement<3, 3>, std::default_delete<dealii::FiniteElement<3, 3> > >, unsigned int> > const&)
#3  ./step-1: dealii::FESystem<3, 3>::FESystem(std::initializer_list<std::pair<std::unique_ptr<dealii::FiniteElement<3, 3>, std::default_delete<dealii::FiniteElement<3, 3> > >, unsigned int> > const&)
#4  ./step-1: dealii::FESystem<3, 3>::FESystem<dealii::FESystem<3, 3>, dealii::FESystem<3, 3>, dealii::FESystem<3, 3>, void>(dealii::FESystem<3, 3>&&, dealii::FESystem<3, 3>&&, dealii::FESystem<3, 3>&&)
#5  ./step-1: Mixed_Mesh::run()
#6  ./step-1: main

Is that what you get as well, and what you are asking about?

[bangerth]

If so, the underlying reason is that in fe_tools.cc (FETools::compute_nonzero_components()) we have the following:

      for (const unsigned int vertex_number :
           fes.front()->reference_cell().vertex_indices())
       ...
                  const unsigned int index_in_base =
                    (fes[base]->n_dofs_per_vertex() * vertex_number +
                     local_index);

                  Assert(comp_start + fes[base]->n_components() <=
                           retval[total_index].size(),
                         ExcInternalError());
                  for (unsigned int c = 0; c < fes[base]->n_components(); ++c)
                    {
                      Assert(c < fes[base]
                                   ->get_nonzero_components(index_in_base)
                                   .size(),
                             ExcInternalError());

Note how in the loop bound we're querying the number of vertices for the reference cell of the first FE in the collection, but then later compute a DoF index based for a concrete FE that may have fewer vertices in the reference cell. That's the core issue here.

[xjtukklaa]

Hello, Thanks for your reply. I upload the solution.vtk file and change the output part code. The test code works fine in my os. I read the "Dealii Simplex Support" and try to solve problem in 3d. In the code, fe is hp::FE_Collection and i don't try to use FESystem. So I don't have this output. By the way, in oder to read mixed msh, the DEAL_II_GMSH_WITH_API must be ON. the problem is that the Dofs are not correct, if i change the order of cell type, and this will lead the linear system crash.

[bangerth]

Yes, I already mentioned what the bug is :-) If you change the order of the elements in the FESystem, you may not get the error message but the result is almost certainly still wrong -- consistent with your observation.

[xjtukklaa]

Yes, this is the problem. Thanks for your help.

[xjtukklaa]

Futher Test In Debug mode(errot order).

Number of active cells: 11728

--------------------------------------------------------
An error occurred in line <657> of file </home/zxl/workspace/dealii/tmp/unpack/deal.II-v9.5.2/source/fe/fe_simplex_p.cc> in function
    dealii::FiniteElementDomination::Domination dealii::FE_SimplexP<dim, spacedim>::compare_for_domination(const dealii::FiniteElement<dim, spacedim>&, unsigned int) const [with int dim = 3; int spacedim = 3]
The violated condition was: 
    false
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).

Stacktrace:
-----------
#0  /home/zxl/workspace/dealii/deal.II-v9.5.2/lib/libdeal_II.g.so.9.5.2: dealii::FE_SimplexP<3, 3>::compare_for_domination(dealii::FiniteElement<3, 3> const&, unsigned int) const
#1  /home/zxl/workspace/dealii/deal.II-v9.5.2/lib/libdeal_II.g.so.9.5.2: dealii::FESystem<3, 3>::compare_for_domination(dealii::FiniteElement<3, 3> const&, unsigned int) const
#2  /home/zxl/workspace/dealii/deal.II-v9.5.2/lib/libdeal_II.g.so.9.5.2: dealii::hp::FECollection<3, 3>::find_dominating_fe(std::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&, unsigned int) const
...

correct order

Number of active cells: 11728
Number of degrees of freedom: 15207

source/fe/fe_pyramid_p.cc<196>

template <int dim, int spacedim>
FiniteElementDomination::Domination
FE_PyramidP<dim, spacedim>::compare_for_domination(
  const FiniteElement<dim, spacedim> &fe_other,
  const unsigned int                  codim) const
{
  Assert(codim <= dim, ExcImpossibleInDim(dim));

  // vertex/line/face domination
  // (if fe_other is derived from FE_SimplexDGP)
  // ------------------------------------
  if (codim > 0)
    if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
        nullptr)
      // there are no requirements between continuous and discontinuous
      // elements
      return FiniteElementDomination::no_requirements;

  // vertex/line/face domination
  // (if fe_other is not derived from FE_SimplexDGP)
  // & cell domination
  // ----------------------------------------
  if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
        dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
    {
      if (this->degree < fe_pp_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_pp_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
             dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
    {
      if (this->degree < fe_p_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_p_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_Q<dim, spacedim> *fe_q_other =
             dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
    {
      if (this->degree < fe_q_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_q_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
             dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
    {
      if (fe_nothing->is_dominating())
        return FiniteElementDomination::other_element_dominates;
      else
        // the FE_Nothing has no degrees of freedom and it is typically used
        // in a context where we don't require any continuity along the
        // interface
        return FiniteElementDomination::no_requirements;
    }

  DEAL_II_NOT_IMPLEMENTED();
  return FiniteElementDomination::neither_element_dominates;
}

source/fe/fe_simplex_p.cc<804>

template <int dim, int spacedim>
FiniteElementDomination::Domination
FE_SimplexP<dim, spacedim>::compare_for_domination(
  const FiniteElement<dim, spacedim> &fe_other,
  const unsigned int                  codim) const
{
  Assert(codim <= dim, ExcImpossibleInDim(dim));

  // vertex/line/face domination
  // (if fe_other is derived from FE_SimplexDGP)
  // ------------------------------------
  if (codim > 0)
    if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
        nullptr)
      // there are no requirements between continuous and discontinuous
      // elements
      return FiniteElementDomination::no_requirements;

  // vertex/line/face domination
  // (if fe_other is not derived from FE_SimplexDGP)
  // & cell domination
  // ----------------------------------------
  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
        dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
    {
      if (this->degree < fe_p_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_p_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_Q<dim, spacedim> *fe_q_other =
             dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
    {
      if (this->degree < fe_q_other->degree)
        return FiniteElementDomination::this_element_dominates;
      else if (this->degree == fe_q_other->degree)
        return FiniteElementDomination::either_element_can_dominate;
      else
        return FiniteElementDomination::other_element_dominates;
    }
  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
             dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
    {
      if (fe_nothing->is_dominating())
        return FiniteElementDomination::other_element_dominates;
      else
        // the FE_Nothing has no degrees of freedom and it is typically used
        // in a context where we don't require any continuity along the
        // interface
        return FiniteElementDomination::no_requirements;
    }

  DEAL_II_NOT_IMPLEMENTED();
  return FiniteElementDomination::neither_element_dominates;
}

In fe_simplex_P.cc, the funciton--compare_for_domination won't check the pyramid cell type. this maybe the one reason for this problem. but i not very sure.

[bangerth]

I thought about this and it occurred to me that what you're doing makes no sense, and that we should produce an error message saying so. FESystem is used to build vector-valued elements where different components of the element correspond to different solution variables. The base elements need to match regarding the reference cell they are defined on. I'll create a patch next week to make that clear.

I suspect that what you were trying to do is define matching elements to work on the cells of a mixed mesh, right? The right class for this is hp::FECollection, not FESystem.

[xjtukklaa]

Yes and i am so sorry for misunderstanding your work. As you said, i just thought about hp::FECollection, not FESystem. I'm very sorry for wasting a lot of time. Thanks for your help.

[bangerth]

No worries, the error message was not very good and this is my opportunity to improve it.