#include <MoFEM.hpp>

Classes
struct	ApproxFunctions

struct	OpCheckValsDiffVals

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

using	EleOp = Ele::UserDataOperator

using	OpDomainMass
	OPerator to integrate mass matrix for least square approximation.

using	OpDomainSource
	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: 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

Initial value:

FormsIntegrators<EleOp>::Assembly<PETSC>::BiLinearForm<

GAUSS>::OpMass<BASE_DIM, SPACE_DIM>

BiLinearForm

MoFEM::GAUSS

@ GAUSS

Definition FormsIntegrators.hpp:136

MoFEM::PETSC

@ PETSC

Definition FormsIntegrators.hpp:105

OpMass

FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass

[Only used with Hooke equation (linear material model)]

Definition seepage.cpp:56

MoFEM::FormsIntegrators

Integrator forms.

Definition FormsIntegrators.hpp:306

OPerator to integrate mass matrix for least square approximation.

Definition at line 25 of file hdiv_check_approx_in_3d.cpp.

◆ OpDomainSource

using OpDomainSource

Initial value:

FormsIntegrators<EleOp>::Assembly<PETSC>::LinearForm<

GAUSS>::OpSource<BASE_DIM, SPACE_DIM>

LinearForm

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_NULL, 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();
}