v0.15.0
Loading...
Searching...
No Matches
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>
12using namespace MoFEM;
13using namespace boost::numeric;
14#include <PostProcOnRefMesh.hpp>
15
16MoFEMErrorCode 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 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79 "data inconsistency %ld!=%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
199MoFEMErrorCode 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
324
325MoFEMErrorCode 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_NULLPTR, PETSC_NULLPTR,
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
500MoFEMErrorCode 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
515
516MoFEMErrorCode 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
539MoFEMErrorCode 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
666MoFEMErrorCode 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
800MoFEMErrorCode 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
885MoFEMErrorCode 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
928template <int RANK>
929MoFEMErrorCode
930PostProcFaceOnRefinedMesh::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
999MoFEMErrorCode 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
1045MoFEMErrorCode 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
1114MoFEMErrorCode 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
1143MoFEMErrorCode 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
1153MoFEMErrorCode 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
1177MoFEMErrorCode 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}
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
FieldSpace
FieldSpace
approximation spaces
Definition definitions.h:82
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition definitions.h:577
order
constexpr int order
Definition dg_projection.cpp:19
N_MBTRI1
#define N_MBTRI1(x, y)
triangle shape function
Definition fem_tools.h:47
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition fem_tools.h:60
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition fem_tools.h:105
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition fem_tools.h:57
N_MBTRI0
#define N_MBTRI0(x, y)
triangle shape function
Definition fem_tools.h:46
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition fem_tools.h:59
N_MBTRI2
#define N_MBTRI2(x, y)
triangle shape function
Definition fem_tools.h:48
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition fem_tools.h:58
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition fem_tools.h:106
eta
double eta
Definition free_surface.cpp:196
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition HODataOperators.hpp:986
zeta
double zeta
Viscous hardening.
Definition plastic.cpp:130
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
g
constexpr double g
Definition shallow_wave.cpp:63
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:152
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition EntitiesFieldData.hpp:1513
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:1486
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
Get DOF values on entity.
Definition EntitiesFieldData.hpp:1422
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
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
Get DOF data structures (const version)
Definition EntitiesFieldData.hpp:1437
MoFEM::FEMethod::getRowDofsPtr
auto getRowDofsPtr() const
Get pointer to row DOFs for the current finite element.
Definition LoopMethods.hpp:791
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition ForcesAndSourcesCore.hpp:1151
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition UserDataOperators.hpp:3020
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition UserDataOperators.hpp:2402
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1490
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor field rank 2, i.e. matrix field.
Definition UserDataOperators.hpp:1084
MoFEM::OpCalculateTensor2SymmetricFieldGradient
Get field gradients at integration pts for symmetric tensorial field rank 2.
Definition UserDataOperators.hpp:1942
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition UserDataOperators.hpp:1173
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1745
MoFEM::OpCalculateVectorFieldValues
Specialization for double precision vector field values calculation.
Definition UserDataOperators.hpp:626
PostProcCommonOnRefMesh::CommonData::fieldMap
std::map< std::string, std::vector< VectorDouble > > fieldMap
Definition PostProcOnRefMesh.hpp:33
PostProcCommonOnRefMesh::OpGetFieldGradientValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition PostProcOnRefMesh.cpp:199
PostProcCommonOnRefMesh::OpGetFieldValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition PostProcOnRefMesh.cpp:16
PostProcCommonOnRefMesh::OpGetFieldValues::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition PostProcOnRefMesh.hpp:53
PostProcEdgeOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Definition PostProcOnRefMesh.cpp:1153
PostProcEdgeOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
Definition PostProcOnRefMesh.cpp:1143
PostProcEdgeOnRefinedMesh::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition PostProcOnRefMesh.cpp:1177
PostProcEdgeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition PostProcOnRefMesh.cpp:1114
PostProcFaceOnRefinedMesh::OpGetFieldValuesOnSkinImpl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition PostProcOnRefMesh.cpp:930
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
PostProcFaceOnRefinedMesh::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition PostProcOnRefMesh.cpp:666
PostProcFaceOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition PostProcOnRefMesh.cpp:539
PostProcFaceOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
Definition PostProcOnRefMesh.cpp:800
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
PostProcFaceOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Definition PostProcOnRefMesh.cpp:885
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:1005
PostProcFatPrismOnRefinedMesh::setGaussPtsTrianglesOnly
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
Definition PostProcOnRefMesh.cpp:325
PostProcFatPrismOnRefinedMesh::preProcess
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
Definition PostProcOnRefMesh.cpp:511
PostProcFatPrismOnRefinedMesh::setGaussPtsThroughThickness
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
Definition PostProcOnRefMesh.cpp:500
PostProcFatPrismOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition PostProcOnRefMesh.cpp:320
PostProcFatPrismOnRefinedMesh::postProcess
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Definition PostProcOnRefMesh.cpp:516
Generated by Doxygen 1.9.8 and hosted at University of Glasgow