Hey,
It might be my fault but I made a simplest case pair of files attached.
In 1D for adaptive meshes, the matrix free operator is not iterating over
the face/point that is between refinement levels.
I attached a code that basically just makes a mesh and refines the last
cell. I use the dof_handler to print out all the coordinates of the faces.
It does iterate over all the faces. However, the matrix free operator in
its face iterator doesn't recognize the point between refinement levels as
a face so I have it print out all the coordinates of faces and its
obviously missing one.
The "mesh" is just a line from [0,16]. First it is divided into 4 cells.
Then the final cell is further refined so there should be faces at:
0,4,8,12,14,16.
Again, feel free to correct me if I made a dumb mistake somewhere, but I
did try the code in 2D and there were no problems there. The code is super
simplistic and will not give you enough information about coordinates if
you step up to 3D.
Thanks for any help you can give,
Sean Johnson
--
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 visit
https://groups.google.com/d/msgid/dealii/4cbe6bb4-961b-40d0-96d5-f686a878cb7bn%40googlegroups.com.
namespace ADM_Problem
{
using namespace dealii;
template <int dim, int degree, int n_points_1d>
class Pi_Operator
{
public:
static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
Pi_Operator();
void reinit(const Mapping<dim> & mapping,
const std::vector<const DoFHandler<dim> *> &dof_handlers,
const std::vector<const AffineConstraints<double> *> &constraints);
void apply(const LinearAlgebra::distributed::Vector<double> &src,
LinearAlgebra::distributed::Vector<double> & dst) const;
void
initialize_vector(LinearAlgebra::distributed::Vector<double> &vector) const;
MatrixFree<dim, double> data;
private:
void local_apply_cell(
const MatrixFree<dim, double> & data,
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> & cell_range) const;
void local_apply_face(
const MatrixFree<dim, double> & data,
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> & face_range) const;
void local_apply_boundary_face(
const MatrixFree<dim, double> & data,
LinearAlgebra::distributed::Vector<double> & dst,
const LinearAlgebra::distributed::Vector<double> &src,
const std::pair<unsigned int, unsigned int> & face_range) const;
};
template <int dim, int degree, int n_points_1d>
Pi_Operator<dim, degree, n_points_1d>::Pi_Operator()
{}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::reinit(
const Mapping<dim> & mapping,
const std::vector<const DoFHandler<dim> *> &dof_handlers,
const std::vector<const AffineConstraints<double> *> &constraints)
{
const std::vector<Quadrature<dim>> quadratures = {QGauss<dim>(n_q_points_1d)};
typename MatrixFree<dim, double>::AdditionalData additional_data;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points |
update_values);
additional_data.mapping_update_flags_inner_faces =
(update_JxW_values | update_quadrature_points | update_normal_vectors |
update_values);
additional_data.mapping_update_flags_boundary_faces =
(update_JxW_values | update_quadrature_points | update_normal_vectors |
update_values);
additional_data.tasks_parallel_scheme =
MatrixFree<dim, double>::AdditionalData::none;
data.reinit(
mapping, dof_handlers, constraints, quadratures, additional_data);
}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::initialize_vector(
LinearAlgebra::distributed::Vector<double> &vector) const
{
data.initialize_dof_vector(vector,0);
}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::local_apply_cell(
const MatrixFree<dim, double> &,
LinearAlgebra::distributed::Vector<double> & /*dst*/,
const LinearAlgebra::distributed::Vector<double> & /*pi_vec*/,
const std::pair<unsigned int, unsigned int> & /*cell_range*/) const
{
}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::local_apply_face(
const MatrixFree<dim, double> &,
LinearAlgebra::distributed::Vector<double> & /*dst*/,
const LinearAlgebra::distributed::Vector<double> &v_vec,
const std::pair<unsigned int, unsigned int> & face_range) const
{
FEFaceEvaluation<dim, degree, n_points_1d> pi_in(data, true);
FEFaceEvaluation<dim, degree, n_points_1d> pi_out(data, false);
std::cout << "#######################################################" << std::endl;
std::cout << "Matrix Free Face Iteration" << std::endl;
std::cout << "#######################################################" << std::endl;
for (unsigned int face = face_range.first; face < face_range.second; ++face)
{
pi_in.reinit(face);
pi_in.gather_evaluate(v_vec, EvaluationFlags::values);
pi_out.reinit(face);
pi_out.gather_evaluate(v_vec, EvaluationFlags::values);
for (unsigned int q = 0; q < pi_in.n_q_points; ++q)
{
VectorizedArray<double> x_val = pi_in.quadrature_point(q)[0];
VectorizedArray<double> y_val;
if (dim == 2) { y_val = pi_in.quadrature_point(q)[1];}
for (unsigned int k = 0; k < x_val.size(); ++k){
if (dim == 1){
std::cout << "x_val: " << x_val[k] << std::endl;
}
else if (dim == 2)
{
std::cout << "(x,y): ( " << x_val[k] << " , " << y_val[k] << " )" << std::endl;
}
}
}
}
}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::local_apply_boundary_face(
const MatrixFree<dim, double> &,
LinearAlgebra::distributed::Vector<double> & /*dst*/,
const LinearAlgebra::distributed::Vector<double> &/*v_vec*/,
const std::pair<unsigned int, unsigned int> & /*face_range*/) const
{
}
template <int dim, int degree, int n_points_1d>
void Pi_Operator<dim, degree, n_points_1d>::apply(
const LinearAlgebra::distributed::Vector<double> &src,
LinearAlgebra::distributed::Vector<double> & dst) const
{
{
data.loop(&Pi_Operator::local_apply_cell,
&Pi_Operator::local_apply_face,
&Pi_Operator::local_apply_boundary_face,
this,
dst,
src,
true,
MatrixFree<dim, double>::DataAccessOnFaces::values,
MatrixFree<dim, double>::DataAccessOnFaces::values);
}
}
} //namespace ADM_Problem
/* ---------------------------------------------------------------------
*
* Copyright (C) 2020 - 2022 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/time_stepping.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_bicgstab.h>
//#include <deal.II/lac/solver_idr.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/hdf5.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/dofs/dof_renumbering.h>
//#include <deal.II/numerics/dof_tools.h>
// Originally neccessary include files that may not be neccessary anymore.
/*
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <deal.II/base/convergence_table.h>
#include <deal.II/distributed/solution_transfer.h>
*/
#include <fstream>
#include <iomanip>
#include <iostream>
// The following file includes the CellwiseInverseMassMatrix data structure
// that we will use for the mass matrix inversion, the only new include
// file for this tutorial program:
#include <deal.II/matrix_free/operators.h>
namespace ADM_Problem
{
using namespace dealii;
constexpr unsigned int dimension = 1;
constexpr unsigned int fe_degree = 3;
constexpr unsigned int n_q_points_1d = 1;
constexpr unsigned int n_global_refinements = 2;
constexpr double outer_boundary = 16.;
constexpr double inner_boundary = 0.;
constexpr double refine_cells_number_above = 2.5;
constexpr bool uniform_refinement = false;
}
#include "pi_operator.h"
namespace ADM_Problem
{
using namespace dealii;
template <int dim>
class AdmProblem
{
public:
AdmProblem();
void run();
private:
void make_grid();
void make_constraints();
void make_operator();
void refinement();
void print_x_of_faces();
ConditionalOStream pcout;
Triangulation<dim> triangulation;
FESystem<dim> fe_DG;
MappingQ<dim> mapping;
DoFHandler<dim> dof_handler_DG;
AffineConstraints<double> constraints_DG;
Pi_Operator<dim, fe_degree, n_q_points_1d> pi_operator;
LinearAlgebra::distributed::Vector<double> pi_solution;
};
// The constructor for this class is unsurprising: We set up a parallel
// triangulation based on the `MPI_COMM_WORLD` communicator, a vector finite
// element with `dim+2` components for density, momentum, and energy, a
// high-order mapping of the same degree as the underlying finite element,
// and initialize the time and time step to zero.
template <int dim>
AdmProblem<dim>::AdmProblem()
: pcout(std::cout)
,triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
, fe_DG(FE_DGQ<dim>(fe_degree))
, mapping(fe_degree)
, dof_handler_DG(triangulation)
, pi_operator()
{}
template <int dim>
void AdmProblem<dim>::make_grid()
{
GridGenerator::hyper_cube(triangulation,
inner_boundary,
outer_boundary);
triangulation.refine_global(n_global_refinements);
dof_handler_DG.distribute_dofs(fe_DG);
}
template <int dim>
void AdmProblem<dim>::make_constraints()
{
const IndexSet locally_relevant_dofs_DG =
DoFTools::extract_locally_relevant_dofs(dof_handler_DG);
constraints_DG.clear();
constraints_DG.reinit(locally_relevant_dofs_DG);
constraints_DG.close();
}
template <int dim>
void AdmProblem<dim>::make_operator()
{
const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler_DG};
const std::vector<const AffineConstraints<double> *> constraints_list = {&constraints_DG};
pi_operator.reinit(mapping, dof_handlers, constraints_list);
pi_operator.initialize_vector(pi_solution);
}
template <int dim>
void AdmProblem<dim>::refinement()
{
Vector<double> cell_numbers(triangulation.n_active_cells());
for (unsigned int i = 0; i < triangulation.n_active_cells(); ++i){
cell_numbers[i] = i;
}
if (!uniform_refinement && triangulation.n_active_cells()-1 > refine_cells_number_above){
GridRefinement::refine(triangulation,cell_numbers,refine_cells_number_above);
triangulation.prepare_coarsening_and_refinement();
triangulation.execute_coarsening_and_refinement();
dof_handler_DG.distribute_dofs(fe_DG);
}
else if (triangulation.n_active_cells()-1 > refine_cells_number_above)
{
triangulation.refine_global(1);
dof_handler_DG.distribute_dofs(fe_DG);
}
}
template<int dim>
void AdmProblem<dim>::print_x_of_faces()
{
std::cout << "############# Printing X values of center of faces ################" << std::endl;
for (const auto &cell: dof_handler_DG.active_cell_iterators()){
for (const auto &f: cell->face_indices()){
const auto center = cell->face(f)->center();
if (dim == 1)
std::cout << "Face at x_val: " << center[0] << std::endl;
else if (dim == 2){
std::cout << "Face at (x,y): ( " << center[0] << " , " << center[1] << " )" << std::endl;
}
}
}
}
template <int dim>
void AdmProblem<dim>::run()
{
make_grid();
refinement();
make_constraints();
make_operator();
print_x_of_faces();
pi_operator.apply(pi_solution,pi_solution);
print_x_of_faces();
}
} // namespace Adm_Problem
int main(int /*argc*/, char** /*argv*/)
{
using namespace ADM_Problem;
using namespace dealii;
try
{
deallog.depth_console(0);
AdmProblem<dimension> adm_problem;
adm_problem.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
