mofem/atom_tests/dg_projection.cpp

Testing Discontinuous Galerkin (DG) projection operators.

Testing Discontinuous Galerkin (DG) projection operatorsThis 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

/**
 * @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>();
 
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getOpDomainRhsPipeline().clear();
 
  auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(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->getOpDomainRhsPipeline().push_back(op_this);
 
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpDGProjectionEvaluation(
      data_ptr, coeffs_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2));
 
  pipeline_mng->getOpDomainRhsPipeline().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));
 
  CHKERR pipeline_mng->loopFiniteElements();
 
  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();
}