v0.9.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 
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  CHKERR postProcMesh.write_file(file_name.c_str(), "MOAB",
236  "PARALLEL=WRITE_PART");
237  // #else
238  // #warning "No parallel HDF5, not most efficient way of writing files"
239  // if(mField.get_comm_rank()==0) {
240  // rval = postProcMesh.write_file(file_name.c_str(),"MOAB","");
241  // CHKERRG(rval);
242  // }
243  // #endif
245  }
246 };
247 
248 template <class VOLUME_ELEMENT>
250  : public PostProcTemplateOnRefineMesh<VOLUME_ELEMENT> {
251 
253 
256 
258  bool ten_nodes_post_proc_tets = true,
259  int nb_ref_levels = -1)
260  : PostProcTemplateOnRefineMesh<VOLUME_ELEMENT>(m_field),
261  tenNodesPostProcTets(ten_nodes_post_proc_tets),
262  nbOfRefLevels(nb_ref_levels) {}
263 
265  moab::Interface &moab = T::coreMesh;
266  ParallelComm *pcomm_post_proc_mesh =
267  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
268  if (pcomm_post_proc_mesh != NULL) {
269  delete pcomm_post_proc_mesh;
270  }
271  }
272 
273  // Gauss pts set on refined mesh
274  int getRule(int order) { return -1; };
275 
278 
280  return commonData;
281  }
282 
283  /** \brief Generate reference mesh on single element
284 
285  Each element is subdivided on smaller elements, i.e. a reference mesh on
286  single element is created. Nodes of such reference mesh are used as
287  integration points at which field values are calculated and to
288  each node a "moab" tag is attached to store those values.
289 
290  */
293 
294  auto get_nb_of_ref_levels_from_options = [this] {
295  if (nbOfRefLevels == -1) {
296  int max_level = 0;
297  PetscBool flg = PETSC_TRUE;
298  PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
299  "-my_max_post_proc_ref_level", &max_level, &flg);
300  return max_level;
301  } else {
302  return nbOfRefLevels;
303  }
304  return 0;
305  };
306  const int max_level = get_nb_of_ref_levels_from_options();
307 
308  moab::Core core_ref;
309  moab::Interface &moab_ref = core_ref;
310 
311  auto create_reference_element = [&moab_ref]() {
313  const double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
314  EntityHandle nodes[4];
315  for (int nn = 0; nn < 4; nn++) {
316  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
317  }
318  EntityHandle tet;
319  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
321  };
322 
323  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
324  MoFEM::Interface &m_field_ref = m_core_ref;
325 
326  auto refine_ref_tetrahedron = [this, &m_field_ref, max_level]() {
328  // seed ref mofem database by setting bit ref level to reference
329  // tetrahedron
330  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
331  0, 3, BitRefLevel().set(0));
332  for (int ll = 0; ll != max_level; ++ll) {
333  PetscPrintf(T::mField.get_comm(), "Refine Level %d\n", ll);
334  Range edges;
335  CHKERR m_field_ref.getInterface<BitRefManager>()
336  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
337  BitRefLevel().set(), MBEDGE, edges);
338  Range tets;
339  CHKERR m_field_ref.getInterface<BitRefManager>()
340  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
341  BitRefLevel(ll).set(), MBTET, tets);
342  // refine mesh
343  MeshRefinement *m_ref;
344  CHKERR m_field_ref.getInterface(m_ref);
345  CHKERR m_ref->add_vertices_in_the_middle_of_edges(
346  edges, BitRefLevel().set(ll + 1));
347  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
348  }
350  };
351 
352  auto get_ref_gauss_pts_and_shape_functions = [this, max_level, &moab_ref,
353  &m_field_ref]() {
355  for (int ll = 0; ll != max_level + 1; ++ll) {
356  Range tets;
357  CHKERR m_field_ref.getInterface<BitRefManager>()
358  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
359  BitRefLevel().set(ll), MBTET, tets);
360  if (tenNodesPostProcTets) {
361  EntityHandle meshset;
362  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
363  CHKERR moab_ref.add_entities(meshset, tets);
364  CHKERR moab_ref.convert_entities(meshset, true, false, false);
365  CHKERR moab_ref.delete_entities(&meshset, 1);
366  }
367  Range elem_nodes;
368  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
369 
370  auto &gauss_pts = levelGaussPtsOnRefMesh[ll];
371  gauss_pts.resize(elem_nodes.size(), 4, false);
372  std::map<EntityHandle, int> little_map;
373  Range::iterator nit = elem_nodes.begin();
374  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
375  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
376  little_map[*nit] = gg;
377  }
378  gauss_pts = trans(gauss_pts);
379 
380  auto &ref_tets = levelRefTets[ll];
381  Range::iterator tit = tets.begin();
382  for (int tt = 0; tit != tets.end(); ++tit, ++tt) {
383  const EntityHandle *conn;
384  int num_nodes;
385  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
386  if (tt == 0) {
387  // Ref tets has number of rows equal to number of tets on element,
388  // columns are number of gauss points
389  ref_tets.resize(tets.size(), num_nodes);
390  }
391  for (int nn = 0; nn != num_nodes; ++nn) {
392  ref_tets(tt, nn) = little_map[conn[nn]];
393  }
394  }
395 
396  auto &shape_functions = levelShapeFunctions[ll];
397  shape_functions.resize(elem_nodes.size(), 4);
398  CHKERR ShapeMBTET(&*shape_functions.data().begin(), &gauss_pts(0, 0),
399  &gauss_pts(1, 0), &gauss_pts(2, 0),
400  elem_nodes.size());
401  }
403  };
404 
405  levelRefTets.resize(max_level + 1);
406  levelGaussPtsOnRefMesh.resize(max_level + 1);
407  levelShapeFunctions.resize(max_level + 1);
408 
409  CHKERR create_reference_element();
410  CHKERR refine_ref_tetrahedron();
411  CHKERR get_ref_gauss_pts_and_shape_functions();
412 
414  }
415 
416  /** \brief Set integration points
417 
418  If reference mesh is generated on single elements. This function maps
419  reference coordinates into physical coordinates and create element
420  on post-processing mesh.
421 
422  */
425 
426  auto get_element_max_dofs_order = [this]() {
427  int max_order = 0;
428  auto &dofs_multi_index = *this->dataPtr;
429  for (auto &dof : dofs_multi_index) {
430  const int dof_order = dof->getDofOrder();
431  max_order = (max_order < dof_order) ? dof_order : max_order;
432  };
433  return max_order;
434  };
435 
436  const int dof_max_order = get_element_max_dofs_order();
437  size_t level = (dof_max_order > 0) ? (dof_max_order - 1) / 2 : 0;
438  if (level > (levelGaussPtsOnRefMesh.size() - 1))
439  level = levelGaussPtsOnRefMesh.size() - 1;
440 
441  auto &level_ref_gauss_pts = levelGaussPtsOnRefMesh[level];
442  auto &level_ref_tets = levelRefTets[level];
443  auto &shape_functions = levelShapeFunctions[level];
444  T::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
445  false);
446  noalias(T::gaussPts) = level_ref_gauss_pts;
447 
448  ReadUtilIface *iface;
449  CHKERR T::postProcMesh.query_interface(iface);
450 
451  const int num_nodes = level_ref_gauss_pts.size2();
452  std::vector<double *> arrays;
453  EntityHandle startv;
454  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
455  T::mapGaussPts.resize(level_ref_gauss_pts.size2());
456  for (int gg = 0; gg != num_nodes; ++gg)
457  T::mapGaussPts[gg] = startv + gg;
458 
459  Tag th;
460  int def_in_the_loop = -1;
461  CHKERR T::postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
462  th, MB_TAG_CREAT | MB_TAG_SPARSE,
463  &def_in_the_loop);
464 
465  commonData.tEts.clear();
466  const int num_el = level_ref_tets.size1();
467  const int num_nodes_on_ele = level_ref_tets.size2();
468  EntityHandle starte;
469  EntityHandle *conn;
470  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
471  starte, conn);
472  for (unsigned int tt = 0; tt != level_ref_tets.size1(); ++tt) {
473  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
474  conn[num_nodes_on_ele * tt + nn] =
475  T::mapGaussPts[level_ref_tets(tt, nn)];
476  }
477  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
478  commonData.tEts = Range(starte, starte + num_el - 1);
479  CHKERR T::postProcMesh.tag_clear_data(th, commonData.tEts,
480  &(T::nInTheLoop));
481 
482  EntityHandle fe_ent = T::numeredEntFiniteElementPtr->getEnt();
483  T::coords.resize(12, false);
484  {
485  const EntityHandle *conn;
486  int num_nodes;
487  T::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true);
488  CHKERR T::mField.get_moab().get_coords(conn, num_nodes, &T::coords[0]);
489  }
490 
493  &*shape_functions.data().begin());
495  arrays[0], arrays[1], arrays[2]);
496  const double *t_coords_ele_x = &T::coords[0];
497  const double *t_coords_ele_y = &T::coords[1];
498  const double *t_coords_ele_z = &T::coords[2];
499  for (int gg = 0; gg != num_nodes; ++gg) {
501  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
502  t_coords(i) = 0;
503  for (int nn = 0; nn != 4; ++nn) {
504  t_coords(i) += t_n * t_ele_coords(i);
505  ++t_ele_coords;
506  ++t_n;
507  }
508  ++t_coords;
509  }
510 
512  }
513 
514  /** \brief Clear operators list
515 
516  Clear operators list, user can use the same mesh instance to post-process
517  different problem or the same problem with different set of post-processed
518  fields.
519 
520  */
523  T::getOpPtrVector().clear();
525  }
526 
529  moab::Interface &moab = T::coreMesh;
530  ParallelComm *pcomm_post_proc_mesh =
531  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
532  if (pcomm_post_proc_mesh != NULL) {
533  delete pcomm_post_proc_mesh;
534  }
535  CHKERR T::postProcMesh.delete_mesh();
537  }
538 
541 
542  moab::Interface &moab = T::coreMesh;
543  ParallelComm *pcomm =
544  ParallelComm::get_pcomm(&T::mField.get_moab(), MYPCOMM_INDEX);
545  ParallelComm *pcomm_post_proc_mesh =
546  ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
547  if (pcomm_post_proc_mesh == NULL) {
548  pcomm_post_proc_mesh = new ParallelComm(&moab, T::mField.get_comm());
549  }
550 
551  Range edges;
552  CHKERR T::postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
553  CHKERR T::postProcMesh.delete_entities(edges);
554  Range tris;
555  CHKERR T::postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
556  CHKERR T::postProcMesh.delete_entities(tris);
557 
558  Range tets;
559  CHKERR T::postProcMesh.get_entities_by_type(0, MBTET, tets, false);
560 
561  int rank = pcomm->rank();
562  Range::iterator tit = tets.begin();
563  for (; tit != tets.end(); tit++) {
564  CHKERR T::postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(),
565  &*tit, 1, &rank);
566  }
567 
568  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
569 
571  }
572 
573  /** \brief Add operator to post-process Hdiv field
574  */
575  MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name) {
577  T::getOpPtrVector().push_back(
578  new OpHdivFunctions(T::postProcMesh, T::mapGaussPts, field_name));
580  }
581 
583 
584  moab::Interface &postProcMesh;
585  std::vector<EntityHandle> &mapGaussPts;
586 
587  OpHdivFunctions(moab::Interface &post_proc_mesh,
588  std::vector<EntityHandle> &map_gauss_pts,
589  const std::string field_name)
590  : VOLUME_ELEMENT::UserDataOperator(field_name,
591  T::UserDataOperator::OPCOL),
592  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
593 
594  MoFEMErrorCode doWork(int side, EntityType type,
597 
598  if (data.getIndices().size() == 0)
600 
601  std::vector<Tag> th;
602  th.resize(data.getFieldData().size());
603 
604  double def_VAL[9] = {0, 0, 0};
605 
606  switch (type) {
607  case MBTRI:
608  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
609  std::ostringstream ss;
610  ss << "HDIV_FACE_" << side << "_" << dd;
611  CHKERR postProcMesh.tag_get_handle(
612  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
613  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
614  }
615  break;
616  case MBTET:
617  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
618  std::ostringstream ss;
619  ss << "HDIV_TET_" << dd;
620  CHKERR postProcMesh.tag_get_handle(
621  ss.str().c_str(), 3, MB_TYPE_DOUBLE, th[dd],
622  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
623  }
624  break;
625  default:
626  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
627  "data inconsistency");
628  }
629 
630  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
631  for (unsigned int dd = 0; dd < data.getN().size2() / 3; dd++) {
632  CHKERR postProcMesh.tag_set_data(th[dd], &mapGaussPts[gg], 1,
633  &data.getVectorN<3>(gg)(dd, 0));
634  }
635  }
636 
638  }
639  };
640 
641 private:
642  std::vector<MatrixDouble> levelShapeFunctions;
643  std::vector<MatrixDouble> levelGaussPtsOnRefMesh;
644  std::vector<ublas::matrix<int>> levelRefTets;
645 };
646 
647 /** \brief Post processing
648  * \ingroup mofem_fs_post_proc
649  */
652  MoFEM::VolumeElementForcesAndSourcesCore> {
653 
655  bool ten_nodes_post_proc_tets = true,
656  int nb_ref_levels = -1)
659 };
660 
661 // /** \deprecated Use PostPocOnRefinedMesh instead
662 // */
663 // DEPRECATED typedef PostProcVolumeOnRefinedMesh PostPocOnRefinedMesh;
664 
665 /**
666  * \brief Postprocess on prism
667  *
668  * \ingroup mofem_fs_post_proc
669  */
672  MoFEM::FatPrismElementForcesAndSourcesCore> {
673 
675 
677  bool ten_nodes_post_proc_tets = true)
680  tenNodesPostProcTets(ten_nodes_post_proc_tets) {}
681 
683  ParallelComm *pcomm_post_proc_mesh =
684  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
685  if (pcomm_post_proc_mesh != NULL) {
686  delete pcomm_post_proc_mesh;
687  }
688  }
689 
690  int getRuleTrianglesOnly(int order) { return -1; };
691  int getRuleThroughThickness(int order) { return -1; };
692 
693  struct PointsMap3D {
694  const int kSi;
695  const int eTa;
696  const int zEta;
697  int nN;
698  PointsMap3D(const int ksi, const int eta, const int zeta, const int nn)
699  : kSi(ksi), eTa(eta), zEta(zeta), nN(nn) {}
700  };
701 
702  typedef multi_index_container<
703  PointsMap3D,
704  indexed_by<ordered_unique<composite_key<
705  PointsMap3D, member<PointsMap3D, const int, &PointsMap3D::kSi>,
706  member<PointsMap3D, const int, &PointsMap3D::eTa>,
707  member<PointsMap3D, const int, &PointsMap3D::zEta>>>>>
709 
711 
712  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
713  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
715 
716  std::map<EntityHandle, EntityHandle> elementsMap;
717 
720 
723 
725  return commonData;
726  }
727 };
728 
729 /**
730  * \brief Postprocess on face
731  *
732  * \ingroup mofem_fs_post_proc
733  */
735  MoFEM::FaceElementForcesAndSourcesCore> {
736 
738 
740  bool six_node_post_proc_tris = true)
742  m_field),
743  sixNodePostProcTris(six_node_post_proc_tris) {}
744 
746  ParallelComm *pcomm_post_proc_mesh =
747  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
748  if (pcomm_post_proc_mesh != NULL) {
749  delete pcomm_post_proc_mesh;
750  }
751  }
752 
753  // Gauss pts set on refined mesh
754  int getRule(int order) { return -1; };
755 
757  MoFEMErrorCode setGaussPts(int order);
758 
759  std::map<EntityHandle, EntityHandle> elementsMap;
760 
763 
766 
768  return commonData;
769  }
770 
773 
774  moab::Interface &postProcMesh;
775  std::vector<EntityHandle> &mapGaussPts;
776  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideOpFe;
777  const std::string feVolName;
778  boost::shared_ptr<MatrixDouble> gradMatPtr;
779  const std::string tagName;
780  const bool saveOnTag;
781 
783  moab::Interface &post_proc_mesh,
784  std::vector<EntityHandle> &map_gauss_pts, const std::string field_name,
785  const std::string tag_name,
786  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe,
787  const std::string vol_fe_name,
788  boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag)
790  field_name, UserDataOperator::OPCOL),
791  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
792  sideOpFe(side_fe), feVolName(vol_fe_name),
793  gradMatPtr(grad_mat_ptr), tagName(tag_name), saveOnTag(save_on_tag) {}
794 
795  MoFEMErrorCode doWork(int side, EntityType type,
797  };
798 
800  const std::string field_name, const std::string vol_fe_name,
801  boost::shared_ptr<MatrixDouble> grad_mat_ptr = nullptr,
802  bool save_on_tag = true) {
804 
805  if (!grad_mat_ptr)
806  grad_mat_ptr = boost::make_shared<MatrixDouble>();
807 
808  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
809  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
810 
811  // check number of coefficients
812  auto field_ptr = mField.get_field_structure(field_name);
813  const int nb_coefficients = field_ptr->getNbOfCoeffs();
814 
815  switch (nb_coefficients) {
816  case 1:
817  my_side_fe->getOpPtrVector().push_back(
818  new OpCalculateScalarFieldGradient<3>(field_name, grad_mat_ptr));
819  break;
820  case 3:
821  my_side_fe->getOpPtrVector().push_back(
822  new OpCalculateVectorFieldGradient<3, 3>(field_name, grad_mat_ptr));
823  break;
824  default:
825  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
826  "field with that number of coefficients is not implemented");
827  }
828 
829  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
831  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
832  my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
833 
835  }
836 
837  private:
840 
841 };
842 
843 #endif //__POSTPROC_ON_REF_MESH_HPP
844 
845 /**
846  * \defgroup mofem_fs_post_proc Post Process
847  * \ingroup user_modules
848  **/
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)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
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
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
MatrixDouble gaussPtsTri
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, const std::string tag_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addHdivFunctionsPostProc(const std::string field_name)
Add operator to post-process Hdiv field.
std::vector< ublas::matrix< int > > levelRefTets
CommonData & commonData
MoFEMErrorCode generateReferenceElementMesh()
bool tenNodesPostProcTets
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
MatrixDouble gaussPtsQuad
moab::Interface & postProcMesh
std::vector< MatrixDouble > levelShapeFunctions
DataForcesAndSourcesCore::EntData EntData
PostProcTemplateVolumeOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true, int nb_ref_levels=-1)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
moab::Core coreMesh
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
int nN
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
PointsMap3D_multiIndex pointsMap
MoFEMErrorCode setGaussPts(int order)
Set integration points.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)