v0.9.0
Classes | Typedefs | Functions | Variables
quad_polynomial_approximation.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ApproxFunction
 
struct  QuadOpCheck
 
struct  QuadOpRhs
 
struct  QuadOpLhs
 

Typedefs

using Ele = FaceElementForcesAndSourcesCore
 
using OpEle = FaceElementForcesAndSourcesCore::UserDataOperator
 
using EntData = DataForcesAndSourcesCore::EntData
 

Functions

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

Variables

static char help [] = "...\n\n"
 
static constexpr int approx_order = 5
 

Typedef Documentation

◆ Ele

◆ EntData

◆ OpEle

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
quad_polynomial_approximation.cpp.

Definition at line 95 of file quad_polynomial_approximation.cpp.

95  {
96 
97  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
98 
99  try {
100 
101  moab::Core mb_instance;
102  moab::Interface &moab = mb_instance;
103 
104  std::array<double, 12> one_quad_coords = {0, 0, 0,
105 
106  1, 0, 0,
107 
108  1, 1, 0,
109 
110  0, 1, 0};
111  std::array<EntityHandle, 4> one_quad_nodes;
112  for (int n = 0; n != 4; ++n)
113  CHKERR moab.create_vertex(&one_quad_coords[3 * n], one_quad_nodes[n]);
114  EntityHandle one_quad;
115  CHKERR moab.create_element(MBQUAD, one_quad_nodes.data(), 4, one_quad);
116  Range one_quad_range;
117  one_quad_range.insert(one_quad);
118  Range one_quad_adj_ents;
119  CHKERR moab.get_adjacencies(one_quad_range, 1, true, one_quad_adj_ents,
120  moab::Interface::UNION);
121 
122  MoFEM::Core core(moab);
123  MoFEM::Interface &m_field = core;
124 
125  BitRefLevel bit_level0 = BitRefLevel().set(0);
126  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
127  0, 2, bit_level0);
128 
129  // Fields
130  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
131  CHKERR m_field.add_ents_to_field_by_type(0, MBQUAD, "FIELD1");
132 
133  CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD1", 1);
134  CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", approx_order);
135  CHKERR m_field.set_field_order(0, MBQUAD, "FIELD1", approx_order);
136  CHKERR m_field.build_fields();
137 
138  // FE
139  CHKERR m_field.add_finite_element("QUAD");
140 
141  // Define rows/cols and element data
142  CHKERR m_field.modify_finite_element_add_field_row("QUAD", "FIELD1");
143  CHKERR m_field.modify_finite_element_add_field_col("QUAD", "FIELD1");
144  CHKERR m_field.modify_finite_element_add_field_data("QUAD", "FIELD1");
145  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBQUAD, "QUAD");
146 
147  // build finite elemnts
148  CHKERR m_field.build_finite_elements();
149  // //build adjacencies
150  CHKERR m_field.build_adjacencies(bit_level0);
151 
152  // Problem
153  CHKERR m_field.add_problem("TEST_PROBLEM");
154 
155  // set finite elements for problem
156  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "QUAD");
157  // set refinement level for problem
158  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
159 
160  // build problem
161  ProblemsManager *prb_mng_ptr;
162  CHKERR m_field.getInterface(prb_mng_ptr);
163  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
164  // partition
165  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
166  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
167  // what are ghost nodes, see Petsc Manual
168  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
169 
170  // Create matrices
173  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("TEST_PROBLEM", A);
175  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
176  ROW, F);
178  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
179  COL, D);
180 
181  auto rule = [&](int, int, int p) { return 2 * (p + 1); };
182 
183  auto assemble_matrices_and_vectors = [&]() {
185  Ele fe(m_field);
186  fe.getRuleHook = rule;
187  MatrixDouble inv_jac;
188  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
189  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
190  fe.getOpPtrVector().push_back(new QuadOpRhs(F));
191  fe.getOpPtrVector().push_back(new QuadOpLhs(A));
192  CHKERR VecZeroEntries(F);
193  CHKERR MatZeroEntries(A);
194  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
195  CHKERR VecAssemblyBegin(F);
196  CHKERR VecAssemblyEnd(F);
197  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
198  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
200  };
201 
202  auto solve_problem = [&] {
204  auto solver = createKSP(PETSC_COMM_WORLD);
205  CHKERR KSPSetOperators(solver, A, A);
206  CHKERR KSPSetFromOptions(solver);
207  CHKERR KSPSetUp(solver);
208  CHKERR KSPSolve(solver, F, D);
209  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
210  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
211  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
212  "TEST_PROBLEM", COL, D, INSERT_VALUES, SCATTER_REVERSE);
214  };
215 
216  auto check_solution = [&] {
218  Ele fe(m_field);
219  fe.getRuleHook = rule;
220  boost::shared_ptr<VectorDouble> field_vals_ptr(new VectorDouble());
221  boost::shared_ptr<MatrixDouble> diff_field_vals_ptr(new MatrixDouble());
222  MatrixDouble inv_jac;
223  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
224  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
225  fe.getOpPtrVector().push_back(
226  new OpCalculateScalarFieldValues("FIELD1", field_vals_ptr));
227  fe.getOpPtrVector().push_back(
228  new OpCalculateScalarFieldGradient<2>("FIELD1", diff_field_vals_ptr));
229  fe.getOpPtrVector().push_back(
230  new QuadOpCheck(field_vals_ptr, diff_field_vals_ptr));
231  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
233  };
234 
235  CHKERR assemble_matrices_and_vectors();
237  CHKERR check_solution();
238  }
239  CATCH_ERRORS;
240 
242  return 0;
243 }
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static char help[]
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
Core (interface) class.
Definition: Core.hpp:50
static constexpr int approx_order
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
PetscErrorCode solve_problem(MoFEM::Interface &m_field, const string &problem_name, const string &fe_name, const string &re_field, const string &im_field, InsertMode mode, FUNEVAL &fun_evaluator, PetscBool is_partitioned)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Calculate inverse of jacobian for face element.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:150
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:36
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:600
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
auto createKSP
Definition: AuxPETSc.hpp:270
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elementsFunction which partition finite elements based on dofs partitioning....
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:437
continuous field
Definition: definitions.h:175
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
Get value at integration points for scalar field.
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
Matrix manager is used to build and partition problems.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

Variable Documentation

◆ approx_order

constexpr int approx_order = 5
static

◆ help

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