discontinuous_fe_constant.h

 
#ifndef PROJECTS_DPG_DISCONTINUOUS_FE_CONSTANT_H
#define PROJECTS_DPG_DISCONTINUOUS_FE_CONSTANT_H
 
#include <lf/fe/fe.h>
#include <lf/uscalfe/uscalfe.h>
 
#include <typeinfo>
 
#include "dpg.h"
 
namespace projects::dpg {
 
template <typename SCALAR>
class FeDiscontinuousO0Tria final
    : public lf::fe::ScalarReferenceFiniteElement<SCALAR> {
 public:
  FeDiscontinuousO0Tria(const FeDiscontinuousO0Tria&) = default;
  FeDiscontinuousO0Tria(FeDiscontinuousO0Tria&&) noexcept = default;
  FeDiscontinuousO0Tria& operator=(const FeDiscontinuousO0Tria&) = default;
  FeDiscontinuousO0Tria& operator=(FeDiscontinuousO0Tria&&) noexcept = default;
  FeDiscontinuousO0Tria() = default;
  ~FeDiscontinuousO0Tria() override = default;
 
  [[nodiscard]] lf::base::RefEl RefEl() const override {
    return lf::base::RefEl::kTria();
  }
 
  [[nodiscard]] unsigned Degree() const override { return 0; }
 
  [[nodiscard]] size_type NumRefShapeFunctions() const override { return 1; }
 
  [[nodiscard]] size_type NumRefShapeFunctions(dim_t codim) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
  [[nodiscard]] size_type NumRefShapeFunctions(
      dim_t codim, sub_idx_t /*subidx*/) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  EvalReferenceShapeFunctions(const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 2,
                  "Reference coordinates must be 2-vectors");
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(n_pts, 1);
    return result;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  GradientsReferenceShapeFunctions(
      const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 2,
                  "Reference coordinates must be 2-vectors");
 
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, 2 * n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(2 * n_pts, 0);
    return result;
  }
 
  [[nodiscard]] Eigen::MatrixXd EvaluationNodes() const override {
    Eigen::MatrixXd nodes(2, 1);
    nodes << 1.0 / 3.0, 1.0 / 3.0;
    return nodes;
  }
 
  [[nodiscard]] size_type NumEvaluationNodes() const override {
    return NumRefShapeFunctions();
  }
};
 
template <typename SCALAR>
class FeDiscontinuousO0Quad final
    : public lf::fe::ScalarReferenceFiniteElement<SCALAR> {
 public:
  FeDiscontinuousO0Quad(const FeDiscontinuousO0Quad&) = default;
  FeDiscontinuousO0Quad(FeDiscontinuousO0Quad&&) noexcept = default;
  FeDiscontinuousO0Quad& operator=(const FeDiscontinuousO0Quad&) = default;
  FeDiscontinuousO0Quad& operator=(FeDiscontinuousO0Quad&&) noexcept = default;
  FeDiscontinuousO0Quad() = default;
  ~FeDiscontinuousO0Quad() override = default;
 
  [[nodiscard]] lf::base::RefEl RefEl() const override {
    return lf::base::RefEl::kQuad();
  }
 
  [[nodiscard]] unsigned Degree() const override { return 0; }
 
  [[nodiscard]] size_type NumRefShapeFunctions() const override { return 1; }
 
  [[nodiscard]] size_type NumRefShapeFunctions(dim_t codim) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
  [[nodiscard]] size_type NumRefShapeFunctions(
      dim_t codim, sub_idx_t /*subidx*/) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  EvalReferenceShapeFunctions(const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 2,
                  "Reference coordinates must be 2-vectors");
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(n_pts, 1);
    return result;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  GradientsReferenceShapeFunctions(
      const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 2,
                  "Reference coordinates must be 2-vectors");
 
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, 2 * n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(2 * n_pts, 0);
    return result;
  }
 
  [[nodiscard]] Eigen::MatrixXd EvaluationNodes() const override {
    Eigen::MatrixXd nodes(2, 1);
    nodes << 0.5, 0.5;
    return nodes;
  }
  [[nodiscard]] size_type NumEvaluationNodes() const override {
    return NumRefShapeFunctions();
  }
};
 
template <typename SCALAR>
class FeDiscontinuousO0Segment final
    : public lf::fe::ScalarReferenceFiniteElement<SCALAR> {
 public:
  FeDiscontinuousO0Segment(const FeDiscontinuousO0Segment&) = default;
  FeDiscontinuousO0Segment(FeDiscontinuousO0Segment&&) noexcept = default;
  FeDiscontinuousO0Segment& operator=(const FeDiscontinuousO0Segment&) =
      default;
  FeDiscontinuousO0Segment& operator=(FeDiscontinuousO0Segment&&) noexcept =
      default;
  FeDiscontinuousO0Segment() = default;
  ~FeDiscontinuousO0Segment() override = default;
 
  [[nodiscard]] lf::base::RefEl RefEl() const override {
    return lf::base::RefEl::kSegment();
  }
 
  [[nodiscard]] unsigned Degree() const override { return 0; }
 
  [[nodiscard]] size_type NumRefShapeFunctions() const override { return 1; }
 
  [[nodiscard]] size_type NumRefShapeFunctions(dim_t codim) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
  [[nodiscard]] size_type NumRefShapeFunctions(
      dim_t codim, sub_idx_t /*subidx*/) const override {
    LF_ASSERT_MSG(codim <= 2, "Illegal codim" << codim);
    return (codim == 0) ? 1 : 0;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  EvalReferenceShapeFunctions(const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 1,
                  "Reference coordinates must be 1-vectors");
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(n_pts, 1);
    return result;
  }
 
  // clang-format off
  // clang-format on
  [[nodiscard]] Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
  GradientsReferenceShapeFunctions(
      const Eigen::MatrixXd& refcoords) const override {
    LF_ASSERT_MSG(refcoords.rows() == 1,
                  "Reference coordinates must be 1-vectors");
 
    const size_type n_pts(refcoords.cols());
    Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic> result(1, n_pts);
    result.row(0) = Eigen::RowVectorXd::Constant(n_pts, 0);
    return result;
  }
 
  [[nodiscard]] Eigen::MatrixXd EvaluationNodes() const override {
    Eigen::MatrixXd nodes(1, 1);
    nodes << 0.5;
    return nodes;
  }
  [[nodiscard]] size_type NumEvaluationNodes() const override {
    return NumRefShapeFunctions();
  }
};
 
}  // namespace projects::dpg
 
#endif  // PROJECTS_DPG_DISCONTINUOUS_FE_CONSTANT_H