#include <MoFEM.hpp>

Classes
struct	ApproxFunctions

struct	OpCheckValsDiffVals

Typedefs
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 = 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

◆ OpDomainMass

using OpDomainMass

Initial value:

FormsIntegrators<FaceEleOp>::Assembly<PETSC>::BiLinearForm<

GAUSS>::OpMass<BASE_DIM, BASE_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 27 of file hcurl_check_approx_in_2d.cpp.

◆ OpDomainSource

using OpDomainSource

Initial value:

FormsIntegrators<FaceEleOp>::Assembly<PETSC>::LinearForm<

GAUSS>::OpSource<BASE_DIM, BASE_DIM>

LinearForm

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