#include <MoFEM.hpp>

Classes
struct	ApproxFunctions

struct	OpCheckValsDiffVals

Typedefs
using	FaceEle = MoFEM::FaceElementForcesAndSourcesCore

using	FaceEleOp = FaceEle::UserDataOperator

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

using	OpDomainSource = FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_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 = 2

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

◆ FaceEle

using FaceEle = MoFEM::FaceElementForcesAndSourcesCore

Definition at line 17 of file hcurl_check_approx_in_2d.cpp.

◆ FaceEleOp

using FaceEleOp = FaceEle::UserDataOperator

Definition at line 18 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainMass

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

OPerator to integrate mass matrix for least square approximation.

Examples: hcurl_check_approx_in_2d.cpp.

Definition at line 28 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainSource

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

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

Examples: hcurl_check_approx_in_2d.cpp.

Definition at line 35 of file hcurl_check_approx_in_2d.cpp.

Function Documentation

◆ main()

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

Examples: hcurl_check_approx_in_2d.cpp.

Definition at line 236 of file hcurl_check_approx_in_2d.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_interface = m_field.getInterface<Simple>();
     PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
     CHKERR simple_interface->getOptions();
     CHKERR simple_interface->loadFile("", "rectangle_tri.h5m");
  
     // 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;
  
     CHKERR simple_interface->addDomainField("FIELD1", HCURL, base, 1);
     constexpr int order = 5;
     CHKERR simple_interface->setFieldOrder("FIELD1", order);
     CHKERR simple_interface->setUp();
     auto dm = simple_interface->getDM();
  
     MatrixDouble vals, diff_vals;
  
     auto assemble_matrices_and_vectors = [&]() {
       MoFEMFunctionBegin;
       auto jac_ptr = boost::make_shared<MatrixDouble>();
       auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
       auto det_ptr = boost::make_shared<VectorDouble>();
  
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpMakeHdivFromHcurl());
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpDomainSource("FIELD1", ApproxFunctions::fUn));
  
       pipeline_mng->getOpDomainLhsPipeline().push_back(
           new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
       pipeline_mng->getOpDomainLhsPipeline().push_back(
           new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
       pipeline_mng->getOpDomainLhsPipeline().push_back(
           new OpMakeHdivFromHcurl());
       pipeline_mng->getOpDomainLhsPipeline().push_back(
           new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
       pipeline_mng->getOpDomainLhsPipeline().push_back(
  
           new OpDomainMass("FIELD1", "FIELD1",
                            [](double, double, double) { return 1; })
  
       );
  
       auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
       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_divergence = boost::make_shared<VectorDouble>();
       auto ptr_grad = boost::make_shared<MatrixDouble>();
       auto ptr_hessian = boost::make_shared<MatrixDouble>();
  
       auto jac_ptr = boost::make_shared<MatrixDouble>();
       auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
       auto det_ptr = boost::make_shared<VectorDouble>();
  
       // Change H-curl to H-div in 2D, and apply Piola transform
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpMakeHdivFromHcurl());
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpSetInvJacHcurlFace(inv_jac_ptr));
  
       // Evaluate base function second derivative
       auto base_mass = boost::make_shared<MatrixDouble>();
       auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
  
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2, base, L2));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpBaseDerivativesSetHOInvJacobian<SPACE_DIM>(data_l2,
                                                            inv_jac_ptr));
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
                                               base_mass, data_l2, base, HCURL));
  
       // 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));
       // Hessian
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpCalculateHdivVectorDivergence<BASE_DIM, SPACE_DIM>(
               "FIELD1", ptr_divergence)); 
       pipeline_mng->getOpDomainRhsPipeline().push_back(
           new OpCalculateHVecVectorHessian<BASE_DIM, SPACE_DIM>("FIELD1",
                                                                 ptr_hessian));
  
       pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
           ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
  
       CHKERR pipeline_mng->loopFiniteElements();
  
       MoFEMFunctionReturn(0);
     };
  
     CHKERR assemble_matrices_and_vectors();
     CHKERR solve_problem();
     CHKERR check_solution();
   }
   CATCH_ERRORS;
  
   CHKERR MoFEM::Core::Finalize();
 }