integrator.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <prismspf/core/system_wide.h>
 
#include <prismspf/config.h>
 
#include <vector>
 
PRISMS_PF_BEGIN_NAMESPACE
 
template <unsigned int dim, unsigned int degree, typename number>
class Integrator
{
public:
  using VectorType = dealii::LinearAlgebra::distributed::Vector<number>;
 
  template <int rank>
  static dealii::Vector<number>
  integrate(const dealii::DoFHandler<dim> &dof_handler, const auto &solution_vector)
  {
    constexpr unsigned int expected_components =
      dealii::Tensor<rank, dim>::n_independent_components;
    Assert(dof_handler.get_fe().n_components() == expected_components,
           dealii::ExcMessage("The provided DoFHandler does not have the same number of "
                              "components as the expected ones. For scalar fields there "
                              "should be 1 component."));
 
    // Update ghosts
    solution_vector.update_ghost_values();
 
    // Set quadrature rule and FEValues to update the JxW values
    dealii::FEValues<dim> fe_values(dof_handler.get_fe(),
                                    SystemWide<dim, degree>::quadrature,
                                    dealii::update_values | dealii::update_JxW_values);
 
    // Get the number of quadrature points
    const unsigned int num_quad_points = SystemWide<dim, degree>::quadrature.size();
 
    // Create a value vector
    std::vector<dealii::Vector<number>> quad_values(num_quad_points,
                                                    dealii::Vector<number>(
                                                      expected_components));
 
    // Loop over the cells provided by the DoFHandler
    dealii::Vector<number> value(expected_components);
    for (const auto &cell : dof_handler.active_cell_iterators())
      {
        if (cell->is_locally_owned())
          {
            // Reinitialize the cell
            fe_values.reinit(cell);
 
            // Get the values
            fe_values.get_function_values(solution_vector, quad_values);
 
            // Sum up the product of the JxW and values at each quadrature point to
            // compute the element integral.
            for (unsigned int q_point = 0; q_point < num_quad_points; ++q_point)
              {
                for (unsigned int component = 0; component < expected_components;
                     component++)
                  {
                    value[component] +=
                      quad_values[q_point][component] * fe_values.JxW(q_point);
                  }
              }
          }
      }
 
    for (unsigned int component = 0; component < expected_components; component++)
      {
        value[component] = dealii::Utilities::MPI::sum(value[component], MPI_COMM_WORLD);
      }
    return value;
  }
};
 
PRISMS_PF_END_NAMESPACE