cut_norm.cpp

#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_signed_distance.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/vectorization.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/q_collection.h>

#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/vector.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <deal.II/non_matching/mapping_info.h>
#include <deal.II/non_matching/mesh_classifier.h>
#include <deal.II/non_matching/quadrature_generator.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools_interpolate.h>

#include <iostream>
#include <memory>
#include <vector>

using namespace dealii;

enum CellCategory
{
  liquid      = 0,
  intersected = 1,
  gas         = 2
};

template <int dim>
void
set_fe_index(const DoFHandler<dim>                  &dof_handler,
             const NonMatching::MeshClassifier<dim> &mesh_classifier)
{
  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      if (!cell->is_locally_owned())
        continue;
      const auto cell_location = mesh_classifier.location_to_level_set(cell);

      if (cell_location == NonMatching::LocationToLevelSet::outside)
        cell->set_active_fe_index(CellCategory::liquid);
      else if (cell_location == NonMatching::LocationToLevelSet::intersected)
        cell->set_active_fe_index(CellCategory::intersected);
      else if (cell_location == NonMatching::LocationToLevelSet::inside)
        cell->set_active_fe_index(CellCategory::gas);
      else
        AssertThrow(false, ExcInternalError());
    }
}

// interpolating this function on the solution vector makes the first phase's
// values 1.0 and the second phase's values 2.0
template <int dim>
class SolutionFunction : public Function<dim>
{
public:
  SolutionFunction()
    : Function<dim>(2)
  {}

  double
  value(const Point<dim> &, const unsigned int component) const override
  {
    return component == 0 ? 1.0 : 2.0;
  }
};

template <int dim, typename number, typename VectorType>
void
compute_intersected_quadrature(
  std::vector<std::shared_ptr<
    NonMatching::MappingInfo<dim, dim, VectorizedArray<number>>>>
                                                          mapping_info_cells,
  const DoFHandler<dim>                                  &level_set_dof_handler,
  const VectorType                                       &level_set,
  const MatrixFree<dim, number, VectorizedArray<number>> &matrix_free,
  const int                                               fe_degree)
{
  hp::QCollection<1> q_collection(QGauss<1>(fe_degree + 1));

  const unsigned int n_lanes = VectorizedArray<number>::size();
  const unsigned int n_cell_batches =
    matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();

  NonMatching::DiscreteQuadratureGenerator<dim> quadrature_generator(
    q_collection, level_set_dof_handler, level_set);

  std::vector<Quadrature<dim>> quad_vec_cells_liquid_outside;
  std::vector<Quadrature<dim>> quad_vec_cells_gas_domain;
  std::vector<typename DoFHandler<dim>::cell_iterator> vector_cell_iterators;
  {
    const unsigned int reserve_size = n_cell_batches * n_lanes;
    quad_vec_cells_liquid_outside.reserve(reserve_size);
    quad_vec_cells_gas_domain.reserve(reserve_size);
    vector_cell_iterators.reserve(reserve_size);
  }

  for (unsigned int cell_batch = 0; cell_batch < n_cell_batches; ++cell_batch)
    for (unsigned int lane = 0; lane < n_lanes; ++lane)
      {
        if (lane < matrix_free.n_active_entries_per_cell_batch(cell_batch))
          {
            vector_cell_iterators.push_back(
              matrix_free.get_cell_iterator(cell_batch, lane));
            quadrature_generator.generate(
              matrix_free.get_cell_iterator(cell_batch, lane));
          }
        else
          {
            // fill empty lanes with dummy data
            vector_cell_iterators.push_back(
              matrix_free.get_cell_iterator(cell_batch, 0));
            quadrature_generator.generate(
              matrix_free.get_cell_iterator(cell_batch, 0));
          }
        quad_vec_cells_liquid_outside.push_back(
          quadrature_generator.get_outside_quadrature());
        quad_vec_cells_gas_domain.push_back(
          quadrature_generator.get_inside_quadrature());
      }

  mapping_info_cells[0]->reinit_cells(vector_cell_iterators,
                                      quad_vec_cells_liquid_outside);
  mapping_info_cells[1]->reinit_cells(vector_cell_iterators,
                                      quad_vec_cells_gas_domain);
}

template <int dim, typename number, int n_components = 1>
inline void
evaluate_intersected_domain(
  FEPointEvaluation<n_components, dim, dim, VectorizedArray<number>>
    &point_eval,
  const FEEvaluation<dim, -1, 0, n_components, number, VectorizedArray<number>>
                                        &cell_eval,
  const EvaluationFlags::EvaluationFlags evaluation_flags,
  const unsigned int                     cell_batch_index,
  const unsigned int                     cell_batch_lane,
  const unsigned int                     n_dofs_per_cell)
{
  static constexpr unsigned int n_lanes = VectorizedArray<number>::size();

  point_eval.reinit(cell_batch_index * n_lanes + cell_batch_lane);
  point_eval.evaluate(StridedArrayView<const number, n_lanes>(
                        &cell_eval.begin_dof_values()[0][cell_batch_lane],
                        n_dofs_per_cell),
                      evaluation_flags);
}

template <int dim, typename number>
number
compute_cut_L1_norm(
  const LinearAlgebra::distributed::Vector<number>       &solution,
  const MatrixFree<dim, number, VectorizedArray<number>> &matrix_free,
  const std::vector<std::shared_ptr<
    NonMatching::MappingInfo<dim, dim, VectorizedArray<number>>>>
                    &mapping_info_cells,
  const FE_Q<dim>   &reference_element,
  const unsigned int dof_idx,
  const unsigned int quad_idx)
{
  if (not solution.has_ghost_elements())
    solution.update_ghost_values();

  number norm_L1_acc = 0.;

  // contributions of the different phases
  number phase_1_contrib             = 0.;
  number phase_2_contrib             = 0.;
  number intersected_phase_1_contrib = 0.;
  number intersected_phase_2_contrib = 0.;

  for (unsigned int cell_batch = 0; cell_batch < matrix_free.n_cell_batches();
       ++cell_batch)
    {
      const auto cell_category = matrix_free.get_cell_category(cell_batch);
      if (cell_category == CellCategory::liquid)
        {
          FEEvaluation<dim, -1, 0, 1, number, VectorizedArray<number>> eval_l(
            matrix_free,
            dof_idx /*dof_no*/,
            quad_idx /*quad_no*/,
            0 /*selected component*/,
            CellCategory::liquid /*active_fe_index*/);

          eval_l.reinit(cell_batch);
          eval_l.read_dof_values_plain(solution);
          eval_l.evaluate(EvaluationFlags::values);

          for (const unsigned int q : eval_l.quadrature_point_indices())
            {
              const auto local_contrib = eval_l.get_value(q) * eval_l.JxW(q);
              for (unsigned int lane = 0;
                   lane <
                   matrix_free.n_active_entries_per_cell_batch(cell_batch);
                   ++lane)
                {
                  phase_1_contrib += local_contrib[lane];
                  norm_L1_acc += local_contrib[lane];
                }
            }
        }
      else if (cell_category == CellCategory::gas)
        {
          FEEvaluation<dim, -1, 0, 1, number, VectorizedArray<number>> eval_g(
            matrix_free,
            dof_idx /*dof_no*/,
            quad_idx /*quad_no*/,
            1 /*selected component*/,
            CellCategory::gas /*active_fe_index*/);

          eval_g.reinit(cell_batch);
          eval_g.read_dof_values_plain(solution);
          eval_g.evaluate(EvaluationFlags::values);

          for (const unsigned int q : eval_g.quadrature_point_indices())
            {
              const auto local_contrib = eval_g.get_value(q) * eval_g.JxW(q);
              for (unsigned int lane = 0;
                   lane <
                   matrix_free.n_active_entries_per_cell_batch(cell_batch);
                   ++lane)
                {
                  phase_2_contrib += local_contrib[lane];
                  norm_L1_acc += local_contrib[lane];
                }
            }
        }
      else if (cell_category == CellCategory::intersected)
        {
          // use FEEvaluation and FEPointEvaluation in combination for
          // intersected cells
          FEEvaluation<dim, -1, 0, 1, number, VectorizedArray<number>>
            eval_cell_l(matrix_free,
                        dof_idx /*dof_no*/,
                        quad_idx /*quad_no*/,
                        0 /*selected component*/,
                        CellCategory::liquid /*active_fe_index*/);
          FEPointEvaluation<1, dim, dim, VectorizedArray<number>>
            eval_subdomain_l(*mapping_info_cells[0], reference_element);
          FEEvaluation<dim, -1, 0, 1, number, VectorizedArray<number>>
            eval_cell_g(matrix_free,
                        dof_idx /*dof_no*/,
                        quad_idx /*quad_no*/,
                        1 /*selected component*/,
                        CellCategory::gas /*active_fe_index*/);
          FEPointEvaluation<1, dim, dim, VectorizedArray<number>>
            eval_subdomain_g(*mapping_info_cells[1], reference_element);

          eval_cell_l.reinit(cell_batch);
          eval_cell_l.read_dof_values_plain(solution);
          eval_cell_g.reinit(cell_batch);
          eval_cell_g.read_dof_values_plain(solution);

          for (unsigned int cell_lane = 0;
               cell_lane <
               matrix_free.n_active_entries_per_cell_batch(cell_batch);
               ++cell_lane)
            {
              evaluate_intersected_domain<dim, number>(
                eval_subdomain_l,
                eval_cell_l,
                EvaluationFlags::values,
                cell_batch,
                cell_lane,
                reference_element.dofs_per_cell);
              evaluate_intersected_domain<dim, number>(
                eval_subdomain_g,
                eval_cell_g,
                EvaluationFlags::values,
                cell_batch,
                cell_lane,
                reference_element.dofs_per_cell);

              for (const unsigned int q_batch :
                   eval_subdomain_l.quadrature_point_indices())
                {
                  const auto local_contrib =
                    eval_subdomain_l.get_value(q_batch) *
                    eval_subdomain_l.JxW(q_batch);
                  for (unsigned int q_lane = 0;
                       q_lane <
                       eval_subdomain_l.n_active_entries_per_quadrature_batch(
                         q_batch);
                       ++q_lane)
                    {
                      intersected_phase_1_contrib += local_contrib[q_lane];
                      norm_L1_acc += local_contrib[q_lane];
                    }
                }

              for (const unsigned int q_batch :
                   eval_subdomain_g.quadrature_point_indices())
                {
                  const auto local_contrib =
                    eval_subdomain_g.get_value(q_batch) *
                    eval_subdomain_g.JxW(q_batch);
                  for (unsigned int q_lane = 0;
                       q_lane <
                       eval_subdomain_g.n_active_entries_per_quadrature_batch(
                         q_batch);
                       ++q_lane)
                    {
                      intersected_phase_2_contrib += local_contrib[q_lane];
                      norm_L1_acc += local_contrib[q_lane];
                    }
                }
            }
        }
    }

  const number L1_norm_phase_1 =
    Utilities::MPI::sum(phase_1_contrib, solution.get_mpi_communicator());
  const number L1_norm_phase_2 =
    Utilities::MPI::sum(phase_2_contrib, solution.get_mpi_communicator());
  const number L1_norm_intersected_phase_1 =
    Utilities::MPI::sum(intersected_phase_1_contrib,
                        solution.get_mpi_communicator());
  const number L1_norm_intersected_phase_2 =
    Utilities::MPI::sum(intersected_phase_2_contrib,
                        solution.get_mpi_communicator());

  if (Utilities::MPI::this_mpi_process(solution.get_mpi_communicator()) == 0)
    {
      std::cout << "L1 norm phase 1:             " << L1_norm_phase_1
                << std::endl;
      std::cout << "L1 norm phase 2:             " << L1_norm_phase_2
                << std::endl;
      std::cout << "L1 norm intersected phase 1: "
                << L1_norm_intersected_phase_1 << std::endl;
      std::cout << "L1 norm intersected phase 2: "
                << L1_norm_intersected_phase_2 << std::endl;
    }

  return Utilities::MPI::sum(norm_L1_acc, solution.get_mpi_communicator());
}

template <int dim>
void
test()
{
  using number = double;

  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
  GridGenerator::subdivided_hyper_cube(tria, 9, -1, 1);

  const unsigned int    degree = 1;
  const FE_Q<dim>       fe_q(degree);
  const FE_Nothing<dim> fe_n;
  hp::FECollection<dim> fe_collection;
  fe_collection.push_back(FESystem<dim, dim>(fe_q, 1, fe_n, 1)); // inside
  fe_collection.push_back(FESystem<dim, dim>(fe_q, 1, fe_q, 1)); // intersected
  fe_collection.push_back(FESystem<dim, dim>(fe_n, 1, fe_q, 1)); // outside
  MappingQ<dim> mapping(degree);

  DoFHandler<dim> dof_handler_cut(tria);
  DoFHandler<dim> dof_handler_ls(tria);

  dof_handler_ls.distribute_dofs(fe_q);
  LinearAlgebra::distributed::Vector<number> level_set;
  level_set.reinit(dof_handler_ls.locally_owned_dofs(),
                   DoFTools::extract_locally_relevant_dofs(dof_handler_ls),
                   tria.get_communicator());
  const Functions::SignedDistance::Plane<dim> signed_distance(
    Point<dim>(), Point<dim>::unit_vector(0));
  VectorTools::interpolate(dof_handler_ls, signed_distance, level_set);
  level_set.update_ghost_values();
  NonMatching::MeshClassifier<dim> mesh_classifier(dof_handler_ls, level_set);
  mesh_classifier.reclassify();
  set_fe_index(dof_handler_cut, mesh_classifier);
  dof_handler_cut.distribute_dofs(fe_collection);

  Quadrature<dim>           quadrature = QGauss<dim>(degree + 1);
  AffineConstraints<number> constraints;
  constraints.close();
  typename MatrixFree<dim, number, VectorizedArray<number>>::AdditionalData
    additional_data;
  additional_data.overlap_communication_computation = false;
  additional_data.hold_all_faces_to_owned_cells     = true;
  additional_data.mapping_update_flags =
    (update_values | update_JxW_values | update_quadrature_points);
  MatrixFree<dim, number, VectorizedArray<number>> matrix_free;
  matrix_free.reinit(
    mapping, dof_handler_cut, constraints, quadrature, additional_data);

  LinearAlgebra::distributed::Vector<number> solution;
  matrix_free.initialize_dof_vector(solution);
  VectorTools::interpolate(dof_handler_cut, SolutionFunction<dim>(), solution);

  std::vector<std::shared_ptr<
    NonMatching::MappingInfo<dim, dim, VectorizedArray<number>>>>
    mapping_info_cells(2);
  mapping_info_cells[0] = std::make_shared<
    NonMatching::MappingInfo<dim, dim, VectorizedArray<number>>>(
    mapping, update_values | update_JxW_values);
  mapping_info_cells[1] = std::make_shared<
    NonMatching::MappingInfo<dim, dim, VectorizedArray<number>>>(
    mapping, update_values | update_JxW_values);
  compute_intersected_quadrature(
    mapping_info_cells, dof_handler_ls, level_set, matrix_free, degree);

  const number norm = compute_cut_L1_norm<dim>(
    solution, matrix_free, mapping_info_cells, fe_q, 0, 0);
  if (Utilities::MPI::this_mpi_process(tria.get_communicator()) == 0)
    std::cout << "L1 norm: " << norm << std::endl;


  // DataOut<dim> data_out;
  // data_out.add_data_vector(dof_handler_cut, solution, "solution");
  // data_out.add_data_vector(dof_handler_ls, level_set, "level_set");
  // Vector<number> mpi_owner(tria.n_active_cells());
  // mpi_owner = Utilities::MPI::this_mpi_process(tria.get_communicator());
  // data_out.add_data_vector(mpi_owner, "owner");
  // data_out.build_patches();
  // data_out.write_vtu_in_parallel("solution.vtu", tria.get_communicator());
}

int
main(int argc, char *argv[])
{
  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
  test<2>();
}