prism_polynomial_approximation.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

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

Variable Documentation

approx_order

constexpr int approx_order = 6

staticconstexpr

Definition at line 14 of file prism_polynomial_approximation.cpp.

help

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

static

Examples

prism_polynomial_approximation.cpp.

Definition at line 12 of file prism_polynomial_approximation.cpp.