v0.9.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  * 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 #include <MoFEM.hpp>
25 using namespace MoFEM;
26 using namespace boost::numeric;
27 #include <PostProcOnRefMesh.hpp>
28 
30  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
32 
33  if (data.getFieldData().size() == 0)
35 
36  if (V) {
37  vAlues.resize(data.getFieldData().size());
38  double *a;
39  CHKERR VecGetArray(V, &a);
40  VectorDofs::iterator it, hi_it;
41  it = data.getFieldDofs().begin();
42  hi_it = data.getFieldDofs().end();
43  for (int ii = 0; it != hi_it; it++, ii++) {
44  int local_idx = getFEMethod()
45  ->rowPtr->find((*it)->getGlobalUniqueId())
46  ->get()
47  ->getPetscLocalDofIdx();
48  vAlues[ii] = a[local_idx];
49  }
50  CHKERR VecRestoreArray(V, &a);
51  vAluesPtr = &vAlues;
52  } else {
53  vAluesPtr = &data.getFieldData();
54  }
55 
56  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
57  int rank = dof_ptr->getNbOfCoeffs();
58 
59  int tag_length = rank;
60  FieldSpace space = dof_ptr->getSpace();
61  switch (space) {
62  case L2:
63  case H1:
64  break;
65  case HCURL:
66  case HDIV:
67  tag_length *= 3;
68  break;
69  default:
70  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
71  "field with that space is not implemented");
72  }
73 
74  if (tag_length > 1 && tag_length < 3)
75  tag_length = 3;
76  else if (tag_length > 3 && tag_length < 9)
77  tag_length = 9;
78 
79  double def_VAL[tag_length];
80  bzero(def_VAL, tag_length * sizeof(double));
81  Tag th;
82  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
83  MB_TYPE_DOUBLE, th,
84  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
85 
86  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
87  // no need for L2
88  const void *tags_ptr[mapGaussPts.size()];
89  int nb_gauss_pts = data.getN().size1();
90  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
91  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
92  "data inconsistency %d!=%d", mapGaussPts.size(), nb_gauss_pts);
93  }
94 
95  switch (space) {
96  case H1:
97  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
98  if (type == MBVERTEX) {
99  for (int gg = 0; gg < nb_gauss_pts; gg++) {
100  rval = postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
101  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
102  (commonData.fieldMap[rowFieldName])[gg].clear();
103  }
104  }
105  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
106  tags_ptr);
107  for (int gg = 0; gg < nb_gauss_pts; gg++) {
108  for (int rr = 0; rr < rank; rr++) {
109  const double val =
110  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
111  &((*vAluesPtr)[rr]), rank);
112  (commonData.fieldMap[rowFieldName])[gg][rr] =
113  ((double *)tags_ptr[gg])[rr] += val;
114  }
115  }
116  break;
117  case L2:
118  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
119  for (int gg = 0; gg < nb_gauss_pts; gg++) {
120  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
121  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
122  (commonData.fieldMap[rowFieldName])[gg].clear();
123  }
124  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
125  tags_ptr);
126  for (int gg = 0; gg < nb_gauss_pts; gg++) {
127  bzero((double *)tags_ptr[gg], sizeof(double) * tag_length);
128  for (int rr = 0; rr < rank; rr++) {
129  const double val =
130  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
131  &((*vAluesPtr)[rr]), rank);
132  (commonData.fieldMap[rowFieldName])[gg][rr] =
133  ((double *)tags_ptr[gg])[rr] = val;
134  }
135  }
136  break;
137  case HCURL:
138  // FIXME: fieldMap not set
139  if (type == MBEDGE && side == 0) {
140  for (int gg = 0; gg < nb_gauss_pts; gg++) {
141  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
142  }
143  }
144  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
145  tags_ptr);
146  {
148  auto t_n_hcurl = data.getFTensor1N<3>();
149  for (int gg = 0; gg != nb_gauss_pts; gg++) {
150  double *ptr = &((double *)tags_ptr[gg])[0];
151  int ll = 0;
152  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
153  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
154  for (int rr = 0; rr != rank; rr++) {
155  const double dof_val = (*vAluesPtr)[ll * rank + rr];
156  t_tag_val(i) += dof_val * t_n_hcurl(i);
157  ++t_tag_val;
158  }
159  ++t_n_hcurl;
160  }
161  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
162  ++t_n_hcurl;
163  }
164  }
165  }
166  break;
167  case HDIV:
168  // FIXME: fieldMap not set
169  if (type == MBTRI && side == 0) {
170  for (int gg = 0; gg < nb_gauss_pts; gg++) {
171  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
172  }
173  }
174  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
175  tags_ptr);
176  {
178  auto t_n_hdiv = data.getFTensor1N<3>();
179  for (int gg = 0; gg != nb_gauss_pts; gg++) {
180  double *ptr = &((double *)tags_ptr[gg])[0];
181  int ll = 0;
182  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
183  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
184  for (int rr = 0; rr != rank; rr++) {
185  const double dof_val = (*vAluesPtr)[ll * rank + rr];
186  t_tag_val(i) += dof_val * t_n_hdiv(i);
187  ++t_tag_val;
188  }
189  ++t_n_hdiv;
190  }
191  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
192  ++t_n_hdiv;
193  }
194  }
195  }
196  break;
197  default:
198  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
199  "field with that space is not implemented");
200  }
201 
203 }
204 
206  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
208 
209  if (data.getFieldData().size() == 0)
211  if (V) {
212  vAlues.resize(data.getFieldData().size());
213  double *a;
214  CHKERR VecGetArray(V, &a);
215  VectorDofs::iterator it, hi_it;
216  it = data.getFieldDofs().begin();
217  hi_it = data.getFieldDofs().end();
218  for (int ii = 0; it != hi_it; it++, ii++) {
219  int local_idx = getFEMethod()
220  ->rowPtr->find((*it)->getGlobalUniqueId())
221  ->get()
222  ->getPetscLocalDofIdx();
223  vAlues[ii] = a[local_idx];
224  }
225  CHKERR VecRestoreArray(V, &a);
226  vAluesPtr = &vAlues;
227  } else {
228  vAluesPtr = &data.getFieldData();
229  }
230 
231  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
232  int rank = dof_ptr->getNbOfCoeffs();
233 
234  int tag_length = rank * 3;
235  FieldSpace space = dof_ptr->getSpace();
236  switch (space) {
237  case L2:
238  case H1:
239  break;
240  case HCURL:
241  case HDIV:
242  tag_length *= 3;
243  break;
244  default:
245  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
246  "field with that space is not implemented");
247  }
248 
249  double def_VAL[tag_length];
250  bzero(def_VAL, tag_length * sizeof(double));
251  Tag th;
252  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
253  MB_TYPE_DOUBLE, th,
254  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
255 
256  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
257  // no need for L2
258  const void *tags_ptr[mapGaussPts.size()];
259  int nb_gauss_pts = data.getN().size1();
260  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
261  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
262  }
263 
264  switch (space) {
265  case H1:
266  commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
267  if (type == MBVERTEX) {
268  for (int gg = 0; gg < nb_gauss_pts; gg++) {
269  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
270  (commonData.gradMap[rowFieldName])[gg].resize(rank, 3);
271  (commonData.gradMap[rowFieldName])[gg].clear();
272  }
273  }
274  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
275  tags_ptr);
276  for (int gg = 0; gg < nb_gauss_pts; gg++) {
277  for (int rr = 0; rr < rank; rr++) {
278  for (int dd = 0; dd < 3; dd++) {
279  for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
280  const double val =
281  data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
282  (commonData.gradMap[rowFieldName])[gg](rr, dd) =
283  ((double *)tags_ptr[gg])[rank * rr + dd] += val;
284  }
285  }
286  }
287  }
288  break;
289  case L2:
290  commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
291  for (int gg = 0; gg < nb_gauss_pts; gg++) {
292  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
293  (commonData.gradMap[rowFieldName])[gg].resize(rank, 3);
294  (commonData.gradMap[rowFieldName])[gg].clear();
295  }
296  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
297  tags_ptr);
298  for (int gg = 0; gg < nb_gauss_pts; gg++) {
299  for (int rr = 0; rr < rank; rr++) {
300  for (int dd = 0; dd < 3; dd++) {
301  for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
302  const double val =
303  data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
304  (commonData.gradMap[rowFieldName])[gg](rr, dd) =
305  ((double *)tags_ptr[gg])[rank * rr + dd] += val;
306  }
307  }
308  }
309  }
310  break;
311  default:
312  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
313  "field with that space is not implemented");
314  }
315 
317 }
318 
322 }
323 
325  int order_triangles_only) {
327  // if(gaussPtsTrianglesOnly.size1()==0 || gaussPtsTrianglesOnly.size2()==0) {
328  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"post-process mesh not
329  // generated");
330  // }
331 
332  // FIXME: Refinement not implement and inefficient implementation, looks ugly.
333 
334  //
335  const EntityHandle *conn;
336  int num_nodes;
337  EntityHandle prism;
338  const int nb_on_triangle = 6;
339  const int nb_through_thickness = 3;
340 
341  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
342  elementsMap.end()) {
343  prism = elementsMap[numeredEntFiniteElementPtr->getEnt()];
344  } else {
345  MatrixDouble gauss_pts_triangles_only, gauss_pts_through_thickness;
346  {
347  gauss_pts_triangles_only.resize(3, 3, false);
348  gauss_pts_triangles_only.clear();
349  gauss_pts_triangles_only(0, 0) = 0;
350  gauss_pts_triangles_only(1, 0) = 0;
351  gauss_pts_triangles_only(0, 1) = 1;
352  gauss_pts_triangles_only(1, 1) = 0;
353  gauss_pts_triangles_only(0, 2) = 0;
354  gauss_pts_triangles_only(1, 2) = 1;
355  gauss_pts_through_thickness.resize(2, 2, false);
356  gauss_pts_through_thickness(0, 0) = 0;
357  gauss_pts_through_thickness(0, 1) = 1;
358  }
359  ublas::vector<EntityHandle> prism_conn(6);
360  VectorDouble coords(3);
361  for (int ggf = 0; ggf != 3; ggf++) {
362  double ksi = gauss_pts_triangles_only(0, ggf);
363  double eta = gauss_pts_triangles_only(1, ggf);
364  coords[0] = ksi;
365  coords[1] = eta;
366  for (int ggt = 0; ggt != 2; ggt++) {
367  double zeta = gauss_pts_through_thickness(0, ggt);
368  coords[2] = zeta;
369  int side = ggt * 3 + ggf;
370  CHKERR postProcMesh.create_vertex(&coords[0], prism_conn[side]);
371  }
372  }
373  CHKERR postProcMesh.create_element(MBPRISM, &prism_conn[0], 6, prism);
374 
375  elementsMap[numeredEntFiniteElementPtr->getEnt()] = prism;
376  // Range faces;
377  // CHKERR postProcMesh.get_adjacencies(&prism,1,2,true,faces);
378  Range edges;
379  CHKERR postProcMesh.get_adjacencies(&prism, 1, 1, true, edges);
380  EntityHandle meshset;
381  CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
382  CHKERR postProcMesh.add_entities(meshset, &prism, 1);
383  // CHKERR postProcMesh.add_entities(meshset,faces);
384  CHKERR postProcMesh.add_entities(meshset, edges);
385  if (tenNodesPostProcTets) {
386  CHKERR postProcMesh.convert_entities(meshset, true, false, false);
387  }
388  CHKERR postProcMesh.delete_entities(&meshset, 1);
389  CHKERR postProcMesh.delete_entities(edges);
390  // CHKERR postProcMesh.delete_entities(faces);
391 
392  CHKERR mField.get_moab().get_connectivity(
393  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
394  MatrixDouble coords_prism_global;
395  coords_prism_global.resize(num_nodes, 3, false);
396  CHKERR mField.get_moab().get_coords(conn, num_nodes,
397  &coords_prism_global(0, 0));
398 
399  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
400  MatrixDouble coords_prism_local;
401  coords_prism_local.resize(num_nodes, 3, false);
402  CHKERR postProcMesh.get_coords(conn, num_nodes, &coords_prism_local(0, 0));
403 
404  if (gaussPtsThroughThickness.size2() != nb_through_thickness) {
405  pointsMap.clear();
406  gaussPtsTrianglesOnly.resize(3, nb_on_triangle, false);
407  gaussPtsTrianglesOnly.clear();
408  gaussPtsThroughThickness.resize(2, nb_through_thickness, false);
409  gaussPtsThroughThickness.clear();
410  const double eps = 1e-6;
411  int ggf = 0, ggt = 0;
412  for (int nn = 0; nn != num_nodes; nn++) {
413  double ksi = coords_prism_local(nn, 0);
414  double eta = coords_prism_local(nn, 1);
415  double zeta = coords_prism_local(nn, 2);
416  pointsMap.insert(PointsMap3D(ksi * 100, eta * 100, zeta * 100, nn));
417  if (fabs(zeta) < eps) {
418  gaussPtsTrianglesOnly(0, ggf) = ksi;
419  gaussPtsTrianglesOnly(1, ggf) = eta;
420  ggf++;
421  }
422  if (fabs(ksi) < eps && fabs(eta) < eps) {
423  gaussPtsThroughThickness(0, ggt) = zeta;
424  ggt++;
425  }
426  }
427  }
428 
429  for (int nn = 0; nn != num_nodes; nn++) {
430  const double ksi = coords_prism_local(nn, 0);
431  const double eta = coords_prism_local(nn, 1);
432  const double zeta = coords_prism_local(nn, 2);
433  const double n0 = N_MBTRI0(ksi, eta);
434  const double n1 = N_MBTRI1(ksi, eta);
435  const double n2 = N_MBTRI2(ksi, eta);
436  const double e0 = N_MBEDGE0(zeta);
437  const double e1 = N_MBEDGE1(zeta);
438  double coords_global[3];
439  for (int dd = 0; dd != 3; dd++) {
440  coords_global[dd] = n0 * e0 * coords_prism_global(0, dd) +
441  n1 * e0 * coords_prism_global(1, dd) +
442  n2 * e0 * coords_prism_global(2, dd) +
443  n0 * e1 * coords_prism_global(3, dd) +
444  n1 * e1 * coords_prism_global(4, dd) +
445  n2 * e1 * coords_prism_global(5, dd);
446  }
447  // cerr << coords_global[0] << " " << coords_global[1] << " " <<
448  // coords_global[2] << endl;
449  CHKERR postProcMesh.set_coords(&conn[nn], 1, coords_global);
450  }
451  }
452 
453  mapGaussPts.clear();
454  mapGaussPts.resize(nb_through_thickness * nb_on_triangle);
455  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
456  fill(mapGaussPts.begin(), mapGaussPts.end(), 0);
457  {
458  int gg = 0;
459  for (unsigned int ggf = 0; ggf != gaussPtsTrianglesOnly.size2(); ggf++) {
460  const double ksi = gaussPtsTrianglesOnly(0, ggf);
461  const double eta = gaussPtsTrianglesOnly(1, ggf);
462  for (unsigned int ggt = 0; ggt != gaussPtsThroughThickness.size2();
463  ggt++, gg++) {
464  const double zeta = gaussPtsThroughThickness(0, ggt);
465  PointsMap3D_multiIndex::iterator it;
466  it =
467  pointsMap.find(boost::make_tuple(ksi * 100, eta * 100, zeta * 100));
468  if (it != pointsMap.end()) {
469  mapGaussPts[gg] = conn[it->nN];
470  }
471  }
472  }
473  }
474 
476 }
477 
479  int order_thickness) {
481  if (gaussPtsThroughThickness.size1() == 0 ||
482  gaussPtsThroughThickness.size2() == 0) {
483  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
484  "post-process mesh not generated");
485  }
487 }
488 
491  // MoAB
492  ParallelComm *pcomm_post_proc_mesh =
493  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
494  if (pcomm_post_proc_mesh != NULL) {
495  delete pcomm_post_proc_mesh;
496  }
497  // CHKERR postProcMesh.delete_mesh();
499 }
500 
503  ParallelComm *pcomm =
504  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
505  ParallelComm *pcomm_post_proc_mesh =
506  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
507  if (pcomm_post_proc_mesh == NULL) {
508  pcomm_post_proc_mesh = new ParallelComm(&postProcMesh, mField.get_comm());
509  }
510  Range prims;
511  CHKERR postProcMesh.get_entities_by_type(0, MBPRISM, prims, false);
512  // std::cerr << "total prims size " << prims.size() << std::endl;
513  int rank = pcomm->rank();
514  Range::iterator pit = prims.begin();
515  for (; pit != prims.end(); pit++) {
516  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
517  &rank);
518  }
519  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
521 }
522 
525 
526  auto generate_tri = [&](auto &gauss_pts) {
528  gauss_pts.resize(3, 3, false);
529  gauss_pts.clear();
530  gauss_pts(0, 0) = 0;
531  gauss_pts(1, 0) = 0;
532  gauss_pts(0, 1) = 1;
533  gauss_pts(1, 1) = 0;
534  gauss_pts(0, 2) = 0;
535  gauss_pts(1, 2) = 1;
536 
537  moab::Core core_ref;
538  moab::Interface &moab_ref = core_ref;
539  const EntityHandle *conn;
540  int num_nodes;
541  EntityHandle tri_conn[3];
542  MatrixDouble coords(6, 3);
543  for (int gg = 0; gg != 3; gg++) {
544  coords(gg, 0) = gauss_pts(0, gg);
545  coords(gg, 1) = gauss_pts(1, gg);
546  coords(gg, 2) = 0;
547  CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
548  }
549 
550  EntityHandle tri;
551  CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
552 
553  if (sixNodePostProcTris) {
554  Range edges;
555  CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges);
556  EntityHandle meshset;
557  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
558  CHKERR moab_ref.add_entities(meshset, &tri, 1);
559  CHKERR moab_ref.add_entities(meshset, edges);
560  CHKERR moab_ref.convert_entities(meshset, true, false, false);
561  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
562  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
563  gauss_pts.resize(3, num_nodes, false);
564  gauss_pts.clear();
565  for (int nn = 0; nn < 3; nn++) {
566  gauss_pts(0, nn) = coords(nn, 0);
567  gauss_pts(1, nn) = coords(nn, 1);
568  gauss_pts(0, 3 + nn) = coords(3 + nn, 0);
569  gauss_pts(1, 3 + nn) = coords(3 + nn, 1);
570  }
571  } else {
572  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
573  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
574  gauss_pts.resize(3, num_nodes, false);
575  gauss_pts.clear();
576  for (int nn = 0; nn < 3; nn++) {
577  gauss_pts(0, nn) = coords(nn, 0);
578  gauss_pts(1, nn) = coords(nn, 1);
579  }
580  }
582  };
583 
584  auto generate_quad = [&](auto &gauss_pts) {
586 
587  gauss_pts.resize(3, 4, false);
588  gauss_pts.clear();
589 
590  gauss_pts(0, 0) = 0;
591  gauss_pts(1, 0) = 0;
592  gauss_pts(0, 1) = 1;
593  gauss_pts(1, 1) = 0;
594  gauss_pts(0, 2) = 1;
595  gauss_pts(1, 2) = 1;
596  gauss_pts(0, 3) = 0;
597  gauss_pts(1, 3) = 1;
598 
599  moab::Core core_ref;
600  moab::Interface &moab_ref = core_ref;
601 
602  const EntityHandle *conn;
603  int num_nodes;
604  std::array<EntityHandle, 4> quad_conn;
605  MatrixDouble coords(8, 3);
606  for (int gg = 0; gg != 4; gg++) {
607  coords(gg, 0) = gauss_pts(0, gg);
608  coords(gg, 1) = gauss_pts(1, gg);
609  coords(gg, 2) = 0;
610  CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
611  }
612 
613  EntityHandle quad;
614  CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
615 
616  if (sixNodePostProcTris) {
617  Range edges;
618  CHKERR moab_ref.get_adjacencies(&quad, 1, 1, true, edges);
619  EntityHandle meshset;
620  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
621  CHKERR moab_ref.add_entities(meshset, &quad, 1);
622  CHKERR moab_ref.add_entities(meshset, edges);
623  CHKERR moab_ref.convert_entities(meshset, true, false, false);
624  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
625  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
626  gauss_pts.resize(3, num_nodes, false);
627  gauss_pts.clear();
628  for (int nn = 0; nn != 4; nn++) {
629  gauss_pts(0, nn) = coords(nn, 0);
630  gauss_pts(1, nn) = coords(nn, 1);
631  gauss_pts(0, 4 + nn) = coords(4 + nn, 0);
632  gauss_pts(1, 4 + nn) = coords(4 + nn, 1);
633  }
634  } else {
635  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
636  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
637  gauss_pts.resize(3, num_nodes, false);
638  gauss_pts.clear();
639  for (int nn = 0; nn != 4; nn++) {
640  gauss_pts(0, nn) = coords(nn, 0);
641  gauss_pts(1, nn) = coords(nn, 1);
642  }
643  }
644 
646  };
647 
648  CHKERR generate_tri(gaussPtsTri);
649  CHKERR generate_quad(gaussPtsQuad);
650 
652 }
653 
656  if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
657  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
658  "post-process mesh not generated");
659 
660  auto create_tri = [&](auto &gauss_pts) {
661  EntityHandle tri;
662  std::array<EntityHandle, 3> tri_conn;
663  MatrixDouble3by3 coords_tri(3, 3);
664  VectorDouble3 coords(3);
665  CHKERR mField.get_moab().get_connectivity(
666  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
667  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_tri(0, 0));
668  for (int gg = 0; gg != 3; gg++) {
669  double ksi = gauss_pts(0, gg);
670  double eta = gauss_pts(1, gg);
671  double n0 = N_MBTRI0(ksi, eta);
672  double n1 = N_MBTRI1(ksi, eta);
673  double n2 = N_MBTRI2(ksi, eta);
674  double x =
675  n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
676  double y =
677  n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
678  double z =
679  n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
680  coords[0] = x;
681  coords[1] = y;
682  coords[2] = z;
683  CHKERR postProcMesh.create_vertex(&coords[0], tri_conn[gg]);
684  }
685  CHKERR postProcMesh.create_element(MBTRI, &tri_conn[0], 3, tri);
686  return tri;
687  };
688 
689  auto create_quad = [&](auto &gauss_pts) {
690  EntityHandle quad;
691  std::array<EntityHandle, 4> quad_conn;
692  MatrixDouble coords_quad(4, 3);
693  VectorDouble coords(4);
694  CHKERR mField.get_moab().get_connectivity(
695  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
696  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_quad(0, 0));
697  for (int gg = 0; gg != 4; gg++) {
698  double ksi = gauss_pts(0, gg);
699  double eta = gauss_pts(1, gg);
700  double n0 = N_MBQUAD0(ksi, eta);
701  double n1 = N_MBQUAD1(ksi, eta);
702  double n2 = N_MBQUAD2(ksi, eta);
703  double n3 = N_MBQUAD3(ksi, eta);
704  double x = n0 * coords_quad(0, 0) + n1 * coords_quad(1, 0) +
705  n2 * coords_quad(2, 0) + n3 * coords_quad(3, 0);
706  double y = n0 * coords_quad(0, 1) + n1 * coords_quad(1, 1) +
707  n2 * coords_quad(2, 1) + n3 * coords_quad(3, 1);
708  double z = n0 * coords_quad(0, 2) + n1 * coords_quad(1, 2) +
709  n2 * coords_quad(2, 2) + n3 * coords_quad(3, 2);
710  coords[0] = x;
711  coords[1] = y;
712  coords[2] = z;
713  CHKERR postProcMesh.create_vertex(&coords[0], quad_conn[gg]);
714  }
715  CHKERR postProcMesh.create_element(MBQUAD, &quad_conn[0], 4, quad);
716  return quad;
717  };
718 
719  auto add_mid_nodes = [&](auto tri) {
721  Range edges;
722  CHKERR postProcMesh.get_adjacencies(&tri, 1, 1, true, edges);
723  EntityHandle meshset;
724  CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
725  CHKERR postProcMesh.add_entities(meshset, &tri, 1);
726  CHKERR postProcMesh.add_entities(meshset, edges);
727  CHKERR postProcMesh.convert_entities(meshset, true, false, false);
728  CHKERR postProcMesh.delete_entities(&meshset, 1);
729  CHKERR postProcMesh.delete_entities(edges);
731  };
732 
733  EntityHandle tri;
734 
735  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
736  elementsMap.end()) {
737  tri = elementsMap[numeredEntFiniteElementPtr->getEnt()];
738  switch (numeredEntFiniteElementPtr->getEntType()) {
739  case MBTRI:
740  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
741  noalias(gaussPts) = gaussPtsTri;
742  break;
743  case MBQUAD:
744  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
745  noalias(gaussPts) = gaussPtsQuad;
746  break;
747  default:
748  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
749  "Not implemented element type");
750  }
751  } else {
752  switch (numeredEntFiniteElementPtr->getEntType()) {
753  case MBTRI:
754  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
755  noalias(gaussPts) = gaussPtsTri;
756  tri = create_tri(gaussPtsTri);
757  break;
758  case MBQUAD:
759  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
760  noalias(gaussPts) = gaussPtsQuad;
761  tri = create_quad(gaussPtsQuad);
762  break;
763  default:
764  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
765  "Not implemented element type");
766  }
767  if (sixNodePostProcTris)
768  CHKERR add_mid_nodes(tri);
769  elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
770  }
771 
772  // Set values which map nodes with integration points on the prism
773  const EntityHandle *conn;
774  int num_nodes;
775 
776  CHKERR postProcMesh.get_connectivity(tri, conn, num_nodes, false);
777  mapGaussPts.resize(num_nodes);
778  for (int nn = 0; nn != num_nodes; nn++)
779  mapGaussPts[nn] = conn[nn];
780 
782 }
783 
786  ParallelComm *pcomm_post_proc_mesh =
787  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
788  if (pcomm_post_proc_mesh != NULL) {
789  delete pcomm_post_proc_mesh;
790  }
792 }
793 
796  ParallelComm *pcomm =
797  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
798  ParallelComm *pcomm_post_proc_mesh =
799  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
800  if (pcomm_post_proc_mesh == NULL) {
801  pcomm_post_proc_mesh = new ParallelComm(&postProcMesh, mField.get_comm());
802  }
803  Range tris;
804  CHKERR postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
805  int rank = pcomm->rank();
806  Range::iterator pit = tris.begin();
807  for (; pit != tris.end(); pit++) {
808  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
809  &rank);
810  }
811  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
813 }
814 
817  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
819 
820  if (type != MBVERTEX)
822 
823  CHKERR loopSideVolumes(feVolName, *sideOpFe);
824 
825  // quit if tag is not needed
826  if (!saveOnTag)
828 
829  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
830  int rank = dof_ptr->getNbOfCoeffs();
831 
832  int tag_length = rank * 3;
833  FieldSpace space = dof_ptr->getSpace();
834 
835  double def_VAL[tag_length];
836  bzero(def_VAL, tag_length * sizeof(double));
837  Tag th;
838  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
839  MB_TYPE_DOUBLE, th,
840  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
841 
842  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
843  // no need for L2
844  const void *tags_ptr[mapGaussPts.size()];
845  int nb_gauss_pts = data.getN().size1();
846  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts ||
847  nb_gauss_pts != gradMatPtr->size2()) {
848  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
849  }
850  switch (space) {
851  case H1:
852  case L2:
853 
854  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
855  tags_ptr);
856  // FIXME: this is not very efficient
857  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
858  for (int rr = 0; rr != rank; ++rr) {
859  for (int dd = 0; dd != 3; ++dd) {
860  const double *my_ptr2 = static_cast<const double *>(tags_ptr[gg]);
861  double *my_ptr = const_cast<double *>(my_ptr2);
862  my_ptr[rank * rr + dd] = (*gradMatPtr)(rank * rr + dd, gg);
863  }
864  }
865  }
866  break;
867  default:
868  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
869  "field with that space is not implemented");
870  }
871 
873 }
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:67
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
field with continuous normal traction
Definition: definitions.h:173
Post-process fields on refined mesh.
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:80
FieldCoefficientsNumber getNbOfCoeffs() const
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode postProcess()
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
keeps information about indexed dofs for the finite element
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode generateReferenceElementMesh()
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:68
FieldSpace getSpace() const
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:79
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:65
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
field with continuous tangents
Definition: definitions.h:172
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:95
FieldSpace
approximation spaces
Definition: definitions.h:168
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
function is run at the end of loop
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
MoFEMErrorCode preProcess()
function is run at the beginning of loop
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:66
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setGaussPts(int order)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:55
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:54
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
field with C-1 continuity
Definition: definitions.h:174