v0.13.0
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;
101  int spaceDim;
102 
104  std::vector<EntityHandle> &map_gauss_pts,
105  const std::string field_name,
106  const std::string tag_name,
107  CommonData &common_data, Vec v = PETSC_NULL,
108  int space_dim = 3)
110  field_name, UserDataOperator::OPCOL),
111  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
112  commonData(common_data), tagName(tag_name), V(v),
113  spaceDim(space_dim) {}
114 
117 
118  MoFEMErrorCode doWork(int side, EntityType type,
120  };
121 };
122 
123 /**
124  * \brief Generic post-processing class
125  *
126  * Generate refined mesh and save data on vertices
127  *
128  * \ingroup mofem_fs_post_proc
129  */
130 template <class ELEMENT> struct PostProcTemplateOnRefineMesh : public ELEMENT {
131 
133 
135  boost::shared_ptr<WrapMPIComm> wrapRefMeshComm;
136 
137  std::vector<EntityHandle> mapGaussPts;
138 
140  : ELEMENT(m_field), postProcMesh(coreMesh) {}
141 
143  ParallelComm *pcomm_post_proc_mesh =
144  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
145  if (pcomm_post_proc_mesh != NULL)
146  delete pcomm_post_proc_mesh;
147  }
148 
150  THROW_MESSAGE("not implemented");
151  }
152 
153  /** \brief Add operator to post-process L2, H1, Hdiv, Hcurl field value
154 
155  \param field_name
156  \param v If vector is given, values from vector are used to set tags on mesh
157 
158  Note:
159  Name of the tag to store values on post-process mesh is the same as field
160  name
161 
162  * \ingroup mofem_fs_post_proc
163 
164  */
165  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
166  Vec v = PETSC_NULL) {
168  ELEMENT::getOpPtrVector().push_back(
170  field_name, field_name,
171  getCommonData(), v));
173  }
174 
175  /** \brief Add operator to post-process L2 or H1 field value
176 
177  \param field_name
178  \param tag_name to store results on post-process mesh
179  \param v If vector is given, values from vector are used to set tags on mesh
180 
181  * \ingroup mofem_fs_post_proc
182 
183  */
184  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
185  const std::string tag_name,
186  Vec v = PETSC_NULL) {
188  ELEMENT::getOpPtrVector().push_back(
190  field_name, tag_name,
191  getCommonData(), v));
193  }
194 
195  /** \brief Add operator to post-process L2 or H1 field gradient
196 
197  \param field_name
198  \param v If vector is given, values from vector are used to set tags on mesh
199 
200  Note:
201  Name of the tag to store values on post-process mesh is the same as field
202  name
203 
204  * \ingroup mofem_fs_post_proc
205 
206  */
207  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
208  Vec v = PETSC_NULL) {
210  ELEMENT::getOpPtrVector().push_back(
212  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
213  getCommonData(), v));
215  }
216 
217  /** \brief Add operator to post-process L2 or H1 field gradient
218 
219  \param field_name
220  \param tag_name to store results on post-process mesh
221  \param v If vector is given, values from vector are used to set tags on mesh
222 
223  * \ingroup mofem_fs_post_proc
224 
225  */
226  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
227  const std::string tag_name,
228  Vec v = PETSC_NULL) {
230  ELEMENT::getOpPtrVector().push_back(
232  postProcMesh, mapGaussPts, field_name, tag_name, getCommonData(),
233  v));
235  }
236 
237  /** \brief Add operator to post-process L2 or H1 field gradient
238 
239  \param field_name
240  \param space_dim the dimension of the problem
241  \param v If vector is given, values from vector are used to set tags on mesh
242 
243  * \ingroup mofem_fs_post_proc
244 
245 */
246  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
247  int space_dim,
248  Vec v = PETSC_NULL) {
250  ELEMENT::getOpPtrVector().push_back(
252  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
253  getCommonData(), v, space_dim));
255  }
256 
257  /**
258  * \brief wrote results in (MOAB) format, use "file_name.h5m"
259  * @param file_name file name (should always end with .h5m)
260  * @return error code
261 
262  * \ingroup mofem_fs_post_proc
263 
264  */
265  MoFEMErrorCode writeFile(const std::string file_name) {
267  ParallelComm *pcomm_post_proc_mesh =
268  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
269  if (pcomm_post_proc_mesh == NULL)
270  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
271  "ParallelComm not allocated");
272  CHKERR postProcMesh.write_file(file_name.c_str(), "MOAB",
273  "PARALLEL=WRITE_PART");
275  }
276 
277 };
278 
279 template <class VOLUME_ELEMENT>
281  : public PostProcTemplateOnRefineMesh<VOLUME_ELEMENT> {
282 
284 
287 
289  bool ten_nodes_post_proc_tets = true,
290  int nb_ref_levels = -1)
291  : PostProcTemplateOnRefineMesh<VOLUME_ELEMENT>(m_field),
292  tenNodesPostProcTets(ten_nodes_post_proc_tets),
294 
295  // Gauss pts set on refined mesh
296  int getRule(int order) { return -1; };
297 
300 
302  return commonData;
303  }
304 
305  /** \brief Generate reference mesh on single element
306 
307  Each element is subdivided on smaller elements, i.e. a reference mesh on
308  single element is created. Nodes of such reference mesh are used as
309  integration points at which field values are calculated and to
310  each node a "moab" tag is attached to store those values.
311 
312  */
315 
316  auto get_nb_of_ref_levels_from_options = [this] {
317  if (nbOfRefLevels == -1) {
318  int max_level = 0;
319  PetscBool flg = PETSC_TRUE;
320  PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
321  "-my_max_post_proc_ref_level", &max_level, &flg);
322  return max_level;
323  } else {
324  return nbOfRefLevels;
325  }
326  return 0;
327  };
328 
329  auto generate_for_hex = [&]() {
331 
332  moab::Core core_ref;
333  moab::Interface &moab_ref = core_ref;
334 
335  auto create_reference_element = [&moab_ref]() {
337  constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
338  0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
339  EntityHandle nodes[8];
340  for (int nn = 0; nn < 8; nn++)
341  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
342  EntityHandle hex;
343  CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
345  };
346 
347  auto add_ho_nodes = [&]() {
349  Range hexes;
350  CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
351  EntityHandle meshset;
352  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
353  CHKERR moab_ref.add_entities(meshset, hexes);
354  CHKERR moab_ref.convert_entities(meshset, true, true, true);
355  CHKERR moab_ref.delete_entities(&meshset, 1);
357  };
358 
359  auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
361  Range hexes;
362  CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
363  Range hexes_nodes;
364  CHKERR moab_ref.get_connectivity(hexes, hexes_nodes, false);
365  auto &gauss_pts = levelGaussPtsOnRefMeshHexes[0];
366  gauss_pts.resize(hexes_nodes.size(), 4, false);
367  size_t gg = 0;
368  for (auto node : hexes_nodes) {
369  CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
370  little_map[node] = gg;
371  ++gg;
372  }
373  gauss_pts = trans(gauss_pts);
375  };
376 
377  auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
379  Range hexes;
380  CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, true);
381  size_t hh = 0;
382  auto &ref_hexes = levelRefHexes[0];
383  for (auto hex : hexes) {
384  const EntityHandle *conn;
385  int num_nodes;
386  CHKERR moab_ref.get_connectivity(hex, conn, num_nodes, false);
387  if (ref_hexes.size2() != num_nodes) {
388  ref_hexes.resize(hexes.size(), num_nodes);
389  }
390  for (int nn = 0; nn != num_nodes; ++nn) {
391  ref_hexes(hh, nn) = little_map[conn[nn]];
392  }
393  ++hh;
394  }
396  };
397 
398  auto set_shape_functions = [&]() {
400  auto &gauss_pts = levelGaussPtsOnRefMeshHexes[0];
401  auto &shape_functions = levelShapeFunctionsHexes[0];
402  const auto nb_gauss_pts = gauss_pts.size2();
403  shape_functions.resize(nb_gauss_pts, 8);
404  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
405  const double ksi = gauss_pts(0, gg);
406  const double zeta = gauss_pts(1, gg);
407  const double eta = gauss_pts(2, gg);
408  shape_functions(gg, 0) = N_MBHEX0(ksi, zeta, eta);
409  shape_functions(gg, 1) = N_MBHEX1(ksi, zeta, eta);
410  shape_functions(gg, 2) = N_MBHEX2(ksi, zeta, eta);
411  shape_functions(gg, 3) = N_MBHEX3(ksi, zeta, eta);
412  shape_functions(gg, 4) = N_MBHEX4(ksi, zeta, eta);
413  shape_functions(gg, 5) = N_MBHEX5(ksi, zeta, eta);
414  shape_functions(gg, 6) = N_MBHEX6(ksi, zeta, eta);
415  shape_functions(gg, 7) = N_MBHEX7(ksi, zeta, eta);
416  }
418  };
419 
420  levelRefHexes.resize(1);
421  levelGaussPtsOnRefMeshHexes.resize(1);
422  levelShapeFunctionsHexes.resize(1);
423 
424  CHKERR create_reference_element();
426  CHKERR add_ho_nodes();
427  std::map<EntityHandle, int> little_map;
428  CHKERR set_gauss_pts(little_map);
429  CHKERR set_ref_hexes(little_map);
430  CHKERR set_shape_functions();
431 
433  };
434 
435  auto generate_for_tet = [&]() {
437 
438  const int max_level = get_nb_of_ref_levels_from_options();
439 
440  moab::Core core_ref;
441  moab::Interface &moab_ref = core_ref;
442 
443  auto create_reference_element = [&moab_ref]() {
445  constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
446  EntityHandle nodes[4];
447  for (int nn = 0; nn < 4; nn++) {
448  CHKERR
449  moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
450  }
451  EntityHandle tet;
452  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
454  };
455 
456  MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
457  MoFEM::Interface &m_field_ref = m_core_ref;
458 
459  auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
461  // seed ref mofem database by setting bit ref level to reference
462  // tetrahedron
463  CHKERR
464  m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
465  0, 3, BitRefLevel().set(0));
466  for (int ll = 0; ll != max_level; ++ll) {
467  MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "PostProc",
468  "Refine Level %d", ll);
469  Range edges;
470  CHKERR m_field_ref.getInterface<BitRefManager>()
471  ->getEntitiesByTypeAndRefLevel(
472  BitRefLevel().set(ll), BitRefLevel().set(), MBEDGE, edges);
473  Range tets;
474  CHKERR m_field_ref.getInterface<BitRefManager>()
475  ->getEntitiesByTypeAndRefLevel(
476  BitRefLevel().set(ll), BitRefLevel(ll).set(), MBTET, tets);
477  // refine mesh
478  MeshRefinement *m_ref;
479  CHKERR m_field_ref.getInterface(m_ref);
480  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
481  edges, BitRefLevel().set(ll + 1));
482  CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
483  }
485  };
486 
487  auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
488  &m_field_ref]() {
490  for (int ll = 0; ll != max_level + 1; ++ll) {
491  Range tets;
492  CHKERR
493  m_field_ref.getInterface<BitRefManager>()
494  ->getEntitiesByTypeAndRefLevel(
495  BitRefLevel().set(ll), BitRefLevel().set(ll), MBTET, tets);
496  if (tenNodesPostProcTets) {
497  EntityHandle meshset;
498  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
499  CHKERR moab_ref.add_entities(meshset, tets);
500  CHKERR moab_ref.convert_entities(meshset, true, false, false);
501  CHKERR moab_ref.delete_entities(&meshset, 1);
502  }
503  Range elem_nodes;
504  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
505 
506  auto &gauss_pts = levelGaussPtsOnRefMeshTets[ll];
507  gauss_pts.resize(elem_nodes.size(), 4, false);
508  std::map<EntityHandle, int> little_map;
509  Range::iterator nit = elem_nodes.begin();
510  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
511  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
512  little_map[*nit] = gg;
513  }
514  gauss_pts = trans(gauss_pts);
515 
516  auto &ref_tets = levelRefTets[ll];
517  Range::iterator tit = tets.begin();
518  for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
519  const EntityHandle *conn;
520  int num_nodes;
521  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
522  if (tt == 0) {
523  ref_tets.resize(tets.size(), num_nodes);
524  }
525  for (int nn = 0; nn != num_nodes; ++nn) {
526  ref_tets(tt, nn) = little_map[conn[nn]];
527  }
528  }
529 
530  auto &shape_functions = levelShapeFunctionsTets[ll];
531  shape_functions.resize(elem_nodes.size(), 4);
532  CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
533  &gauss_pts(1, 0), &gauss_pts(2, 0),
534  elem_nodes.size());
535  }
537  };
538 
539  levelRefTets.resize(max_level + 1);
540  levelGaussPtsOnRefMeshTets.resize(max_level + 1);
541  levelShapeFunctionsTets.resize(max_level + 1);
542 
543  CHKERR create_reference_element();
544  CHKERR refine_ref_tetrahedron();
545  CHKERR get_ref_gauss_pts_and_shape_functions();
546 
548  };
549 
550  CHKERR generate_for_hex();
551  CHKERR generate_for_tet();
552 
554  }
555 
556  size_t getMaxLevel() const {
557  auto get_element_max_dofs_order = [&]() {
558  int max_order = 0;
559  auto dofs_vec = this->getDataVectorDofsPtr();
560  for (auto &dof : *dofs_vec) {
561  const int dof_order = dof->getDofOrder();
562  max_order = (max_order < dof_order) ? dof_order : max_order;
563  };
564  return max_order;
565  };
566  const auto dof_max_order = get_element_max_dofs_order();
567  return (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
568  };
569 
570  /** \brief Set integration points
571 
572  If reference mesh is generated on single elements. This function maps
573  reference coordinates into physical coordinates and create element
574  on post-processing mesh.
575 
576  */
579 
580  auto type = type_from_handle(this->getFEEntityHandle());
581 
582  auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
583  auto &level_shape_functions,
584 
585  auto start_vert_handle,
586  auto start_ele_handle,
587  auto &verts_array,
588  auto &conn,
589  auto &ver_count,
590  auto &ele_count
591 
592  ) {
594 
595  auto level =
596  std::min(getMaxLevel(), level_gauss_pts_on_ref_mesh.size() - 1);
597 
598  auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
599  auto &level_ref_ele = level_ref[level];
600  auto &shape_functions = level_shape_functions[level];
601  T::gaussPts.resize(level_ref_gauss_pts.size1(),
602  level_ref_gauss_pts.size2(), false);
603  noalias(T::gaussPts) = level_ref_gauss_pts;
604 
605  EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
606  {
607  const EntityHandle *conn;
608  int num_nodes;
609  T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
610  T::coords.resize(3 * num_nodes, false);
611  CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
612  }
613 
614  const int num_nodes = level_ref_gauss_pts.size2();
615  T::mapGaussPts.resize(level_ref_gauss_pts.size2());
616 
619  &*shape_functions.data().begin());
621  &verts_array[0][ver_count], &verts_array[1][ver_count],
622  &verts_array[2][ver_count]);
623  for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
624 
625  T::mapGaussPts[gg] = start_vert_handle + ver_count;
626 
627  auto set_float_precision = [](const double x) {
628  if (std::abs(x) < std::numeric_limits<float>::epsilon())
629  return 0.;
630  else
631  return x;
632  };
633 
634  t_coords(i) = 0;
635  auto t_ele_coords = getFTensor1FromArray<3, 3>(T::coords);
636  for (int nn = 0; nn != CN::VerticesPerEntity(type); ++nn) {
637  t_coords(i) += t_n * t_ele_coords(i);
638  ++t_ele_coords;
639  ++t_n;
640  }
641 
642  for (auto ii : {0, 1, 2})
643  t_coords(ii) = set_float_precision(t_coords(ii));
644 
645  ++t_coords;
646  }
647 
648  Tag th;
649  int def_in_the_loop = -1;
650  CHKERR T::postProcMesh.tag_get_handle(
651  "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th,
652  MB_TAG_CREAT | MB_TAG_SPARSE, &def_in_the_loop);
653 
654  commonData.tEts.clear();
655  const int num_el = level_ref_ele.size1();
656  const int num_nodes_on_ele = level_ref_ele.size2();
657  auto start_e = start_ele_handle + ele_count;
658  commonData.tEts = Range(start_e, start_e + num_el - 1);
659  for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
660  for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
661  conn[num_nodes_on_ele * ele_count + nn] =
662  T::mapGaussPts[level_ref_ele(tt, nn)];
663  }
664  }
665 
666  const int n_in_the_loop = T::nInTheLoop;
667  CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
668  &n_in_the_loop);
669 
671  };
672 
673  switch (type) {
674  case MBTET:
675  return set_gauss_pts(levelGaussPtsOnRefMeshTets, levelRefTets,
677 
681  tetConn,
682  countVertTet,
683  countTet
684 
685 
686 
687 
688  );
689  case MBHEX:
690  return set_gauss_pts(levelGaussPtsOnRefMeshHexes, levelRefHexes,
692 
696  hexConn,
697  countVertHex,
698  countHex
699 
700  );
701  default:
702  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
703  "Element type not implemented");
704  }
705 
707  }
708 
709  /** \brief Clear operators list
710 
711  Clear operators list, user can use the same mesh instance to post-process
712  different problem or the same problem with different set of post-processed
713  fields.
714 
715  */
718  T::getOpPtrVector().clear();
720  }
721 
725  ParallelComm *pcomm_post_proc_mesh =
726  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
727  if (pcomm_post_proc_mesh != NULL)
728  delete pcomm_post_proc_mesh;
729 
730  CHKERR T::postProcMesh.delete_mesh();
731 
732  auto alloc_vertives_and_elements_on_post_proc_mesh = [&]() {
734 
735  auto fe_name = this->feName;
736  auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
737 
738  auto miit =
739  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
740  boost::make_tuple(fe_name, this->getLoFERank()));
741  auto hi_miit =
742  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
743  boost::make_tuple(fe_name, this->getHiFERank()));
744 
745  const int number_of_ents_in_the_loop = this->getLoopSize();
746  if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
747  SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
748  "Wrong size of indicices. Inconsistent size number of iterated "
749  "elements iterated by problem and from range.");
750  }
751 
752  int nb_tet_vertices = 0;
753  int nb_tets = 0;
754  int nb_hex_vertices = 0;
755  int nb_hexes = 0;
756 
757  for (; miit != hi_miit; ++miit) {
758  auto type = (*miit)->getEntType();
759 
760  // Set pointer to element. So that getDataVectorDofsPtr in getMaxLevel
761  // can work
762  this->numeredEntFiniteElementPtr = *miit;
763  auto level = getMaxLevel();
764 
765  switch (type) {
766  case MBTET:
767  level = std::min(level, levelGaussPtsOnRefMeshTets.size() - 1);
768  nb_tet_vertices += levelGaussPtsOnRefMeshTets[level].size2();
769  nb_tets += levelRefTets[level].size1();
770  break;
771  case MBHEX:
772  level = std::min(level, levelGaussPtsOnRefMeshHexes.size() - 1);
773  nb_hex_vertices += levelGaussPtsOnRefMeshHexes[level].size2();
774  nb_hexes += levelRefHexes[level].size1();
775  break;
776  default:
777  SETERRQ(this->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
778  "Element type not implemented");
779  break;
780  }
781  }
782 
783  ReadUtilIface *iface;
784  CHKERR this->postProcMesh.query_interface(iface);
785 
786  if (nb_tets) {
787  CHKERR iface->get_node_coords(3, nb_tet_vertices, 0,
790  CHKERR iface->get_element_connect(nb_tets, levelRefTets[0].size2(),
791  MBTET, 0, startingEleTetHandle,
792  tetConn);
793  }
794 
795  if (nb_hexes) {
796  CHKERR iface->get_node_coords(
797  3, nb_hex_vertices, 0, startingVertHexHandle, verticesOnHexArrays);
798  CHKERR iface->get_element_connect(nb_hexes, levelRefHexes[0].size2(),
799  MBHEX, 0, startingEleHexHandle,
800  hexConn);
801  }
802 
803  countTet = 0;
804  countVertTet = 0;
805  countHex = 0;
806  countVertHex = 0;
807 
809  };
810 
811  CHKERR alloc_vertives_and_elements_on_post_proc_mesh();
812 
814  }
815 
818 
819  auto update_elements = [&]() {
821  ReadUtilIface *iface;
822  CHKERR this->postProcMesh.query_interface(iface);
823 
824  if (countTet) {
825  MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
826  << "Update tets " << countTet;
827 
828  MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
829  << "Nb nodes on tets " << levelRefTets[0].size2();
830 
831  CHKERR iface->update_adjacencies(startingEleTetHandle, countTet,
832  levelRefTets[0].size2(), tetConn);
833 
834  }
835  if (countHex) {
836  MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "PostProc")
837  << "Update hexes " << countHex;
838  CHKERR iface->update_adjacencies(startingEleHexHandle, countHex,
839  levelRefHexes[0].size2(), hexConn);
840  }
842  };
843 
844  auto resolve_shared_ents = [&]() {
846  ParallelComm *pcomm_post_proc_mesh =
847  ParallelComm::get_pcomm(&(T::postProcMesh), MYPCOMM_INDEX);
848  if (pcomm_post_proc_mesh == NULL) {
850  boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
851  pcomm_post_proc_mesh = new ParallelComm(
852  &(T::postProcMesh), (T::wrapRefMeshComm)->get_comm());
853  }
854 
855  Range edges;
856  CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
857  CHKERR T::postProcMesh.delete_entities(edges);
858  Range faces;
859  CHKERR T::postProcMesh.get_entities_by_dimension(0, 2, faces, false);
860  CHKERR T::postProcMesh.delete_entities(faces);
861 
862  Range ents;
863  CHKERR T::postProcMesh.get_entities_by_dimension(0, 3, ents, false);
864 
865  int rank = T::mField.get_comm_rank();
866  CHKERR T::postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(),
867  ents, &rank);
868 
869  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
870 
872  };
873 
874  CHKERR resolve_shared_ents();
875  CHKERR update_elements();
876 
878  }
879 
880  /** \brief Add operator to post-process Hdiv field
881  */
882  MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
884  T::getOpPtrVector().push_back(
887  }
888 
890 
892  std::vector<EntityHandle> &mapGaussPts;
893 
895  std::vector<EntityHandle> &map_gauss_pts,
896  const std::string field_name)
897  : VOLUME_ELEMENT::UserDataOperator(field_name,
898  T::UserDataOperator::OPCOL),
899  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
900 
901  MoFEMErrorCode doWork(int side, EntityType type,
904 
905  if (data.getIndices().size() == 0)
907 
908  std::vector<Tag> th;
909  th.resize(data.getFieldData().size());
910 
911  double def_VAL[9] = {0, 0, 0};
912 
913  switch (type) {
914  case MBTRI:
915  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
916  std::ostringstream ss;
917  ss << "HDIV_FACE_" << side << "_" << dd;
918  CHKERR postProcMesh.tag_get_handle(
919  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
920  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
921  }
922  break;
923  case MBTET:
924  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
925  std::ostringstream ss;
926  ss << "HDIV_TET_" << dd;
927  CHKERR postProcMesh.tag_get_handle(
928  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
929  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
930  }
931  break;
932  default:
933  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
934  "data inconsistency");
935  }
936 
937  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
938  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
939  CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
940  &data.getVectorN<3>(gg)(dd, 0));
941  }
942  }
943 
945  }
946  };
947 
948 private:
949  std::vector<MatrixDouble> levelShapeFunctionsTets;
950  std::vector<MatrixDouble> levelGaussPtsOnRefMeshTets;
951  std::vector<ublas::matrix<int>> levelRefTets;
952  std::vector<MatrixDouble> levelShapeFunctionsHexes;
953  std::vector<MatrixDouble> levelGaussPtsOnRefMeshHexes;
954  std::vector<ublas::matrix<int>> levelRefHexes;
955 
957  std::vector<double *> verticesOnTetArrays;
960  int countTet;
962 
964  std::vector<double *> verticesOnHexArrays;
967  int countHex;
969 
970 };
971 
972 /** \brief Post processing
973  * \ingroup mofem_fs_post_proc
974  */
977  MoFEM::VolumeElementForcesAndSourcesCore> {
978 
982 };
983 
984 // /** \deprecated Use PostPocOnRefinedMesh instead
985 // */
986 // DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
987 
988 /**
989  * \brief Postprocess on prism
990  *
991  * \ingroup mofem_fs_post_proc
992  */
995  MoFEM::FatPrismElementForcesAndSourcesCore> {
996 
998 
1000  bool ten_nodes_post_proc_tets = true)
1003  tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
1004 
1005  int getRuleTrianglesOnly(int order) { return -1; };
1006  int getRuleThroughThickness(int order) { return -1; };
1007 
1008  struct PointsMap3D {
1009  const int kSi;
1010  const int eTa;
1011  const int zEta;
1012  int nN;
1013  PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
1014  : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
1015  };
1016 
1017  typedef multi_index_container<
1018  PointsMap3D,
1019  indexed_by<ordered_unique<composite_key<
1020  PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
1021  member<PointsMap3D, const int, &PointsMap3D::eTa>,
1022  member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
1024 
1025  // PointsMap3D_multiIndex pointsMap;
1026 
1027  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
1028  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
1030 
1031  // std::map<EntityHandle, EntityHandle> elementsMap;
1032  std::map<EntityHandle, std::vector<EntityHandle>> elementsMap;
1033  std::map<EntityHandle, std::vector<PointsMap3D_multiIndex>>
1035 
1038 
1041 
1043  return commonData;
1044  }
1045 };
1046 
1047 /**
1048  * \brief Postprocess on face
1049  *
1050  * \ingroup mofem_fs_post_proc
1051  */
1053  MoFEM::FaceElementForcesAndSourcesCore> {
1054 
1056 
1058  bool six_node_post_proc_tris = true)
1060  m_field),
1061  sixNodePostProcTris(six_node_post_proc_tris), counterTris(0),
1062  counterQuads(0) {}
1063 
1064  // Gauss pts set on refined mesh
1065  int getRule(int order) { return -1; };
1066 
1069 
1070  std::map<EntityHandle, EntityHandle> elementsMap;
1071 
1074 
1077 
1079  return commonData;
1080  }
1081 
1082  template <int RANK>
1085 
1087  std::vector<EntityHandle> &mapGaussPts;
1088  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
1089  const std::string feVolName;
1090  boost::shared_ptr<MatrixDouble> matPtr;
1091  const std::string tagName;
1092  const std::string fieldName;
1093  const bool saveOnTag;
1094 
1096  moab::Interface &post_proc_mesh,
1097  std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
1098  const std::string tag_name,
1099  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
1100  const std::string vol_fe_name, boost::shared_ptr<MatrixDouble> mat_ptr,
1101  bool save_on_tag)
1103  field_name, UserDataOperator::OPCOL),
1104  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1105  sideOpFe(side_fe), feVolName(vol_fe_name), matPtr(mat_ptr),
1106  tagName(tag_name), fieldName(field_name), saveOnTag(save_on_tag) {}
1107 
1108  MoFEMErrorCode doWork(int side, EntityType type,
1110  };
1111 
1112  typedef struct OpGetFieldValuesOnSkinImpl<3> OpGetFieldGradientValuesOnSkin;
1113  typedef struct OpGetFieldValuesOnSkinImpl<1> OpGetFieldValuesOnSkin;
1114 
1116  const std::string field_name, const std::string vol_fe_name = "dFE",
1117  boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
1118  bool save_on_tag = true);
1119 
1121  const std::string field_name, const std::string vol_fe_name = "dFE",
1122  boost::shared_ptr<MatrixDouble> mat_ptr = nullptr,
1123  bool save_on_tag = true);
1124 
1125 private:
1126  MatrixDouble gaussPtsTri; ///< Gauss points coordinates on reference triangle
1127  MatrixDouble gaussPtsQuad; ///< Gauss points coordinates on reference quad
1128 
1129  EntityHandle *triConn; ///< Connectivity for created tri elements
1130  EntityHandle *quadConn; ///< Connectivity for created quad elements
1131  EntityHandle
1132  startingVertTriHandle; ///< Starting handle for vertices on triangles
1133  EntityHandle
1134  startingVertQuadHandle; ///< Starting handle for vertices on quads
1135  std::vector<double *>
1136  verticesOnTriArrays; /// pointers to memory allocated by MoAB for
1137  /// storing X, Y, and Z coordinates
1138  std::vector<double *>
1139  verticesOnQuadArrays; /// pointers to memory allocated by MoAB for
1140  /// storing X, Y, and Z coordinates
1141 
1142  EntityHandle
1143  startingEleTriHandle; ///< Starting handle for triangles post proc
1144  EntityHandle startingEleQuadHandle; ///< Starting handle for quads post proc
1145 
1146  int numberOfTriangles; ///< Number of triangles to create
1147  int numberOfQuads; ///< NUmber of quads to create
1150 };
1151 
1152 /**
1153  * @brief Postprocess 2d face elements
1154  *
1155  */
1157 
1159 
1161 };
1162 
1164 
1165 /**
1166  * \brief Postprocess on edge
1167  *
1168  * \ingroup mofem_fs_post_proc
1169  */
1171  : public PostProcTemplateOnRefineMesh<EdgeEleBasePostProc> {
1172 
1174 
1176  bool six_node_post_proc_tris = true)
1178  sixNodePostProcTris(six_node_post_proc_tris) {}
1179 
1180  // Gauss pts set on refined mesh
1181  int getRule(int order) { return -1; };
1182 
1185 
1186  std::map<EntityHandle, EntityHandle> elementsMap;
1187 
1190 
1193 
1195  return commonData;
1196  }
1197 };
1198 
1199 #endif //__POSTPROC_ON_REF_MESH_HPP
1200 
1201 /**
1202  * \defgroup mofem_fs_post_proc Post Process
1203  * \ingroup user_modules
1204  **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:363
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:355
EntitiesFieldData::EntData EntData
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:88
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:84
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:86
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:85
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
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:81
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:87
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:83
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:82
constexpr double eta
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
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.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, int space_dim, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
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.
constexpr int nb_ref_levels
Three levels of refinement.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto type_from_handle
get type from entity handle
Definition: Templates.hpp:1441
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Range tEts
std::map< std::string, std::vector< VectorDouble > > fieldMap
std::map< std::string, std::vector< MatrixDouble > > gradMap
operator to post-process (save gradients on refined post-processing mesh) field gradient
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
const std::string tagName
VectorDouble vAlues
Vec V
CommonData & commonData
VectorDouble * vAluesPtr
int spaceDim
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, int space_dim=3)
operator to post-process (save gradients on refined post-processing mesh) field gradient
const std::string tagName
VectorDouble vAlues
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
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)
std::vector< EntityHandle > & mapGaussPts
VectorDouble * vAluesPtr
moab::Interface & postProcMesh
Vec V
CommonData & commonData
Set of operators and data structures used for post-processing.
Postprocess on edge.
MoFEMErrorCode postProcess()
function is run at the end of loop
bool sixNodePostProcTris
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPts(int order)
int getRule(int order)
PostProcEdgeOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
CommonData commonData
std::map< EntityHandle, EntityHandle > elementsMap
MoFEMErrorCode generateReferenceElementMesh()
const bool saveOnTag
const std::string feVolName
boost::shared_ptr< MatrixDouble > matPtr
moab::Interface & postProcMesh
OpGetFieldValuesOnSkinImpl(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 > mat_ptr, bool save_on_tag)
const std::string tagName
const std::string fieldName
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideOpFe
Postprocess 2d face elements.
MoFEMErrorCode operator()()
Postprocess on face.
struct OpGetFieldValuesOnSkinImpl< 3 > OpGetFieldGradientValuesOnSkin
EntityHandle * quadConn
Connectivity for created quad elements.
int numberOfTriangles
Number of triangles to create.
MoFEMErrorCode addFieldValuesPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > mat_ptr=nullptr, bool save_on_tag=true)
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode generateReferenceElementMesh()
std::vector< double * > verticesOnTriArrays
MoFEMErrorCode preProcess()
int numberOfQuads
NUmber of quads to create.
std::map< EntityHandle, EntityHandle > elementsMap
std::vector< double * > verticesOnQuadArrays
EntityHandle startingVertQuadHandle
Starting handle for vertices on quads.
struct OpGetFieldValuesOnSkinImpl< 1 > OpGetFieldValuesOnSkin
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
CommonData commonData
PostProcFaceOnRefinedMesh(MoFEM::Interface &m_field, bool six_node_post_proc_tris=true)
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
MatrixDouble gaussPtsQuad
Gauss points coordinates on reference quad.
EntityHandle startingVertTriHandle
Starting handle for vertices on triangles.
bool sixNodePostProcTris
EntityHandle startingEleQuadHandle
Starting handle for quads post proc.
int counterTris
MoFEMErrorCode postProcess()
int getRule(int order)
int counterQuads
MatrixDouble gaussPtsTri
Gauss points coordinates on reference triangle.
EntityHandle startingEleTriHandle
Starting handle for triangles post proc.
EntityHandle * triConn
Connectivity for created tri elements.
PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
const int zEta
const int eTa
const int kSi
int nN
Postprocess on prism.
int getRuleThroughThickness(int order)
CommonData commonData
std::map< EntityHandle, std::vector< EntityHandle > > elementsMap
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
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
PostProcFatPrismOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true)
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
int getRuleTrianglesOnly(int order)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
std::map< EntityHandle, std::vector< PointsMap3D_multiIndex > > pointsMapVectorMap
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode postProcess()
function is run at the end of loop
bool tenNodesPostProcTets
Generic post-processing class.
std::vector< EntityHandle > mapGaussPts
boost::shared_ptr< WrapMPIComm > wrapRefMeshComm
virtual ~PostProcTemplateOnRefineMesh()
moab::Interface & postProcMesh
moab::Core coreMesh
PostProcTemplateOnRefineMesh(MoFEM::Interface &m_field)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
moab::Interface & postProcMesh
OpHdivFunctions(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
int countVertHex
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
std::vector< MatrixDouble > levelGaussPtsOnRefMeshTets
std::vector< ublas::matrix< int > > levelRefTets
size_t getMaxLevel() const
EntityHandle startingEleHexHandle
EntityHandle startingVertHexHandle
std::vector< MatrixDouble > levelShapeFunctionsHexes
int countVertTet
PostProcTemplateVolumeOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true, int nb_ref_levels=-1)
std::vector< double * > verticesOnHexArrays
MoFEMErrorCode setGaussPts(int order)
Set integration points.
CommonData commonData
int nbOfRefLevels
std::vector< double * > verticesOnTetArrays
bool tenNodesPostProcTets
std::vector< MatrixDouble > levelGaussPtsOnRefMeshHexes
MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name)
Add operator to post-process Hdiv field.
EntityHandle startingVertTetHandle
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
MoFEMErrorCode preProcess()
std::vector< MatrixDouble > levelShapeFunctionsTets
int countHex
EntityHandle * hexConn
std::vector< ublas::matrix< int > > levelRefHexes
MoFEMErrorCode postProcess()
EntityHandle * tetConn
int countTet
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
MoFEMErrorCode clearOperators()
Clear operators list.
EntityHandle startingEleTetHandle
int getRule(int order)
Post processing.