variable_coefficients.h

/*  ______________________________________________________________________
 *
 *  ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
 *
 *  Copyright (C) 2023 by the ExaDG authors
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 *  ______________________________________________________________________
 */
#ifndef INCLUDE_EXADG_OPERATORS_VARIABLE_COEFFICIENTS_H_
#define INCLUDE_EXADG_OPERATORS_VARIABLE_COEFFICIENTS_H_
 
namespace ExaDG
{
template<typename coefficient_type>
class VariableCoefficients
{
public:
  template<int dim, typename Number>
  void
  initialize(dealii::MatrixFree<dim, Number> const & matrix_free,
             unsigned int const                      quad_index,
             bool const                              store_face_data_in,
             bool const                              store_cell_based_face_data_in)
  {
    if(not store_face_data_in)
      AssertThrow(not store_cell_based_face_data_in,
                  dealii::ExcMessage("Storing only cell-based face data does not make sense"
                                     " if storing face data is disabled."));
 
    store_face_data            = store_face_data_in;
    store_cell_based_face_data = store_cell_based_face_data_in;
 
    reinit(matrix_free, quad_index);
  }
 
  void
  set_coefficients(coefficient_type const & constant_coefficient)
  {
    fill(constant_coefficient);
  }
 
  coefficient_type
  get_coefficient_cell(unsigned int const cell, unsigned int const q) const
  {
    return coefficients_cell[cell][q];
  }
 
  void
  set_coefficient_cell(unsigned int const       cell,
                       unsigned int const       q,
                       coefficient_type const & coefficient)
  {
    coefficients_cell[cell][q] = coefficient;
  }
 
  coefficient_type
  get_coefficient_face(unsigned int const face, unsigned int const q) const
  {
    return coefficients_face[face][q];
  }
 
  void
  set_coefficient_face(unsigned int const       face,
                       unsigned int const       q,
                       coefficient_type const & coefficient)
  {
    coefficients_face[face][q] = coefficient;
  }
 
  coefficient_type
  get_coefficient_face_neighbor(unsigned int const face, unsigned int const q) const
  {
    return coefficients_face_neighbor[face][q];
  }
 
  void
  set_coefficient_face_neighbor(unsigned int const       face,
                                unsigned int const       q,
                                coefficient_type const & coefficient)
  {
    coefficients_face_neighbor[face][q] = coefficient;
  }
 
  coefficient_type
  get_coefficient_face_cell_based(unsigned int const cell_based_face, unsigned int const q) const
  {
    return coefficients_face_cell_based[cell_based_face][q];
  }
 
  void
  set_coefficient_face_cell_based(unsigned int const       cell_based_face,
                                  unsigned int const       q,
                                  coefficient_type const & coefficient)
  {
    coefficients_face_cell_based[cell_based_face][q] = coefficient;
  }
 
private:
  template<int dim, typename Number>
  void
  reinit(dealii::MatrixFree<dim, Number> const & matrix_free, unsigned int const quad_index)
  {
    coefficients_cell.reinit(matrix_free.n_cell_batches(), matrix_free.get_n_q_points(quad_index));
 
    if(store_face_data)
    {
      coefficients_face.reinit(matrix_free.n_inner_face_batches() +
                                 matrix_free.n_boundary_face_batches(),
                               matrix_free.get_n_q_points_face(quad_index));
      coefficients_face_neighbor.reinit(matrix_free.n_inner_face_batches(),
                                        matrix_free.get_n_q_points_face(quad_index));
 
      if(store_cell_based_face_data)
      {
        unsigned int const n_faces_per_cell =
          matrix_free.get_dof_handler().get_triangulation().get_reference_cells()[0].n_faces();
 
        coefficients_face_cell_based.reinit(matrix_free.n_cell_batches() * n_faces_per_cell,
                                            matrix_free.get_n_q_points_face(quad_index));
      }
    }
  }
 
  void
  fill(coefficient_type const & constant_coefficient)
  {
    coefficients_cell.fill(constant_coefficient);
 
    if(store_face_data)
    {
      coefficients_face.fill(constant_coefficient);
      coefficients_face_neighbor.fill(constant_coefficient);
 
      if(store_cell_based_face_data)
        coefficients_face_cell_based.fill(constant_coefficient);
    }
  }
 
  dealii::Table<2, coefficient_type> coefficients_cell;
 
  dealii::Table<2, coefficient_type> coefficients_face;
 
  dealii::Table<2, coefficient_type> coefficients_face_neighbor;
 
  dealii::Table<2, coefficient_type> coefficients_face_cell_based;
 
  bool store_face_data{false};
 
  bool store_cell_based_face_data{false};
};
 
} // namespace ExaDG
 
#endif /* INCLUDE_EXADG_OPERATORS_VARIABLE_COEFFICIENTS_H_ */