v0.8.23
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 
64  moab::Interface &postProcMesh;
65  std::vector<EntityHandle> &mapGaussPts;
67  const std::string tagName;
68  Vec V;
69 
70  OpGetFieldValues(moab::Interface &post_proc_mesh,
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)
74  : MoFEM::ForcesAndSourcesCore::UserDataOperator(
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 
96  moab::Interface &postProcMesh;
97  std::vector<EntityHandle> &mapGaussPts;
99  const std::string tagName;
100  Vec V;
101 
102  OpGetFieldGradientValues(moab::Interface &post_proc_mesh,
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)
107  : MoFEM::ForcesAndSourcesCore::UserDataOperator(
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;
130  moab::Interface &postProcMesh;
131  std::vector<EntityHandle> mapGaussPts;
132 
134  : ELEMENT(m_field), postProcMesh(coreMesh) {}
135 
137  THROW_MESSAGE("not implemented");
138  }
139 
140  /** \brief Add operator to post-process L2, H1, Hdiv, Hcurl field value
141 
142  \param field_name
143  \param v If vector is given, values from vector are used to set tags on mesh
144 
145  Note:
146  Name of the tag to store values on post-process mesh is the same as field
147  name
148 
149  * \ingroup mofem_fs_post_proc
150 
151  */
152  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
153  Vec v = PETSC_NULL) {
155  ELEMENT::getOpPtrVector().push_back(
157  field_name, field_name,
158  getCommonData(), v));
160  }
161 
162  /** \brief Add operator to post-process L2 or H1 field value
163 
164  \param field_name
165  \param tag_name to store results on post-process mesh
166  \param v If vector is given, values from vector are used to set tags on mesh
167 
168  * \ingroup mofem_fs_post_proc
169 
170  */
171  MoFEMErrorCode addFieldValuesPostProc(const std::string field_name,
172  const std::string tag_name,
173  Vec v = PETSC_NULL) {
175  ELEMENT::getOpPtrVector().push_back(
177  field_name, tag_name,
178  getCommonData(), v));
180  }
181 
182  /** \brief Add operator to post-process L2 or H1 field gradient
183 
184  \param field_name
185  \param v If vector is given, values from vector are used to set tags on mesh
186 
187  Note:
188  Name of the tag to store values on post-process mesh is the same as field
189  name
190 
191  * \ingroup mofem_fs_post_proc
192 
193  */
194  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
195  Vec v = PETSC_NULL) {
197  ELEMENT::getOpPtrVector().push_back(
199  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
200  getCommonData(), v));
202  }
203 
204  /** \brief Add operator to post-process L2 or H1 field gradient
205 
206  \param field_name
207  \param tag_name to store results on post-process mesh
208  \param v If vector is given, values from vector are used to set tags on mesh
209 
210  * \ingroup mofem_fs_post_proc
211 
212  */
213  MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name,
214  const std::string tag_name,
215  Vec v = PETSC_NULL) {
217  ELEMENT::getOpPtrVector().push_back(
219  postProcMesh, mapGaussPts, field_name, tag_name, getCommonData(),
220  v));
222  }
223 
224  /**
225  * \brief wrote results in (MOAB) format, use "file_name.h5m"
226  * @param file_name file name (should always end with .h5m)
227  * @return error code
228 
229  * \ingroup mofem_fs_post_proc
230 
231  */
232  MoFEMErrorCode writeFile(const std::string file_name) {
234  // #ifdef MOAB_HDF5_PARALLEL
235  rval = postProcMesh.write_file(file_name.c_str(), "MOAB",
236  "PARALLEL=WRITE_PART");
237  CHKERRG(rval);
238  // #else
239  // #warning "No parallel HDF5, not most efficient way of writing files"
240  // if(mField.get_comm_rank()==0) {
241  // rval = postProcMesh.write_file(file_name.c_str(),"MOAB","");
242  // CHKERRG(rval);
243  // }
244  // #endif
246  }
247 };
248 
249 template <class VOLUME_ELEMENT>
251  : public PostProcTemplateOnRefineMesh<VOLUME_ELEMENT> {
252 
254 
257 
259  bool ten_nodes_post_proc_tets = true,
260  int nb_ref_levels = -1)
261  : PostProcTemplateOnRefineMesh<VOLUME_ELEMENT>(m_field),
262  tenNodesPostProcTets(ten_nodes_post_proc_tets),
263  nbOfRefLevels(nb_ref_levels) {}
264 
266  moab::Interface &moab = T::coreMesh;
267  ParallelComm *pcomm_post_proc_mesh =
268  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
269  if (pcomm_post_proc_mesh != NULL) {
270  delete pcomm_post_proc_mesh;
271  }
272  }
273 
274  // Gauss pts set on refined mesh
275  int getRule(int order) { return -1; };
276 
279 
281  return commonData;
282  }
283 
284  /** \brief Generate reference mesh on single element
285 
286  Each element is subdivided on smaller elements, i.e. a reference mesh on
287  single element is created. Nodes of such reference mesh are used as
288  integration points at which field values are calculated and to
289  each node a "moab" tag is attached to store those values.
290 
291  */
294 
295  auto get_nb_of_ref_levels_from_options = [this] {
296  if (nbOfRefLevels == -1) {
297  int max_level = 0;
298  PetscBool flg = PETSC_TRUE;
299  PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
300  "-my_max_post_proc_ref_level", &max_level, &flg);
301  return max_level;
302  } else {
303  return nbOfRefLevels;
304  }
305  return 0;
306  };
307  const int max_level = get_nb_of_ref_levels_from_options();
308 
309  moab::Core core_ref;
310  moab::Interface &moab_ref = core_ref;
311 
312  auto create_reference_element = [&moab_ref]() {
314  const double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
315  EntityHandle nodes[4];
316  for (int nn = 0; nn < 4; nn++) {
317  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
318  }
319  EntityHandle tet;
320  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
322  };
323 
324  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
325  MoFEM::Interface &m_field_ref = m_core_ref;
326 
327  auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
329  // seed ref mofem database by setting bit ref level to reference
330  // tetrahedron
331  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
332  0, 3, BitRefLevel().set(0));
333  for (int ll = 0; ll != max_level; ++ll) {
334  PetscPrintf(T::mField.get_comm(), "Refine Level %d\n", ll);
335  Range edges;
336  CHKERR m_field_ref.getInterface<BitRefManager>()
337  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
338  BitRefLevel().set(), MBEDGE, edges);
339  Range tets;
340  CHKERR m_field_ref.getInterface<BitRefManager>()
341  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
342  BitRefLevel(ll).set(), MBTET, tets);
343  // refine mesh
344  MeshRefinement *m_ref;
345  CHKERR m_field_ref.getInterface(m_ref);
346  CHKERR m_ref->add_vertices_in_the_middel_of_edges(
347  edges, BitRefLevel().set(ll + 1));
348  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
349  }
351  };
352 
353  auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
354  &m_field_ref]() {
356  for (int ll = 0; ll != max_level + 1; ++ll) {
357  Range tets;
358  CHKERR m_field_ref.getInterface<BitRefManager>()
359  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
360  BitRefLevel().set(ll), MBTET, tets);
361  if (tenNodesPostProcTets) {
362  EntityHandle meshset;
363  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
364  CHKERR moab_ref.add_entities(meshset, tets);
365  CHKERR moab_ref.convert_entities(meshset, true, false, false);
366  CHKERR moab_ref.delete_entities(&meshset, 1);
367  }
368  Range elem_nodes;
369  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
370 
371  auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
372  gauss_pts.resize(elem_nodes.size(), 4, false);
373  std::map<EntityHandle, int> little_map;
374  Range::iterator nit = elem_nodes.begin();
375  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
376  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
377  little_map[*nit] = gg;
378  }
379  gauss_pts = trans(gauss_pts);
380 
381  auto &ref_tets = levelRefTets[ll];
382  Range::iterator tit = tets.begin();
383  for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
384  const EntityHandle *conn;
385  int num_nodes;
386  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
387  if (tt == 0) {
388  // Ref tets has number of rows equal to number of tets on element,
389  // columns are number of gauss points
390  ref_tets.resize(tets.size(), num_nodes);
391  }
392  for (int nn = 0; nn != num_nodes; ++nn) {
393  ref_tets(tt, nn) = little_map[conn[nn]];
394  }
395  }
396 
397  auto &shape_functions = levelShapeFunctions[ll];
398  shape_functions.resize(elem_nodes.size(), 4);
399  CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
400  &gauss_pts(1, 0), &gauss_pts(2, 0),
401  elem_nodes.size());
402  }
404  };
405 
406  levelRefTets.resize(max_level + 1);
407  levelGaussPtsOnRefMesh.resize(max_level + 1);
408  levelShapeFunctions.resize(max_level + 1);
409 
410  CHKERR create_reference_element();
411  CHKERR refine_ref_tetrahedron();
412  CHKERR get_ref_gauss_pts_and_shape_functions();
413 
415  }
416 
417  /** \brief Set integration points
418 
419  If reference mesh is generated on single elements. This function maps
420  reference coordinates into physical coordinates and create element
421  on post-processing mesh.
422 
423  */
426 
427  auto get_element_max_dofs_order = [this]() {
428  int max_order = 0;
429  auto &dofs_multi_index = *this->dataPtr;
430  for (auto &dof : dofs_multi_index) {
431  const int dof_order = dof->getDofOrder();
432  max_order = (max_order < dof_order) ? dof_order : max_order;
433  };
434  return max_order;
435  };
436 
437  const int dof_max_order = get_element_max_dofs_order();
438  int level = (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
439  if (level > levelGaussPtsOnRefMesh.size() - 1) {
440  level = levelGaussPtsOnRefMesh.size() - 1;
441  }
442 
443  auto &level_ref_gauss_pts = levelGaussPtsOnRefMesh[level];
444  auto &level_ref_tets = levelRefTets[level];
445  auto &shape_functions = levelShapeFunctions[level];
446  T::gaussPts.resize(level_ref_gauss_pts.size1(),
447  level_ref_gauss_pts.size2(), false);
448  noalias(T::gaussPts) = level_ref_gauss_pts;
449 
450  ReadUtilIface *iface;
451  CHKERR T::postProcMesh.query_interface(iface);
452 
453  const int num_nodes = level_ref_gauss_pts.size2();
454  std::vector<double *> arrays;
455  EntityHandle startv;
456  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
457  T::mapGaussPts.resize(level_ref_gauss_pts.size2());
458  for (int gg = 0; gg != num_nodes; ++gg)
459  T::mapGaussPts[gg] = startv + gg;
460 
461 
462  Tag th;
463  int def_in_the_loop = -1;
464  CHKERR T::postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
465  th, MB_TAG_CREAT | MB_TAG_SPARSE,
466  &def_in_the_loop);
467 
468  commonData.tEts.clear();
469  const int num_el = level_ref_tets.size1();
470  const int num_nodes_on_ele = level_ref_tets.size2();
471  EntityHandle starte;
472  EntityHandle *conn;
473  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
474  starte, conn);
475  for (unsigned int tt = 0; tt != level_ref_tets.size1(); ++tt) {
476  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
477  conn[num_nodes_on_ele * tt + nn] =
478  T::mapGaussPts[level_ref_tets(tt, nn)];
479  }
480  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
481  commonData.tEts = Range(starte, starte + num_el - 1);
482  CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
483  &(T::nInTheLoop));
484 
485  EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
486  T::coords.resize(12, false);
487  {
488  const EntityHandle *conn;
489  int num_nodes;
490  T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
491  CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
492  }
493 
496  &*shape_functions.data().begin());
498  arrays[0], arrays[1], arrays[2]);
499  const double *t_coords_ele_x = &T::coords[0];
500  const double *t_coords_ele_y = &T::coords[1];
501  const double *t_coords_ele_z = &T::coords[2];
502  for (unsigned int gg = 0; gg != num_nodes; ++gg) {
504  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
505  t_coords(i) = 0;
506  for (int nn = 0; nn != 4; ++nn) {
507  t_coords(i) += t_n * t_ele_coords(i);
508  ++t_ele_coords;
509  ++t_n;
510  }
511  ++t_coords;
512  }
513 
515  }
516 
517  /** \brief Clear operators list
518 
519  Clear operators list, user can use the same mesh instance to post-process
520  different problem or the same problem with different set of post-processed
521  fields.
522 
523  */
526  T::getOpPtrVector().clear();
528  }
529 
532  moab::Interface &moab = T::coreMesh;
533  ParallelComm *pcomm_post_proc_mesh =
534  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
535  if (pcomm_post_proc_mesh != NULL) {
536  delete pcomm_post_proc_mesh;
537  }
538  CHKERR T::postProcMesh.delete_mesh();
540  }
541 
544 
545  moab::Interface &moab = T::coreMesh;
546  ParallelComm *pcomm =
547  ParallelComm::get_pcomm(&T::mField.get_moab(), MYPCOMM_INDEX);
548  ParallelComm *pcomm_post_proc_mesh =
549  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
550  if (pcomm_post_proc_mesh == NULL) {
551  pcomm_post_proc_mesh = new ParallelComm(&moab, T::mField.get_comm());
552  }
553 
554  Range edges;
555  CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
556  CHKERR T::postProcMesh.delete_entities(edges);
557  Range tris;
558  CHKERR T::postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
559  CHKERR T::postProcMesh.delete_entities(tris);
560 
561  Range tets;
562  CHKERR T::postProcMesh.get_entities_by_type(0, MBTET, tets, false);
563 
564  int rank = pcomm->rank();
565  Range::iterator tit = tets.begin();
566  for (; tit != tets.end(); tit++) {
567  CHKERR T::postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(),
568  &*tit, 1, &rank);
569  }
570 
571  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
572 
574  }
575 
576  /** \brief Add operator to post-process Hdiv field
577  */
578  MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
580  T::getOpPtrVector().push_back(
581  new OpHdivFunctions(T::postProcMesh, T::mapGaussPts, field_name));
583  }
584 
586 
587  moab::Interface &postProcMesh;
588  std::vector<EntityHandle> &mapGaussPts;
589 
590  OpHdivFunctions(moab::Interface &post_proc_mesh,
591  std::vector<EntityHandle> &map_gauss_pts,
592  const std::string field_name)
593  : VOLUME_ELEMENT::UserDataOperator(field_name,
594  T::UserDataOperator::OPCOL),
595  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
596 
597  MoFEMErrorCode doWork(int side, EntityType type,
600 
601  if (data.getIndices().size() == 0)
603 
604  std::vector<Tag> th;
605  th.resize(data.getFieldData().size());
606 
607  double def_VAL[9] = {0, 0, 0};
608 
609  switch (type) {
610  case MBTRI:
611  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
612  std::ostringstream ss;
613  ss << "HDIV_FACE_" << side << "_" << dd;
614  CHKERR postProcMesh.tag_get_handle(
615  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
616  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
617  }
618  break;
619  case MBTET:
620  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
621  std::ostringstream ss;
622  ss << "HDIV_TET_" << dd;
623  CHKERR postProcMesh.tag_get_handle(
624  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
625  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
626  }
627  break;
628  default:
629  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
630  "data inconsistency");
631  }
632 
633  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
634  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
635  CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
636  &data.getVectorN<3>(gg)(dd, 0));
637  }
638  }
639 
641  }
642  };
643 
644  private:
645  std::vector<MatrixDouble> levelShapeFunctions;
646  std::vector<MatrixDouble> levelGaussPtsOnRefMesh;
647  std::vector<ublas::matrix<int>> levelRefTets;
648  };
649 
650  /** \brief Post processing
651  * \ingroup mofem_fs_post_proc
652  */
655  MoFEM::VolumeElementForcesAndSourcesCore> {
656 
658  bool ten_nodes_post_proc_tets = true,
659  int nb_ref_levels = -1)
662 };
663 
664 // /** \deprecated Use PostPocOnRefinedMesh instead
665 // */
666 // DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
667 
668 /**
669  * \brief Postprocess on prism
670  *
671  * \ingroup mofem_fs_post_proc
672  */
675  MoFEM::FatPrismElementForcesAndSourcesCore> {
676 
678 
680  bool ten_nodes_post_proc_tets = true)
683  tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
684 
686  ParallelComm *pcomm_post_proc_mesh =
687  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
688  if (pcomm_post_proc_mesh != NULL) {
689  delete pcomm_post_proc_mesh;
690  }
691  }
692 
693  int getRuleTrianglesOnly(int order) { return -1; };
694  int getRuleThroughThickness(int order) { return -1; };
695 
696  struct PointsMap3D {
697  const int kSi;
698  const int eTa;
699  const int zEta;
700  int nN;
701  PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
702  : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
703  };
704 
705  typedef multi_index_container<
706  PointsMap3D,
707  indexed_by<ordered_unique<composite_key<
708  PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
709  member<PointsMap3D, const int, &PointsMap3D::eTa>,
710  member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
712 
714 
715  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
716  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
718 
719  std::map<EntityHandle, EntityHandle> elementsMap;
720 
723 
726 
728  return commonData;
729  }
730 };
731 
732 /**
733  * \brief Postprocess on face
734  *
735  * \ingroup mofem_fs_post_proc
736  */
738  MoFEM::FaceElementForcesAndSourcesCore> {
739 
741 
743  bool six_node_post_proc_tris = true)
745  m_field),
746  sixNodePostProcTris(six_node_post_proc_tris) {}
747 
749  ParallelComm *pcomm_post_proc_mesh =
750  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
751  if (pcomm_post_proc_mesh != NULL) {
752  delete pcomm_post_proc_mesh;
753  }
754  }
755 
756  // Gauss pts set on refined mesh
757  int getRule(int order) { return -1; };
758 
761 
762  std::map<EntityHandle, EntityHandle> elementsMap;
763 
766 
769 
771  return commonData;
772  }
773 };
774 
775 #endif //__POSTPROC_ON_REF_MESH_HPP
776 
777 /***************************************************************************/ /**
778  * \defgroup mofem_fs_post_proc Post Process
779  * \ingroup user_modules
780 ******************************************************************************/
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)
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
MoFEMErrorCode preProcess()
const int kSi
PostProcTemplateOnRefineMesh< VOLUME_ELEMENT > T
std::vector< EntityHandle > & mapGaussPts
CommonData commonData
int getRule(int order)
VectorDouble * vAluesPtr
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
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:500
PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
Vec V
PostProcTemplateOnRefineMesh(MoFEM::Interface &m_field)
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:476
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:543
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
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:619
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:507
CommonData & commonData
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
OpHdivFunctions(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
MoFEMErrorCode postProcess()
function is run at the end of loop
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Postprocess on prism.
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
Vec V
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
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
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
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
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
moab::Interface & postProcMesh
std::vector< MatrixDouble > levelShapeFunctions
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:66
moab::Core coreMesh
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
int nN
DataForcesAndSourcesCore::EntData EntData
VectorDouble * vAluesPtr
virtual ~PostProcTemplateVolumeOnRefinedMesh()
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.
virtual ~PostProcFaceOnRefinedMesh()
bool sixNodePostProcTris
#define CHKERR
Inline error check.
Definition: definitions.h:595
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
std::vector< EntityHandle > mapGaussPts
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRuleTrianglesOnly(int order)
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
ForcesAndSourcesCore::UserDataOperator UserDataOperator
virtual ~PostProcFatPrismOnRefinedMesh()
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:284
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Postprocess on face.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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:406
CommonData commonData
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
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
Generic post-processing class.
const int order
approximation order
PointsMap3D_multiIndex pointsMap
MoFEMErrorCode setGaussPts(int order)
Set integration points.