Dear David,
Thank you for your quick reply.
Yes, I need the neighbour information of each cell to compute the boundary
integration in the DG formulation. Therefore, I get the neighbouring cell
in the loop of the active cell by using "const auto neighbor =
cell->neighbor(face_n)"; (similar to step 30).
I am quite sure that the problem is what you are describing, thus it occurs
only when I use the adaptative refinement.
The problem appears here:
#0 ./FEM-DGdeal: dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>,
false>::get_dof_indices(std::vector<unsigned int, std::allocator<unsigned
int> >&) const
#1 ./FEM-DGdeal: FEMwDG<2>::assemble_system(unsigned int)
#2 ./FEM-DGdeal: FEMwDG<2>::run(unsigned int, unsigned int, unsigned int,
unsigned int)
#3 ./FEM-DGdeal: main
And the part of the code where I call the get_dof_indices is this:
void FEMwDG<dim>::assemble_system(const unsigned int typecase)
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<types::global_dof_index>
neighbor_dof_indices(dofs_per_cell);
const UpdateFlags update_flags = update_values | update_gradients |
update_quadrature_points |
update_JxW_values;
const UpdateFlags face_update_flags = update_values |
update_quadrature_points | update_JxW_values |
update_normal_vectors;
const UpdateFlags neighbor_face_update_flags = update_values;
FEValues<dim> fe_values(mapping, fe, quadrature, update_flags);
FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature,
face_update_flags);
FESubfaceValues<dim> fe_subface_values(mapping, fe, face_quadrature,
face_update_flags);
FEFaceValues<dim> fe_face_neighbor(mapping, fe,
face_quadrature,neighbor_face_update_flags);
FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
//matrix dominio
FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_vector(dofs_per_cell);
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell); ui_vi_matrix = 0; cell_vector = 0;
//--------------------------------------DOMAIN
INTEGRATOR--------------------------------------------------------------------------
adv.assemble_cell_term(fe_values, ui_vi_matrix, cell_vector,
verror, variable,cell, typecase);
cell->get_dof_indices(local_dof_indices);
for (unsigned int face_n = 0; face_n <
GeometryInfo<dim>::faces_per_cell; ++face_n)
{
//-------------------------------------FACE AT THE
EXTERNAL
BOUNDARY-----------------------------------------------------------------------------
const auto face = cell->face(face_n);
if (face->at_boundary())
{
fe_face_values.reinit(cell, face_n);
adv.assemble_boundary_term(fe_face_values,
ui_vi_matrix, cell_vector, verror, variable);
}
//-------------------------------------FACE AT THE
INTERNAL
BOUNDARY-----------------------------------------------------------------------------
else
{
const auto neighbor = cell->neighbor(face_n);
if ((dim > 1) && (cell->index() <
neighbor->index()))
{
const unsigned int neighbor2
=cell->neighbor_face_no(face_n);
ue_vi_matrix = 0; ui_ve_matrix = 0; ue_ve_matrix
= 0;
fe_face_values.reinit(cell, face_n);
fe_face_neighbor.reinit(neighbor, neighbor2);
adv.assemble_face_term(fe_face_values,
fe_face_neighbor, ui_vi_matrix, ue_vi_matrix, ui_ve_matrix, ue_ve_matrix,
verror, variable, typecase);
neighbor->get_dof_indices(neighbor_dof_indices);
// ---------------LOCAL ASSEMBLE
----------------------------------------------------------
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
system_matrix.add(local_dof_indices[i],
neighbor_dof_indices[j], ue_vi_matrix(i, j));
system_matrix.add(neighbor_dof_indices[i], local_dof_indices[j],
ui_ve_matrix(i, j));
system_matrix.add(neighbor_dof_indices[i],
neighbor_dof_indices[j], ue_ve_matrix(i, j));
}
}
}
}
// ----------------------------------GLOBAL ASSEMBLE
----------------------------------------------------------
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{for (unsigned int j = 0; j < dofs_per_cell; ++j)
{system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
ui_vi_matrix(i, j));}
system_rhs(local_dof_indices[i]) += cell_vector(i);}
}
}
Thanks :),
Juan Giraldo
El lunes, 9 de diciembre de 2019, 3:15:33 (UTC+8), David Wells escribió:
>
> Dear Juan,
>
> My best guess is that, at some point in the loop over all active
> cells, you are also doing some computations involving the neighbor of
> each cell. TriaAccessor::neighbor() returns the neighboring cell which
> has, at most, the same level number as the current cell: i.e., if the
> neighbor cell has been refined but the current cell has not, then
> trying to get the degrees of freedom from the neighbor will throw
> exactly this exception since the cell returned by the call to
> neighbor() is inactive.
>
> This is only a guess. Could you show us the part of your code that
> triggers this exception? That will help us give you a better answer :)
>
> Thanks,
> David
>
> On Fri, Dec 6, 2019 at 3:53 AM Juan Felipe Giraldo <[email protected]
> <javascript:>> wrote:
> >
> >
> > Hi all!
> >
> > I am currently working on the linear advection problem by using a DG
> formulation, similar to the problem in step 12
> https://www.dealii.org/current/doxygen/deal.II/step_12.html, but using a
> mixed formulation as it is shown in step 20
> https://www.dealii.org/current/doxygen/deal.II/step_20.html. I can solve
> the problem without problems using a global refinement, however, when I
> implement an adaptative refinement I have the following problem after the
> 1st refinement cycle
> >
> > An error occurred in line <3792> of file
> </home/arsen/Documents/instaladores/dealii/include/deal.II/dofs/dof_accessor.templates.h>
>
> in function
> > void dealii::DoFCellAccessor<DoFHandlerType,
> lda>::get_dof_indices(std::vector<unsigned int>&) const [with
> DoFHandlerType = dealii::DoFHandler<2, 2>; bool level_dof_access = false]
> > The violated condition was:
> > this->active()
> > Additional information:
> > get_dof_indices() only works on active cells.
> >
> >
> > I have seen that some people had this problem a couple of years ago, but
> I have not seen a proper solution to solve it.
> >
> > Anyone knows how can I solve it?
> >
> > Thank you so much,
> >
> > Regards,
> >
> > Juan Felipe Giraldo
> >
> > --
> > The deal.II project is located at http://www.dealii.org/
> > For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> > ---
> > You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> > To unsubscribe from this group and stop receiving emails from it, send
> an email to [email protected] <javascript:>.
> > To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/38d23dcb-3028-45eb-a539-737e72c597b0%40googlegroups.com.
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/3d98145f-d1ed-4b4d-ae62-22d61df5ebe9%40googlegroups.com.
