#include <MoFEM.hpp>

Classes
struct	ApproxFunctions

struct	OpCheckValsDiffVals

Typedefs
using	Ele = MoFEM::PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle

using	EleOp = Ele::UserDataOperator

using	OpDomainMass = FormsIntegrators< EleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, SPACE_DIM >
	OPerator to integrate mass matrix for least square approximation.

using	OpDomainSource = FormsIntegrators< EleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, SPACE_DIM >
	Operator to integrate the right hand side matrix for the problem.

Functions
int	main (int argc, char *argv[])

Variables
static char	help [] = "...\n\n"

constexpr int	BASE_DIM = 3

constexpr int	SPACE_DIM = 3

constexpr double	a0 = 0.0

constexpr double	a1 = 2.0

constexpr double	a2 = -15.0 * a0

constexpr double	a3 = -20.0 / 6 * a1

constexpr double	a4 = 15.0 * a0

constexpr double	a5 = a1

constexpr double	a6 = -a0

Typedef Documentation

◆ Ele

using Ele = MoFEM::PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle

Examples: mofem/tutorials/cor-9/reaction_diffusion.cpp.

Definition at line 18 of file hdiv_check_approx_in_3d.cpp.

◆ EleOp

using EleOp = Ele::UserDataOperator

Definition at line 19 of file hdiv_check_approx_in_3d.cpp.

◆ OpDomainMass

using OpDomainMass = FormsIntegrators<EleOp>::Assembly<PETSC>::BiLinearForm< GAUSS>::OpMass<BASE_DIM, SPACE_DIM>

OPerator to integrate mass matrix for least square approximation.

Definition at line 25 of file hdiv_check_approx_in_3d.cpp.

◆ OpDomainSource

using OpDomainSource = FormsIntegrators<EleOp>::Assembly<PETSC>::LinearForm< GAUSS>::OpSource<BASE_DIM, SPACE_DIM>

Operator to integrate the right hand side matrix for the problem.

Definition at line 32 of file hdiv_check_approx_in_3d.cpp.

Function Documentation

◆ main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 129 of file hdiv_check_approx_in_3d.cpp.

                                 {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    // Create MoFEM instance
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    Simple *simple = m_field.getInterface<Simple>();
    PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
    CHKERR simple->getOptions();
    CHKERR simple->loadFile();
 
    // Declare elements
    enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
    const char *list_bases[] = {"ainsworth", "demkowicz"};
    PetscBool flg;
    PetscInt choice_base_value = AINSWORTH;
    CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                                LASBASETOP, &choice_base_value, &flg);
 
    if (flg != PETSC_TRUE)
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
    FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
    if (choice_base_value == AINSWORTH)
      base = AINSWORTH_LEGENDRE_BASE;
    else if (choice_base_value == DEMKOWICZ)
      base = DEMKOWICZ_JACOBI_BASE;
 
    AinsworthOrderHooks::broken_nbfacetri_edge_hdiv = [](int p) { return p; };
    AinsworthOrderHooks::broken_nbfacetri_face_hdiv = [](int p) { return p; };
    AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv = [](int p) { return p; };
    AinsworthOrderHooks::broken_nbvolumetet_face_hdiv = [](int p) { return p; };
    AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv = [](int p) {
      return p;
    };
 
    CHKERR simple->addDomainField("FIELD1", HDIV, base, 1);
    constexpr int order = 4;
    CHKERR simple->setFieldOrder("FIELD1", order);
    CHKERR simple->setUp();
    auto dm = simple->getDM();
 
    MatrixDouble vals, diff_vals;
 
    auto assemble_matrices_and_vectors = [&]() {
      MoFEMFunctionBegin;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          pipeline_mng->getOpDomainRhsPipeline(), {HDIV});
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpDomainSource("FIELD1", ApproxFunctions::fUn));
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          pipeline_mng->getOpDomainLhsPipeline(), {HDIV});
      pipeline_mng->getOpDomainLhsPipeline().push_back(
 
          new OpDomainMass("FIELD1", "FIELD1",
                           [](double x, double, double) { return 1; })
 
      );
 
      auto integration_rule = [](int, int, int p_data) {
        return 2 * p_data + 1;
      };
      CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
      CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
 
      MoFEMFunctionReturn(0);
    };
 
    auto solve_problem = [&] {
      MoFEMFunctionBegin;
      auto solver = pipeline_mng->createKSP();
      CHKERR KSPSetFromOptions(solver);
      CHKERR KSPSetUp(solver);
 
      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);
    };
 
    auto check_solution = [&] {
      MoFEMFunctionBegin;
 
      pipeline_mng->getOpDomainLhsPipeline().clear();
      pipeline_mng->getOpDomainRhsPipeline().clear();
 
      auto ptr_values = boost::make_shared<MatrixDouble>();
      auto ptr_grad = boost::make_shared<MatrixDouble>();
 
      // Change H-curl to H-div in 2D, and apply Piola transform
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          pipeline_mng->getOpDomainRhsPipeline(), {HDIV});
 
      // Calculate field values at integration points
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpCalculateHVecVectorField<BASE_DIM>("FIELD1", ptr_values));
      // Gradient
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpCalculateHVecVectorGradient<BASE_DIM, SPACE_DIM>("FIELD1",
                                                                 ptr_grad));
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpCheckValsDiffVals(ptr_values, ptr_grad));
 
      CHKERR pipeline_mng->loopFiniteElements();
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR assemble_matrices_and_vectors();
    CHKERR solve_problem();
    CHKERR check_solution();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}