v0.9.1
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 
112  std::array<EntityHandle, 4> one_quad_nodes;
113  for (int n = 0; n != 4; ++n)
114  CHKERR moab.create_vertex(&one_quad_coords[3 * n], one_quad_nodes[n]);
115  EntityHandle one_quad;
116  CHKERR moab.create_element(MBQUAD, one_quad_nodes.data(), 4, one_quad);
117  Range one_quad_range;
118  one_quad_range.insert(one_quad);
119  Range one_quad_adj_ents;
120  CHKERR moab.get_adjacencies(one_quad_range, 1, true, one_quad_adj_ents,
121  moab::Interface::UNION);
122 
123  MoFEM::Core core(moab);
124  MoFEM::Interface &m_field = core;
125 
126  BitRefLevel bit_level0 = BitRefLevel().set(0);
127  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
128  0, 2, bit_level0);
129 
130  // Fields
131  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
132  CHKERR m_field.add_ents_to_field_by_type(0, MBQUAD, "FIELD1");
133 
134  CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD1", 1);
135  CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", approx_order);
136  CHKERR m_field.set_field_order(0, MBQUAD, "FIELD1", approx_order);
137  CHKERR m_field.build_fields();
138 
139  // FE
140  CHKERR m_field.add_finite_element("QUAD");
141 
142  // Define rows/cols and element data
143  CHKERR m_field.modify_finite_element_add_field_row("QUAD", "FIELD1");
144  CHKERR m_field.modify_finite_element_add_field_col("QUAD", "FIELD1");
145  CHKERR m_field.modify_finite_element_add_field_data("QUAD", "FIELD1");
146  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBQUAD, "QUAD");
147 
148  // build finite elemnts
149  CHKERR m_field.build_finite_elements();
150  // //build adjacencies
151  CHKERR m_field.build_adjacencies(bit_level0);
152 
153  // Problem
154  CHKERR m_field.add_problem("TEST_PROBLEM");
155 
156  // set finite elements for problem
157  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "QUAD");
158  // set refinement level for problem
159  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
160 
161  // build problem
162  ProblemsManager *prb_mng_ptr;
163  CHKERR m_field.getInterface(prb_mng_ptr);
164  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
165  // partition
166  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
167  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
168  // what are ghost nodes, see Petsc Manual
169  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
170 
171  // Create matrices
174  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("TEST_PROBLEM", A);
176  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
177  ROW, F);
179  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
180  COL, D);
181 
182  auto rule = [&](int, int, int p) { return 2 * (p + 1); };
183 
184  auto assemble_matrices_and_vectors = [&]() {
186  Ele fe(m_field);
187  fe.getRuleHook = rule;
188  MatrixDouble inv_jac;
189  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
190  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
191  fe.getOpPtrVector().push_back(new QuadOpRhs(F));
192  fe.getOpPtrVector().push_back(new QuadOpLhs(A));
193  CHKERR VecZeroEntries(F);
194  CHKERR MatZeroEntries(A);
195  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
196  CHKERR VecAssemblyBegin(F);
197  CHKERR VecAssemblyEnd(F);
198  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
199  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
201  };
202 
203  auto solve_problem = [&] {
205  auto solver = createKSP(PETSC_COMM_WORLD);
206  CHKERR KSPSetOperators(solver, A, A);
207  CHKERR KSPSetFromOptions(solver);
208  CHKERR KSPSetUp(solver);
209  CHKERR KSPSolve(solver, F, D);
210  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
211  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
212  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
213  "TEST_PROBLEM", COL, D, INSERT_VALUES, SCATTER_REVERSE);
215  };
216 
217  auto check_solution = [&] {
219  Ele fe(m_field);
220  fe.getRuleHook = rule;
221  boost::shared_ptr<VectorDouble> field_vals_ptr(new VectorDouble());
222  boost::shared_ptr<MatrixDouble> diff_field_vals_ptr(new MatrixDouble());
223  MatrixDouble inv_jac;
224  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
225  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
226  fe.getOpPtrVector().push_back(
227  new OpCalculateScalarFieldValues("FIELD1", field_vals_ptr));
228  fe.getOpPtrVector().push_back(
229  new OpCalculateScalarFieldGradient<2>("FIELD1", diff_field_vals_ptr));
230  fe.getOpPtrVector().push_back(
231  new QuadOpCheck(field_vals_ptr, diff_field_vals_ptr));
232  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
234  };
235 
236  CHKERR assemble_matrices_and_vectors();
238  CHKERR check_solution();
239  }
240  CATCH_ERRORS;
241 
243  return 0;
244 }
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
Deprecated interface functions.
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:75
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
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:151
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:601
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:50
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
auto createKSP
Definition: AuxPETSc.hpp:289
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.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
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:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
continuous field
Definition: definitions.h:176
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
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