12 using namespace MoFEM;
13 using namespace boost::numeric;
27 VectorDofs::iterator it, hi_it;
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())
34 ->getPetscLocalDofIdx();
35 vAlues[ii] =
a[local_idx];
44 int rank = dof_ptr->getNbOfCoeffs();
46 int tag_length = rank;
58 "field with that space is not implemented");
61 if (tag_length > 1 && tag_length < 3)
63 else if (tag_length > 3 && tag_length < 9)
66 double def_VAL[tag_length];
67 bzero(def_VAL, tag_length *
sizeof(
double));
69 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
71 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
75 const void *tags_ptr[mapGaussPts.size()];
76 int nb_gauss_pts = data.
getN().size1();
77 if (mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
79 "data inconsistency %d!=%d", mapGaussPts.size(), nb_gauss_pts);
82 auto set_float_precision = [](
const double x) {
83 if (std::abs(x) < std::numeric_limits<float>::epsilon())
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();
99 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
101 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
102 for (
int rr = 0; rr < rank; rr++) {
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);
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();
118 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
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++) {
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);
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);
138 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
143 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
144 double *ptr = &((
double *)tags_ptr[gg])[0];
146 for (; ll !=
static_cast<int>(vAluesPtr->size() / rank); ll++) {
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);
155 for (; ll !=
static_cast<int>(data.
getN().size2() / 3); ll++) {
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);
168 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
173 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
174 double *ptr = &((
double *)tags_ptr[gg])[0];
176 for (; ll !=
static_cast<int>(vAluesPtr->size() / rank); ll++) {
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);
185 for (; ll !=
static_cast<int>(data.
getN().size2() / 3); ll++) {
193 "field with that space is not implemented");
209 VectorDofs::iterator it, hi_it;
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())
216 ->getPetscLocalDofIdx();
217 vAlues[ii] =
a[local_idx];
226 int rank = dof_ptr->getNbOfCoeffs();
229 int space_dim = spaceDim;
233 int tag_length = rank * space_dim;
241 tag_length *= space_dim;
245 "field with that space is not implemented");
248 double def_VAL[tag_length];
249 bzero(def_VAL, tag_length *
sizeof(
double));
251 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
253 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
257 const void *tags_ptr[mapGaussPts.size()];
258 int nb_gauss_pts = data.
getN().size1();
259 if (mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
263 auto set_float_precision = [](
const double x) {
264 if (std::abs(x) < std::numeric_limits<float>::epsilon())
270 auto clear_vals = [&]() {
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);
280 auto set_vals = [&]() {
282 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
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++) {
289 data.
getDiffN(gg)(dof,
dd) * (*vAluesPtr)[rank * dof + rr];
290 (commonData.gradMap[rowFieldName])[gg](rr,
dd) += val;
292 ((
double *)tags_ptr[gg])[rank * rr +
dd] =
293 (commonData.gradMap[rowFieldName])[gg](rr,
dd);
300 commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
304 if (
type == MBVERTEX)
314 "field with that space is not implemented");
326 int order_triangles_only) {
333 int max_ref_level_thick = 0;
334 PetscBool flg = PETSC_TRUE;
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;
342 const int nb_on_triangle = 6;
343 const int nb_through_thickness = 3;
345 std::vector<PointsMap3D_multiIndex> pointsMapVector;
348 if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
350 pointsMapVector = pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()];
352 pointsMapVector.resize(0);
353 gaussPtsTrianglesOnly.resize(3, nb_on_triangle,
false);
354 gaussPtsTrianglesOnly.clear();
355 gaussPtsThroughThickness.resize(2, nb_through_thickness * level_thick,
357 gaussPtsThroughThickness.clear();
359 MatrixDouble gauss_pts_triangles_only, gauss_pts_through_thickness;
361 double inc_thick = 1.0 / (
double)level_thick;
362 for (
int ll = 0; ll < level_thick; ll++) {
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;
376 ublas::vector<EntityHandle> prism_conn(6);
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);
383 for (
int ggt = 0; ggt != 2; ggt++) {
384 double zeta = gauss_pts_through_thickness(0, ggt);
386 int side = ggt * 3 + ggf;
387 CHKERR postProcMesh.create_vertex(&coords[0], prism_conn[side]);
390 CHKERR postProcMesh.create_element(MBPRISM, &prism_conn[0], 6, prism);
392 elementsMap[numeredEntFiniteElementPtr->getEnt()].push_back(prism);
396 CHKERR postProcMesh.get_adjacencies(&prism, 1, 1,
true, edges);
398 CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
399 CHKERR postProcMesh.add_entities(meshset, &prism, 1);
401 CHKERR postProcMesh.add_entities(meshset, edges);
402 if (tenNodesPostProcTets) {
403 CHKERR postProcMesh.convert_entities(meshset,
true,
false,
false);
405 CHKERR postProcMesh.delete_entities(&meshset, 1);
406 CHKERR postProcMesh.delete_entities(edges);
409 CHKERR mField.get_moab().get_connectivity(
410 numeredEntFiniteElementPtr->getEnt(), conn, num_nodes,
true);
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));
416 CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes,
false);
418 coords_prism_local.resize(num_nodes, 3,
false);
419 CHKERR postProcMesh.get_coords(conn, num_nodes,
420 &coords_prism_local(0, 0));
422 const double eps = 1e-6;
423 int ggf = 0, ggt = 0;
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);
433 gaussPtsTrianglesOnly(0, ggf) = ksi;
434 gaussPtsTrianglesOnly(1, ggf) =
eta;
437 if (fabs(ksi) <
eps && fabs(
eta) <
eps) {
438 gaussPtsThroughThickness(0, ggt + nb_through_thickness * ll) =
zeta;
443 pointsMapVector.push_back(pointsMap);
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);
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);
463 CHKERR postProcMesh.set_coords(&conn[nn], 1, coords_global);
467 pointsMapVectorMap[numeredEntFiniteElementPtr->getEnt()] = pointsMapVector;
471 mapGaussPts.resize(nb_through_thickness * nb_on_triangle * level_thick);
472 fill(mapGaussPts.begin(), mapGaussPts.end(), 0);
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();
482 prism = elementsMap[numeredEntFiniteElementPtr->getEnt()][ll];
483 pointsMap = pointsMapVector[ll];
484 CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes,
false);
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];
501 int order_thickness) {
503 if (gaussPtsThroughThickness.size1() == 0 ||
504 gaussPtsThroughThickness.size2() == 0) {
506 "post-process mesh not generated");
518 ParallelComm *pcomm_post_proc_mesh =
520 if (pcomm_post_proc_mesh == NULL) {
522 pcomm_post_proc_mesh =
523 new ParallelComm(&postProcMesh, PETSC_COMM_WORLD);
527 CHKERR postProcMesh.get_entities_by_type(0, MBPRISM, prims,
false);
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,
535 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
542 auto generate_tri = [&](
auto &gauss_pts) {
544 gauss_pts.resize(3, 3,
false);
559 for (
int gg = 0; gg != 3; gg++) {
560 coords(gg, 0) = gauss_pts(0, gg);
561 coords(gg, 1) = gauss_pts(1, gg);
563 CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
567 CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
569 if (sixNodePostProcTris) {
571 CHKERR moab_ref.get_adjacencies(&tri, 1, 1,
true, edges);
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);
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);
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);
590 for (
int nn = 0; nn < 3; nn++) {
591 gauss_pts(0, nn) = coords(nn, 0);
592 gauss_pts(1, nn) = coords(nn, 1);
598 auto generate_quad = [&](
auto &gauss_pts) {
601 gauss_pts.resize(3, 4,
false);
618 std::array<EntityHandle, 4> quad_conn;
620 for (
int gg = 0; gg != 4; gg++) {
621 coords(gg, 0) = gauss_pts(0, gg);
622 coords(gg, 1) = gauss_pts(1, gg);
624 CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
628 CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
630 if (sixNodePostProcTris) {
632 CHKERR moab_ref.get_adjacencies(&quad, 1, 1,
true, edges);
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);
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);
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);
651 for (
int nn = 0; nn != 4; nn++) {
652 gauss_pts(0, nn) = coords(nn, 0);
653 gauss_pts(1, nn) = coords(nn, 1);
660 CHKERR generate_tri(gaussPtsTri);
661 CHKERR generate_quad(gaussPtsQuad);
668 if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
670 "post-process mesh not generated");
672 auto create_tri = [&](
auto &gauss_pts) {
673 std::array<EntityHandle, 3> tri_conn;
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();
680 for (
int gg = 0; gg != num_nodes_on_ele; gg++) {
682 const double ksi = gauss_pts(0, gg);
683 const double eta = gauss_pts(1, gg);
688 n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
690 n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
692 n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
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;
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];
706 const auto tri = startingEleTriHandle + counterTris;
711 auto create_quad = [&](
auto &gauss_pts) {
712 std::array<EntityHandle, 4> quad_conn;
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);
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;
736 mapGaussPts.resize(num_nodes_on_ele);
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];
744 const auto quad = startingEleQuadHandle + counterQuads;
750 if (elementsMap.size() == getLoopSize()) {
752 tri = elementsMap.at(numeredEntFiniteElementPtr->getEnt());
753 switch (numeredEntFiniteElementPtr->getEntType()) {
755 gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(),
false);
756 noalias(gaussPts) = gaussPtsTri;
759 gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(),
false);
760 noalias(gaussPts) = gaussPtsQuad;
764 "Not implemented element type");
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];
776 switch (numeredEntFiniteElementPtr->getEntType()) {
778 gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(),
false);
779 noalias(gaussPts) = gaussPtsTri;
780 tri = create_tri(gaussPtsTri);
784 gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(),
false);
785 noalias(gaussPts) = gaussPtsQuad;
786 tri = create_quad(gaussPtsQuad);
791 "Not implemented element type");
794 elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
803 ReadUtilIface *iface;
804 CHKERR postProcMesh.query_interface(iface);
806 const int number_of_ents_in_the_loop = getLoopSize();
807 if (elementsMap.size() != number_of_ents_in_the_loop) {
810 postProcMesh.delete_mesh();
812 auto get_number_of_computational_elements = [&]() {
813 auto fe_ptr = this->problemPtr->numeredFiniteElementsPtr;
816 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().lower_bound(
817 boost::make_tuple(this->getFEName(), this->getLoFERank()));
819 fe_ptr->template get<Composite_Name_And_Part_mi_tag>().upper_bound(
820 boost::make_tuple(this->getFEName(), this->getHiFERank()));
822 const int number_of_ents_in_the_loop = this->getLoopSize();
823 if (std::distance(miit, hi_miit) != number_of_ents_in_the_loop) {
825 "Wrong size of indicices. Inconsistent size number of iterated "
826 "elements iterated by problem and from range.");
829 std::array<int, MBMAXTYPE> nb_elemms_by_type;
830 std::fill(nb_elemms_by_type.begin(), nb_elemms_by_type.end(), 0);
832 for (; miit != hi_miit; ++miit) {
836 numeredEntFiniteElementPtr = *miit;
837 add = exeTestHook(
this);
841 auto type = (*miit)->getEntType();
842 ++nb_elemms_by_type[
type];
846 return nb_elemms_by_type;
849 auto nb_computational_elements_by_type =
850 get_number_of_computational_elements();
852 const int numberOfTriangles = nb_computational_elements_by_type[MBTRI];
853 const int numberOfQuads = nb_computational_elements_by_type[MBQUAD];
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();
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,
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,
888 auto update_elements = [&]() {
890 ReadUtilIface *iface;
891 CHKERR postProcMesh.query_interface(iface);
894 CHKERR iface->update_adjacencies(startingEleTriHandle, counterTris,
895 gaussPtsTri.size2(), triConn);
897 CHKERR iface->update_adjacencies(startingEleQuadHandle, counterQuads,
898 gaussPtsQuad.size2(), quadConn);
902 auto resolve_shared = [&]() {
904 ParallelComm *pcomm_post_proc_mesh =
906 if (pcomm_post_proc_mesh == NULL) {
909 pcomm_post_proc_mesh =
new ParallelComm(
910 &postProcMesh, PETSC_COMM_WORLD );
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,
918 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
934 if (
type != MBVERTEX)
937 CHKERR loopSideVolumes(feVolName, *sideOpFe);
943 auto &m_field = getPtrFE()->mField;
944 auto field_ptr = m_field.get_field_structure(fieldName);
945 const int rank = field_ptr->getNbOfCoeffs();
952 int full_size = rank * RANK;
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));
960 CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
962 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
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()) {
977 CHKERR postProcMesh.tag_get_by_ptr(
th, &mapGaussPts[0], mapGaussPts.size(),
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);
993 "field with that space is not implemented");
1000 const std::string
field_name,
const std::string vol_fe_name,
1001 boost::shared_ptr<MatrixDouble> grad_mat_ptr,
bool save_on_tag) {
1005 grad_mat_ptr = boost::make_shared<MatrixDouble>();
1007 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
1008 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1010 if (mField.check_field(
"MESH_NODE_POSITIONS"))
1015 auto field_ptr = mField.get_field_structure(
field_name);
1016 const int nb_coefficients = field_ptr->getNbOfCoeffs();
1018 switch (nb_coefficients) {
1020 my_side_fe->getOpPtrVector().push_back(
1024 my_side_fe->getOpPtrVector().push_back(
1028 my_side_fe->getOpPtrVector().push_back(
1034 "field with that number of coefficients is not implemented");
1040 my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
1046 const std::string
field_name,
const std::string vol_fe_name,
1047 boost::shared_ptr<MatrixDouble> mat_ptr,
bool save_on_tag) {
1051 mat_ptr = boost::make_shared<MatrixDouble>();
1054 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
1057 auto field_ptr = mField.get_field_structure(
field_name);
1058 const int nb_coefficients = field_ptr->getNbOfCoeffs();
1064 switch (nb_coefficients) {
1066 my_side_fe->getOpPtrVector().push_back(
1070 my_side_fe->getOpPtrVector().push_back(
1074 my_side_fe->getOpPtrVector().push_back(
1078 my_side_fe->getOpPtrVector().push_back(
1083 "field with that number of coefficients is not implemented");
1087 switch (nb_coefficients) {
1089 my_side_fe->getOpPtrVector().push_back(
1093 my_side_fe->getOpPtrVector().push_back(
1098 "field with that number of coefficients is not implemented");
1103 "field with that space is not implemented.");
1109 mat_ptr, save_on_tag));
1117 gaussPts.resize(3, 2,
false);
1123 mapGaussPts.resize(gaussPts.size2());
1131 for (
int gg = 0; gg != 2; gg++) {
1132 coords(gg, 0) = gaussPts(0, gg);
1133 coords(gg, 1) = gaussPts(1, gg);
1135 CHKERR moab_ref.create_vertex(&coords(gg, 0), edge_conn[gg]);
1139 CHKERR moab_ref.create_element(MBEDGE, edge_conn, 2, edge);
1145 ParallelComm *pcomm_post_proc_mesh =
1147 if (pcomm_post_proc_mesh != NULL) {
1148 delete pcomm_post_proc_mesh;
1155 ParallelComm *pcomm_post_proc_mesh =
1157 if (pcomm_post_proc_mesh == NULL) {
1159 pcomm_post_proc_mesh =
new ParallelComm(
1160 &postProcMesh, PETSC_COMM_WORLD );
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,
1171 set_edges_rank(rank, edges);
1173 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
1179 if (gaussPts.size1() == 0) {
1181 "post-process mesh not generated");
1188 if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
1189 elementsMap.end()) {
1190 edge = elementsMap[numeredEntFiniteElementPtr->getEnt()];
1192 ublas::vector<EntityHandle> edge_conn(2);
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);
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);
1212 CHKERR postProcMesh.create_vertex(&coords[0], edge_conn[gg]);
1214 CHKERR postProcMesh.create_element(MBEDGE, &edge_conn[0], 2, edge);
1215 elementsMap[numeredEntFiniteElementPtr->getEnt()] = edge;
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];