v0.15.5
Loading...
Searching...
No Matches
Static Public Member Functions | List of all members
EdgeSlidingConstrains::CalculateEdgeBase Struct Reference

#include "tools/src/SurfaceSlidingConstrains.hpp"

Static Public Member Functions

static MoFEMErrorCode createTag (moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
 
static MoFEMErrorCode numberSurfaces (moab::Interface &moab, Range edges, Range tris)
 
static MoFEMErrorCode setTags (moab::Interface &moab, Range edges, Range tris, bool number_pathes=true, boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > surface_orientation=nullptr, MoFEM::Interface *m_field_ptr=nullptr)
 
static MoFEMErrorCode saveEdges (moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
 
static MoFEMErrorCode createTag (moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
 
static MoFEMErrorCode numberSurfaces (moab::Interface &moab, Range edges, Range tris)
 
static MoFEMErrorCode setTags (moab::Interface &moab, Range edges, Range tris, bool number_pathes=true, boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > surface_orientation=nullptr, MoFEM::Interface *m_field_ptr=nullptr)
 
static MoFEMErrorCode saveEdges (moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
 

Detailed Description

Definition at line 632 of file SurfaceSlidingConstrains.hpp.

Member Function Documentation

◆ createTag() [1/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::createTag ( moab::Interface &  moab,
Tag th0,
Tag th1,
Tag th2,
Tag th3 
)
inlinestatic

Definition at line 634 of file SurfaceSlidingConstrains.hpp.

635 {
637 double def_val[] = {0, 0, 0};
638 CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
639 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
640 CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
641 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
642 int def_numb_val[] = {-1};
643 CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
644 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
645 int def_orientation_val[] = {1};
646 CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
647 MB_TAG_CREAT | MB_TAG_SPARSE,
648 def_orientation_val);
649
651 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.

◆ createTag() [2/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::createTag ( moab::Interface &  moab,
Tag th0,
Tag th1,
Tag th2,
Tag th3 
)
inlinestatic

Definition at line 638 of file SurfaceSlidingConstrains.hpp.

639 {
641 double def_val[] = {0, 0, 0};
642 CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
643 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
644 CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
645 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
646 int def_numb_val[] = {-1};
647 CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
648 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
649 int def_orientation_val[] = {1};
650 CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
651 MB_TAG_CREAT | MB_TAG_SPARSE,
652 def_orientation_val);
653
655 }

◆ numberSurfaces() [1/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::numberSurfaces ( moab::Interface &  moab,
Range  edges,
Range  tris 
)
inlinestatic

Definition at line 653 of file SurfaceSlidingConstrains.hpp.

654 {
656
657 auto get_edges = [&](const Range &ents) {
658 Range edges;
659 CHKERR moab.get_adjacencies(ents, 1, false, edges,
660 moab::Interface::UNION);
661 return edges;
662 };
663
664 auto get_face_adj = [&, get_edges](const Range &faces) {
665 Range adj_faces;
666 CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
667 adj_faces, moab::Interface::UNION);
668 return intersect(adj_faces, tris);
669 };
670
671 auto get_patch = [&, get_face_adj](const EntityHandle face) {
672 Range patch_ents;
673 patch_ents.insert(face);
674 unsigned int nb0;
675 do {
676 nb0 = patch_ents.size();
677 patch_ents.merge(get_face_adj(patch_ents));
678 } while (nb0 != patch_ents.size());
679 return patch_ents;
680 };
681
682 auto get_patches = [&]() {
683 std::vector<Range> patches;
684 while (!tris.empty()) {
685 patches.push_back(get_patch(tris[0]));
686 tris = subtract(tris, patches.back());
687 }
688 return patches;
689 };
690
691 Tag th0, th1, th2, th3;
692 CHKERR createTag(moab, th0, th1, th2, th3);
693
694 auto patches = get_patches();
695 int pp = 0;
696 for (auto patch : patches) {
697 std::vector<int> tags_vals(patch.size(), pp);
698 CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
699 ++pp;
700 }
701
703 }
static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)

◆ numberSurfaces() [2/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::numberSurfaces ( moab::Interface &  moab,
Range  edges,
Range  tris 
)
inlinestatic

Definition at line 657 of file SurfaceSlidingConstrains.hpp.

658 {
660
661 auto get_edges = [&](const Range &ents) {
662 Range edges;
663 CHKERR moab.get_adjacencies(ents, 1, false, edges,
664 moab::Interface::UNION);
665 return edges;
666 };
667
668 auto get_face_adj = [&, get_edges](const Range &faces) {
669 Range adj_faces;
670 CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
671 adj_faces, moab::Interface::UNION);
672 return intersect(adj_faces, tris);
673 };
674
675 auto get_patch = [&, get_face_adj](const EntityHandle face) {
676 Range patch_ents;
677 patch_ents.insert(face);
678 unsigned int nb0;
679 do {
680 nb0 = patch_ents.size();
681 patch_ents.merge(get_face_adj(patch_ents));
682 } while (nb0 != patch_ents.size());
683 return patch_ents;
684 };
685
686 auto get_patches = [&]() {
687 std::vector<Range> patches;
688 while (!tris.empty()) {
689 patches.push_back(get_patch(tris[0]));
690 tris = subtract(tris, patches.back());
691 }
692 return patches;
693 };
694
695 Tag th0, th1, th2, th3;
696 CHKERR createTag(moab, th0, th1, th2, th3);
697
698 auto patches = get_patches();
699 int pp = 0;
700 for (auto patch : patches) {
701 std::vector<int> tags_vals(patch.size(), pp);
702 CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
703 ++pp;
704 }
705
707 }

◆ saveEdges() [1/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::saveEdges ( moab::Interface &  moab,
std::string  name,
Range  edges,
Range faces = nullptr 
)
inlinestatic

Definition at line 846 of file SurfaceSlidingConstrains.hpp.

847 {
849 EntityHandle meshset;
850 Tag ths[4];
851 CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
852 CHKERR moab.create_meshset(MESHSET_SET, meshset);
853 CHKERR moab.add_entities(meshset, edges);
854 if (faces != nullptr) {
855 CHKERR moab.add_entities(meshset, *faces);
856 }
857 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
858 CHKERR moab.delete_entities(&meshset, 1);
860 }

◆ saveEdges() [2/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::saveEdges ( moab::Interface &  moab,
std::string  name,
Range  edges,
Range faces = nullptr 
)
inlinestatic

Definition at line 850 of file SurfaceSlidingConstrains.hpp.

851 {
853 EntityHandle meshset;
854 Tag ths[4];
855 CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
856 CHKERR moab.create_meshset(MESHSET_SET, meshset);
857 CHKERR moab.add_entities(meshset, edges);
858 if (faces != nullptr) {
859 CHKERR moab.add_entities(meshset, *faces);
860 }
861 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
862 CHKERR moab.delete_entities(&meshset, 1);
864 }

◆ setTags() [1/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::setTags ( moab::Interface &  moab,
Range  edges,
Range  tris,
bool  number_pathes = true,
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation surface_orientation = nullptr,
MoFEM::Interface m_field_ptr = nullptr 
)
inlinestatic

Definition at line 705 of file SurfaceSlidingConstrains.hpp.

710 {
712 Tag th0, th1, th2, th3;
713 CHKERR createTag(moab, th0, th1, th2, th3);
714 if (number_pathes) {
715 CHKERR numberSurfaces(moab, edges, tris);
716 }
717 for (auto edge : edges) {
718 Range adj_faces;
719 CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
720 adj_faces = intersect(adj_faces, tris);
721 if (adj_faces.size() != 2) {
722 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
723 "Should be 2 faces adjacent to edge but is %ld",
724 adj_faces.size());
725 }
727 auto get_tensor_from_vec = [](VectorDouble3 &v) {
729 &v[0], &v[1], &v[2]);
730 };
731
732 FTensor::Index<'i', 3> i;
733 FTensor::Index<'j', 3> j;
734 FTensor::Index<'k', 3> k;
735
736 std::array<double, 3> areas;
737 auto calculate_normals = [&, get_tensor_from_vec]() {
738 int ff = 0;
739 for (auto face : adj_faces) {
740 double &x = (v[ff])[0];
741 double &y = (v[ff])[1];
742 double &z = (v[ff])[2];
743 moab::Util::normal(&moab, face, x, y, z);
744 auto t_n = get_tensor_from_vec(v[ff]);
745 areas[ff] = sqrt(t_n(i) * t_n(i));
746 t_n(i) /= areas[ff];
747 int orientation;
748 // FIXME: this tag is empty
749 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
750 if (orientation == -1) {
751 t_n(i) *= -1;
752 }
753 if (surface_orientation && m_field_ptr) {
754 CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
755 face);
756 int eo = surface_orientation->elementOrientation;
757 if (eo == -1) {
758 t_n(i) *= -1;
759 }
760 }
761 ++ff;
762 }
763 };
764 calculate_normals();
765
766 auto get_patch_number = [&]() {
767 std::vector<int> p = {0, 0};
768 CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
769 return p;
770 };
771
772 std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
773 adj_faces[1]};
774 auto order_normals = [&, get_patch_number]() {
775 auto p = get_patch_number();
776 if (p[0] < p[1]) {
777 v[0].swap(v[1]);
778 faces_handles[0] = adj_faces[1];
779 faces_handles[1] = adj_faces[0];
780 }
781 };
782 order_normals();
783
785 auto t_n0 = get_tensor_from_vec(v[0]);
786 auto t_n1 = get_tensor_from_vec(v[1]);
787 t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
788
789 auto get_base_for_coplanar_faces = [&]() {
791
792 std::vector<EntityHandle> face_conn;
793 CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
794 std::vector<EntityHandle> edge_conn;
795 CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
796 std::sort(face_conn.begin(), face_conn.end());
797 std::sort(edge_conn.begin(), edge_conn.end());
798 int n = 0;
799 for (; n != 2; ++n)
800 if (face_conn[n] != edge_conn[n])
801 break;
802 VectorDouble3 coords_face(3);
803 CHKERR moab.get_coords(&face_conn[n], 1,
804 &*coords_face.data().begin());
805 VectorDouble6 coords_edge(6);
806 CHKERR moab.get_coords(&edge_conn[0], 2,
807 &*coords_edge.data().begin());
808 auto t_edge_n0 = FTensor::Tensor1<double, 3>(
809 coords_edge[0], coords_edge[1], coords_edge[2]);
810 auto t_edge_n1 = FTensor::Tensor1<double, 3>(
811 coords_edge[3], coords_edge[4], coords_edge[5]);
812 auto t_face_n = get_tensor_from_vec(coords_face);
813 t_face_n(i) -= t_edge_n0(i);
815 t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
816 t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
817 t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
818 if (t_n1(i) * t_face_n(i) < 0)
819 t_n1(i) *= -1;
820
822 };
823
824 auto get_base_for_not_planar_faces = [&]() {
826 t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
827 t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
829 };
830
831 // This a case when faces adjacent to edge are coplanar !!!
832 constexpr double tol = 1e-6;
833 if ((t_cross(i) * t_cross(i)) < tol)
834 CHKERR get_base_for_coplanar_faces();
835 else
836 CHKERR get_base_for_not_planar_faces();
837
838 VectorDouble3 &v0 = v[0];
839 VectorDouble3 &v1 = v[1];
840 CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
841 CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
842 }
844 }
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto cross(const Tensor1_Expr< A, T, 3, i > &a, const Tensor1_Expr< B, U, 3, j > &b, const Index< k, 3 > &)
Definition cross.hpp:10
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
VectorBoundedArray< double, 6 > VectorDouble6
Definition Types.hpp:95
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)
double tol

◆ setTags() [2/2]

static MoFEMErrorCode EdgeSlidingConstrains::CalculateEdgeBase::setTags ( moab::Interface &  moab,
Range  edges,
Range  tris,
bool  number_pathes = true,
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation surface_orientation = nullptr,
MoFEM::Interface m_field_ptr = nullptr 
)
inlinestatic

Definition at line 709 of file SurfaceSlidingConstrains.hpp.

714 {
716 Tag th0, th1, th2, th3;
717 CHKERR createTag(moab, th0, th1, th2, th3);
718 if (number_pathes) {
719 CHKERR numberSurfaces(moab, edges, tris);
720 }
721 for (auto edge : edges) {
722 Range adj_faces;
723 CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
724 adj_faces = intersect(adj_faces, tris);
725 if (adj_faces.size() != 2) {
726 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
727 "Should be 2 faces adjacent to edge but is %ld",
728 adj_faces.size());
729 }
731 auto get_tensor_from_vec = [](VectorDouble3 &v) {
733 &v[0], &v[1], &v[2]);
734 };
735
736 FTensor::Index<'i', 3> i;
737 FTensor::Index<'j', 3> j;
738 FTensor::Index<'k', 3> k;
739
740 std::array<double, 3> areas;
741 auto calculate_normals = [&, get_tensor_from_vec]() {
742 int ff = 0;
743 for (auto face : adj_faces) {
744 double &x = (v[ff])[0];
745 double &y = (v[ff])[1];
746 double &z = (v[ff])[2];
747 moab::Util::normal(&moab, face, x, y, z);
748 auto t_n = get_tensor_from_vec(v[ff]);
749 areas[ff] = sqrt(t_n(i) * t_n(i));
750 t_n(i) /= areas[ff];
751 int orientation;
752 // FIXME: this tag is empty
753 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
754 if (orientation == -1) {
755 t_n(i) *= -1;
756 }
757 if (surface_orientation && m_field_ptr) {
758 CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
759 face);
760 int eo = surface_orientation->elementOrientation;
761 if (eo == -1) {
762 t_n(i) *= -1;
763 }
764 }
765 ++ff;
766 }
767 };
768 calculate_normals();
769
770 auto get_patch_number = [&]() {
771 std::vector<int> p = {0, 0};
772 CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
773 return p;
774 };
775
776 std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
777 adj_faces[1]};
778 auto order_normals = [&, get_patch_number]() {
779 auto p = get_patch_number();
780 if (p[0] < p[1]) {
781 v[0].swap(v[1]);
782 faces_handles[0] = adj_faces[1];
783 faces_handles[1] = adj_faces[0];
784 }
785 };
786 order_normals();
787
789 auto t_n0 = get_tensor_from_vec(v[0]);
790 auto t_n1 = get_tensor_from_vec(v[1]);
791 t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
792
793 auto get_base_for_coplanar_faces = [&]() {
795
796 std::vector<EntityHandle> face_conn;
797 CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
798 std::vector<EntityHandle> edge_conn;
799 CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
800 std::sort(face_conn.begin(), face_conn.end());
801 std::sort(edge_conn.begin(), edge_conn.end());
802 int n = 0;
803 for (; n != 2; ++n)
804 if (face_conn[n] != edge_conn[n])
805 break;
806 VectorDouble3 coords_face(3);
807 CHKERR moab.get_coords(&face_conn[n], 1,
808 &*coords_face.data().begin());
809 VectorDouble6 coords_edge(6);
810 CHKERR moab.get_coords(&edge_conn[0], 2,
811 &*coords_edge.data().begin());
812 auto t_edge_n0 = FTensor::Tensor1<double, 3>(
813 coords_edge[0], coords_edge[1], coords_edge[2]);
814 auto t_edge_n1 = FTensor::Tensor1<double, 3>(
815 coords_edge[3], coords_edge[4], coords_edge[5]);
816 auto t_face_n = get_tensor_from_vec(coords_face);
817 t_face_n(i) -= t_edge_n0(i);
819 t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
820 t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
821 t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
822 if (t_n1(i) * t_face_n(i) < 0)
823 t_n1(i) *= -1;
824
826 };
827
828 auto get_base_for_not_planar_faces = [&]() {
830 t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
831 t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
833 };
834
835 // This a case when faces adjacent to edge are coplanar !!!
836 constexpr double tol = 1e-6;
837 if ((t_cross(i) * t_cross(i)) < tol)
838 CHKERR get_base_for_coplanar_faces();
839 else
840 CHKERR get_base_for_not_planar_faces();
841
842 VectorDouble3 &v0 = v[0];
843 VectorDouble3 &v1 = v[1];
844 CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
845 CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
846 }
848 }

The documentation for this struct was generated from the following files: