v0.15.0
Loading...
Searching...
No Matches
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 638 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 641 of file SurfaceSlidingConstrains.hpp.

642 {
644 double def_val[] = {0, 0, 0};
645 CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
646 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
647 CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
648 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
649 int def_numb_val[] = {-1};
650 CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
651 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
652 int def_orientation_val[] = {1};
653 CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
654 MB_TAG_CREAT | MB_TAG_SPARSE,
655 def_orientation_val);
656
658 }
#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 640 of file SurfaceSlidingConstrains.hpp.

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

◆ numberSurfaces() [1/2]

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

Definition at line 660 of file SurfaceSlidingConstrains.hpp.

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

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

◆ saveEdges() [1/2]

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

Definition at line 855 of file SurfaceSlidingConstrains.hpp.

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

◆ saveEdges() [2/2]

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

Definition at line 854 of file SurfaceSlidingConstrains.hpp.

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

◆ 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 714 of file SurfaceSlidingConstrains.hpp.

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

◆ 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 713 of file SurfaceSlidingConstrains.hpp.

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

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