v0.14.0
Static Public Member Functions | List of all members
EdgeSlidingConstrains::CalculateEdgeBase Struct Reference

#include <users_modules/basic_finite_elements/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)
 

Detailed Description

Definition at line 641 of file SurfaceSlidingConstrains.hpp.

Member Function Documentation

◆ createTag()

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

Definition at line 643 of file SurfaceSlidingConstrains.hpp.

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

◆ numberSurfaces()

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

Definition at line 662 of file SurfaceSlidingConstrains.hpp.

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

◆ saveEdges()

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

Definition at line 857 of file SurfaceSlidingConstrains.hpp.

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

◆ setTags()

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

721  {
722  MoFEMFunctionBegin;
723  Tag th0, th1, th2, th3;
724  CHKERR createTag(moab, th0, th1, th2, th3);
725  if (number_pathes) {
726  CHKERR numberSurfaces(moab, edges, tris);
727  }
728  for (auto edge : edges) {
729  Range adj_faces;
730  CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
731  adj_faces = intersect(adj_faces, tris);
732  if (adj_faces.size() != 2) {
733  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
734  "Should be 2 faces adjacent to edge but is %d",
735  adj_faces.size());
736  }
737  VectorDouble3 v[2] = {VectorDouble3(3), VectorDouble3(3)};
738  auto get_tensor_from_vec = [](VectorDouble3 &v) {
739  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
740  &v[0], &v[1], &v[2]);
741  };
742 
743  FTensor::Index<'i', 3> i;
744  FTensor::Index<'j', 3> j;
745  FTensor::Index<'k', 3> k;
746 
747  std::array<double, 3> areas;
748  auto calculate_normals = [&, get_tensor_from_vec]() {
749  int ff = 0;
750  for (auto face : adj_faces) {
751  double &x = (v[ff])[0];
752  double &y = (v[ff])[1];
753  double &z = (v[ff])[2];
754  moab::Util::normal(&moab, face, x, y, z);
755  auto t_n = get_tensor_from_vec(v[ff]);
756  areas[ff] = sqrt(t_n(i) * t_n(i));
757  t_n(i) /= areas[ff];
758  int orientation;
759  // FIXME: this tag is empty
760  CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
761  if (orientation == -1) {
762  t_n(i) *= -1;
763  }
764  if (surface_orientation && m_field_ptr) {
765  CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
766  face);
767  int eo = surface_orientation->elementOrientation;
768  if (eo == -1) {
769  t_n(i) *= -1;
770  }
771  }
772  ++ff;
773  }
774  };
775  calculate_normals();
776 
777  auto get_patch_number = [&]() {
778  std::vector<int> p = {0, 0};
779  CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
780  return p;
781  };
782 
783  std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
784  adj_faces[1]};
785  auto order_normals = [&, get_patch_number]() {
786  auto p = get_patch_number();
787  if (p[0] < p[1]) {
788  v[0].swap(v[1]);
789  faces_handles[0] = adj_faces[1];
790  faces_handles[1] = adj_faces[0];
791  }
792  };
793  order_normals();
794 
795  FTensor::Tensor1<double, 3> t_cross;
796  auto t_n0 = get_tensor_from_vec(v[0]);
797  auto t_n1 = get_tensor_from_vec(v[1]);
798  t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
799 
800  auto get_base_for_coplanar_faces = [&]() {
801  MoFEMFunctionBegin;
802 
803  std::vector<EntityHandle> face_conn;
804  CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
805  std::vector<EntityHandle> edge_conn;
806  CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
807  std::sort(face_conn.begin(), face_conn.end());
808  std::sort(edge_conn.begin(), edge_conn.end());
809  int n = 0;
810  for (; n != 2; ++n)
811  if (face_conn[n] != edge_conn[n])
812  break;
813  VectorDouble3 coords_face(3);
814  CHKERR moab.get_coords(&face_conn[n], 1,
815  &*coords_face.data().begin());
816  VectorDouble6 coords_edge(6);
817  CHKERR moab.get_coords(&edge_conn[0], 2,
818  &*coords_edge.data().begin());
819  auto t_edge_n0 = FTensor::Tensor1<double, 3>(
820  coords_edge[0], coords_edge[1], coords_edge[2]);
821  auto t_edge_n1 = FTensor::Tensor1<double, 3>(
822  coords_edge[3], coords_edge[4], coords_edge[5]);
823  auto t_face_n = get_tensor_from_vec(coords_face);
824  t_face_n(i) -= t_edge_n0(i);
825  FTensor::Tensor1<double, 3> t_delta;
826  t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
827  t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
828  t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
829  if (t_n1(i) * t_face_n(i) < 0)
830  t_n1(i) *= -1;
831 
832  MoFEMFunctionReturn(0);
833  };
834 
835  auto get_base_for_not_planar_faces = [&]() {
836  MoFEMFunctionBeginHot;
837  t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
838  t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
839  MoFEMFunctionReturnHot(0);
840  };
841 
842  // This a case when faces adjacent to edge are coplanar !!!
843  constexpr double tol = 1e-6;
844  if ((t_cross(i) * t_cross(i)) < tol)
845  CHKERR get_base_for_coplanar_faces();
846  else
847  CHKERR get_base_for_not_planar_faces();
848 
849  VectorDouble3 &v0 = v[0];
850  VectorDouble3 &v1 = v[1];
851  CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
852  CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
853  }
854  MoFEMFunctionReturn(0);
855  }
The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::Types::VectorDouble6
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
FTensor::cross
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EdgeSlidingConstrains::CalculateEdgeBase::createTag
static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
Definition: SurfaceSlidingConstrains.hpp:643
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
EdgeSlidingConstrains::CalculateEdgeBase::numberSurfaces
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)
Definition: SurfaceSlidingConstrains.hpp:662
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
tol
double tol
Definition: mesh_smoothing.cpp:27
Generated by Doxygen 1.8.17 and hosted at University of Glasgow