Hi dealii community,
I think I found where the problem lies. In order to reproduce the assertion that gets thrown, i.e.,
An error occurred
in line <1690> of file
</calculate/temp/iwtm009/spack-stage-dealii-9.3.3-h7vrbckbqhxjllhyiwpmxjdp3242wjvv/spack-src/include/deal.II/dofs/dof_accessor.templates.h>
in function
const dealii::FiniteElement<dim,
spacedim>& dealii::DoFAccessor<structdim, dim, spacedim,
lda>::get_fe(unsigned int) const [with int structdim = 1; int dim =
2; int spacedim = 2; bool level_dof_access = false]
The violated condition was:
fe_index_is_active(fe_index) == true
Additional information:
This function can only be called for active FE indices
which comes from calls face_x->get_fe(face_x->nth_active_fe_index(0)) in make_periodicity_constraints and its internal methods, I wrote a MWE on the unit cube with the left and right faces (face_1 and face_2, respectively) being periodic and with subdomains_id (when run with mpirun -np 3) and active_fe_index as seen below
and printed out the face's height, mpi rank, the number of active fe_indices and the first active fe index of the face pair when the assertion's condition is not met. The terminal output is
face_1->center()[1]= 0.875, rank = 0, face_1->n_active_fe_indices() = 0, face_1->nth_active_fe_index(0) = 1, face_2->n_active_fe_indices() = 0, face_2->nth_active_fe_index(0) = 0
face_2->center()[1]= 0.875, rank = 0, face_1->n_active_fe_indices() = 0, face_1->nth_active_fe_index(0) = 1, face_2->n_active_fe_indices() = 0, face_2->nth_active_fe_index(0) = 0
face_1->center()[1]= 0.125, rank = 2, face_1->n_active_fe_indices() = 0, face_1->nth_active_fe_index(0) = 0, face_2->n_active_fe_indices() = 0, face_2->nth_active_fe_index(0) = 0
face_2->center()[1]= 0.125, rank = 2, face_1->n_active_fe_indices() = 0, face_1->nth_active_fe_index(0) = 0, face_2->n_active_fe_indices() = 0, face_2->nth_active_fe_index(0) = 0
As can be seen, both n_active_fe_indices() and nth_active_fe_index(0) print out erroneous values. What seems to happen is that a processor tries to access hp-information of a cell which is neither locally owned nor a ghost cell; thus, triggering the assertion. If that is the case, a probable solution would be to store the pertinent hp-information of the faces at the periodic boundaries in a container which is shared by all processors and relay all the calls made to obtain the hp-information from the face_iterator to the container. Said container should be passed to all the internal methods of make_periodicity_constraints. I will post an update when I have tested this.
Attached is the minimum working example with which the .vtu and terminal output was obtained.
Cheers,
Jose