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)
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)
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_middle_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(), level_ref_gauss_pts.size2(),
447  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  Tag th;
462  int def_in_the_loop = -1;
463  CHKERR T::postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
464  th, MB_TAG_CREAT | MB_TAG_SPARSE,
465  &def_in_the_loop);
466 
467  commonData.tEts.clear();
468  const int num_el = level_ref_tets.size1();
469  const int num_nodes_on_ele = level_ref_tets.size2();
470  EntityHandle starte;
471  EntityHandle *conn;
472  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
473  starte, conn);
474  for (unsigned int tt = 0; tt != level_ref_tets.size1(); ++tt) {
475  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
476  conn[num_nodes_on_ele * tt + nn] =
477  T::mapGaussPts[level_ref_tets(tt, nn)];
478  }
479  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
480  commonData.tEts = Range(starte, starte + num_el - 1);
481  CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
482  &(T::nInTheLoop));
483 
484  EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
485  T::coords.resize(12, false);
486  {
487  const EntityHandle *conn;
488  int num_nodes;
489  T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
490  CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
491  }
492 
495  &*shape_functions.data().begin());
497  arrays[0], arrays[1], arrays[2]);
498  const double *t_coords_ele_x = &T::coords[0];
499  const double *t_coords_ele_y = &T::coords[1];
500  const double *t_coords_ele_z = &T::coords[2];
501  for (unsigned int gg = 0; gg != num_nodes; ++gg) {
503  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
504  t_coords(i) = 0;
505  for (int nn = 0; nn != 4; ++nn) {
506  t_coords(i) += t_n * t_ele_coords(i);
507  ++t_ele_coords;
508  ++t_n;
509  }
510  ++t_coords;
511  }
512 
514  }
515 
516  /** \brief Clear operators list
517 
518  Clear operators list, user can use the same mesh instance to post-process
519  different problem or the same problem with different set of post-processed
520  fields.
521 
522  */
525  T::getOpPtrVector().clear();
527  }
528 
531  moab::Interface &moab = T::coreMesh;
532  ParallelComm *pcomm_post_proc_mesh =
533  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
534  if (pcomm_post_proc_mesh != NULL) {
535  delete pcomm_post_proc_mesh;
536  }
537  CHKERR T::postProcMesh.delete_mesh();
539  }
540 
543 
544  moab::Interface &moab = T::coreMesh;
545  ParallelComm *pcomm =
546  ParallelComm::get_pcomm(&T::mField.get_moab(), MYPCOMM_INDEX);
547  ParallelComm *pcomm_post_proc_mesh =
548  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
549  if (pcomm_post_proc_mesh == NULL) {
550  pcomm_post_proc_mesh = new ParallelComm(&moab, T::mField.get_comm());
551  }
552 
553  Range edges;
554  CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
555  CHKERR T::postProcMesh.delete_entities(edges);
556  Range tris;
557  CHKERR T::postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
558  CHKERR T::postProcMesh.delete_entities(tris);
559 
560  Range tets;
561  CHKERR T::postProcMesh.get_entities_by_type(0, MBTET, tets, false);
562 
563  int rank = pcomm->rank();
564  Range::iterator tit = tets.begin();
565  for (; tit != tets.end(); tit++) {
566  CHKERR T::postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(),
567  &*tit, 1, &rank);
568  }
569 
570  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
571 
573  }
574 
575  /** \brief Add operator to post-process Hdiv field
576  */
577  MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
579  T::getOpPtrVector().push_back(
580  new OpHdivFunctions(T::postProcMesh, T::mapGaussPts, field_name));
582  }
583 
585 
586  moab::Interface &postProcMesh;
587  std::vector<EntityHandle> &mapGaussPts;
588 
589  OpHdivFunctions(moab::Interface &post_proc_mesh,
590  std::vector<EntityHandle> &map_gauss_pts,
591  const std::string field_name)
592  : VOLUME_ELEMENT::UserDataOperator(field_name,
593  T::UserDataOperator::OPCOL),
594  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
595 
596  MoFEMErrorCode doWork(int side, EntityType type,
599 
600  if (data.getIndices().size() == 0)
602 
603  std::vector<Tag> th;
604  th.resize(data.getFieldData().size());
605 
606  double def_VAL[9] = {0, 0, 0};
607 
608  switch (type) {
609  case MBTRI:
610  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
611  std::ostringstream ss;
612  ss << "HDIV_FACE_" << side << "_" << dd;
613  CHKERR postProcMesh.tag_get_handle(
614  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
615  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
616  }
617  break;
618  case MBTET:
619  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
620  std::ostringstream ss;
621  ss << "HDIV_TET_" << dd;
622  CHKERR postProcMesh.tag_get_handle(
623  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
624  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
625  }
626  break;
627  default:
628  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
629  "data inconsistency");
630  }
631 
632  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
633  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
634  CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
635  &data.getVectorN<3>(gg)(dd, 0));
636  }
637  }
638 
640  }
641  };
642 
643 private:
644  std::vector<MatrixDouble> levelShapeFunctions;
645  std::vector<MatrixDouble> levelGaussPtsOnRefMesh;
646  std::vector<ublas::matrix<int>> levelRefTets;
647 };
648 
649 /** \brief Post processing
650  * \ingroup mofem_fs_post_proc
651  */
654  MoFEM::VolumeElementForcesAndSourcesCore> {
655 
657  bool ten_nodes_post_proc_tets = true,
658  int nb_ref_levels = -1)
661 };
662 
663 // /** \deprecated Use PostPocOnRefinedMesh instead
664 // */
665 // DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
666 
667 /**
668  * \brief Postprocess on prism
669  *
670  * \ingroup mofem_fs_post_proc
671  */
674  MoFEM::FatPrismElementForcesAndSourcesCore> {
675 
677 
679  bool ten_nodes_post_proc_tets = true)
682  tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
683 
685  ParallelComm *pcomm_post_proc_mesh =
686  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
687  if (pcomm_post_proc_mesh != NULL) {
688  delete pcomm_post_proc_mesh;
689  }
690  }
691 
692  int getRuleTrianglesOnly(int order) { return -1; };
693  int getRuleThroughThickness(int order) { return -1; };
694 
695  struct PointsMap3D {
696  const int kSi;
697  const int eTa;
698  const int zEta;
699  int nN;
700  PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
701  : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
702  };
703 
704  typedef multi_index_container<
705  PointsMap3D,
706  indexed_by<ordered_unique<composite_key<
707  PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
708  member<PointsMap3D, const int, &PointsMap3D::eTa>,
709  member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
711 
713 
714  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
715  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
717 
718  std::map<EntityHandle, EntityHandle> elementsMap;
719 
722 
725 
727  return commonData;
728  }
729 };
730 
731 /**
732  * \brief Postprocess on face
733  *
734  * \ingroup mofem_fs_post_proc
735  */
737  MoFEM::FaceElementForcesAndSourcesCore> {
738 
740 
742  bool six_node_post_proc_tris = true)
744  m_field),
745  sixNodePostProcTris(six_node_post_proc_tris) {}
746 
748  ParallelComm *pcomm_post_proc_mesh =
749  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
750  if (pcomm_post_proc_mesh != NULL) {
751  delete pcomm_post_proc_mesh;
752  }
753  }
754 
755  // Gauss pts set on refined mesh
756  int getRule(int order) { return -1; };
757 
760 
761  std::map<EntityHandle, EntityHandle> elementsMap;
762 
765 
768 
770  return commonData;
771  }
772 
775 
776  moab::Interface &postProcMesh;
777  std::vector<EntityHandle> &mapGaussPts;
778  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
779  const std::string tagName;
780  const bool saveOnTag;
781 
782  const std::string feVolName;
783  boost::shared_ptr<MatrixDouble> gradMatPtr;
784 
786  moab::Interface &post_proc_mesh,
787  std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
788  const std::string tag_name,
789  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
790  const std::string vol_fe_name,
791  boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag)
793  field_name, UserDataOperator::OPCOL),
794  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
795  tagName(tag_name), sideOpFe(side_fe), feVolName(vol_fe_name),
796  gradMatPtr(grad_mat_ptr), saveOnTag(save_on_tag) {}
797 
798  MoFEMErrorCode doWork(int side, EntityType type,
800  };
801 
803  const std::string field_name, const std::string vol_fe_name,
804  boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
805  bool save_on_tag = true) {
807 
808  if (!grad_mat_ptr)
809  grad_mat_ptr = boost::make_shared<MatrixDouble>();
810 
811  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
812  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
813 
814  // check number of coefficients
815  auto field_ptr = mField.get_field_structure(field_name);
816  const int nb_coefficients = field_ptr->getNbOfCoeffs();
817 
818  switch (nb_coefficients) {
819  case 1:
820  my_side_fe->getOpPtrVector().push_back(
821  new OpCalculateScalarFieldGradient<3>(field_name, grad_mat_ptr));
822  break;
823  case 3:
824  my_side_fe->getOpPtrVector().push_back(
825  new OpCalculateVectorFieldGradient<3, 3>(field_name, grad_mat_ptr));
826  break;
827  default:
828  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
829  "field with that number of coefficients is not implemented");
830  }
831 
832  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
834  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
835  my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
836 
838  }
839 };
840 
841 #endif //__POSTPROC_ON_REF_MESH_HPP
842 
843 /**
844  * \defgroup mofem_fs_post_proc Post Process
845  * \ingroup user_modules
846  **/
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
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
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:501
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:477
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
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:620
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:508
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()
PostProcCommonOnRefMesh::CommonDataForVolume CommonData
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Postprocess on prism.
const std::string feVolName
virtual PostProcCommonOnRefMesh::CommonData & getCommonData()
Vec V
boost::shared_ptr< MatrixDouble > gradMatPtr
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
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
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
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideOpFe
#define CHKERR
Inline error check.
Definition: definitions.h:596
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode preProcess()
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:285
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.
std::vector< EntityHandle > & mapGaussPts
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:407
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:72
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 bool saveOnTag
const int order
approximation order
PointsMap3D_multiIndex pointsMap
MoFEMErrorCode setGaussPts(int order)
Set integration points.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)