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        }
  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;
  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 
  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, 3 > VectorDouble3
VectorBoundedArray< double, 6 > VectorDouble6
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)