v0.11.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 = 6
 

Typedef Documentation

◆ Ele

Examples
reaction_diffusion.cpp.

Definition at line 25 of file quad_polynomial_approximation.cpp.

◆ EntData

◆ OpEle

Examples
reaction_diffusion.cpp.

Definition at line 26 of file quad_polynomial_approximation.cpp.

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  // Declare elements
102  enum bases {
103  AINSWORTH,
104  AINSWORTH_LOBATTO,
105  DEMKOWICZ,
106  BERNSTEIN,
107  LASBASETOP
108  };
109  const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
110  "bernstein"};
111  PetscBool flg;
112  PetscInt choice_base_value = AINSWORTH;
113  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
114  LASBASETOP, &choice_base_value, &flg);
115 
116  if (flg != PETSC_TRUE)
117  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
119  if (choice_base_value == AINSWORTH)
121  if (choice_base_value == AINSWORTH_LOBATTO)
122  base = AINSWORTH_LOBATTO_BASE;
123  else if (choice_base_value == DEMKOWICZ)
124  base = DEMKOWICZ_JACOBI_BASE;
125  else if (choice_base_value == BERNSTEIN)
127 
128  enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
129  const char *list_spaces[] = {"h1", "l2"};
130  PetscInt choice_space_value = H1SPACE;
131  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
132  LASBASETSPACE, &choice_space_value, &flg);
133  if (flg != PETSC_TRUE)
134  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
135  FieldSpace space = H1;
136  if (choice_space_value == H1SPACE)
137  space = H1;
138  else if (choice_space_value == L2SPACE)
139  space = L2;
140 
141  moab::Core mb_instance;
142  moab::Interface &moab = mb_instance;
143 
144  std::array<double, 12> one_quad_coords = {0, 0, 0,
145 
146  2, 0, 0,
147 
148  1, 1, 0,
149 
150  0, 1, 0};
151 
152  std::array<EntityHandle, 4> one_quad_nodes;
153  for (int n = 0; n != 4; ++n)
154  CHKERR moab.create_vertex(&one_quad_coords[3 * n], one_quad_nodes[n]);
155  EntityHandle one_quad;
156  CHKERR moab.create_element(MBQUAD, one_quad_nodes.data(), 4, one_quad);
157  Range one_quad_range;
158  one_quad_range.insert(one_quad);
159  Range one_quad_adj_ents;
160  CHKERR moab.get_adjacencies(one_quad_range, 1, true, one_quad_adj_ents,
161  moab::Interface::UNION);
162 
163  MoFEM::Core core(moab);
164  MoFEM::Interface &m_field = core;
165 
166  BitRefLevel bit_level0 = BitRefLevel().set(0);
167  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
168  0, 2, bit_level0);
169 
170  // Fields
171  CHKERR m_field.add_field("FIELD1", space, base, 1);
172  CHKERR m_field.add_ents_to_field_by_type(0, MBQUAD, "FIELD1");
173 
174  CHKERR m_field.set_field_order(0, MBVERTEX, "FIELD1", 1);
175  CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", approx_order + 1);
176  CHKERR m_field.set_field_order(0, MBQUAD, "FIELD1", approx_order + 1);
177  CHKERR m_field.build_fields();
178 
179  // FE
180  CHKERR m_field.add_finite_element("QUAD");
181 
182  // Define rows/cols and element data
183  CHKERR m_field.modify_finite_element_add_field_row("QUAD", "FIELD1");
184  CHKERR m_field.modify_finite_element_add_field_col("QUAD", "FIELD1");
185  CHKERR m_field.modify_finite_element_add_field_data("QUAD", "FIELD1");
186  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBQUAD, "QUAD");
187 
188  // build finite elemnts
189  CHKERR m_field.build_finite_elements();
190  // //build adjacencies
191  CHKERR m_field.build_adjacencies(bit_level0);
192 
193  // Problem
194  CHKERR m_field.add_problem("TEST_PROBLEM");
195 
196  // set finite elements for problem
197  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "QUAD");
198  // set refinement level for problem
199  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
200 
201  // build problem
202  ProblemsManager *prb_mng_ptr;
203  CHKERR m_field.getInterface(prb_mng_ptr);
204  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
205  // partition
206  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
207  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
208  // what are ghost nodes, see Petsc Manual
209  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
210 
211  // Create matrices
214  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("TEST_PROBLEM", A);
216  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
217  ROW, F);
219  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
220  COL, D);
221 
222  auto rule = [&](int, int, int p) { return 2 * (p + 1); };
223 
224  auto assemble_matrices_and_vectors = [&]() {
226  Ele fe(m_field);
227  fe.getRuleHook = rule;
228  MatrixDouble inv_jac;
229  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
230  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
231  fe.getOpPtrVector().push_back(new OpSetInvJacL2ForFace(inv_jac));
232  fe.getOpPtrVector().push_back(new OpMakeHighOrderGeometryWeightsOnFace());
233  fe.getOpPtrVector().push_back(new QuadOpRhs(F));
234  fe.getOpPtrVector().push_back(new QuadOpLhs(A));
235  CHKERR VecZeroEntries(F);
236  CHKERR MatZeroEntries(A);
237  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
238  CHKERR VecAssemblyBegin(F);
239  CHKERR VecAssemblyEnd(F);
240  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
241  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
243  };
244 
245  auto solve_problem = [&] {
247  auto solver = createKSP(PETSC_COMM_WORLD);
248  CHKERR KSPSetOperators(solver, A, A);
249  CHKERR KSPSetFromOptions(solver);
250  CHKERR KSPSetUp(solver);
251  CHKERR KSPSolve(solver, F, D);
252  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
253  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
254  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
255  "TEST_PROBLEM", COL, D, INSERT_VALUES, SCATTER_REVERSE);
257  };
258 
259  auto check_solution = [&] {
261  Ele fe(m_field);
262  fe.getRuleHook = rule;
263  boost::shared_ptr<VectorDouble> field_vals_ptr(new VectorDouble());
264  boost::shared_ptr<MatrixDouble> diff_field_vals_ptr(new MatrixDouble());
265  MatrixDouble inv_jac;
266  fe.getOpPtrVector().push_back(
267  new OpCalculateScalarFieldValues("FIELD1", field_vals_ptr));
268  fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
269  fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
270  fe.getOpPtrVector().push_back(new OpSetInvJacL2ForFace(inv_jac));
271  fe.getOpPtrVector().push_back(new OpMakeHighOrderGeometryWeightsOnFace());
272  fe.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<2>(
273  "FIELD1", diff_field_vals_ptr, space == L2 ? MBQUAD : MBVERTEX));
274  fe.getOpPtrVector().push_back(
275  new QuadOpCheck(field_vals_ptr, diff_field_vals_ptr));
276  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe);
278  };
279 
280  CHKERR assemble_matrices_and_vectors();
281  CHKERR solve_problem();
282  CHKERR check_solution();
283  }
284  CATCH_ERRORS;
285 
287  return 0;
288 }
static Index< 'n', 3 > n
static Index< 'p', 3 > p
@ COL
Definition: definitions.h:192
@ ROW
Definition: definitions.h:192
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
FieldApproximationBase
approximation base
Definition: definitions.h:150
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:154
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:156
FieldSpace
approximation spaces
Definition: definitions.h:174
@ L2
field with C-1 continuity
Definition: definitions.h:180
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:127
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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 build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
auto createKSP
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1953
const double D
diffusivity
static char help[]
static constexpr int approx_order
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
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)
Matrix manager is used to build and partition problems.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Modify integration weights on face to take in account higher-order geometry.
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:36

Variable Documentation

◆ approx_order

constexpr int approx_order = 6
staticconstexpr

◆ help

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