dg_projection.cpp

Go to the documentation of this file.

/**
 * @file dg_projection.cpp
 * @example mofem/atom_tests/dg_projection.cpp 
 * 
 * @brief Testing Discontinuous Galerkin (DG) projection operators
 *
 * This test validates the accuracy of DG projection algorithms used to
 * represent functions on finite element meshes. The projection process
 * involves solving local least-squares problems to find the best
 * approximation of a given function in the finite element space.
 *
 * Mathematical background:
 * - DG projection finds coefficients c_i such that ∑c_i φ_i minimizes
 *   ||u - ∑c_i φ_i||_{L2} where u is the target function and φ_i are basis functions
 * - This leads to the normal equation: M c = f, where M is the mass matrix
 *   and f is the load vector from projecting the target function
 * - The test verifies that the projection accurately represents a polynomial
 *   function that should be exactly representable in the chosen finite element space
 *
 * Test workflow:
 * 1. Setup finite element space with polynomial basis functions
 * 2. Assemble mass matrix and load vector for the projection problem
 * 3. Solve the projection system to get approximation coefficients
 * 4. Evaluate the projected function and compare with analytical values
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
static char help[] = "DG Projection Test - validates discontinuous Galerkin "
                     "projection accuracy\n\n";
 
constexpr char FIELD_NAME[] = "U";
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
constexpr int order = 2;
 
template <int DIM> struct ElementsAndOps {};
 
template <> struct ElementsAndOps<2> {
  using DomainEle = PipelineManager::FaceEle;
  using DomainEleOp = DomainEle::UserDataOperator;
};
 
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using EntData = EntitiesFieldData::EntData;
 
auto fun = [](const double x, const double y, const double z) {
  return x + y + x * x + y * y;
};
 
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpMass<BASE_DIM, FIELD_DIM>;
 
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
 
struct AtomTest {
 
  AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
  Simple *simpleInterface;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode checkResults();
 
  struct CommonData {
    boost::shared_ptr<MatrixDouble> invJacPtr;
    boost::shared_ptr<VectorDouble> approxVals;
    boost::shared_ptr<MatrixDouble> approxGradVals;
    boost::shared_ptr<MatrixDouble> approxHessianVals;
    SmartPetscObj<Vec> L2Vec;
  };
 
  struct OpError;
};
 
struct AtomTest::OpError : public DomainEleOp {
  boost::shared_ptr<CommonData> commonDataPtr;
 
  OpError(boost::shared_ptr<MatrixDouble> data_ptr)
      : DomainEleOp(NOSPACE, OPSPACE), dataPtr(data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
 
    const int nb_integration_pts = getGaussPts().size2();
 
    auto t_val = getFTensor1FromMat<1>(*(dataPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      double projected_value = t_val(0);
      double analytical_value = fun(t_coords(0), t_coords(1), t_coords(2));
      double error = projected_value - analytical_value;
 
      constexpr double eps = 1e-8;
      if (std::abs(error) > eps) {
        MOFEM_LOG("SELF", Sev::error)
            << "Projection error too large: " << error << " at point ("
            << t_coords(0) << ", " << t_coords(1) << ")"
            << " projected=" << projected_value
            << " analytical=" << analytical_value;
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "DG projection failed accuracy test");
      }
 
      ++t_val;
      ++t_coords;
    }
 
    MOFEM_LOG("SELF", Sev::noisy) << "DG projection accuracy validation passed";
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dataPtr;
};
 
//! [Run programme]
MoFEMErrorCode AtomTest::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run programme]
 
//! [Read mesh]
MoFEMErrorCode AtomTest::readMesh() {
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(simpleInterface);
 
  CHKERR simpleInterface->getOptions();
 
  CHKERR simpleInterface->loadFile();
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode AtomTest::setupProblem() {
  MoFEMFunctionBegin;
 
  CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
                                         AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
 
  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
 
  CHKERR simpleInterface->setUp();
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Push operators to pipeline]
MoFEMErrorCode AtomTest::assembleSystem() {
  MoFEMFunctionBegin;
 
  auto rule = [](int, int, int p) -> int { return 2 * p; };
 
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
 
  auto beta = [](const double, const double, const double) { return 1; };
 
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
 
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource(FIELD_NAME, fun));
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode AtomTest::solveSystem() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  MOFEM_LOG("WORLD", Sev::inform) << "Solving DG projection system";
 
  auto solver = pipeline_mng->createKSP();
  CHKERR KSPSetFromOptions(solver);
  CHKERR KSPSetUp(solver);
 
  auto dm = simpleInterface->getDM();
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  CHKERR KSPSolve(solver, F, D);
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Check results]
MoFEMErrorCode AtomTest::checkResults() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
  CHKERR pipeline_mng->setEvaluationIntegrationRule(rule);
 
  auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
 
  auto mass_ptr = boost::make_shared<MatrixDouble>();
  auto coeffs_ptr = boost::make_shared<MatrixDouble>();
  auto data_ptr = boost::make_shared<MatrixDouble>();
 
  auto op_this =
      new OpLoopThis<DomainEle>(mField, simple->getDomainFEName(), Sev::noisy);
  pipeline_mng->getOpEvaluationPipeline().push_back(op_this);
 
  pipeline_mng->getOpEvaluationPipeline().push_back(
      new OpDGProjectionEvaluation(data_ptr, coeffs_ptr, entity_data_l2,
                                   AINSWORTH_LEGENDRE_BASE, L2));
 
  pipeline_mng->getOpEvaluationPipeline().push_back(new OpError(data_ptr));
 
  auto fe_physics_ptr = op_this->getThisFEPtr();
  fe_physics_ptr->getRuleHook = [](int, int, int p) { return 2 * p; };
 
  fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
      order, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2));
 
  fe_physics_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME, data_ptr));
 
  fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
      data_ptr, coeffs_ptr, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE,
      L2, Sev::noisy));
 
  CHKERR pipeline_mng->loopFiniteElementsEvaluation();
 
  MoFEMFunctionReturn(0);
}
//! [Check results]
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, NULL, help);
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc]
 
    //! [Create MoAB]
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
    //! [Create MoFEM]
 
    //! [Execute DG Projection Test]
    AtomTest ex(m_field);
    CHKERR ex.runProblem();
    //! [Execute DG Projection Test]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}