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

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

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

Definition at line 34 of file hcurl_check_approx_in_2d.cpp.

Function Documentation

◆ main()

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

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_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;
 
    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;
 
 
      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->getOpEvaluationPipeline().push_back(
          new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpMakeHdivFromHcurl());
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
      pipeline_mng->getOpEvaluationPipeline().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->getOpEvaluationPipeline().push_back(
          new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2, base, L2));
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpBaseDerivativesSetHOInvJacobian<SPACE_DIM>(data_l2,
                                                           inv_jac_ptr));
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
                                              base_mass, data_l2, base, HCURL));
 
      // Calculate field values at integration points
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpCalculateHVecVectorField<BASE_DIM>("FIELD1", ptr_values));
      // Gradient
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpCalculateHVecVectorGradient<BASE_DIM, SPACE_DIM>("FIELD1",
                                                                 ptr_grad));
      // Hessian
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpCalculateHdivVectorDivergence<BASE_DIM, SPACE_DIM>(
              "FIELD1", ptr_divergence));
      pipeline_mng->getOpEvaluationPipeline().push_back(
          new OpCalculateHVecVectorHessian<BASE_DIM, SPACE_DIM>("FIELD1",
                                                                ptr_hessian));
 
      pipeline_mng->getOpEvaluationPipeline().push_back(new OpCheckValsDiffVals(
          ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
 
      CHKERR pipeline_mng->loopFiniteElementsEvaluation();
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR assemble_matrices_and_vectors();
    CHKERR solve_problem();
    CHKERR check_solution();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}