710 {
712 Tag th0, th1, th2, th3;
714 if (number_pathes) {
716 }
717 for (auto edge : edges) {
719 CHKERR moab.get_adjacencies(&edge, 1, 2,
false, adj_faces);
720 adj_faces = intersect(adj_faces, tris);
721 if (adj_faces.size() != 2) {
723 "Should be 2 faces adjacent to edge but is %ld",
724 adj_faces.size());
725 }
729 &
v[0], &
v[1], &
v[2]);
730 };
731
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));
747 int orientation;
748
749 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
750 if (orientation == -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) {
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]) {
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]);
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());
800 if (face_conn[
n] != edge_conn[
n])
801 break;
803 CHKERR moab.get_coords(&face_conn[
n], 1,
804 &*coords_face.data().begin());
806 CHKERR moab.get_coords(&edge_conn[0], 2,
807 &*coords_edge.data().begin());
809 coords_edge[0], coords_edge[1], coords_edge[2]);
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));
818 if (t_n1(
i) * t_face_n(
i) < 0)
820
822 };
823
824 auto get_base_for_not_planar_faces = [&]() {
827 t_n1(
i) /= sqrt(t_n1(
i) * t_n1(
i));
829 };
830
831
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
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
#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, 3 > VectorDouble3
VectorBoundedArray< double, 6 > VectorDouble6
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)