Program Listing for File tensor_product_space.h

Return to documentation for file (include/tensor_product_space.h)

// ---------------------------------------------------------------------
//
// Copyright (C) 2024 by Luca Heltai
//
// This file is part of the ImmersX application, based on
// the deal.II library.
//
// The ImmersX application is free software; you can use
// it, redistribute it, and/or modify it under the terms of the Apache-2.0
// License WITH LLVM-exception as published by the Free Software Foundation;
// either version 3.0 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
// of the ImmersX distribution.
//
// ---------------------------------------------------------------------

#ifndef tensor_product_space_h
#define tensor_product_space_h

#include <deal.II/base/mpi.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>
#include <deal.II/base/parameter_acceptor.h>

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

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

#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q1.h>

#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>

#include <deal.II/numerics/vector_tools.h>

#include <deal.II/particles/particle_handler.h>
#include <deal.II/particles/utilities.h>

#include "reference_cross_section.h"

using namespace dealii;

template <int reduced_dim, int dim, int spacedim, int n_components>
struct TensorProductSpaceParameters : public ParameterAcceptor
{
  TensorProductSpaceParameters();

  static constexpr int cross_section_dim = dim - reduced_dim;

  ReferenceCrossSectionParameters<cross_section_dim, spacedim, n_components>
    section;

  unsigned int fe_degree = 1;

  std::string quadrature_type = "gauss";

  unsigned int n_q_points = 0;

  unsigned int n_quadrature_repetitions = 1;

  double thickness = 0.01;

  std::string thickness_field_name = "";

  std::string reduced_grid_name = "";
};


template <int reduced_dim, int dim, int spacedim, int n_components>
class TensorProductSpace
{
public:
  TensorProductSpace(
    const TensorProductSpaceParameters<reduced_dim, dim, spacedim, n_components>
            &par,
    MPI_Comm mpi_communicator = MPI_COMM_WORLD);

  static constexpr int cross_section_dim = dim - reduced_dim;

  std::function<void(Triangulation<reduced_dim, spacedim> &)>
    preprocess_serial_triangulation = [](auto &) {};

  std::function<void(
    parallel::fullydistributed::Triangulation<reduced_dim, spacedim> &)>
    set_partitioner = [](auto &) {};

  void
  initialize();

  const ReferenceCrossSection<dim - reduced_dim, spacedim, n_components> &
  get_reference_cross_section() const;

  void
  make_reduced_grid_and_properties();

  const DoFHandler<reduced_dim, spacedim> &
  get_dof_handler() const;

  const std::vector<Point<spacedim>> &
  get_locally_owned_qpoints() const;

  const std::vector<std::vector<double>> &
  get_locally_owned_weights() const;

  const std::vector<Point<spacedim>> &
  get_locally_owned_reduced_qpoints() const;

  const std::vector<std::vector<double>> &
  get_locally_owned_reduced_weights() const;

  const std::vector<std::vector<double>> &
  get_locally_owned_section_measure() const;

  void
  update_local_dof_indices(
    const std::map<unsigned int, IndexSet> &remote_q_point_indices);


  const std::vector<types::global_dof_index> &
  get_dof_indices(const types::global_cell_index cell_index) const;


  std::tuple<unsigned int, unsigned int, unsigned int>
  particle_id_to_cell_and_qpoint_indices(const unsigned int qpoint_index) const;


  IndexSet
  locally_owned_qpoints() const;


  IndexSet
  locally_relevant_indices() const;

  auto
  get_quadrature() const -> const Quadrature<reduced_dim> &;

  void
  compute_points_and_weights();

  const parallel::fullydistributed::Triangulation<reduced_dim, spacedim> &
  get_triangulation() const;

  double
  get_scaling(const unsigned int) const;

  const LinearAlgebra::distributed::Vector<double> &
  get_properties() const;

  const DoFHandler<reduced_dim, spacedim> &
  get_properties_dh() const;

  DoFHandler<reduced_dim, spacedim> &
  get_properties_dh();

  const std::vector<std::string> &
  get_properties_names() const;

  std::vector<std::string> &
  get_properties_names();

protected:
  void
  setup_dofs();

  std::map<unsigned int, IndexSet>
  local_q_point_indices_to_global_cell_indices(
    const std::map<unsigned int, IndexSet> &local_q_point_indices) const;

  MPI_Comm mpi_communicator;

  const TensorProductSpaceParameters<reduced_dim, dim, spacedim, n_components>
    &par;

  ReferenceCrossSection<cross_section_dim, spacedim, n_components>
    reference_cross_section;

  parallel::fullydistributed::Triangulation<reduced_dim, spacedim>
    triangulation;

  FESystem<reduced_dim, spacedim> fe;

  Quadrature<reduced_dim> quadrature_formula;

  DoFHandler<reduced_dim, spacedim> dof_handler;

  std::map<types::global_cell_index, std::vector<types::global_dof_index>>
    global_cell_to_dof_indices;

  std::vector<Point<spacedim>> all_qpoints;
  std::vector<std::vector<double>> all_weights;

  std::vector<Point<spacedim>> reduced_qpoints;
  std::vector<std::vector<double>> reduced_weights;

  LinearAlgebra::distributed::Vector<double> properties;

  DoFHandler<reduced_dim, spacedim> properties_dh;

  std::vector<std::string> properties_names;
};


// Template specializations for the TensorProductSpaceParameters



#endif // tensor_product_space_h