v0.14.0
scalar_check_approximation.cpp
Go to the documentation of this file.
1 /**
2  * \file scalar_check_approximation.cpp
3  * \example scalar_check_approximation.cpp
4  *
5  * Approximate polynomial in 2D and check derivatives
6  *
7  */
8 
9 #include <MoFEM.hpp>
10 
11 using namespace MoFEM;
12 
13 static char help[] = "...\n\n";
14 
15 static int approx_order = 4;
16 
17 template <int DIM> struct ApproxFunctionsImpl {};
18 
19 template <int DIM> struct ElementsAndOps {};
20 
21 constexpr int SPACE_DIM =
22  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
23 
25 using BoundaryEle =
28 
31 
32 template <> struct ApproxFunctionsImpl<2> {
33  static double fUn(const double x, const double y, double z) {
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 FTensor::Tensor1<double, 2> diffFun(const double x, const double y,
46  double z) {
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 template <> struct ApproxFunctionsImpl<3> {
62  static double fUn(const double x, const double y, double z) {
63  double r = 1;
64  for (int o = 1; o <= approx_order; ++o) {
65  for (int i = 0; i <= o; ++i) {
66  for (int j = 0; j <= o - i; j++) {
67  int k = o - i - j;
68  if (k >= 0) {
69  r += pow(x, i) * pow(y, j) * pow(z, k);
70  }
71  }
72  }
73  }
74  return r;
75  }
76 
77  static FTensor::Tensor1<double, 3> diffFun(const double x, const double y,
78  double z) {
79  FTensor::Tensor1<double, 3> r{0., 0., 0.};
80  for (int o = 1; o <= approx_order; ++o) {
81  for (int i = 0; i <= o; ++i) {
82  for (int j = 0; j <= o - i; j++) {
83  int k = o - i - j;
84  if (k >= 0) {
85  r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) * pow(z, k) : 0;
86  r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) * pow(z, k) : 0;
87  r(2) += k > 0 ? k * pow(x, i) * pow(y, j) * pow(z, k - 1) : 0;
88  }
89  }
90  }
91  }
92  return r;
93  }
94 };
95 
97 
98 struct OpValsDiffVals : public DomainEleOp {
101  const bool checkGradients;
102  OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
103  : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
104  checkGradients(check_grads) {}
105 
107 
108  MoFEMErrorCode doWork(int side, EntityType type,
111  const int nb_gauss_pts = getGaussPts().size2();
112  if (type == MBVERTEX) {
113  vAls.resize(nb_gauss_pts, false);
114  diffVals.resize(SPACE_DIM, nb_gauss_pts, false);
115  vAls.clear();
116  diffVals.clear();
117  }
118 
119  const int nb_dofs = data.getIndices().size();
120  if (nb_dofs) {
121 
122  MOFEM_LOG("AT", Sev::noisy)
123  << "Type " << moab::CN::EntityTypeName(type) << " side " << side;
124  MOFEM_LOG("AT", Sev::noisy) << data.getN();
125  MOFEM_LOG("AT", Sev::noisy) << data.getDiffN();
126 
127  auto t_vals = getFTensor0FromVec(vAls);
128  auto t_base_fun = data.getFTensor0N();
129  for (int gg = 0; gg != nb_gauss_pts; gg++) {
130  auto t_data = data.getFTensor0FieldData();
131  for (int bb = 0; bb != nb_dofs; bb++) {
132  t_vals += t_base_fun * t_data;
133  ++t_base_fun;
134  ++t_data;
135  }
136  const double v = t_vals;
137  if (!std::isnormal(v))
138  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
139  ++t_vals;
140  }
141 
142  if (checkGradients) {
143  auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
144  auto t_diff_base_fun = data.getFTensor1DiffN<SPACE_DIM>();
145  for (int gg = 0; gg != nb_gauss_pts; gg++) {
146  auto t_data = data.getFTensor0FieldData();
147  for (int bb = 0; bb != nb_dofs; bb++) {
148  t_diff_vals(i) += t_diff_base_fun(i) * t_data;
149  ++t_diff_base_fun;
150  ++t_data;
151  }
152  for (int d = 0; d != SPACE_DIM; ++d)
153  if (!std::isnormal(t_diff_vals(d)))
154  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
155  ++t_diff_vals;
156  }
157  }
158  }
160  }
161 };
162 
163 struct OpCheckValsDiffVals : public DomainEleOp {
166  boost::shared_ptr<VectorDouble> ptrVals;
167  boost::shared_ptr<MatrixDouble> ptrDiffVals;
168  const bool checkGradients;
169 
171  boost::shared_ptr<VectorDouble> &ptr_vals,
172  boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
173  bool check_grads)
174  : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
175  ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
176  checkGradients(check_grads) {
177  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
178  }
179 
181 
182  MoFEMErrorCode doWork(int side, EntityType type,
185  const double eps = 1e-6;
186  const int nb_gauss_pts = data.getN().size1();
187 
188  auto t_vals = getFTensor0FromVec(vAls);
189 
190  auto t_ptr_vals = getFTensor0FromVec(*ptrVals);
191 
192  for (int gg = 0; gg != nb_gauss_pts; gg++) {
193 
194  double err_val;
195 
196  // Check user data operators
197  err_val = std::abs(t_vals - t_ptr_vals);
198  MOFEM_LOG("AT", Sev::noisy) << "Val op error " << err_val;
199 
200  if (err_val > eps)
201  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
202  "Wrong value from operator %4.3e", err_val);
203 
204  const double x = getCoordsAtGaussPts()(gg, 0);
205  const double y = getCoordsAtGaussPts()(gg, 1);
206  const double z = getCoordsAtGaussPts()(gg, 2);
207 
208  // Check approximation
209  const double delta_val = t_vals - ApproxFunctions::fUn(x, y, z);
210 
211  err_val = std::fabs(delta_val * delta_val);
212  MOFEM_LOG("AT", Sev::verbose) << err_val << " : " << t_vals;
213  if (err_val > eps)
214  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong value %4.3e",
215  err_val);
216 
217  ++t_vals;
218  ++t_ptr_vals;
219  }
220 
221  if (checkGradients) {
222 
223  auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
224  auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
225 
226  for (int gg = 0; gg != nb_gauss_pts; gg++) {
227 
228  FTensor::Tensor1<double, SPACE_DIM> t_delta_diff_val;
229  double err_diff_val;
230 
231  t_delta_diff_val(i) = t_diff_vals(i) - t_ptr_diff_vals(i);
232  err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
233  MOFEM_LOG("AT", Sev::noisy) << "Diff val op error " << err_diff_val;
234 
235  if (err_diff_val > eps)
236  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
237  "Wrong derivatives from operator %4.3e", err_diff_val);
238 
239  const double x = getCoordsAtGaussPts()(gg, 0);
240  const double y = getCoordsAtGaussPts()(gg, 1);
241  const double z = getCoordsAtGaussPts()(gg, 2);
242 
243  // Check approximation
244  auto t_diff_anal = ApproxFunctions::diffFun(x, y, z);
245  t_delta_diff_val(i) = t_diff_vals(i) - t_diff_anal(i);
246 
247  err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
248  if (SPACE_DIM == 3)
249  MOFEM_LOG("AT", Sev::noisy)
250  << "Diff val " << err_diff_val << " : "
251  << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
252  << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
253  << t_diff_vals(1) << " (" << t_diff_anal(1) << ") "
254  << t_diff_vals(2) << " (" << t_diff_anal(2) << ")";
255  else
256  MOFEM_LOG("AT", Sev::noisy)
257  << "Diff val " << err_diff_val << " : "
258  << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
259  << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
260  << t_diff_vals(1) << " (" << t_diff_anal(1) << ")";
261 
262  MOFEM_LOG("AT", Sev::verbose)
263  << getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
264  MOFEM_LOG("AT", Sev::verbose)
265  << getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
266  MOFEM_LOG("AT", Sev::verbose)
267  << getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
268 
269  MOFEM_LOG("AT", Sev::verbose) << "Diff val error " << err_diff_val;
270  if (err_diff_val > eps)
271  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
272  "Wrong derivative of value %4.3e %4.3e", err_diff_val,
273  t_diff_anal.l2());
274 
275  ++t_diff_vals;
276  ++t_ptr_diff_vals;
277  }
278  }
279 
281  }
282 };
283 
284 int main(int argc, char *argv[]) {
285 
286  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
287 
288  try {
289 
290  DMType dm_name = "DMMOFEM";
291  CHKERR DMRegister_MoFEM(dm_name);
292 
293  moab::Core mb_instance;
294  moab::Interface &moab = mb_instance;
295 
296  // Add logging channel for example
297  auto core_log = logging::core::get();
298  core_log->add_sink(
300  LogManager::setLog("AT");
301  MOFEM_LOG_TAG("AT", "atom_test");
302 
303  // Create MoFEM instance
304  MoFEM::Core core(moab);
305  MoFEM::Interface &m_field = core;
306 
307  Simple *simple = m_field.getInterface<Simple>();
308  PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
309  CHKERR simple->getOptions();
310 
311  simple->getAddBoundaryFE() = true;
312 
313  CHKERR simple->loadFile("", "");
314 
315  // Declare elements
316  enum bases {
317  AINSWORTH,
318  AINSWORTH_LOBATTO,
319  DEMKOWICZ,
320  BERNSTEIN,
321  LASBASETOP
322  };
323  const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
324  "bernstein"};
325  PetscBool flg;
326  PetscInt choice_base_value = AINSWORTH;
327  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
328  LASBASETOP, &choice_base_value, &flg);
329 
330  if (flg != PETSC_TRUE)
331  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
333  if (choice_base_value == AINSWORTH)
335  if (choice_base_value == AINSWORTH_LOBATTO)
336  base = AINSWORTH_LOBATTO_BASE;
337  else if (choice_base_value == DEMKOWICZ)
338  base = DEMKOWICZ_JACOBI_BASE;
339  else if (choice_base_value == BERNSTEIN)
341 
342  enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
343  const char *list_spaces[] = {"h1", "l2"};
344  PetscInt choice_space_value = H1SPACE;
345  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
346  LASBASETSPACE, &choice_space_value, &flg);
347  if (flg != PETSC_TRUE)
348  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
349  FieldSpace space = H1;
350  if (choice_space_value == H1SPACE)
351  space = H1;
352  else if (choice_space_value == L2SPACE)
353  space = L2;
354 
355  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
356  PETSC_NULL);
357 
358  CHKERR simple->addDomainField("FIELD1", space, base, 1);
359  CHKERR simple->setFieldOrder("FIELD1", approx_order);
360 
361  Range edges, faces;
362  CHKERR moab.get_entities_by_dimension(0, 1, edges);
363  CHKERR moab.get_entities_by_dimension(0, 2, faces);
364 
365  if (choice_base_value != BERNSTEIN) {
366  Range rise_order;
367 
368  int i = 0;
369  for (auto e : edges) {
370  if (!(i % 2)) {
371  rise_order.insert(e);
372  }
373  ++i;
374  }
375 
376  for (auto f : faces) {
377  if (!(i % 3)) {
378  rise_order.insert(f);
379  }
380  ++i;
381  }
382 
383  Range rise_twice;
384  for (auto e : rise_order) {
385  if (!(i % 2)) {
386  rise_twice.insert(e);
387  }
388  ++i;
389  }
390 
391  CHKERR simple->setFieldOrder("FIELD1", approx_order + 1, &rise_order);
392 
393  CHKERR simple->setFieldOrder("FIELD1", approx_order + 2, &rise_twice);
394  }
395 
396  CHKERR simple->defineFiniteElements();
397 
398  auto volume_adj = [](moab::Interface &moab, const Field &field,
399  const EntFiniteElement &fe,
400  std::vector<EntityHandle> &adjacency) {
402  EntityHandle fe_ent = fe.getEnt();
403  switch (field.getSpace()) {
404  case H1:
405  CHKERR moab.get_connectivity(&fe_ent, 1, adjacency, true);
406  case HCURL:
407  CHKERR moab.get_adjacencies(&fe_ent, 1, 1, false, adjacency,
408  moab::Interface::UNION);
409  case HDIV:
410  case L2:
411  CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, adjacency,
412  moab::Interface::UNION);
413  adjacency.push_back(fe_ent);
414  // build side table
415  for (auto ent : adjacency)
416  fe.getSideNumberPtr(ent);
417  break;
418  case NOFIELD: {
419  CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
420  false);
421  for (auto ent : adjacency) {
422  const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
423  .insert(
424  boost::shared_ptr<SideNumber>(new SideNumber(ent, -1, 0, 0)));
425  }
426  } break;
427  default:
428  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
429  "this field is not implemented for TRI finite element");
430  }
431 
433  };
434 
436  simple->getDomainFEName(), MBTET, volume_adj);
438  simple->getDomainFEName(), MBHEX, volume_adj);
439 
440  CHKERR simple->defineProblem(PETSC_TRUE);
441  CHKERR simple->buildFields();
442  CHKERR simple->buildFiniteElements();
443  CHKERR simple->buildProblem();
444 
445  auto dm = simple->getDM();
446 
447  VectorDouble vals;
448  auto jac_ptr = boost::make_shared<MatrixDouble>();
449  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
450  auto det_ptr = boost::make_shared<VectorDouble>();
451  MatrixDouble diff_vals;
452 
453  auto assemble_matrices_and_vectors = [&]() {
455 
457  PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
458 
460  pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
462  pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
463 
466 
467  pipeline_mng->getOpDomainLhsPipeline().push_back(
469  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
470  "FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
471  pipeline_mng->getOpDomainLhsPipeline().push_back(
472  createOpSchurAssembleEnd({}, {}, {}, {}, {}, true));
473 
474  pipeline_mng->getOpDomainRhsPipeline().push_back(
475  new OpSource("FIELD1", ApproxFunctions::fUn));
476 
477  auto integration_rule = [](int, int, int p_data) {
478  return 2 * p_data + 1;
479  };
480  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
481  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
482 
484  };
485 
486  auto solve_problem = [&] {
488  auto solver = pipeline_mng->createKSP();
489  CHKERR KSPSetFromOptions(solver);
490  CHKERR KSPSetUp(solver);
491 
492  auto dm = simple->getDM();
493  auto D = createDMVector(dm);
494  auto F = vectorDuplicate(D);
495 
496  CHKERR KSPSolve(solver, F, D);
497  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
498  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
499  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
500 
502  };
503 
504  auto check_solution = [&] {
506 
507  auto ptr_values = boost::make_shared<VectorDouble>();
508  auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
509 
510  pipeline_mng->getDomainLhsFE().reset();
511  pipeline_mng->getOpDomainRhsPipeline().clear();
512 
514  pipeline_mng->getOpDomainRhsPipeline(), {space});
515 
516  pipeline_mng->getOpDomainRhsPipeline().push_back(
517  new OpValsDiffVals(vals, diff_vals, true));
518  pipeline_mng->getOpDomainRhsPipeline().push_back(
519  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
520  pipeline_mng->getOpDomainRhsPipeline().push_back(
522  ptr_diff_vals));
523  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
524  vals, diff_vals, ptr_values, ptr_diff_vals, true));
525 
526  CHKERR pipeline_mng->loopFiniteElements();
527 
529  };
530 
531  auto post_proc = [&] {
533 
535 
536  auto get_domain_post_proc_fe = [&](auto post_proc_mesh_ptr) {
537  auto post_proc_fe =
538  boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
539  m_field, post_proc_mesh_ptr);
540 
542  post_proc_fe->getOpPtrVector(), {space});
543 
544  auto ptr_values = boost::make_shared<VectorDouble>();
545  auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
546 
547  post_proc_fe->getOpPtrVector().push_back(
548  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
549  post_proc_fe->getOpPtrVector().push_back(
551  ptr_diff_vals));
552 
553  post_proc_fe->getOpPtrVector().push_back(
554 
555  new OpPPMap(
556 
557  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
558 
559  {{"FIELD1", ptr_values}},
560 
561  {{"FIELD1_GRAD", ptr_diff_vals}},
562 
563  {}, {})
564 
565  );
566  return post_proc_fe;
567  };
568 
569  auto get_bdy_post_proc_fe = [&](auto post_proc_mesh_ptr) {
570  auto bdy_post_proc_fe =
571  boost::make_shared<PostProcBrokenMeshInMoabBaseCont<BoundaryEle>>(
572  m_field, post_proc_mesh_ptr);
573 
574  auto op_loop_side = new OpLoopSide<EleOnSide>(
575  m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy,
576  boost::make_shared<
578 
579  auto ptr_values = boost::make_shared<VectorDouble>();
580  auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
581 
582  // push operators to side element
584  op_loop_side->getOpPtrVector(), {space});
585  op_loop_side->getOpPtrVector().push_back(
586  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
587  op_loop_side->getOpPtrVector().push_back(
589  ptr_diff_vals));
590  // push op to boundary element
591  bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
592 
593  bdy_post_proc_fe->getOpPtrVector().push_back(
594 
595  new OpPPMap(
596 
597  bdy_post_proc_fe->getPostProcMesh(),
598  bdy_post_proc_fe->getMapGaussPts(),
599 
600  {{"FIELD1", ptr_values}},
601 
602  {{"FIELD1_GRAD", ptr_diff_vals}},
603 
604  {}, {})
605 
606  );
607 
608  return bdy_post_proc_fe;
609  };
610 
611  auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
612  auto post_proc_begin =
613  boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
614  m_field, post_proc_mesh_ptr);
615  auto post_proc_end =
616  boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
617  m_field, post_proc_mesh_ptr);
618  auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
619  auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
620 
621  CHKERR DMoFEMPreProcessFiniteElements(dm, post_proc_begin->getFEMethod());
622  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
623  domain_post_proc_fe);
624  CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
625  bdy_post_proc_fe);
626  CHKERR DMoFEMPostProcessFiniteElements(dm, post_proc_end->getFEMethod());
627  CHKERR post_proc_end->writeFile("out_post_proc.h5m");
628 
630  };
631 
632  CHKERR assemble_matrices_and_vectors();
633  CHKERR solve_problem();
634  CHKERR check_solution();
635  CHKERR post_proc();
636  }
637  CATCH_ERRORS;
638 
640 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreInterface::modify_finite_element_adjacency_table
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
help
static char help[]
Definition: scalar_check_approximation.cpp:13
MoFEM::ForcesAndSourcesCore::UserDataOperator::AdjCache
std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > >> AdjCache
Definition: ForcesAndSourcesCore.hpp:904
OpCheckValsDiffVals::diffVals
MatrixDouble & diffVals
Definition: scalar_check_approximation.cpp:165
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
ApproxFunctionsImpl< 3 >::fUn
static double fUn(const double x, const double y, double z)
Definition: scalar_check_approximation.cpp:62
ApproxFunctions
Definition: hcurl_check_approx_in_2d.cpp:44
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpCheckValsDiffVals::ptrDiffVals
boost::shared_ptr< MatrixDouble > ptrDiffVals
Definition: scalar_check_approximation.cpp:167
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
SPACE_DIM
constexpr int SPACE_DIM
Definition: scalar_check_approximation.cpp:21
OpValsDiffVals
Definition: scalar_check_approximation.cpp:98
EleOnSide
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
Definition: scalar_check_approximation.cpp:27
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
OpValsDiffVals::checkGradients
const bool checkGradients
Definition: scalar_check_approximation.cpp:101
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
ApproxFunctionsImpl< 2 >::diffFun
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y, double z)
Definition: scalar_check_approximation.cpp:45
MoFEM::EntFiniteElement
Finite element data for entity.
Definition: FEMultiIndices.hpp:501
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1316
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
ApproxFunctions::fUn
static FTensor::Tensor1< double, BASE_DIM > fUn(const double x, const double y, double z)
Definition: hcurl_check_approx_in_2d.cpp:46
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
sdf.r
int r
Definition: sdf.py:8
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
main
int main(int argc, char *argv[])
Definition: scalar_check_approximation.cpp:284
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
OpCheckValsDiffVals::OpCheckValsDiffVals
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals, bool check_grads)
Definition: scalar_check_approximation.cpp:170
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, std::vector< SmartPetscObj< AO >> sequence_of_aos, std::vector< SmartPetscObj< Mat >> sequence_of_mats, std::vector< bool > sym_schur, bool symm_op, boost::shared_ptr< BlockStructure > diag_blocks)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2442
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
approx_order
static int approx_order
Definition: scalar_check_approximation.cpp:15
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
SideNumber_multiIndex
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
Definition: RefEntsMultiIndices.hpp:101
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2437
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
ApproxFunctionsImpl< 3 >::diffFun
static FTensor::Tensor1< double, 3 > diffFun(const double x, const double y, double z)
Definition: scalar_check_approximation.cpp:77
FTensor::Index< 'i', SPACE_DIM >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
OpCheckValsDiffVals::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: scalar_check_approximation.cpp:182
MoFEM::EntitiesFieldData::EntData::getFTensor0FieldData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Definition: EntitiesFieldData.cpp:506
DomainEleOp
MoFEM::SideNumber
keeps information about side number for the finite element
Definition: RefEntsMultiIndices.hpp:57
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
ApproxFunctionsImpl
Definition: scalar_check_approximation.cpp:17
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
OpCheckValsDiffVals::i
FTensor::Index< 'i', SPACE_DIM > i
Definition: scalar_check_approximation.cpp:180
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
OpCheckValsDiffVals::checkGradients
const bool checkGradients
Definition: scalar_check_approximation.cpp:168
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
ApproxFunctions::diffFun
static FTensor::Tensor2< double, BASE_DIM, SPACE_DIM > diffFun(const double x, const double y)
Definition: hcurl_check_approx_in_2d.cpp:68
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
EntData
EntitiesFieldData::EntData EntData
Definition: scalar_check_approximation.cpp:29
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
OpValsDiffVals::vAls
VectorDouble & vAls
Definition: scalar_check_approximation.cpp:99
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ApproxFunctionsImpl< 2 >::fUn
static double fUn(const double x, const double y, double z)
Definition: scalar_check_approximation.cpp:33
OpCheckValsDiffVals::vAls
VectorDouble & vAls
Definition: scalar_check_approximation.cpp:164
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
OpCheckValsDiffVals::ptrVals
boost::shared_ptr< VectorDouble > ptrVals
Definition: scalar_check_approximation.cpp:166
OpCheckValsDiffVals
Definition: hcurl_check_approx_in_2d.cpp:161
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
OpValsDiffVals::diffVals
MatrixDouble & diffVals
Definition: scalar_check_approximation.cpp:100
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
OpValsDiffVals::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: scalar_check_approximation.cpp:108
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
OpValsDiffVals::i
FTensor::Index< 'i', SPACE_DIM > i
Definition: scalar_check_approximation.cpp:106
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
F
@ F
Definition: free_surface.cpp:394
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
OpValsDiffVals::OpValsDiffVals
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
Definition: scalar_check_approximation.cpp:102