Re: [deal.II] get_dof_indices() only works on active cells - blocks vector
Dear Luca,
Thank you very much for your reply, your answer was very precise and,
therefore, I could solve my problem. Effectively, I was not considering the
subfaces when I was assembling, and of course the "get_dof_indices"
function was not recognizing the nodes neighbour nodes those cells.
Note: As I am working with mixed formulation, I need to compute the cell,
boundary and face terms for two different variables (in blocks). I had some
problems to declare them by separate using the meshworker. So, it was
easier for me to declare the functions independently.
Regards,
Juan Giraldo
El lunes, 9 de diciembre de 2019, 17:31:24 (UTC+8), luca.heltai escribió:
>
> Dear Juan,
>
> The logic you should use is: when your neighbour is finer, you should not
> assemble the full face term from the coarse neighbour. Instead, you have to
> assemble your DG terms using a FESubfaceValues initialized with the correct
> arguments on the coarse face, to match *each* of the fine FEFaceValues.
> This logic is tricky to get right.
>
> May I suggest that you use the MeshWorker::mesh_loop function?
>
> That function takes care of these issues for you, by calling a user
> provided assembly function for cells, faces and boundary faces, and passing
> to that function the right arguments with which you can initialise your
> FEFaceValues and FESubFaceValues.
>
> You can see it in action in step-12 (development version of deal.II), or
> in the tests/meshworker/mesh_loop_*.cc files and
> tests/meshworker/scratch_data_*.cc
>
> Best,
> Luca.
>
> > On 9 Dec 2019, at 4:55, Juan Felipe Giraldo > wrote:
> >
> > 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,
> false>::get_dof_indices(std::vector 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::assemble_system(const unsigned int typecase)
> > {
> > const unsigned int dofs_per_cell =
> dof_handler.get_fe().dofs_per_cell;
> > std::vector
> local_dof_indices(dofs_per_cell);
> > std::vector
> 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 fe_values(mapping, fe, quadrature,
> update_flags);
> > FEFaceValues fe_face_values(mapping, fe, face_quadrature,
> face_update_flags);
> > FESubfaceValues fe_subface_values(mapping, fe,
> face_quadrature, face_update_flags);
> > FEFaceValues fe_face_neighbor(mapping, fe,
> face_quadrature,neighbor_face_update_flags);
> >
> > FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell);
> //matrix dominio
> > FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell);
> > FullMatrix ui_ve_matrix(dofs_per_cell, dofs_per_cell);
> > FullMatrix ue_ve_matrix(dofs_per_cell, dofs_per_cell);
> >
> > Vector 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::faces_per_cell; ++face_n)
> > {
> > //-FACE AT THE
> EXTERNAL
> BOUNDARY-
>
>
> > const auto face = cell->face(face_n);
> >
Re: [deal.II] get_dof_indices() only works on active cells - blocks vector
Dear Juan,
The logic you should use is: when your neighbour is finer, you should not
assemble the full face term from the coarse neighbour. Instead, you have to
assemble your DG terms using a FESubfaceValues initialized with the correct
arguments on the coarse face, to match *each* of the fine FEFaceValues. This
logic is tricky to get right.
May I suggest that you use the MeshWorker::mesh_loop function?
That function takes care of these issues for you, by calling a user provided
assembly function for cells, faces and boundary faces, and passing to that
function the right arguments with which you can initialise your FEFaceValues
and FESubFaceValues.
You can see it in action in step-12 (development version of deal.II), or in the
tests/meshworker/mesh_loop_*.cc files and tests/meshworker/scratch_data_*.cc
Best,
Luca.
> On 9 Dec 2019, at 4:55, Juan Felipe Giraldo wrote:
>
> 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,
> false>::get_dof_indices(std::vector 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::assemble_system(const unsigned int typecase)
> {
> const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
> std::vector local_dof_indices(dofs_per_cell);
> std::vector
> 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 fe_values(mapping, fe, quadrature, update_flags);
> FEFaceValues fe_face_values(mapping, fe, face_quadrature,
> face_update_flags);
> FESubfaceValues fe_subface_values(mapping, fe, face_quadrature,
> face_update_flags);
> FEFaceValues fe_face_neighbor(mapping, fe,
> face_quadrature,neighbor_face_update_flags);
>
> FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell); //matrix
> dominio
> FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell);
> FullMatrix ui_ve_matrix(dofs_per_cell, dofs_per_cell);
> FullMatrix ue_ve_matrix(dofs_per_cell, dofs_per_cell);
>
> Vector 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::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,
Re: [deal.II] get_dof_indices() only works on active cells - blocks vector
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,
false>::get_dof_indices(std::vector >&) 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::assemble_system(const unsigned int typecase)
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
std::vector local_dof_indices(dofs_per_cell);
std::vector
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 fe_values(mapping, fe, quadrature, update_flags);
FEFaceValues fe_face_values(mapping, fe, face_quadrature,
face_update_flags);
FESubfaceValues fe_subface_values(mapping, fe, face_quadrature,
face_update_flags);
FEFaceValues fe_face_neighbor(mapping, fe,
face_quadrature,neighbor_face_update_flags);
FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell);
//matrix dominio
FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix ui_ve_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix ue_ve_matrix(dofs_per_cell, dofs_per_cell);
Vector 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::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));
}
}
}
Re: [deal.II] get_dof_indices() only works on active cells - blocks vector
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 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 > > in function > void dealii::DoFCellAccessor lda>::get_dof_indices(std::vector&) 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 dealii+unsubscr...@googlegroups.com. > 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 dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CABrTbYS7_D7cz8BVjObSSy3a6zAHMSCrwn75cj79XBVnOicCpw%40mail.gmail.com.
[deal.II] get_dof_indices() only works on active cells - blocks vector
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 in function void dealii::DoFCellAccessor::get_dof_indices(std::vector&) 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 dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/38d23dcb-3028-45eb-a539-737e72c597b0%40googlegroups.com.
