#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. More...

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

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

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.

Examples: hdiv_check_approx_in_3d.cpp.

Definition at line 26 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.

Examples: hdiv_check_approx_in_3d.cpp.

Definition at line 33 of file hdiv_check_approx_in_3d.cpp.

Function Documentation

◆ main()

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

Examples: hdiv_check_approx_in_3d.cpp.

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