v0.14.0
PostProcOnRefMesh.cpp
Go to the documentation of this file.
1 /** \file PostProcOnRefMesh.cpp
2  * \brief Postprocess fields on refined mesh made for 10 Node tets
3  *
4  * Create refined mesh, without enforcing continuity between element. Calculate
5  * field values on nodes of that mesh.
6  * \ingroup mofem_fs_post_proc
7  */
8 
9 
10 
11 #include <MoFEM.hpp>
12 using namespace MoFEM;
13 using namespace boost::numeric;
14 #include <PostProcOnRefMesh.hpp>
15 
16 MoFEMErrorCode PostProcCommonOnRefMesh::OpGetFieldValues::doWork(
17  int side, EntityType type, EntitiesFieldData::EntData &data) {
18  MoFEMFunctionBegin;
19 
20  if (data.getFieldData().size() == 0)
21  MoFEMFunctionReturnHot(0);
22 
23  if (V) {
24  vAlues.resize(data.getFieldData().size());
25  double *a;
26  CHKERR VecGetArray(V, &a);
27  VectorDofs::iterator it, hi_it;
28  it = data.getFieldDofs().begin();
29  hi_it = data.getFieldDofs().end();
30  auto row_dofs = getFEMethod()->getRowDofsPtr();
31  for (int ii = 0; it != hi_it; it++, ii++) {
32  int local_idx = row_dofs->find((*it)->getLocalUniqueId())
33  ->get()
34  ->getPetscLocalDofIdx();
35  vAlues[ii] = a[local_idx];
36  }
37  CHKERR VecRestoreArray(V, &a);
38  vAluesPtr = &vAlues;
39  } else {
40  vAluesPtr = &data.getFieldData();
41  }
42 
43  auto dof_ptr = data.getFieldDofs()[0];
44  int rank = dof_ptr->getNbOfCoeffs();
45 
46  int tag_length = rank;
47  FieldSpace space = dof_ptr->getSpace();
48  switch (space) {
49  case L2:
50  case H1:
51  break;
52  case HCURL:
53  case HDIV:
54  tag_length *= 3;
55  break;
56  default:
57  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
58  "field with that space is not implemented");
59  }
60 
61  if (tag_length > 1 && tag_length < 3)
62  tag_length = 3;
63  else if (tag_length > 3 && tag_length < 9)
64  tag_length = 9;
65 
66  double def_VAL[tag_length];
67  bzero(def_VAL, tag_length * sizeof(double));
68  Tag th;
69  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
70  MB_TYPE_DOUBLE, th,
71  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
72 
73  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
74  // no need for L2
75  const void *tags_ptr[mapGaussPts.size()];
76  int nb_gauss_pts = data.getN().size1();
77  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
78  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79  "data inconsistency %d!=%d", mapGaussPts.size(), nb_gauss_pts);
80  }
81 
82  auto set_float_precision = [](const double x) {
83  if (std::abs(x) < std::numeric_limits<float>::epsilon())
84  return 0.;
85  else
86  return x;
87  };
88 
89  switch (space) {
90  case H1:
91  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
92  if (type == MBVERTEX) {
93  for (int gg = 0; gg < nb_gauss_pts; gg++) {
94  rval = postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
95  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
96  (commonData.fieldMap[rowFieldName])[gg].clear();
97  }
98  }
99  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
100  tags_ptr);
101  for (int gg = 0; gg < nb_gauss_pts; gg++) {
102  for (int rr = 0; rr < rank; rr++) {
103  const double val =
104  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
105  &((*vAluesPtr)[rr]), rank);
106  (commonData.fieldMap[rowFieldName])[gg][rr] =
107  ((double *)tags_ptr[gg])[rr] += set_float_precision(val);
108  }
109  }
110  break;
111  case L2:
112  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
113  for (int gg = 0; gg < nb_gauss_pts; gg++) {
114  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
115  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
116  (commonData.fieldMap[rowFieldName])[gg].clear();
117  }
118  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
119  tags_ptr);
120  for (int gg = 0; gg < nb_gauss_pts; gg++) {
121  bzero((double *)tags_ptr[gg], sizeof(double) * tag_length);
122  for (int rr = 0; rr < rank; rr++) {
123  const double val =
124  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
125  &((*vAluesPtr)[rr]), rank);
126  (commonData.fieldMap[rowFieldName])[gg][rr] =
127  ((double *)tags_ptr[gg])[rr] = set_float_precision(val);
128  }
129  }
130  break;
131  case HCURL:
132  // FIXME: fieldMap not set
133  if (type == MBEDGE && side == 0) {
134  for (int gg = 0; gg < nb_gauss_pts; gg++) {
135  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
136  }
137  }
138  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
139  tags_ptr);
140  {
141  FTensor::Index<'i', 3> i;
142  auto t_n_hcurl = data.getFTensor1N<3>();
143  for (int gg = 0; gg != nb_gauss_pts; gg++) {
144  double *ptr = &((double *)tags_ptr[gg])[0];
145  int ll = 0;
146  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
147  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
148  for (int rr = 0; rr != rank; rr++) {
149  const double dof_val = (*vAluesPtr)[ll * rank + rr];
150  t_tag_val(i) += set_float_precision(dof_val) * t_n_hcurl(i);
151  ++t_tag_val;
152  }
153  ++t_n_hcurl;
154  }
155  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
156  ++t_n_hcurl;
157  }
158  }
159  }
160  break;
161  case HDIV:
162  // FIXME: fieldMap not set
163  if (moab::CN::Dimension(type) == 2 && side == 0) {
164  for (int gg = 0; gg < nb_gauss_pts; gg++) {
165  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
166  }
167  }
168  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
169  tags_ptr);
170  {
171  FTensor::Index<'i', 3> i;
172  auto t_n_hdiv = data.getFTensor1N<3>();
173  for (int gg = 0; gg != nb_gauss_pts; gg++) {
174  double *ptr = &((double *)tags_ptr[gg])[0];
175  int ll = 0;
176  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
177  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
178  for (int rr = 0; rr != rank; rr++) {
179  const double dof_val = (*vAluesPtr)[ll * rank + rr];
180  t_tag_val(i) += set_float_precision(dof_val) * t_n_hdiv(i);
181  ++t_tag_val;
182  }
183  ++t_n_hdiv;
184  }
185  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
186  ++t_n_hdiv;
187  }
188  }
189  }
190  break;
191  default:
192  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
193  "field with that space is not implemented");
194  }
195 
196  MoFEMFunctionReturn(0);
197 }
198 
199 MoFEMErrorCode PostProcCommonOnRefMesh::OpGetFieldGradientValues::doWork(
200  int side, EntityType type, EntitiesFieldData::EntData &data) {
201  MoFEMFunctionBegin;
202 
203  if (data.getFieldData().size() == 0)
204  MoFEMFunctionReturnHot(0);
205  if (V) {
206  vAlues.resize(data.getFieldData().size());
207  double *a;
208  CHKERR VecGetArray(V, &a);
209  VectorDofs::iterator it, hi_it;
210  it = data.getFieldDofs().begin();
211  hi_it = data.getFieldDofs().end();
212  auto row_dofs = getFEMethod()->getRowDofsPtr();
213  for (int ii = 0; it != hi_it; it++, ii++) {
214  int local_idx = row_dofs->find((*it)->getLocalUniqueId())
215  ->get()
216  ->getPetscLocalDofIdx();
217  vAlues[ii] = a[local_idx];
218  }
219  CHKERR VecRestoreArray(V, &a);
220  vAluesPtr = &vAlues;
221  } else {
222  vAluesPtr = &data.getFieldData();
223  }
224 
225  auto dof_ptr = data.getFieldDofs()[0];
226  int rank = dof_ptr->getNbOfCoeffs();
227 
228  FieldSpace space = dof_ptr->getSpace();
229  int space_dim = spaceDim;
230  if (space == HCURL || space == HDIV)
231  space_dim = 3;
232 
233  int tag_length = rank * space_dim;
234  // FieldSpace space = dof_ptr->getSpace();
235  switch (space) {
236  case L2:
237  case H1:
238  break;
239  case HCURL:
240  case HDIV:
241  tag_length *= space_dim;
242  break;
243  default:
244  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
245  "field with that space is not implemented");
246  }
247 
248  double def_VAL[tag_length];
249  bzero(def_VAL, tag_length * sizeof(double));
250  Tag th;
251  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
252  MB_TYPE_DOUBLE, th,
253  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
254 
255  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
256  // no need for L2
257  const void *tags_ptr[mapGaussPts.size()];
258  int nb_gauss_pts = data.getN().size1();
259  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
260  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
261  }
262 
263  auto set_float_precision = [](const double x) {
264  if (std::abs(x) < std::numeric_limits<float>::epsilon())
265  return 0.;
266  else
267  return x;
268  };
269 
270  auto clear_vals = [&]() {
271  MoFEMFunctionBeginHot;
272  for (int gg = 0; gg < nb_gauss_pts; gg++) {
273  (commonData.gradMap[rowFieldName])[gg].resize(rank, space_dim);
274  (commonData.gradMap[rowFieldName])[gg].clear();
275  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
276  }
277  MoFEMFunctionReturnHot(0);
278  };
279 
280  auto set_vals = [&]() {
281  MoFEMFunctionBeginHot;
282  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
283  tags_ptr);
284  for (int gg = 0; gg < nb_gauss_pts; gg++) {
285  for (int rr = 0; rr < rank; rr++) {
286  for (int dd = 0; dd < space_dim; dd++) {
287  for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
288  const double val =
289  data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
290  (commonData.gradMap[rowFieldName])[gg](rr, dd) += val;
291  }
292  ((double *)tags_ptr[gg])[rank * rr + dd] =
293  (commonData.gradMap[rowFieldName])[gg](rr, dd);
294  }
295  }
296  }
297  MoFEMFunctionReturnHot(0);
298  };
299 
300  commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
301 
302  switch (space) {
303  case H1:
304  if (type == MBVERTEX)
305  clear_vals();
306  CHKERR set_vals();
307  break;
308  case L2:
309  clear_vals();
310  CHKERR set_vals();
311  break;
312  default:
313  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
314  "field with that space is not implemented");
315  }
316 
317  MoFEMFunctionReturn(0);
318 }
319 
320 MoFEMErrorCode PostProcFatPrismOnRefinedMesh::generateReferenceElementMesh() {
321  MoFEMFunctionBeginHot;
322  MoFEMFunctionReturnHot(0);
323 }
324 
325 MoFEMErrorCode PostProcFatPrismOnRefinedMesh::setGaussPtsTrianglesOnly(
326  int order_triangles_only) {
327  MoFEMFunctionBegin;
328  // if(gaussPtsTrianglesOnly.size1()==0 || gaussPtsTrianglesOnly.size2()==0) {
329  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"post-process mesh not
330  // generated");
331  // }
332 
333  int max_ref_level_thick = 0;
334  PetscBool flg = PETSC_TRUE;
335  PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
336  "-my_max_post_proc_ref_level_prism_thick",
337  &max_ref_level_thick, &flg);
338  int const level_thick = max_ref_level_thick + 1;
339  const EntityHandle *conn;
340  int num_nodes;
341  EntityHandle prism;
342  const int nb_on_triangle = 6;
343  const int nb_through_thickness = 3;
344 
345  std::vector<PointsMap3D_multiIndex> pointsMapVector;
346  PointsMap3D_multiIndex pointsMap;
347 
348  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
349  elementsMap.end()) {
350  pointsMapVector = pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()];
351  } else {
352  pointsMapVector.resize(0);
353  gaussPtsTrianglesOnly.resize(3, nb_on_triangle, false);
354  gaussPtsTrianglesOnly.clear();
355  gaussPtsThroughThickness.resize(2, nb_through_thickness * level_thick,
356  false);
357  gaussPtsThroughThickness.clear();
358 
359  MatrixDouble gauss_pts_triangles_only, gauss_pts_through_thickness;
360 
361  double inc_thick = 1.0 / (double)level_thick;
362  for (int ll = 0; ll < level_thick; ll++) {
363 
364  gauss_pts_triangles_only.resize(3, 3, false);
365  gauss_pts_triangles_only.clear();
366  gauss_pts_triangles_only(0, 0) = 0;
367  gauss_pts_triangles_only(1, 0) = 0;
368  gauss_pts_triangles_only(0, 1) = 1;
369  gauss_pts_triangles_only(1, 1) = 0;
370  gauss_pts_triangles_only(0, 2) = 0;
371  gauss_pts_triangles_only(1, 2) = 1;
372  gauss_pts_through_thickness.resize(2, 2, false);
373  gauss_pts_through_thickness(0, 0) = ll * inc_thick;
374  gauss_pts_through_thickness(0, 1) = (ll + 1) * inc_thick;
375 
376  ublas::vector<EntityHandle> prism_conn(6);
377  VectorDouble coords(3);
378  for (int ggf = 0; ggf != 3; ggf++) {
379  double ksi = gauss_pts_triangles_only(0, ggf);
380  double eta = gauss_pts_triangles_only(1, ggf);
381  coords[0] = ksi;
382  coords[1] = eta;
383  for (int ggt = 0; ggt != 2; ggt++) {
384  double zeta = gauss_pts_through_thickness(0, ggt);
385  coords[2] = zeta;
386  int side = ggt * 3 + ggf;
387  CHKERR postProcMesh.create_vertex(&coords[0], prism_conn[side]);
388  }
389  }
390  CHKERR postProcMesh.create_element(MBPRISM, &prism_conn[0], 6, prism);
391 
392  elementsMap[numeredEntFiniteElementPtr->getEnt()].push_back(prism);
393  // Range faces;
394  // CHKERR postProcMesh.get_adjacencies(&prism,1,2,true,faces);
395  Range edges;
396  CHKERR postProcMesh.get_adjacencies(&prism, 1, 1, true, edges);
397  EntityHandle meshset;
398  CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
399  CHKERR postProcMesh.add_entities(meshset, &prism, 1);
400  // CHKERR postProcMesh.add_entities(meshset,faces);
401  CHKERR postProcMesh.add_entities(meshset, edges);
402  if (tenNodesPostProcTets) {
403  CHKERR postProcMesh.convert_entities(meshset, true, false, false);
404  }
405  CHKERR postProcMesh.delete_entities(&meshset, 1);
406  CHKERR postProcMesh.delete_entities(edges);
407  // CHKERR postProcMesh.delete_entities(faces);
408 
409  CHKERR mField.get_moab().get_connectivity(
410  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
411  MatrixDouble coords_prism_global;
412  coords_prism_global.resize(num_nodes, 3, false);
413  CHKERR mField.get_moab().get_coords(conn, num_nodes,
414  &coords_prism_global(0, 0));
415 
416  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
417  MatrixDouble coords_prism_local;
418  coords_prism_local.resize(num_nodes, 3, false);
419  CHKERR postProcMesh.get_coords(conn, num_nodes,
420  &coords_prism_local(0, 0));
421 
422  const double eps = 1e-6;
423  int ggf = 0, ggt = 0;
424 
425  pointsMap.clear();
426 
427  for (int nn = 0; nn != num_nodes; nn++) {
428  double ksi = coords_prism_local(nn, 0);
429  double eta = coords_prism_local(nn, 1);
430  double zeta = coords_prism_local(nn, 2);
431  pointsMap.insert(PointsMap3D(ksi * 100, eta * 100, zeta * 100, nn));
432  if (fabs(zeta) < eps) {
433  gaussPtsTrianglesOnly(0, ggf) = ksi;
434  gaussPtsTrianglesOnly(1, ggf) = eta;
435  ggf++;
436  }
437  if (fabs(ksi) < eps && fabs(eta) < eps) {
438  gaussPtsThroughThickness(0, ggt + nb_through_thickness * ll) = zeta;
439  ggt++;
440  }
441  }
442 
443  pointsMapVector.push_back(pointsMap);
444 
445  for (int nn = 0; nn != num_nodes; nn++) {
446  const double ksi = coords_prism_local(nn, 0);
447  const double eta = coords_prism_local(nn, 1);
448  const double zeta = coords_prism_local(nn, 2);
449  const double n0 = N_MBTRI0(ksi, eta);
450  const double n1 = N_MBTRI1(ksi, eta);
451  const double n2 = N_MBTRI2(ksi, eta);
452  const double e0 = N_MBEDGE0(zeta);
453  const double e1 = N_MBEDGE1(zeta);
454  double coords_global[3];
455  for (int dd = 0; dd != 3; dd++) {
456  coords_global[dd] = n0 * e0 * coords_prism_global(0, dd) +
457  n1 * e0 * coords_prism_global(1, dd) +
458  n2 * e0 * coords_prism_global(2, dd) +
459  n0 * e1 * coords_prism_global(3, dd) +
460  n1 * e1 * coords_prism_global(4, dd) +
461  n2 * e1 * coords_prism_global(5, dd);
462  }
463  CHKERR postProcMesh.set_coords(&conn[nn], 1, coords_global);
464  }
465  }
466 
467  pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()] = pointsMapVector;
468  }
469 
470  mapGaussPts.clear();
471  mapGaussPts.resize(nb_through_thickness * nb_on_triangle * level_thick);
472  fill(mapGaussPts.begin(), mapGaussPts.end(), 0);
473  int gg = 0;
474  int ll;
475 
476  for (unsigned int ggf = 0; ggf != gaussPtsTrianglesOnly.size2(); ggf++) {
477  const double ksi = gaussPtsTrianglesOnly(0, ggf);
478  const double eta = gaussPtsTrianglesOnly(1, ggf);
479  for (unsigned int ggt = 0; ggt != gaussPtsThroughThickness.size2();
480  ggt++, gg++) {
481  ll = ggt / 3;
482  prism = elementsMap[numeredEntFiniteElementPtr->getEnt()][ll];
483  pointsMap = pointsMapVector[ll];
484  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
485 
486  const double zeta = gaussPtsThroughThickness(0, ggt);
487  PointsMap3D_multiIndex::iterator it;
488  it = pointsMap.find(boost::make_tuple(ksi * 100, eta * 100, zeta * 100));
489  if (it != pointsMap.end()) {
490  mapGaussPts[gg] = conn[it->nN];
491  }
492  }
493  }
494 
495  int g = 0;
496 
497  MoFEMFunctionReturn(0);
498 }
499 
500 MoFEMErrorCode PostProcFatPrismOnRefinedMesh::setGaussPtsThroughThickness(
501  int order_thickness) {
502  MoFEMFunctionBeginHot;
503  if (gaussPtsThroughThickness.size1() == 0 ||
504  gaussPtsThroughThickness.size2() == 0) {
505  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
506  "post-process mesh not generated");
507  }
508  MoFEMFunctionReturnHot(0);
509 }
510 
511 MoFEMErrorCode PostProcFatPrismOnRefinedMesh::preProcess() {
512  MoFEMFunctionBegin;
513  MoFEMFunctionReturn(0);
514 }
515 
516 MoFEMErrorCode PostProcFatPrismOnRefinedMesh::postProcess() {
517  MoFEMFunctionBegin;
518  ParallelComm *pcomm_post_proc_mesh =
519  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
520  if (pcomm_post_proc_mesh == NULL) {
521  // wrapRefMeshComm = boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
522  pcomm_post_proc_mesh =
523  new ParallelComm(&postProcMesh, PETSC_COMM_WORLD/*wrapRefMeshComm->get_comm()*/);
524  }
525 
526  Range prims;
527  CHKERR postProcMesh.get_entities_by_type(0, MBPRISM, prims, false);
528  // std::cerr << "total prims size " << prims.size() << std::endl;
529  int rank = mField.get_comm_rank();
530  Range::iterator pit = prims.begin();
531  for (; pit != prims.end(); pit++) {
532  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
533  &rank);
534  }
535  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
536  MoFEMFunctionReturn(0);
537 }
538 
539 MoFEMErrorCode PostProcFaceOnRefinedMesh::generateReferenceElementMesh() {
540  MoFEMFunctionBegin;
541 
542  auto generate_tri = [&](auto &gauss_pts) {
543  MoFEMFunctionBegin;
544  gauss_pts.resize(3, 3, false);
545  gauss_pts.clear();
546  gauss_pts(0, 0) = 0;
547  gauss_pts(1, 0) = 0;
548  gauss_pts(0, 1) = 1;
549  gauss_pts(1, 1) = 0;
550  gauss_pts(0, 2) = 0;
551  gauss_pts(1, 2) = 1;
552 
553  moab::Core core_ref;
554  moab::Interface &moab_ref = core_ref;
555  const EntityHandle *conn;
556  int num_nodes;
557  EntityHandle tri_conn[3];
558  MatrixDouble coords(6, 3);
559  for (int gg = 0; gg != 3; gg++) {
560  coords(gg, 0) = gauss_pts(0, gg);
561  coords(gg, 1) = gauss_pts(1, gg);
562  coords(gg, 2) = 0;
563  CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
564  }
565 
566  EntityHandle tri;
567  CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
568 
569  if (sixNodePostProcTris) {
570  Range edges;
571  CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges);
572  EntityHandle meshset;
573  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
574  CHKERR moab_ref.add_entities(meshset, &tri, 1);
575  CHKERR moab_ref.add_entities(meshset, edges);
576  CHKERR moab_ref.convert_entities(meshset, true, false, false);
577  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
578  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
579  gauss_pts.resize(3, num_nodes, false);
580  gauss_pts.clear();
581  for (int nn = 0; nn < num_nodes; nn++) {
582  gauss_pts(0, nn) = coords(nn, 0);
583  gauss_pts(1, nn) = coords(nn, 1);
584  }
585  } else {
586  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
587  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
588  gauss_pts.resize(3, num_nodes, false);
589  gauss_pts.clear();
590  for (int nn = 0; nn < 3; nn++) {
591  gauss_pts(0, nn) = coords(nn, 0);
592  gauss_pts(1, nn) = coords(nn, 1);
593  }
594  }
595  MoFEMFunctionReturn(0);
596  };
597 
598  auto generate_quad = [&](auto &gauss_pts) {
599  MoFEMFunctionBegin;
600 
601  gauss_pts.resize(3, 4, false);
602  gauss_pts.clear();
603 
604  gauss_pts(0, 0) = 0;
605  gauss_pts(1, 0) = 0;
606  gauss_pts(0, 1) = 1;
607  gauss_pts(1, 1) = 0;
608  gauss_pts(0, 2) = 1;
609  gauss_pts(1, 2) = 1;
610  gauss_pts(0, 3) = 0;
611  gauss_pts(1, 3) = 1;
612 
613  moab::Core core_ref;
614  moab::Interface &moab_ref = core_ref;
615 
616  const EntityHandle *conn;
617  int num_nodes;
618  std::array<EntityHandle, 4> quad_conn;
619  MatrixDouble coords(8, 3);
620  for (int gg = 0; gg != 4; gg++) {
621  coords(gg, 0) = gauss_pts(0, gg);
622  coords(gg, 1) = gauss_pts(1, gg);
623  coords(gg, 2) = 0;
624  CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
625  }
626 
627  EntityHandle quad;
628  CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
629 
630  if (sixNodePostProcTris) {
631  Range edges;
632  CHKERR moab_ref.get_adjacencies(&quad, 1, 1, true, edges);
633  EntityHandle meshset;
634  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
635  CHKERR moab_ref.add_entities(meshset, &quad, 1);
636  CHKERR moab_ref.add_entities(meshset, edges);
637  CHKERR moab_ref.convert_entities(meshset, true, false, false);
638  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
639  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
640  gauss_pts.resize(3, num_nodes, false);
641  gauss_pts.clear();
642  for (int nn = 0; nn != num_nodes; nn++) {
643  gauss_pts(0, nn) = coords(nn, 0);
644  gauss_pts(1, nn) = coords(nn, 1);
645  }
646  } else {
647  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
648  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
649  gauss_pts.resize(3, num_nodes, false);
650  gauss_pts.clear();
651  for (int nn = 0; nn != 4; nn++) {
652  gauss_pts(0, nn) = coords(nn, 0);
653  gauss_pts(1, nn) = coords(nn, 1);
654  }
655  }
656 
657  MoFEMFunctionReturn(0);
658  };
659 
660  CHKERR generate_tri(gaussPtsTri);
661  CHKERR generate_quad(gaussPtsQuad);
662 
663  MoFEMFunctionReturn(0);
664 }
665 
666 MoFEMErrorCode PostProcFaceOnRefinedMesh::setGaussPts(int order) {
667  MoFEMFunctionBegin;
668  if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
669  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
670  "post-process mesh not generated");
671 
672  auto create_tri = [&](auto &gauss_pts) {
673  std::array<EntityHandle, 3> tri_conn;
674  MatrixDouble3by3 coords_tri(3, 3);
675  CHKERR mField.get_moab().get_connectivity(
676  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
677  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_tri(0, 0));
678  const int num_nodes_on_ele = gauss_pts.size2();
679 
680  for (int gg = 0; gg != num_nodes_on_ele; gg++) {
681 
682  const double ksi = gauss_pts(0, gg);
683  const double eta = gauss_pts(1, gg);
684  const double n0 = N_MBTRI0(ksi, eta);
685  const double n1 = N_MBTRI1(ksi, eta);
686  const double n2 = N_MBTRI2(ksi, eta);
687  const double x =
688  n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
689  const double y =
690  n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
691  const double z =
692  n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
693 
694  verticesOnTriArrays[0][counterTris * num_nodes_on_ele + gg] = x;
695  verticesOnTriArrays[1][counterTris * num_nodes_on_ele + gg] = y;
696  verticesOnTriArrays[2][counterTris * num_nodes_on_ele + gg] = z;
697  }
698 
699  mapGaussPts.resize(num_nodes_on_ele);
700  for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
701  triConn[num_nodes_on_ele * counterTris + nn] =
702  num_nodes_on_ele * counterTris + nn + startingVertTriHandle;
703  mapGaussPts[nn] = triConn[num_nodes_on_ele * counterTris + nn];
704  }
705 
706  const auto tri = startingEleTriHandle + counterTris;
707 
708  return tri;
709  };
710 
711  auto create_quad = [&](auto &gauss_pts) {
712  std::array<EntityHandle, 4> quad_conn;
713  MatrixDouble coords_quad(4, 3);
714  CHKERR mField.get_moab().get_connectivity(
715  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
716  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_quad(0, 0));
717  const int num_nodes_on_ele = gauss_pts.size2();
718  for (int gg = 0; gg != num_nodes_on_ele; ++gg) {
719  double ksi = gauss_pts(0, gg);
720  double eta = gauss_pts(1, gg);
721  double n0 = N_MBQUAD0(ksi, eta);
722  double n1 = N_MBQUAD1(ksi, eta);
723  double n2 = N_MBQUAD2(ksi, eta);
724  double n3 = N_MBQUAD3(ksi, eta);
725  double x = n0 * coords_quad(0, 0) + n1 * coords_quad(1, 0) +
726  n2 * coords_quad(2, 0) + n3 * coords_quad(3, 0);
727  double y = n0 * coords_quad(0, 1) + n1 * coords_quad(1, 1) +
728  n2 * coords_quad(2, 1) + n3 * coords_quad(3, 1);
729  double z = n0 * coords_quad(0, 2) + n1 * coords_quad(1, 2) +
730  n2 * coords_quad(2, 2) + n3 * coords_quad(3, 2);
731  verticesOnQuadArrays[0][counterQuads * num_nodes_on_ele + gg] = x;
732  verticesOnQuadArrays[1][counterQuads * num_nodes_on_ele + gg] = y;
733  verticesOnQuadArrays[2][counterQuads * num_nodes_on_ele + gg] = z;
734  }
735 
736  mapGaussPts.resize(num_nodes_on_ele);
737 
738  for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
739  quadConn[num_nodes_on_ele * counterQuads + nn] =
740  num_nodes_on_ele * counterQuads + nn + startingVertQuadHandle;
741  mapGaussPts[nn] = quadConn[num_nodes_on_ele * counterQuads + nn];
742  }
743 
744  const auto quad = startingEleQuadHandle + counterQuads;
745  return quad;
746  };
747 
748  EntityHandle tri;
749 
750  if (elementsMap.size() == getLoopSize()) {
751  // Note "at" that will trigger error if element is not there.
752  tri = elementsMap.at(numeredEntFiniteElementPtr->getEnt());
753  switch (numeredEntFiniteElementPtr->getEntType()) {
754  case MBTRI:
755  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
756  noalias(gaussPts) = gaussPtsTri;
757  break;
758  case MBQUAD:
759  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
760  noalias(gaussPts) = gaussPtsQuad;
761  break;
762  default:
763  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
764  "Not implemented element type");
765  }
766 
767  // Set values which map nodes with integration points on the prism
768  const EntityHandle *conn;
769  int num_nodes;
770  CHKERR postProcMesh.get_connectivity(tri, conn, num_nodes, false);
771  mapGaussPts.resize(num_nodes);
772  for (int nn = 0; nn != num_nodes; nn++)
773  mapGaussPts[nn] = conn[nn];
774 
775  } else {
776  switch (numeredEntFiniteElementPtr->getEntType()) {
777  case MBTRI:
778  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
779  noalias(gaussPts) = gaussPtsTri;
780  tri = create_tri(gaussPtsTri);
781  ++counterTris;
782  break;
783  case MBQUAD:
784  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
785  noalias(gaussPts) = gaussPtsQuad;
786  tri = create_quad(gaussPtsQuad);
787  ++counterQuads;
788  break;
789  default:
790  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
791  "Not implemented element type");
792  }
793 
794  elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
795  }
796 
797  MoFEMFunctionReturn(0);
798 }
799 
800 MoFEMErrorCode PostProcFaceOnRefinedMesh::preProcess() {
801  MoFEMFunctionBegin;
802 
803  ReadUtilIface *iface;
804  CHKERR postProcMesh.query_interface(iface);
805 
806  const int number_of_ents_in_the_loop = getLoopSize();
807  if (elementsMap.size() != number_of_ents_in_the_loop) {
808 
809  elementsMap.clear();
810  postProcMesh.delete_mesh();
811 
812  auto get_number_of_computational_elements = [&]() {
813  auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
814 
815  auto miit =
816  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
817  boost::make_tuple(this->getFEName(), this->getLoFERank()));
818  auto hi_miit =
819  fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
820  boost::make_tuple(this->getFEName(), this->getHiFERank()));
821 
822  const int number_of_ents_in_the_loop = this->getLoopSize();
823  if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
824  THROW_MESSAGE(
825  "Wrong size of indicices. Inconsistent size number of iterated "
826  "elements iterated by problem and from range.");
827  }
828 
829  std::array<int, MBMAXTYPE> nb_elemms_by_type;
830  std::fill(nb_elemms_by_type.begin(), nb_elemms_by_type.end(), 0);
831 
832  for (; miit != hi_miit; ++miit) {
833 
834  bool add = true;
835  if (exeTestHook) {
836  numeredEntFiniteElementPtr = *miit;
837  add = exeTestHook(this);
838  }
839 
840  if (add) {
841  auto type = (*miit)->getEntType();
842  ++nb_elemms_by_type[type];
843  }
844  }
845 
846  return nb_elemms_by_type;
847  };
848 
849  auto nb_computational_elements_by_type =
850  get_number_of_computational_elements();
851 
852  const int numberOfTriangles = nb_computational_elements_by_type[MBTRI];
853  const int numberOfQuads = nb_computational_elements_by_type[MBQUAD];
854 
855  // Here we create vertices using ReadUtilface
856  const int total_number_of_nodes_on_tri =
857  numberOfTriangles * gaussPtsTri.size2();
858  const int total_number_of_nodes_on_quad =
859  numberOfQuads * gaussPtsQuad.size2();
860 
861  if (total_number_of_nodes_on_tri) {
862  CHKERR iface->get_node_coords(3, total_number_of_nodes_on_tri, 0,
863  startingVertTriHandle, verticesOnTriArrays);
864  CHKERR iface->get_element_connect(numberOfTriangles, gaussPtsTri.size2(),
865  MBTRI, 0, startingEleTriHandle,
866  triConn);
867  }
868 
869  if (total_number_of_nodes_on_quad) {
870  CHKERR iface->get_node_coords(3, total_number_of_nodes_on_quad, 0,
871  startingVertQuadHandle,
872  verticesOnQuadArrays);
873  CHKERR iface->get_element_connect(numberOfQuads, gaussPtsQuad.size2(),
874  MBQUAD, 0, startingEleQuadHandle,
875  quadConn);
876  }
877  }
878 
879  counterTris = 0;
880  counterQuads = 0;
881 
882  MoFEMFunctionReturn(0);
883 }
884 
885 MoFEMErrorCode PostProcFaceOnRefinedMesh::postProcess() {
886  MoFEMFunctionBegin;
887 
888  auto update_elements = [&]() {
889  MoFEMFunctionBegin;
890  ReadUtilIface *iface;
891  CHKERR postProcMesh.query_interface(iface);
892 
893  if (counterTris)
894  CHKERR iface->update_adjacencies(startingEleTriHandle, counterTris,
895  gaussPtsTri.size2(), triConn);
896  if (counterQuads)
897  CHKERR iface->update_adjacencies(startingEleQuadHandle, counterQuads,
898  gaussPtsQuad.size2(), quadConn);
899  MoFEMFunctionReturn(0);
900  };
901 
902  auto resolve_shared = [&]() {
903  MoFEMFunctionBegin;
904  ParallelComm *pcomm_post_proc_mesh =
905  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
906  if (pcomm_post_proc_mesh == NULL) {
907  // wrapRefMeshComm =
908  // boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
909  pcomm_post_proc_mesh = new ParallelComm(
910  &postProcMesh, PETSC_COMM_WORLD /* wrapRefMeshComm->get_comm()*/);
911  }
912 
913  Range faces;
914  CHKERR postProcMesh.get_entities_by_dimension(0, 2, faces, false);
915  int rank = mField.get_comm_rank();
916  CHKERR postProcMesh.tag_clear_data(pcomm_post_proc_mesh->part_tag(), faces,
917  &rank);
918  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
919  MoFEMFunctionReturn(0);
920  };
921 
922  CHKERR update_elements();
923  CHKERR resolve_shared();
924 
925  MoFEMFunctionReturn(0);
926 }
927 
928 template <int RANK>
929 MoFEMErrorCode
930 PostProcFaceOnRefinedMesh::OpGetFieldValuesOnSkinImpl<RANK>::doWork(
931  int side, EntityType type, EntitiesFieldData::EntData &data) {
932  MoFEMFunctionBegin;
933 
934  if (type != MBVERTEX)
935  MoFEMFunctionReturnHot(0);
936 
937  CHKERR loopSideVolumes(feVolName, *sideOpFe);
938 
939  // quit if tag is not needed
940  if (!saveOnTag)
941  MoFEMFunctionReturnHot(0);
942 
943  auto &m_field = getPtrFE()->mField;
944  auto field_ptr = m_field.get_field_structure(fieldName);
945  const int rank = field_ptr->getNbOfCoeffs();
946  FieldSpace space = field_ptr->getSpace();
947 
948  // auto dof_ptr = data.getFieldDofs()[0];
949  // int rank = dof_ptr->getNbOfCoeffs();
950  // FieldSpace space = dof_ptr->getSpace();
951 
952  int full_size = rank * RANK;
953  if (space == HDIV)
954  full_size *= 3;
955  // for paraview
956  int tag_length = full_size > 3 && full_size < 9 ? 9 : full_size;
957  double def_VAL[tag_length];
958  bzero(def_VAL, tag_length * sizeof(double));
959  Tag th;
960  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
961  MB_TYPE_DOUBLE, th,
962  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
963 
964  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
965  // no need for L2
966  const void *tags_ptr[mapGaussPts.size()];
967  int nb_gauss_pts = data.getN().size1();
968  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts ||
969  nb_gauss_pts != matPtr->size2()) {
970  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
971  }
972  switch (space) {
973  case H1:
974  case L2:
975  case HDIV:
976 
977  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
978  tags_ptr);
979 
980  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
981  const double *my_ptr2 = static_cast<const double *>(tags_ptr[gg]);
982  double *my_ptr = const_cast<double *>(my_ptr2);
983  for (int rr = 0; rr != RANK; ++rr) {
984  for (int dd = 0; dd != rank; ++dd) {
985  my_ptr[rank * rr + dd] = (*matPtr)(rank * rr + dd, gg);
986  }
987  }
988  }
989 
990  break;
991  default:
992  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
993  "field with that space is not implemented");
994  }
995 
996  MoFEMFunctionReturn(0);
997 }
998 
999 MoFEMErrorCode PostProcFaceOnRefinedMesh::addFieldValuesGradientPostProcOnSkin(
1000  const std::string field_name, const std::string vol_fe_name,
1001  boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag) {
1002  MoFEMFunctionBegin;
1003 
1004  if (!grad_mat_ptr)
1005  grad_mat_ptr = boost::make_shared<MatrixDouble>();
1006 
1007  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
1008  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1009 
1010  if (mField.check_field("MESH_NODE_POSITIONS"))
1011  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *my_side_fe, true, false, false,
1012  false);
1013 
1014  // check number of coefficients
1015  auto field_ptr = mField.get_field_structure(field_name);
1016  const int nb_coefficients = field_ptr->getNbOfCoeffs();
1017 
1018  switch (nb_coefficients) {
1019  case 1:
1020  my_side_fe->getOpPtrVector().push_back(
1021  new OpCalculateScalarFieldGradient<3>(field_name, grad_mat_ptr));
1022  break;
1023  case 3:
1024  my_side_fe->getOpPtrVector().push_back(
1025  new OpCalculateVectorFieldGradient<3, 3>(field_name, grad_mat_ptr));
1026  break;
1027  case 6:
1028  my_side_fe->getOpPtrVector().push_back(
1029  new OpCalculateTensor2SymmetricFieldGradient<3, 3>(field_name,
1030  grad_mat_ptr));
1031  break;
1032  default:
1033  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1034  "field with that number of coefficients is not implemented");
1035  }
1036 
1037  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1038  new OpGetFieldValuesOnSkinImpl<3>(
1039  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
1040  my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
1041 
1042  MoFEMFunctionReturn(0);
1043 }
1044 
1045 MoFEMErrorCode PostProcFaceOnRefinedMesh::addFieldValuesPostProcOnSkin(
1046  const std::string field_name, const std::string vol_fe_name,
1047  boost::shared_ptr<MatrixDouble> mat_ptr, bool save_on_tag) {
1048  MoFEMFunctionBegin;
1049 
1050  if (!mat_ptr)
1051  mat_ptr = boost::make_shared<MatrixDouble>();
1052 
1053  auto my_side_fe =
1054  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1055 
1056  // check number of coefficients
1057  auto field_ptr = mField.get_field_structure(field_name);
1058  const int nb_coefficients = field_ptr->getNbOfCoeffs();
1059  FieldSpace space = field_ptr->getSpace();
1060 
1061  switch (space) {
1062  case L2:
1063  case H1:
1064  switch (nb_coefficients) {
1065  case 1:
1066  my_side_fe->getOpPtrVector().push_back(
1067  new OpCalculateVectorFieldValues<1>(field_name, mat_ptr));
1068  break;
1069  case 3:
1070  my_side_fe->getOpPtrVector().push_back(
1071  new OpCalculateVectorFieldValues<3>(field_name, mat_ptr));
1072  break;
1073  case 6:
1074  my_side_fe->getOpPtrVector().push_back(
1075  new OpCalculateTensor2SymmetricFieldValues<3>(field_name, mat_ptr));
1076  break;
1077  case 9:
1078  my_side_fe->getOpPtrVector().push_back(
1079  new OpCalculateTensor2FieldValues<3, 3>(field_name, mat_ptr));
1080  break;
1081  default:
1082  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1083  "field with that number of coefficients is not implemented");
1084  }
1085  break;
1086  case HDIV:
1087  switch (nb_coefficients) {
1088  case 1:
1089  my_side_fe->getOpPtrVector().push_back(
1090  new OpCalculateHVecVectorField<3>(field_name, mat_ptr));
1091  break;
1092  case 3:
1093  my_side_fe->getOpPtrVector().push_back(
1094  new OpCalculateHVecTensorField<3, 3>(field_name, mat_ptr));
1095  break;
1096  default:
1097  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1098  "field with that number of coefficients is not implemented");
1099  }
1100  break;
1101  default:
1102  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1103  "field with that space is not implemented.");
1104  }
1105 
1106  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1107  new OpGetFieldValuesOnSkinImpl<1>(postProcMesh, mapGaussPts, field_name,
1108  field_name, my_side_fe, vol_fe_name,
1109  mat_ptr, save_on_tag));
1110 
1111  MoFEMFunctionReturn(0);
1112 }
1113 
1114 MoFEMErrorCode PostProcEdgeOnRefinedMesh::generateReferenceElementMesh() {
1115  MoFEMFunctionBegin;
1116 
1117  gaussPts.resize(3, 2, false);
1118  gaussPts.clear();
1119  gaussPts(0, 0) = 0;
1120  gaussPts(1, 0) = 0;
1121  gaussPts(0, 1) = 1;
1122  gaussPts(1, 1) = 0;
1123  mapGaussPts.resize(gaussPts.size2());
1124 
1125  moab::Core core_ref;
1126  moab::Interface &moab_ref = core_ref;
1127  const EntityHandle *conn;
1128  int num_nodes;
1129  EntityHandle edge_conn[2];
1130  MatrixDouble coords(2, 3);
1131  for (int gg = 0; gg != 2; gg++) {
1132  coords(gg, 0) = gaussPts(0, gg);
1133  coords(gg, 1) = gaussPts(1, gg);
1134  coords(gg, 2) = 0;
1135  CHKERR moab_ref.create_vertex(&coords(gg, 0), edge_conn[gg]);
1136  }
1137 
1138  EntityHandle edge;
1139  CHKERR moab_ref.create_element(MBEDGE, edge_conn, 2, edge);
1140  MoFEMFunctionReturn(0);
1141 }
1142 
1143 MoFEMErrorCode PostProcEdgeOnRefinedMesh::preProcess() {
1144  MoFEMFunctionBeginHot;
1145  ParallelComm *pcomm_post_proc_mesh =
1146  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
1147  if (pcomm_post_proc_mesh != NULL) {
1148  delete pcomm_post_proc_mesh;
1149  }
1150  MoFEMFunctionReturnHot(0);
1151 }
1152 
1153 MoFEMErrorCode PostProcEdgeOnRefinedMesh::postProcess() {
1154  MoFEMFunctionBegin;
1155  ParallelComm *pcomm_post_proc_mesh =
1156  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
1157  if (pcomm_post_proc_mesh == NULL) {
1158  // wrapRefMeshComm = boost::make_shared<WrapMPIComm>(mField.get_comm(), false);
1159  pcomm_post_proc_mesh = new ParallelComm(
1160  &postProcMesh, PETSC_COMM_WORLD /*wrapRefMeshComm->get_comm()*/);
1161  }
1162 
1163  Range edges;
1164  CHKERR postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
1165  int rank = mField.get_comm_rank();
1166  auto set_edges_rank = [&](const auto rank, const auto &edges) {
1167  std::vector<EntityHandle> ranks(edges.size(), rank);
1168  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), edges,
1169  &*ranks.begin());
1170  };
1171  set_edges_rank(rank, edges);
1172 
1173  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
1174  MoFEMFunctionReturn(0);
1175 }
1176 
1177 MoFEMErrorCode PostProcEdgeOnRefinedMesh::setGaussPts(int order) {
1178  MoFEMFunctionBegin;
1179  if (gaussPts.size1() == 0) {
1180  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1181  "post-process mesh not generated");
1182  }
1183 
1184  const EntityHandle *conn;
1185  int num_nodes;
1186  EntityHandle edge;
1187 
1188  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
1189  elementsMap.end()) {
1190  edge = elementsMap[numeredEntFiniteElementPtr->getEnt()];
1191  } else {
1192  ublas::vector<EntityHandle> edge_conn(2);
1193  MatrixDouble coords_edge(2, 3);
1194  VectorDouble coords(3);
1195  CHKERR mField.get_moab().get_connectivity(
1196  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
1197  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_edge(0, 0));
1198  for (int gg = 0; gg != 2; gg++) {
1199  double ksi = gaussPts(0, gg);
1200  // double eta = gaussPts(1, gg);
1201  double n0 = N_MBEDGE0(ksi);
1202  double n1 = N_MBEDGE1(ksi);
1203 
1204  double x = n0 * coords_edge(0, 0) + n1 * coords_edge(1, 0);
1205  double y = n0 * coords_edge(0, 1) + n1 * coords_edge(1, 1);
1206  double z = n0 * coords_edge(0, 2) + n1 * coords_edge(1, 2);
1207 
1208  coords[0] = x;
1209  coords[1] = y;
1210  coords[2] = z;
1211 
1212  CHKERR postProcMesh.create_vertex(&coords[0], edge_conn[gg]);
1213  }
1214  CHKERR postProcMesh.create_element(MBEDGE, &edge_conn[0], 2, edge);
1215  elementsMap[numeredEntFiniteElementPtr->getEnt()] = edge;
1216  }
1217 
1218  {
1219  CHKERR postProcMesh.get_connectivity(edge, conn, num_nodes, false);
1220  mapGaussPts.resize(num_nodes);
1221  for (int nn = 0; nn < num_nodes; nn++) {
1222  mapGaussPts[nn] = conn[nn];
1223  }
1224  }
1225 
1226  MoFEMFunctionReturn(0);
1227 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
PostProcFaceOnRefinedMesh::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition: PostProcOnRefMesh.cpp:666
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
g
constexpr double g
Definition: shallow_wave.cpp:63
PostProcCommonOnRefMesh::OpGetFieldGradientValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: PostProcOnRefMesh.cpp:199
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
PostProcFatPrismOnRefinedMesh::setGaussPtsTrianglesOnly
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
Definition: PostProcOnRefMesh.cpp:325
H1
@ H1
continuous field
Definition: definitions.h:85
PostProcFaceOnRefinedMesh::OpGetFieldValuesOnSkinImpl
Definition: PostProcOnRefMesh.hpp:1063
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:887
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
PostProcFatPrismOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: PostProcOnRefMesh.cpp:320
PostProcFatPrismOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.cpp:516
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
PostProcFaceOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.cpp:800
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
PostProcFaceOnRefinedMesh::OpGetFieldValuesOnSkinImpl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: PostProcOnRefMesh.cpp:930
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
MoFEM.hpp
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1316
PostProcEdgeOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.cpp:1153
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::OpCalculateTensor2SymmetricFieldGradient
Get field gradients at integration pts for symmetric tensorial field rank 2.
Definition: UserDataOperators.hpp:1734
PostProcEdgeOnRefinedMesh::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition: PostProcOnRefMesh.cpp:1177
PostProcFatPrismOnRefinedMesh::PointsMap3D_multiIndex
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
Definition: PostProcOnRefMesh.hpp:1003
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2116
PostProcFatPrismOnRefinedMesh::PointsMap3D
Definition: PostProcOnRefMesh.hpp:988
double
convert.type
type
Definition: convert.py:64
eta
double eta
Definition: free_surface.cpp:170
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2576
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
PostProcEdgeOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.cpp:1143
PostProcEdgeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: PostProcOnRefMesh.cpp:1114
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
N_MBTRI0
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:46
PostProcFaceOnRefinedMesh::addFieldValuesGradientPostProcOnSkin
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
Definition: PostProcOnRefMesh.cpp:999
Range
FTensor::dd
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
PostProcCommonOnRefMesh::OpGetFieldValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: PostProcOnRefMesh.cpp:16
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
N_MBTRI1
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:47
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
PostProcFaceOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.cpp:885
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:978
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
PostProcFatPrismOnRefinedMesh::setGaussPtsThroughThickness
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
Definition: PostProcOnRefMesh.cpp:500
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
PostProcFaceOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: PostProcOnRefMesh.cpp:539
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
PostProcFaceOnRefinedMesh::addFieldValuesPostProcOnSkin
MoFEMErrorCode addFieldValuesPostProcOnSkin(const std::string field_name, const std::string vol_fe_name="dFE", boost::shared_ptr< MatrixDouble > mat_ptr=nullptr, bool save_on_tag=true)
Definition: PostProcOnRefMesh.cpp:1045
N_MBTRI2
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:48
PostProcFatPrismOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.cpp:511
Generated by Doxygen 1.8.17 and hosted at University of Glasgow