dpg_tools.h

#ifndef PROJECTS_DPG_DPG_TOOLS_H
#define PROJECTS_DPG_DPG_TOOLS_H
 
#include <lf/base/base.h>
#include <lf/fe/fe.h>
#include <lf/geometry/geometry.h>
#include <lf/mesh/mesh.h>
#include <lf/mesh/utils/utils.h>
#include <lf/quad/quad.h>
 
#include <algorithm>
#include <cmath>
#include <vector>
 
#include "product_dofhandler.h"
#include "product_element_matrix_provider.h"
#include "product_element_vector_provider.h"
#include "product_fe_space.h"
 
namespace projects::dpg {
 
std::vector<lf::quad::QuadRule> BoundaryQuadRule(lf::base::RefEl ref_el,
                                                 const lf::quad::QuadRule& qr);
 
Eigen::MatrixXd OuterNormals(const lf::geometry::Geometry& geometry);
 
class PrescribedSignProvider {
 public:
  PrescribedSignProvider() = default;
  PrescribedSignProvider(const PrescribedSignProvider&) = delete;
  PrescribedSignProvider(PrescribedSignProvider&&) noexcept = default;
  PrescribedSignProvider& operator=(const PrescribedSignProvider&) = delete;
  PrescribedSignProvider& operator=(PrescribedSignProvider&&) noexcept =
      default;
  virtual ~PrescribedSignProvider() = default;
 
  explicit PrescribedSignProvider(
      std::shared_ptr<const lf::mesh::Mesh> mesh_ptr)
      : mesh_ptr_(std::move(mesh_ptr)),
        maxElement_ptr_(std::move(
            lf::mesh::utils::make_CodimMeshDataSet(mesh_ptr_, 1, -1))) {
    init();
  }
 
  [[nodiscard]] int PrescribedSign(const lf::mesh::Entity& element,
                                   const lf::mesh::Entity& edge) const;
 
 private:
  void init();
 
  std::shared_ptr<const lf::mesh::Mesh> mesh_ptr_;
 
  std::shared_ptr<lf::mesh::utils::CodimMeshDataSet<int>> maxElement_ptr_;
};
 
template <typename CONVECTION_COEFF>
lf::mesh::utils::CodimMeshDataSet<bool> flagEntitiesOnInflowBoundary(
    const std::shared_ptr<const lf::mesh::Mesh>& mesh_p,
    CONVECTION_COEFF beta) {
  // flag edges on the boundary
  auto bd_flags = lf::mesh::utils::flagEntitiesOnBoundary(mesh_p, 1);
 
  // initialize inflow boundary flags to fals
  lf::mesh::utils::CodimMeshDataSet<bool> inflow_boundary{mesh_p, 1, false};
 
  // iterate over all cells:
  for (const auto* const cell : mesh_p->Entities(0)) {
    // query geometry
    auto geo_ptr = cell->Geometry();
    auto sub_entities = cell->SubEntities(1);
 
    // calcualte normals
    Eigen::MatrixXd normals = OuterNormals(*geo_ptr);
 
    // iterate over all edges of the cell
    for (int i = 0; i < cell->RefEl().NumSubEntities(1); ++i) {
      const auto* const edge = sub_entities[i];
 
      // check the inflow condition at the center of the edge
      bool inflow = normals.col(i).transpose() *
                        beta(*(sub_entities[i]),
                             (Eigen::MatrixXd(1, 1) << 0.5).finished())[0] <
                    0.0;
      inflow_boundary(*edge) = bd_flags(*edge) && inflow;
    }
  }
  return inflow_boundary;
}
 
template <typename CONVECTION_COEFF>
lf::mesh::utils::CodimMeshDataSet<bool> flagEntitiesOnOutflowBoundary(
    const std::shared_ptr<const lf::mesh::Mesh>& mesh_p,
    CONVECTION_COEFF beta) {
  // flag all entities on the boundary:
  auto bd_flags = lf::mesh::utils::flagEntitiesOnBoundary(mesh_p, 1);
 
  // flag all entities on the inflow boundary:
  auto inflow_boundary = flagEntitiesOnInflowBoundary(mesh_p, beta);
 
  // initialize outflow flags:
  lf::mesh::utils::CodimMeshDataSet<bool> outflow_boundary{mesh_p, 1, false};
 
  // outflow boundary dofs are those, that lie on the boundary, but not on the
  // inflow boundary This guarantues a partition of the boundary edges into
  // inflow and outflow edges.
  for (const auto* const edge : mesh_p->Entities(1)) {
    outflow_boundary(*edge) = bd_flags(*edge) && (!inflow_boundary(*edge));
  }
  return outflow_boundary;
}
 
template <typename SCALAR, typename EDGESELECTOR_DIRICHLET,
          typename EDGESELECTOR_NEUMANN, typename FUNCTION_G,
          typename FUNCTION_H>
std::vector<std::pair<bool, SCALAR>> InitEssentialConditionsFromFunctions(
    const ProductUniformFESpace<SCALAR>& fes, size_type trace_component,
    size_type flux_component, EDGESELECTOR_DIRICHLET&& dirichletcondflag,
    EDGESELECTOR_NEUMANN&& neumanncondflag, FUNCTION_G&& g, FUNCTION_H&& h) {
  // check, that the two components actually differ (necessary?)
  LF_ASSERT_MSG(trace_component != flux_component,
                "Boundary conditions imposed on same component");
 
  // retrive the dofhandlers for the corresponding  trace and flux components
  const auto& trace_fes = fes.ComponentFESpace(trace_component);
  const auto& flux_fes = fes.ComponentFESpace(flux_component);
 
  // flag dofs on the component spaces
  std::vector<std::pair<bool, SCALAR>> componentTraceConditions =
      lf::fe::InitEssentialConditionFromFunction(*trace_fes, dirichletcondflag,
                                                 g);
  std::vector<std::pair<bool, SCALAR>> componentFluxConditions =
      lf::fe::InitEssentialConditionFromFunction(*flux_fes, neumanncondflag, h);
 
  // calculate offsets to transform flags on component spaces
  // to the full product space
  size_type trace_offset = fes.LocGlobMap().Offset(trace_component);
  size_type flux_offset = fes.LocGlobMap().Offset(flux_component);
 
  std::vector<std::pair<bool, SCALAR>> EssentialConditions(
      fes.LocGlobMap().NumDofs(), {false, SCALAR{}});
 
  // flag dofs on product space corresponding to fixed trace dofs
  for (size_type dof = 0; dof < trace_fes->LocGlobMap().NumDofs(); dof++) {
    EssentialConditions[trace_offset + dof] = componentTraceConditions[dof];
  }
 
  // flag dofs on product space corresponding to fixed flux dofs
  for (size_type dof = 0; dof < flux_fes->LocGlobMap().NumDofs(); dof++) {
    EssentialConditions[flux_offset + dof] = componentFluxConditions[dof];
  }
 
  return EssentialConditions;
}
 
template <typename SCALAR>
lf::mesh::utils::CodimMeshDataSet<SCALAR> ElementErrorEstimators(
    const ProductUniformFEDofHandler& trial_dofh,
    // const ProductUniformFEDofHandler& test_dofh,
    std::shared_ptr<ProductElementMatrixProvider<SCALAR>> stiffness_provider,
    std::shared_ptr<ProductElementMatrixProvider<SCALAR>> gramian_provider,
    std::shared_ptr<ProductElementVectorProvider<SCALAR>> load_provider,
    const Eigen::Matrix<SCALAR, Eigen::Dynamic, 1>& sol_vec) {
  // retrive the underlying mesh:
  auto mesh_p = trial_dofh.Mesh();
 
  // retrive the number of cells, initialize the return vector
  lf::mesh::utils::CodimMeshDataSet<SCALAR> element_errors{mesh_p, 0, 0.0};
 
  // compute the element wise errors
  for (const auto* const cell : mesh_p->Entities(0)) {
    // evaluate the providers on the current cell
    auto G = gramian_provider->Eval(*cell);
    auto B = stiffness_provider->Eval(*cell);
    auto l = load_provider->Eval(*cell);
 
    // retrive local solution
    Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> local_solution(
        trial_dofh.NumLocalDofs(*cell));
    auto global_dofs = trial_dofh.GlobalDofIndices(*cell);
    for (size_type i = 0; i < trial_dofh.NumLocalDofs(*cell); ++i) {
      local_solution(i) = sol_vec(global_dofs[i]);
    }
 
    // evaluate posterior error on current element:
    Eigen::Matrix<SCALAR, Eigen::Dynamic, 1> r = l - B * local_solution;
 
    SCALAR e = std::sqrt(r.transpose() * G.ldlt().solve(r));
    element_errors(*cell) = e;
  }
  return element_errors;
}
 
template <typename SCALAR>
SCALAR EvalPosteriorError(
    const std::shared_ptr<const lf::mesh::Mesh>& mesh_p,
    const lf::mesh::utils::CodimMeshDataSet<SCALAR>& local_errors) {
  SCALAR error = 0.0;
  for (const auto* const cell : mesh_p->Entities(0)) {
    error += local_errors(*cell) * local_errors(*cell);
  }
  return std::sqrt(error);
}
 
}  // namespace projects::dpg
#endif  // PROJECTS_DPG_DPG_TOOLS_H