13using namespace boost::numeric;
27 VectorDofs::iterator it, hi_it;
31 for (
int ii = 0; it != hi_it; it++, ii++) {
32 int local_idx = row_dofs->find((*it)->getLocalUniqueId())
34 ->getPetscLocalDofIdx();
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));
71 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
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())
92 if (type == MBVERTEX) {
93 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
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,
107 ((
double *)tags_ptr[gg])[rr] += set_float_precision(val);
113 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
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,
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++) {
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++) {
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);
554 moab::Interface &moab_ref = core_ref;
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);
614 moab::Interface &moab_ref = core_ref;
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");
1037 FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
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.");
1106 FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
1109 mat_ptr, save_on_tag));
1117 gaussPts.resize(3, 2,
false);
1123 mapGaussPts.resize(gaussPts.size2());
1125 moab::Core core_ref;
1126 moab::Interface &moab_ref = core_ref;
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];
Post-process fields on refined mesh.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr auto field_name
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
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
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
auto getRowDofsPtr() const
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get vector field for H-div approximation.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Get field gradients at integration pts for symmetric tensorial field rank 2.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
std::map< std::string, std::vector< VectorDouble > > fieldMap
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
const std::string tagName
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
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)
MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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)
MoFEMErrorCode postProcess()
function is run at the end of loop
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
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode postProcess()
function is run at the end of loop