v0.9.2
PostProcOnRefMesh.hpp
Go to the documentation of this file.
1 /** \file PostProcOnRefMesh.hpp
2  * \brief Post-process fields on refined mesh
3  *
4  * Create refined mesh, without enforcing continuity between element. Calculate
5  * field values on nodes of that mesh.
6  *
7  * \ingroup mofem_fs_post_proc
8  */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #ifndef __POSTPROC_ON_REF_MESH_HPP
25 #define __POSTPROC_ON_REF_MESH_HPP
26 
27 /** \brief Set of operators and data structures used for post-processing
28 
29 This set of functions works that for given problem a new MoAB instance is crated
30 only used for post-processing. For each post-processed element in the problem
31 an element entity in post-processing mesh is created. Post-processed elements do
32 not share nodes and any others entities between them, such that discontinuities
33 between element could be shown.
34 
35 Post-processed entities could be represented by ho-elements, for example 10 node
36 tetrahedrons. Moreover each element could be refined such that higher order
37 polynomials could be well represented.
38 
39 * \ingroup mofem_fs_post_proc
40 
41 */
43 
44  struct CommonData {
45  std::map<std::string, std::vector<VectorDouble>> fieldMap;
46  std::map<std::string, std::vector<MatrixDouble>> gradMap;
47  };
48 
50  Range tEts;
51  };
52 
53  /**
54  * \brief operator to post-process (save gradients on refined post-processing
55  * mesh) field gradient \ingroup mofem_fs_post_proc
56  *
57  * \todo Implamentation of setting values to fieldMap for Hcurl and Hdiv not
58  * implemented
59  *
60  */
63 
65  std::vector<EntityHandle> &mapGaussPts;
67  const std::string tagName;
68  Vec V;
69 
71  std::vector<EntityHandle> &map_gauss_pts,
72  const std::string field_name, const std::string tag_name,
73  CommonData &common_data, Vec v = PETSC_NULL)
75  field_name, UserDataOperator::OPCOL),
76  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
77  commonData(common_data), tagName(tag_name), V(v) {}
78 
81 
82  MoFEMErrorCode doWork(int side, EntityType type,
84  };
85 
86  /**
87  * \brief operator to post-process (save gradients on refined post-processing
88  * mesh) field gradient \ingroup mofem_fs_post_proc
89  *
90  * \todo Implementation for Hdiv and Hcurl to be implemented
91  *
92  */
95 
97  std::vector<EntityHandle> &mapGaussPts;
99  const std::string tagName;
100  Vec V;
101 
103  std::vector<EntityHandle> &map_gauss_pts,
104  const std::string field_name,
105  const std::string tag_name,
106  CommonData &common_data, Vec v = PETSC_NULL)
108  field_name, UserDataOperator::OPCOL),
109  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
110  commonData(common_data), tagName(tag_name), V(v) {}
111 
114 
115  MoFEMErrorCode doWork(int side, EntityType type,
117  };
118 };
119 
120 /**
121  * \brief Generic post-processing class
122  *
123  * Generate refined mesh and save data on vertices
124  *
125  * \ingroup mofem_fs_post_proc
126  */
127 template <class ELEMENT> struct PostProcTemplateOnRefineMesh : public ELEMENT {
128 
129  moab::Core coreMesh;
131  std::vector<EntityHandle> mapGaussPts;
132 
134  : ELEMENT(m_field), postProcMesh(coreMesh) {}
135 
137  ParallelComm *pcomm_post_proc_mesh =
138  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
139  if (pcomm_post_proc_mesh != NULL)
140  delete pcomm_post_proc_mesh;
141  }
142 
144  THROW_MESSAGE("not implemented");
145  }
146 
147  /** \brief Add operator to post-process L2, H1, Hdiv, Hcurl field value
148 
149  \param field_name
150  \param v If vector is given, values from vector are used to set tags on mesh
151 
152  Note:
153  Name of the tag to store values on post-process mesh is the same as field
154  name
155 
156  * \ingroup mofem_fs_post_proc
157 
158  */
159  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
160  Vec v = PETSC_NULL) {
162  ELEMENT::getOpPtrVector().push_back(
164  field_name, field_name,
165  getCommonData(), v));
167  }
168 
169  /** \brief Add operator to post-process L2 or H1 field value
170 
171  \param field_name
172  \param tag_name to store results on post-process mesh
173  \param v If vector is given, values from vector are used to set tags on mesh
174 
175  * \ingroup mofem_fs_post_proc
176 
177  */
178  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
179  const std::string tag_name,
180  Vec v = PETSC_NULL) {
182  ELEMENT::getOpPtrVector().push_back(
184  field_name, tag_name,
185  getCommonData(), v));
187  }
188 
189  /** \brief Add operator to post-process L2 or H1 field gradient
190 
191  \param field_name
192  \param v If vector is given, values from vector are used to set tags on mesh
193 
194  Note:
195  Name of the tag to store values on post-process mesh is the same as field
196  name
197 
198  * \ingroup mofem_fs_post_proc
199 
200  */
201  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
202  Vec v = PETSC_NULL) {
204  ELEMENT::getOpPtrVector().push_back(
206  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
207  getCommonData(), v));
209  }
210 
211  /** \brief Add operator to post-process L2 or H1 field gradient
212 
213  \param field_name
214  \param tag_name to store results on post-process mesh
215  \param v If vector is given, values from vector are used to set tags on mesh
216 
217  * \ingroup mofem_fs_post_proc
218 
219  */
220  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
221  const std::string tag_name,
222  Vec v = PETSC_NULL) {
224  ELEMENT::getOpPtrVector().push_back(
226  postProcMesh, mapGaussPts, field_name, tag_name, getCommonData(),
227  v));
229  }
230 
231  /**
232  * \brief wrote results in (MOAB) format, use "file_name.h5m"
233  * @param file_name file name (should always end with .h5m)
234  * @return error code
235 
236  * \ingroup mofem_fs_post_proc
237 
238  */
239  MoFEMErrorCode writeFile(const std::string file_name) {
241  CHKERR postProcMesh.write_file(file_name.c_str(), "MOAB",
242  "PARALLEL=WRITE_PART");
244  }
245 };
246 
247 template <class VOLUME_ELEMENT>
249  : public PostProcTemplateOnRefineMesh<VOLUME_ELEMENT> {
250 
252 
255 
257  bool ten_nodes_post_proc_tets = true,
258  int nb_ref_levels = -1)
259  : PostProcTemplateOnRefineMesh<VOLUME_ELEMENT>(m_field),
260  tenNodesPostProcTets(ten_nodes_post_proc_tets),
261  nbOfRefLevels(nb_ref_levels) {}
262 
263  // Gauss pts set on refined mesh
264  int getRule(int order) { return -1; };
265 
268 
270  return commonData;
271  }
272 
273  /** \brief Generate reference mesh on single element
274 
275  Each element is subdivided on smaller elements, i.e. a reference mesh on
276  single element is created. Nodes of such reference mesh are used as
277  integration points at which field values are calculated and to
278  each node a "moab" tag is attached to store those values.
279 
280  */
283 
284  auto get_nb_of_ref_levels_from_options = [this] {
285  if (nbOfRefLevels == -1) {
286  int max_level = 0;
287  PetscBool flg = PETSC_TRUE;
288  PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
289  "-my_max_post_proc_ref_level", &max_level, &flg);
290  return max_level;
291  } else {
292  return nbOfRefLevels;
293  }
294  return 0;
295  };
296  const int max_level = get_nb_of_ref_levels_from_options();
297 
298  moab::Core core_ref;
299  moab::Interface &moab_ref = core_ref;
300 
301  auto create_reference_element = [&moab_ref]() {
303  const double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
304  EntityHandle nodes[4];
305  for (int nn = 0; nn < 4; nn++) {
306  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
307  }
308  EntityHandle tet;
309  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
311  };
312 
313  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
314  MoFEM::Interface &m_field_ref = m_core_ref;
315 
316  auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
318  // seed ref mofem database by setting bit ref level to reference
319  // tetrahedron
320  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
321  0, 3, BitRefLevel().set(0));
322  for (int ll = 0; ll != max_level; ++ll) {
323  PetscPrintf(T::mField.get_comm(), "Refine Level %d\n", ll);
324  Range edges;
325  CHKERR m_field_ref.getInterface<BitRefManager>()
326  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
327  BitRefLevel().set(), MBEDGE, edges);
328  Range tets;
329  CHKERR m_field_ref.getInterface<BitRefManager>()
330  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
331  BitRefLevel(ll).set(), MBTET, tets);
332  // refine mesh
333  MeshRefinement *m_ref;
334  CHKERR m_field_ref.getInterface(m_ref);
335  CHKERR m_ref->add_vertices_in_the_middle_of_edges(
336  edges, BitRefLevel().set(ll + 1));
337  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
338  }
340  };
341 
342  auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
343  &m_field_ref]() {
345  for (int ll = 0; ll != max_level + 1; ++ll) {
346  Range tets;
347  CHKERR m_field_ref.getInterface<BitRefManager>()
348  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
349  BitRefLevel().set(ll), MBTET, tets);
350  if (tenNodesPostProcTets) {
351  EntityHandle meshset;
352  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
353  CHKERR moab_ref.add_entities(meshset, tets);
354  CHKERR moab_ref.convert_entities(meshset, true, false, false);
355  CHKERR moab_ref.delete_entities(&meshset, 1);
356  }
357  Range elem_nodes;
358  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
359 
360  auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
361  gauss_pts.resize(elem_nodes.size(), 4, false);
362  std::map<EntityHandle, int> little_map;
363  Range::iterator nit = elem_nodes.begin();
364  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
365  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
366  little_map[*nit] = gg;
367  }
368  gauss_pts = trans(gauss_pts);
369 
370  auto &ref_tets = levelRefTets[ll];
371  Range::iterator tit = tets.begin();
372  for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
373  const EntityHandle *conn;
374  int num_nodes;
375  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
376  if (tt == 0) {
377  // Ref tets has number of rows equal to number of tets on element,
378  // columns are number of gauss points
379  ref_tets.resize(tets.size(), num_nodes);
380  }
381  for (int nn = 0; nn != num_nodes; ++nn) {
382  ref_tets(tt, nn) = little_map[conn[nn]];
383  }
384  }
385 
386  auto &shape_functions = levelShapeFunctions[ll];
387  shape_functions.resize(elem_nodes.size(), 4);
388  CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
389  &gauss_pts(1, 0), &gauss_pts(2, 0),
390  elem_nodes.size());
391  }
393  };
394 
395  levelRefTets.resize(max_level + 1);
396  levelGaussPtsOnRefMesh.resize(max_level + 1);
397  levelShapeFunctions.resize(max_level + 1);
398 
399  CHKERR create_reference_element();
400  CHKERR refine_ref_tetrahedron();
401  CHKERR get_ref_gauss_pts_and_shape_functions();
402 
404  }
405 
406  /** \brief Set integration points
407 
408  If reference mesh is generated on single elements. This function maps
409  reference coordinates into physical coordinates and create element
410  on post-processing mesh.
411 
412  */
415 
416  auto get_element_max_dofs_order = [this]() {
417  int max_order = 0;
418  auto &dofs_multi_index = *this->dataPtr;
419  for (auto &dof : dofs_multi_index) {
420  const int dof_order = dof->getDofOrder();
421  max_order = (max_order < dof_order) ? dof_order : max_order;
422  };
423  return max_order;
424  };
425 
426  const int dof_max_order = get_element_max_dofs_order();
427  size_t level = (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
428  if (level > (levelGaussPtsOnRefMesh.size() - 1))
429  level = levelGaussPtsOnRefMesh.size() - 1;
430 
431  auto &level_ref_gauss_pts = levelGaussPtsOnRefMesh[level];
432  auto &level_ref_tets = levelRefTets[level];
433  auto &shape_functions = levelShapeFunctions[level];
434  T::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
435  false);
436  noalias(T::gaussPts) = level_ref_gauss_pts;
437 
438  ReadUtilIface *iface;
439  CHKERR T::postProcMesh.query_interface(iface);
440 
441  const int num_nodes = level_ref_gauss_pts.size2();
442  std::vector<double *> arrays;
443  EntityHandle startv;
444  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
445  T::mapGaussPts.resize(level_ref_gauss_pts.size2());
446  for (int gg = 0; gg != num_nodes; ++gg)
447  T::mapGaussPts[gg] = startv + gg;
448 
449  Tag th;
450  int def_in_the_loop = -1;
451  CHKERR T::postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
452  th, MB_TAG_CREAT | MB_TAG_SPARSE,
453  &def_in_the_loop);
454 
455  commonData.tEts.clear();
456  const int num_el = level_ref_tets.size1();
457  const int num_nodes_on_ele = level_ref_tets.size2();
458  EntityHandle starte;
459  EntityHandle *conn;
460  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
461  starte, conn);
462  for (unsigned int tt = 0; tt != level_ref_tets.size1(); ++tt) {
463  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
464  conn[num_nodes_on_ele * tt + nn] =
465  T::mapGaussPts[level_ref_tets(tt, nn)];
466  }
467  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
468  commonData.tEts = Range(starte, starte + num_el - 1);
469  CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
470  &(T::nInTheLoop));
471 
472  EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
473  T::coords.resize(12, false);
474  {
475  const EntityHandle *conn;
476  int num_nodes;
477  T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
478  CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
479  }
480 
483  &*shape_functions.data().begin());
485  arrays[0], arrays[1], arrays[2]);
486  const double *t_coords_ele_x = &T::coords[0];
487  const double *t_coords_ele_y = &T::coords[1];
488  const double *t_coords_ele_z = &T::coords[2];
489  for (int gg = 0; gg != num_nodes; ++gg) {
491  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
492  t_coords(i) = 0;
493  for (int nn = 0; nn != 4; ++nn) {
494  t_coords(i) += t_n * t_ele_coords(i);
495  ++t_ele_coords;
496  ++t_n;
497  }
498  ++t_coords;
499  }
500 
502  }
503 
504  /** \brief Clear operators list
505 
506  Clear operators list, user can use the same mesh instance to post-process
507  different problem or the same problem with different set of post-processed
508  fields.
509 
510  */
513  T::getOpPtrVector().clear();
515  }
516 
520  ParallelComm *pcomm_post_proc_mesh =
521  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
522  if (pcomm_post_proc_mesh != NULL)
523  delete pcomm_post_proc_mesh;
524 
525  CHKERR T::postProcMesh.delete_mesh();
527  }
528 
531 
533  ParallelComm *pcomm =
534  ParallelComm::get_pcomm(&T::mField.get_moab(), MYPCOMM_INDEX);
535  ParallelComm *pcomm_post_proc_mesh =
536  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
537  if (pcomm_post_proc_mesh == NULL) {
538  pcomm_post_proc_mesh = new ParallelComm(&moab, T::mField.get_comm());
539  }
540 
541  Range edges;
542  CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
543  CHKERR T::postProcMesh.delete_entities(edges);
544  Range tris;
545  CHKERR T::postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
546  CHKERR T::postProcMesh.delete_entities(tris);
547 
548  Range tets;
549  CHKERR T::postProcMesh.get_entities_by_type(0, MBTET, tets, false);
550 
551  int rank = pcomm->rank();
552  Range::iterator tit = tets.begin();
553  for (; tit != tets.end(); tit++) {
554  CHKERR T::postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(),
555  &*tit, 1, &rank);
556  }
557 
558  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
559 
561  }
562 
563  /** \brief Add operator to post-process Hdiv field
564  */
565  MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
567  T::getOpPtrVector().push_back(
568  new OpHdivFunctions(T::postProcMesh, T::mapGaussPts, field_name));
570  }
571 
573 
575  std::vector<EntityHandle> &mapGaussPts;
576 
578  std::vector<EntityHandle> &map_gauss_pts,
579  const std::string field_name)
580  : VOLUME_ELEMENT::UserDataOperator(field_name,
581  T::UserDataOperator::OPCOL),
582  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
583 
584  MoFEMErrorCode doWork(int side, EntityType type,
587 
588  if (data.getIndices().size() == 0)
590 
591  std::vector<Tag> th;
592  th.resize(data.getFieldData().size());
593 
594  double def_VAL[9] = {0, 0, 0};
595 
596  switch (type) {
597  case MBTRI:
598  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
599  std::ostringstream ss;
600  ss << "HDIV_FACE_" << side << "_" << dd;
601  CHKERR postProcMesh.tag_get_handle(
602  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
603  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
604  }
605  break;
606  case MBTET:
607  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
608  std::ostringstream ss;
609  ss << "HDIV_TET_" << dd;
610  CHKERR postProcMesh.tag_get_handle(
611  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
612  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
613  }
614  break;
615  default:
616  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
617  "data inconsistency");
618  }
619 
620  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
621  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
622  CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
623  &data.getVectorN<3>(gg)(dd, 0));
624  }
625  }
626 
628  }
629  };
630 
631 private:
632  std::vector<MatrixDouble> levelShapeFunctions;
633  std::vector<MatrixDouble> levelGaussPtsOnRefMesh;
634  std::vector<ublas::matrix<int>> levelRefTets;
635 };
636 
637 /** \brief Post processing
638  * \ingroup mofem_fs_post_proc
639  */
642  MoFEM::VolumeElementForcesAndSourcesCore> {
643 
645  bool ten_nodes_post_proc_tets = true,
646  int nb_ref_levels = -1)
649 };
650 
651 // /** \deprecated Use PostPocOnRefinedMesh instead
652 // */
653 // DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
654 
655 /**
656  * \brief Postprocess on prism
657  *
658  * \ingroup mofem_fs_post_proc
659  */
662  MoFEM::FatPrismElementForcesAndSourcesCore> {
663 
665 
667  bool ten_nodes_post_proc_tets = true)
670  tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
671 
672  int getRuleTrianglesOnly(int order) { return -1; };
673  int getRuleThroughThickness(int order) { return -1; };
674 
675  struct PointsMap3D {
676  const int kSi;
677  const int eTa;
678  const int zEta;
679  int nN;
680  PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
681  : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
682  };
683 
684  typedef multi_index_container<
685  PointsMap3D,
686  indexed_by<ordered_unique<composite_key<
687  PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
688  member<PointsMap3D, const int, &PointsMap3D::eTa>,
689  member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
691 
693 
694  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
695  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
697 
698  std::map<EntityHandle, EntityHandle> elementsMap;
699 
702 
705 
707  return commonData;
708  }
709 };
710 
711 /**
712  * \brief Postprocess on face
713  *
714  * \ingroup mofem_fs_post_proc
715  */
717  MoFEM::FaceElementForcesAndSourcesCore> {
718 
720 
722  bool six_node_post_proc_tris = true)
724  m_field),
725  sixNodePostProcTris(six_node_post_proc_tris) {}
726 
727  // Gauss pts set on refined mesh
728  int getRule(int order) { return -1; };
729 
732 
733  std::map<EntityHandle, EntityHandle> elementsMap;
734 
737 
740 
742  return commonData;
743  }
744 
747 
749  std::vector<EntityHandle> &mapGaussPts;
750  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
751  const std::string feVolName;
752  boost::shared_ptr<MatrixDouble> gradMatPtr;
753  const std::string tagName;
754  const bool saveOnTag;
755 
757  moab::Interface &post_proc_mesh,
758  std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
759  const std::string tag_name,
760  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
761  const std::string vol_fe_name,
762  boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag)
764  field_name, UserDataOperator::OPCOL),
765  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
766  sideOpFe(side_fe), feVolName(vol_fe_name), gradMatPtr(grad_mat_ptr),
767  tagName(tag_name), saveOnTag(save_on_tag) {}
768 
769  MoFEMErrorCode doWork(int side, EntityType type,
771  };
772 
774  const std::string field_name, const std::string vol_fe_name,
775  boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
776  bool save_on_tag = true);
777 
778 private:
781 };
782 
783 /**
784  * @brief Postprocess 2d face elements
785  *
786  */
788 
790 
792 };
793 
794 /**
795  * \brief Postprocess on edge
796  *
797  * \ingroup mofem_fs_post_proc
798  */
800  MoFEM::EdgeElementForcesAndSourcesCore> {
801 
803 
805  bool six_node_post_proc_tris = true)
807  m_field),
808  sixNodePostProcTris(six_node_post_proc_tris) {}
809 
810  // Gauss pts set on refined mesh
811  int getRule(int order) { return -1; };
812 
815 
816  std::map<EntityHandle, EntityHandle> elementsMap;
817 
820 
823 
825  return commonData;
826  }
827 
828 };
829 
830 #endif //__POSTPROC_ON_REF_MESH_HPP
831 
832 /**
833  * \defgroup mofem_fs_post_proc Post Process
834  * \ingroup user_modules
835  **/
PostProcFaceOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
PostProcVolumeOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true, int nb_ref_levels=-1)
int getRule(int order)
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
operator to post-process (save gradients on refined post-processing mesh) field gradient
operator to post-process (save gradients on refined post-processing mesh) field gradient
CommonData commonData
Deprecated interface functions.
MoFEMErrorCode preProcess()
const int kSi
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
std::vector< EntityHandle > & mapGaussPts
CommonData commonData
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
int getRule(int order)
VectorDouble * vAluesPtr
moab::Interface & postProcMesh
bool tenNodesPostProcTets
int getRuleThroughThickness(int order)
Set of operators and data structures used for post-processing.
VectorDouble vAlues
std::map< std::string, std::vector< MatrixDouble > > gradMap
moab::Interface & postProcMesh
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:318
bool sixNodePostProcTris
OpGetFieldGradientValues(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, CommonData &common_data, Vec v=PETSC_NULL)
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
std::map< EntityHandle, EntityHandle > elementsMap
Range tEts
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
std::map< EntityHandle, EntityHandle > elementsMap
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
Vec V
PostProcTemplateOnRefineMesh(MoFEM::Interface &m_field)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
const std::string tagName
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field value.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode postProcess()
function is run at the end of loop
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode postProcess()
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode clearOperators()
Clear operators list.
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:626
const std::string tagName
moab::Interface & postProcMesh
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
CommonData & commonData
OpHdivFunctions(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
MoFEMErrorCode operator()()
MoFEMErrorCode postProcess()
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Postprocess on prism.
const std::string feVolName
PostProcEdgeOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
Vec V
boost::shared_ptr< MatrixDouble > gradMatPtr
MatrixDouble gaussPtsTri
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name)
Add operator to post-process Hdiv field.
std::vector< ublas::matrix< int > > levelRefTets
CommonData & commonData
MoFEMErrorCode generateReferenceElementMesh()
bool tenNodesPostProcTets
virtual ~PostProcTemplateOnRefineMesh()
MoFEMErrorCode generateReferenceElementMesh()
OpGetFieldGradientValuesOnSkin(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > side_fe, const std::string vol_fe_name, boost::shared_ptr< MatrixDouble > grad_mat_ptr, bool save_on_tag)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
int nbOfRefLevels
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
std::map< std::string, std::vector< VectorDouble > > fieldMap
Postprocess 2d face elements.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
const int zEta
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
const int eTa
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
MatrixDouble gaussPtsQuad
moab::Interface & postProcMesh
std::vector< MatrixDouble > levelShapeFunctions
DataForcesAndSourcesCore::EntData EntData
PostProcTemplateVolumeOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true, int nb_ref_levels=-1)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
moab::Core coreMesh
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
int nN
VectorDouble * vAluesPtr
OpGetFieldValues(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, const std::string tag_name, CommonData &common_data, Vec v=PETSC_NULL)
Post processing.
bool sixNodePostProcTris
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideOpFe
#define CHKERR
Inline error check.
Definition: definitions.h:602
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
std::map< EntityHandle, EntityHandle > elementsMap
MoFEMErrorCode preProcess()
Postprocess on edge.
std::vector< EntityHandle > mapGaussPts
constexpr int order
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRuleTrianglesOnly(int order)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:291
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Postprocess on face.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode preProcess()
function is run at the beginning of loop
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
MoFEMErrorCode setGaussPts(int order)
multi_index_container< PointsMap3D, indexed_by< ordered_unique< composite_key< PointsMap3D, member< PointsMap3D, const int, &PointsMap3D::kSi >, member< PointsMap3D, const int, &PointsMap3D::eTa >, member< PointsMap3D, const int, &PointsMap3D::zEta > > > > > PointsMap3D_multiIndex
CommonData commonData
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
CommonData commonData
const std::string tagName
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin(const std::string field_name, const std::string vol_fe_name, boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
PostProcFatPrismOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true)
int getRule(int order)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
VectorDouble vAlues
MoFEMErrorCode setGaussPts(int order)
Generic post-processing class.
const bool saveOnTag
PointsMap3D_multiIndex pointsMap
MoFEMErrorCode setGaussPts(int order)
Set integration points.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26