Hello everyone,
I'm solving a 1D explicit DG-FEM problem and I've encountered behavior that
I don't understand using MeshWorker::mesh_loop.
- I'm using the MeshWorker::assemble_own_interior_faces_both flag since
I want to do work on each face corresponding to the same cell interior
- If the grid is created using typical GridGenerator calls, the order is
exactly as I'd expect: first, work is done on the cell, then the faces
(boundary or interior), followed by the next cell (further, all mass
matrices are the same)
- However, if I manually adapt the mesh by refining some cells, the
order seems to change; the face work is done without respect to the cell
interior last worked on
I've attached a stripped-down program that illustrates this behavior on a
toy mesh in 1D, as well as the output.
- Ultimately, my goal is to assemble an inverse mass matrix on each
cell, and apply it to a residual vector containing interior and face
contributions (which can be done element-wise since the elements are
FE_DGQ)
- I was attempting to store the inverse via CopyData and then apply it
in the copy worker.
- However, I'm finding that due to the order of execution, I can't rely
on the face work being done immediately after the cell work, and the
inverse mass matrix stored to CopyData is often from another cell than the
faces being worked on.
- I could do the assembly and application of the inverse mass matrices
separately, or store the inverse mass matrices in a map to the cell
iterators, but I'm curious why this ordering occurs.
Am I misunderstanding how mesh_loop is supposed to work? Any guidance would
be appreciated!
Corbin
--
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/ac4d4550-b014-4440-b2c2-227dba71fffen%40googlegroups.com.
// mesh_loop_order.cc
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/meshworker/copy_data.h>
#include <deal.II/meshworker/mesh_loop.h>
#include <deal.II/meshworker/scratch_data.h>
using namespace dealii;
struct CopyData
{
FullMatrix<double> inverse_mass_matrix;
Vector<double> residual_contribution;
std::vector<types::global_dof_index> local_dof_indices;
template<class Iterator>
void reinit(const Iterator &cell, unsigned int dofs_per_cell)
{
inverse_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
residual_contribution.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
}
};
template<int dim>
class Wrapper
{
public:
Wrapper(DoFHandler<dim>& dof_handler,
FE_DGQ<dim>& fe)
: dofh(&dof_handler)
, fe(&fe)
, cell_quadrature(fe.tensor_degree() + 1)
, face_quadrature(fe.tensor_degree() + 1)
{ }
void do_mesh_loop();
protected:
const QGauss<dim> cell_quadrature;
const QGauss<dim-1> face_quadrature;
SmartPointer<DoFHandler<dim>> dofh;
const SmartPointer<FE_DGQ<dim>> fe;
};
template<int dim>
void Wrapper<dim>::do_mesh_loop()
{
using Iterator = typename DoFHandler<dim>::active_cell_iterator;
using ScratchData = MeshWorker::ScratchData<dim>;
const auto cell_worker = [&](const Iterator &cell,
ScratchData &sd,
CopyData ©_data) {
std::cout << "cell worker entered on cell " << cell << std::endl;
const FEValues<dim> &fev = sd.reinit(cell);
const unsigned int n_cell_dofs = fev.get_fe().n_dofs_per_cell();
std::cout << "\tcopy data reinitialized" << std::endl;
copy_data.reinit(cell, n_cell_dofs);
// cell work here (assemble inverse mass matrix)
};
const auto face_worker = [&](const Iterator &cell,
const unsigned int &f,
const unsigned int &sf,
const Iterator &ncell,
const unsigned int &nf,
const unsigned int &nsf,
ScratchData &sd,
CopyData ©_data) {
std::cout << "face worker called on cell " << cell << std::endl;
const FEInterfaceValues<dim> &fe_iv = sd.reinit(cell, f, sf, ncell, nf, nsf);
};
const auto boundary_worker = [&](const Iterator &cell,
const unsigned int &face_no,
ScratchData &sd,
CopyData ©_data) {
std::cout << "boundary worker entered on cell " << cell << std::endl;
const FEFaceValues<dim> &fe_fv = sd.reinit(cell, face_no);
};
AffineConstraints<double> constraints;
constraints.close();
const auto copier = [&](const CopyData &c) {
std::cout << "copier entered: " << std::endl;
};
const MappingQ1<dim> mapping;
const UpdateFlags cell_flags(update_values
| update_gradients
| update_JxW_values
| update_quadrature_points);
const UpdateFlags face_flags(update_values
| update_normal_vectors
| update_JxW_values
| update_quadrature_points);
ScratchData scratch_data(mapping,
*fe,
cell_quadrature, cell_flags,
face_quadrature, face_flags);
CopyData copy_data;
MeshWorker::mesh_loop(dofh->begin_active(),
dofh->end(),
cell_worker,
copier,
scratch_data,
copy_data,
MeshWorker::assemble_own_cells |
MeshWorker::assemble_boundary_faces |
MeshWorker::assemble_own_interior_faces_both,
boundary_worker,
face_worker,
/* serial */ 1,1);
}
template<int dim>
void print_mesh_loop_order()
{
// declare a normal triangulation
Triangulation<dim> tri;
GridGenerator::subdivided_hyper_cube(
tri,
/* nx */ 3,
/*left*/ -4,
/*right*/ 4);
const unsigned int p_order = 2;
FE_DGQ<dim> fe(p_order);
DoFHandler<dim> dofh(tri);
{ // non-adaptive case
dofh.distribute_dofs(fe);
std::cout << "standard mesh test: " << dofh.n_dofs() << " dofs" << std::endl;
Wrapper<dim> wrapper(dofh, fe);
wrapper.do_mesh_loop();
}
{ // make selective refinement (every other cell)
unsigned int counter = 0;
for (const auto &cell : dofh.active_cell_iterators())
{
if (counter % 2 == 0)
cell->set_refine_flag();
counter++;
}
tri.execute_coarsening_and_refinement();
dofh.reinit(tri);
dofh.distribute_dofs(fe);
std::cout << "\n\nadaptive mesh test: " << dofh.n_dofs() << " dofs" << std::endl;
Wrapper<dim> wrapper(dofh, fe);
wrapper.do_mesh_loop();
}
}
int main()
{
print_mesh_loop_order<1>();
return 0;
}
standard mesh test: 9 dofs
cell worker entered on cell 0.0
copy data reinitialized
boundary worker entered on cell 0.0
face worker called on cell 0.0
cell worker entered on cell 0.1
copy data reinitialized
face worker called on cell 0.1
face worker called on cell 0.1
cell worker entered on cell 0.2
copy data reinitialized
face worker called on cell 0.2
boundary worker entered on cell 0.2
copier entered:
copier entered:
copier entered:
adaptive mesh test: 15 dofs
cell worker entered on cell 0.1
copy data reinitialized
cell worker entered on cell 1.0
copy data reinitialized
boundary worker entered on cell 1.0
face worker called on cell 1.0
cell worker entered on cell 1.1
copy data reinitialized
face worker called on cell 1.1
face worker called on cell 1.1
face worker called on cell 0.1
cell worker entered on cell 1.2
copy data reinitialized
face worker called on cell 1.2
face worker called on cell 0.1
face worker called on cell 1.2
cell worker entered on cell 1.3
copy data reinitialized
face worker called on cell 1.3
boundary worker entered on cell 1.3
copier entered:
copier entered:
copier entered:
copier entered:
copier entered:
