#include <MoFEM.hpp>

Classes
struct	ApproxFunction

struct	PrismOpCheck

struct	PrismOpRhs

struct	PrismOpLhs

struct	PrismFE

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

Variables
static char	help [] = "...\n\n"

static constexpr int	approx_order = 6

Function Documentation

◆ main()

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

Examples: prism_polynomial_approximation.cpp.

Definition at line 96 of file prism_polynomial_approximation.cpp.

                                  {
  
   MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
  
   try {
  
     moab::Core mb_instance;
     moab::Interface &moab = mb_instance;
  
     ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
     auto moab_comm_wrap =
         boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
     if (pcomm == NULL)
       pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
  
     std::array<double, 18> one_prism_coords = {0, 0, 0, 1, 0, 0, 0, 1, 0,
                                                0, 0, 1, 1, 0, 1, 0, 1, 1};
     std::array<EntityHandle, 6> one_prism_nodes;
     for (int n = 0; n != 6; ++n)
       CHKERR moab.create_vertex(&one_prism_coords[3 * n], one_prism_nodes[n]);
     EntityHandle one_prism;
     CHKERR moab.create_element(MBPRISM, one_prism_nodes.data(), 6, one_prism);
     Range one_prism_range;
     one_prism_range.insert(one_prism);
     Range one_prism_adj_ents;
     for (int d = 1; d != 3; ++d)
       CHKERR moab.get_adjacencies(one_prism_range, d, true, one_prism_adj_ents,
                                   moab::Interface::UNION);
  
     MoFEM::Core core(moab);
     MoFEM::Interface &m_field = core;
  
     PrismsFromSurfaceInterface *prisms_from_surface_interface;
     CHKERR m_field.getInterface(prisms_from_surface_interface);
     BitRefLevel bit_level0;
     bit_level0.set(0);
     CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
         one_prism_range, bit_level0);
     CHKERR prisms_from_surface_interface->seedPrismsEntities(one_prism_range,
                                                              bit_level0);
  
     // Fields
     CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
     CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "FIELD1");
  
     CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD1", 1);
     CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", approx_order);
     CHKERR m_field.set_field_order(0, MBTRI, "FIELD1", approx_order);
     CHKERR m_field.set_field_order(0, MBQUAD, "FIELD1", approx_order);
     CHKERR m_field.set_field_order(0, MBPRISM, "FIELD1", approx_order);
     CHKERR m_field.build_fields();
  
     // FE
     CHKERR m_field.add_finite_element("PRISM");
  
     // Define rows/cols and element data
     CHKERR m_field.modify_finite_element_add_field_row("PRISM", "FIELD1");
     CHKERR m_field.modify_finite_element_add_field_col("PRISM", "FIELD1");
     CHKERR m_field.modify_finite_element_add_field_data("PRISM", "FIELD1");
     CHKERR m_field.add_ents_to_finite_element_by_type(0, MBPRISM, "PRISM");
  
     // build finite elemnts
     CHKERR m_field.build_finite_elements();
     // //build adjacencies
     CHKERR m_field.build_adjacencies(bit_level0);
  
     // Problem
     CHKERR m_field.add_problem("TEST_PROBLEM");
  
     // set finite elements for problem
     CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "PRISM");
     // set refinement level for problem
     CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
  
     // build problem
     ProblemsManager *prb_mng_ptr;
     CHKERR m_field.getInterface(prb_mng_ptr);
     CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
     // partition
     CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
     CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
     // what are ghost nodes, see Petsc Manual
     CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
  
     // Create matrices
     SmartPetscObj<Mat> A;
     CHKERR m_field.getInterface<MatrixManager>()
         ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("TEST_PROBLEM", A);
     SmartPetscObj<Vec> F;
     CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
                                                               ROW, F);
     SmartPetscObj<Vec> D;
     CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
                                                               COL, D);
  
     auto assemble_matrices_and_vectors = [&]() {
       MoFEMFunctionBegin;
       PrismFE fe(m_field);
       MatrixDouble inv_jac;
       fe.getOpPtrVector().push_back(
           new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
       fe.getOpPtrVector().push_back(
           new MoFEM::OpCalculateInvJacForFatPrism(inv_jac));
       fe.getOpPtrVector().push_back(
           new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
       fe.getOpPtrVector().push_back(new PrismOpRhs(F));
       fe.getOpPtrVector().push_back(new PrismOpLhs(A));
       CHKERR VecZeroEntries(F);
       CHKERR MatZeroEntries(A);
       CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe);
       CHKERR VecAssemblyBegin(F);
       CHKERR VecAssemblyEnd(F);
       CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
       CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
       MoFEMFunctionReturn(0);
     };
  
     auto solve_problem = [&] {
       MoFEMFunctionBegin;
       auto solver = createKSP(PETSC_COMM_WORLD);
       CHKERR KSPSetOperators(solver, A, A);
       CHKERR KSPSetFromOptions(solver);
       CHKERR KSPSetUp(solver);
       CHKERR KSPSolve(solver, F, D);
       CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
       CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
       CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
           "TEST_PROBLEM", COL, D, INSERT_VALUES, SCATTER_REVERSE);
       MoFEMFunctionReturn(0);
     };
  
     auto check_solution = [&] {
       MoFEMFunctionBegin;
       PrismFE fe(m_field);
       boost::shared_ptr<VectorDouble> field_vals_ptr(new VectorDouble());
       boost::shared_ptr<MatrixDouble> diff_field_vals_ptr(new MatrixDouble());
       MatrixDouble inv_jac;
       fe.getOpPtrVector().push_back(
           new OpCalculateInvJacForFatPrism(inv_jac));
       fe.getOpPtrVector().push_back(
           new OpSetInvJacH1ForFatPrism(inv_jac));
       fe.getOpPtrVector().push_back(
           new OpCalculateScalarFieldValues("FIELD1", field_vals_ptr));
       fe.getOpPtrVector().push_back(
           new OpCalculateScalarFieldGradient<3>("FIELD1", diff_field_vals_ptr));
       fe.getOpPtrVector().push_back(
           new PrismOpCheck(field_vals_ptr, diff_field_vals_ptr));
       CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe);
       MoFEMFunctionReturn(0);
     };
  
     CHKERR assemble_matrices_and_vectors();
     CHKERR solve_problem();
     CHKERR check_solution();
   }
   CATCH_ERRORS;
  
   MoFEM::Core::Finalize();
   return 0;
 }