719 {
721 Tag th0, th1, th2, th3;
723 if (number_pathes) {
725 }
726 for (auto edge : edges) {
728 CHKERR moab.get_adjacencies(&edge, 1, 2,
false, adj_faces);
729 adj_faces = intersect(adj_faces, tris);
730 if (adj_faces.size() != 2) {
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
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));
756 int orientation;
757
758 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
759 if (orientation == -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) {
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]) {
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]);
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());
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());
815 CHKERR moab.get_coords(&edge_conn[0], 2,
816 &*coords_edge.data().begin());
818 coords_edge[0], coords_edge[1], coords_edge[2]);
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));
827 if (t_n1(
i) * t_face_n(
i) < 0)
829
831 };
832
833 auto get_base_for_not_planar_faces = [&]() {
836 t_n1(
i) /= sqrt(t_n1(
i) * t_n1(
i));
838 };
839
840
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
#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 > &)
VectorBoundedArray< double, 6 > VectorDouble6
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)