v0.9.1
quad_polynomial_approximation.cpp
Go to the documentation of this file.
1 /** \file quad_polynomial_approximation.cpp
2  \example quad_polynomial_approximation.cpp
3  \brief Checking approximation functions for quad
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #include <MoFEM.hpp>
22 
23 using namespace MoFEM;
24 
28 
29 static char help[] = "...\n\n";
30 
31 static constexpr int approx_order = 5;
32 struct ApproxFunction {
33  static inline double fun(double x, double y) {
34  double r = 1;
35  for (int o = 1; o <= approx_order; ++o) {
36  for (int i = 0; i <= o; ++i) {
37  int j = o - i;
38  if (j >= 0)
39  r += pow(x, i) * pow(y, j);
40  }
41  }
42  return r;
43  }
44 
45  static inline VectorDouble3 diff_fun(double x, double y) {
46  VectorDouble3 r(2);
47  r.clear();
48  for (int o = 1; o <= approx_order; ++o) {
49  for (int i = 0; i <= o; ++i) {
50  int j = o - i;
51  if (j >= 0) {
52  r[0] += i > 0 ? i * pow(x, i - 1) * pow(y, j) : 0;
53  r[1] += j > 0 ? j * pow(x, i) * pow(y, j - 1) : 0;
54  }
55  }
56  }
57  return r;
58  }
59 };
60 
61 struct QuadOpCheck : public OpEle {
62 
63  QuadOpCheck(boost::shared_ptr<VectorDouble> &field_vals,
64  boost::shared_ptr<MatrixDouble> &diff_field_vals);
65  MoFEMErrorCode doWork(int side, EntityType type,
67 
68 private:
69  boost::shared_ptr<VectorDouble> fieldVals;
70  boost::shared_ptr<MatrixDouble> diffFieldVals;
71 };
72 
73 struct QuadOpRhs : public OpEle {
74 
76  MoFEMErrorCode doWork(int side, EntityType type,
78 
79 private:
81 };
82 
83 struct QuadOpLhs : public OpEle {
84 
86  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
87  EntityType col_type,
90 
91 private:
93 };
94 
95 int main(int argc, char *argv[]) {
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 }
245 
246 QuadOpCheck::QuadOpCheck(boost::shared_ptr<VectorDouble> &field_vals,
247  boost::shared_ptr<MatrixDouble> &diff_field_vals)
248  : OpEle("FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
249  fieldVals(field_vals), diffFieldVals(diff_field_vals) {}
250 
254 
255  if (type == MBVERTEX) {
256  const int nb_gauss_pts = data.getN().size1();
257  auto t_coords = getFTensor1CoordsAtGaussPts();
258  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
259  double f = ApproxFunction::fun(t_coords(0), t_coords(1));
260  constexpr double eps = 1e-6;
261  if (std::abs(f - (*fieldVals)[gg]) > eps)
262  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
263  "Wrong value %6.4e != %6.4e", f, (*fieldVals)[gg]);
264  VectorDouble3 diff_f = ApproxFunction::diff_fun(t_coords(0), t_coords(1));
265 
266  std::cout << f - (*fieldVals)[gg] << " : ";
267  for (auto d : {0, 1})
268  std::cout << diff_f[d] - (*diffFieldVals)(d, gg) << " ";
269  std::cout << std::endl;
270 
271  for (auto d : {0, 1})
272  if (std::abs(diff_f[d] - (*diffFieldVals)(d, gg)) > eps)
273  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
274  "Wrong directive value (%d) %6.4e != %6.4e", diff_f[d],
275  (*diffFieldVals)(d, gg));
276 
277  ++t_coords;
278  }
279  }
281 }
282 
284  : OpEle("FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
285  F(f) {}
286 
287 MoFEMErrorCode QuadOpRhs::doWork(int side, EntityType type,
291  const int nb_dofs = data.getN().size2();
292  if (nb_dofs) {
293  const int nb_gauss_pts = data.getN().size1();
294  VectorDouble nf(nb_dofs);
295  nf.clear();
296  auto t_base = data.getFTensor0N();
297  auto t_coords = getFTensor1CoordsAtGaussPts();
298  auto t_w = getFTensor0IntegrationWeight();
299  auto t_normal = getFTensor1NormalsAtGaussPts();
300  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
301  double f = ApproxFunction::fun(t_coords(0), t_coords(1));
302  double v = t_w * f * sqrt(t_normal(i) * t_normal(i));
303  double *val = &*nf.begin();
304  for (int bb = 0; bb != nb_dofs; ++bb) {
305  *val += v * t_base;
306  ++t_base;
307  ++val;
308  }
309  ++t_coords;
310  ++t_w;
311  ++t_normal;
312  }
313  CHKERR VecSetValues(F, data, &*nf.data().begin(), ADD_VALUES);
314  }
316 }
317 
319  : OpEle("FIELD1", "FIELD1",
321  A(a) {
322  // FIXME: Can be symmetric, is not for simplicity
323  sYmm = false;
324 }
325 
326 MoFEMErrorCode QuadOpLhs::doWork(int row_side, int col_side,
327  EntityType row_type, EntityType col_type,
332  const int row_nb_dofs = row_data.getN().size2();
333  const int col_nb_dofs = col_data.getN().size2();
334  if (row_nb_dofs && col_nb_dofs) {
335  const int nb_gauss_pts = row_data.getN().size1();
336  MatrixDouble m(row_nb_dofs, col_nb_dofs);
337  m.clear();
338  auto t_w = getFTensor0IntegrationWeight();
339  auto t_normal = getFTensor1NormalsAtGaussPts();
340 
341  double *row_base_ptr = &*row_data.getN().data().begin();
342  double *col_base_ptr = &*col_data.getN().data().begin();
343  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
344 
345  double v = t_w * sqrt(t_normal(i) * t_normal(i));
346  cblas_dger(CblasRowMajor, row_nb_dofs, col_nb_dofs, v, row_base_ptr, 1,
347  col_base_ptr, 1, &*m.data().begin(), col_nb_dofs);
348 
349  row_base_ptr += row_nb_dofs;
350  col_base_ptr += col_nb_dofs;
351  ++t_w;
352  ++t_normal;
353  }
354  CHKERR MatSetValues(A, row_data, col_data, &*m.data().begin(), ADD_VALUES);
355  }
357 }
QuadOpLhs(SmartPetscObj< Mat > &a)
static double fun(double x, double y)
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
QuadOpCheck(boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &diff_field_vals)
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[]
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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.
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda)
Definition: cblas_dger.c:12
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 double fun(double x, double y, double z)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static constexpr int approx_order
SmartPetscObj< Vec > F
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
static VectorDouble3 diff_fun(double x, double y)
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.
bool sYmm
If true assume that matrix is symmetric structure.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Calculate inverse of jacobian for face element.
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
QuadOpRhs(SmartPetscObj< Vec > &f)
auto getFTensor0IntegrationWeight()
Get integration weights.
static VectorDouble3 diff_fun(double x, double y, double z)
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
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
boost::shared_ptr< VectorDouble > fieldVals
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
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
SmartPetscObj< Mat > A
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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
Data on single entity (This is passed as argument to DataOperator::doWork)
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....
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
structure to get information form mofem into DataForcesAndSourcesCore
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
static const double eps
int main(int argc, char *argv[])
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
boost::shared_ptr< MatrixDouble > diffFieldVals
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
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.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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)
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26