elem/Tetra4.cpp

#include "DAPTA.hpp"

namespace DAPTA { // Define namespace DAPTA

/*
Tetra4<V>::Tetra4(VertexHandle i, VertexHandle j, VertexHandle k, VertexHandle l)
Create a Tetra element from 4 Vertex definitions.
*/
template <typename V>
Tetra4<V>::Tetra4(VertexHandle i, VertexHandle j, VertexHandle k, VertexHandle l) {
   assert(i != NULL);
   assert(j != NULL);
   assert(k != NULL);
   assert(l != NULL);

   for(int x=0; x<4; x++)
      vertices.at(x) = NULL;

   vertices.at(0) = i;
   vertices.at(1) = j;
   vertices.at(2) = k;
   vertices.at(3) = l;
      
   for(int x=0; x<num_adjacents; x++) {
      adjacents.at(x) = NULL;
      coh_adjacents.at(x) = NULL;
   }
   
   globalID = 0;
}

/*
Tetra4<V>::Tetra4(boost::array<VertexHandle, num_vertices> nodes)
Create a Tetra element from a vector of Vertex definitions.
*/
template <typename V>
Tetra4<V>::Tetra4(boost::array<VertexHandle, num_vertices> nodes) {
   for(int i=0; i<num_vertices; i++) {
      assert(nodes.at(i) != NULL);
      vertices.at(i) = nodes.at(i);
   }
   
   for(int x=0; x<num_adjacents; x++) {
      adjacents.at(x) = NULL;
      coh_adjacents.at(x) = NULL;
   }
   
   globalID = 0;
}

/*
int Tetra4<V>::FindVertex(VertexHandle v)
Returns the index in the vertices vector where v is located.
*/
template <typename V>
int Tetra4<V>::FindVertex(VertexHandle v) {
   assert(v != NULL);
   
   int j=0;
   while(j<vertices.size() && vertices.at(j)!=v)
      j++;
   return j;
}

/*
bool Tetra4<V>::IsBoundaryEdge(VertexHandle v1, VertexHandle v2)
*/
template <typename V>
bool Tetra4<V>::IsBoundaryEdge(VertexHandle v1, VertexHandle v2) {
   return this->IsBoundaryEdge(FindVertex(v1), FindVertex(v2));
}

/*
bool Tetra4<V>::IsBoundaryEdge(int v1, int v2)
*/
template <typename V>
bool Tetra4<V>::IsBoundaryEdge(int v1, int v2) {
   return true;
/* assert(v1 < vertices.size());
   assert(v2 < vertices.size());
   
   ElementHandle cur_tet = this;
   
   boost::array<VertexHandle,2> vvs;
   vvs[0] = vertices.at(v1);
   vvs[1] = vertices.at(v2);
   
   EdgeElementType edge(vvs);
   
   do {
      v1 = cur_tet->FindVertex(vvs[0]);
      v2 = cur_tet->FindVertex(vvs[1]);
      int test1=-1, test2=-1;
      for(int i=0; i<4; i++) {
         if((i != v1) && (i != v2)) {
            (test1 == -1) ? test1 = i : test2 = i;
         }
      }
      if((cur_tet->adjacents.at(test1) == NULL) || (cur_tet->coh_adjacents.at(test1) != NULL) || 
         (cur_tet->adjacents.at(test2) == NULL) || (cur_tet->coh_adjacents.at(test2) != NULL) )
         // If we've found a boundary face then return true
         // Boundary is if there's a cohesive element there, or if there isn't any element there
         return true;
      cur_tet = cur_tet->NextOnEdge(cur_tet->GetEdgeIndex(edge));
      std::cout << this << ":" << cur_tet << std::endl;
   } while(cur_tet != this);
   
   // We didn't find a boundary face, so this isn't a boundary edge
   return false;
*/
   
/* assert(v1 < vertices.size());
   assert(v2 < vertices.size());

   // Start iterating round with this tetra
   ElementHandle cur_tet = this;

   // Choose a direction to go round in
   int next_face=0;
   while((next_face == v1) || (next_face == v2))
      next_face++;

   do {
      //std::cout << "1: " << v1 << ", 2: " << v2 << ", f: " << next_face << std::endl;
      //std::cout << cur_tet->adjacents.at(next_face) << ", " << cur_tet->coh_adjacents.at(next_face) << std::endl;
      if((cur_tet->adjacents.at(next_face) == NULL) || (cur_tet->coh_adjacents.at(next_face) != NULL))
         // If we've found a boundary face then return true
         // Boundary is if there's a cohesive element there, or if there isn't any element there
         return true;
      cur_tet = cur_tet->NextTetOnEdge(v1, v2, next_face);
   } while(cur_tet != this);
   
   // We didn't find a boundary face, so this isn't a boundary edge
   return false; */
}

/*
typename Tetra4<V>::ElementHandle Tetra4<V>::NextTetOnEdge(int &v1, int &v2, int &f)
Update v1, v2 and f indexes and return pointer to the next tetra along when iterating
around an edge defined by v1 and v2 going in the direction of f.
*/
template <typename V>
typename Tetra4<V>::ElementHandle Tetra4<V>::NextTetOnEdge(int &v1, int &v2, int &f) {
   ElementHandle new_tetra = adjacents.at(f);

   assert(new_tetra != NULL);

   int new_f=this->FaceToVertex(v1, v2, f);
   //std::cout << "v1: " << v1 << ", v2: " << v2 << ", f: " << f << std::endl;
   //std::cout << "new_f: " << new_f << std::endl;

   VertexHandle v1_vertex, v2_vertex, new_f_vertex;
   v1_vertex = this->vertices.at(v1);
   v2_vertex = this->vertices.at(v2);
   new_f_vertex = this->vertices.at(new_f);

   v1 = new_tetra->FindVertex(v1_vertex);
   v2 = new_tetra->FindVertex(v2_vertex);
   f = new_tetra->FindVertex(new_f_vertex);
   
   //std::cout << "v1: " << v1 << ", v2: " << v2 << ", f: " << f << std::endl;
   //std::cout << "old: " << this << ", new: " << new_tetra << std::endl;
   
   assert(v1 < new_tetra->vertices.size());
   assert(v2 < new_tetra->vertices.size());
   assert(f < new_tetra->vertices.size());

   return new_tetra;
}

/*
int Tetra4<V>::FaceToVertex(VertexHandle v1, VertexHandle v2, VertexHandle v3)
Return the index of the vertex defining the face defined by v1, v2 and v3
*/
template <typename V>
int Tetra4<V>::FaceToVertex(VertexHandle v1, VertexHandle v2, VertexHandle v3) { 
   return this->FaceToVertex(FindVertex(v1), FindVertex(v2), FindVertex(v3));
}

/*
int Tetra4<V>::FaceToVertex(int v1, int v2, int v3)
Return the index of the vertex defining the face defined by v1, v2 and v3
*/
template <typename V>
int Tetra4<V>::FaceToVertex(int v1, int v2, int v3) {
   for(int i=0; i<vertices.size(); i++)
      if((i != v1) && (i != v2) && (i != v3))
         return i;
}

/*
int Tetra4<V>::FaceToVertex(boost::array<VertexHandle, 3> vs)
Return the index of the vertex defining the face defined by v1, v2 and v3
*/
template <typename V>
int Tetra4<V>::FaceToVertex(boost::array<VertexHandle, 3> vs) {
   return this->FaceToVertex(vs.at(0), vs.at(1), vs.at(2));
}

/*
Tetra4<V>::FaceElementType Tetra4<V>::GetFace(int ind)
Return the edge data structure for the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::FaceElementType Tetra4<V>::GetFace(int ind) {
   int f[4][3] = { {1,2,3}, {0,3,2}, {0,1,3}, {0,2,1} };
   
   boost::array<VertexHandle, FaceElementType::num_vertices> vs;
   
   for(int i=0; i<3; i++) {
      vs.at(i) = vertices.at(f[ind][i]);
   }
   return FaceElementType(vs);
}

/*
int Tetra4<V>::GetFaceIndex(Tetra4<V>::FaceElementType face)
Return the index for the face as defined
*/
template <typename V>
int Tetra4<V>::GetFaceIndex(Tetra4<V>::FaceElementType face) {
   return this->FaceToVertex(face.vertices);
   
   for(int i=0; i<num_adjacents; i++) {
      if(face == this->GetFace(i)) {
         return i;
      }
   }
   return num_adjacents;
}

/*
Tetra4<V>::EdgeElementType Tetra4<V>::GetEdge(int ind)
Return the edge data structure for the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::EdgeElementType Tetra4<V>::GetEdge(int ind) {
   int e[6][2] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
   boost::array<VertexHandle, EdgeElementType::num_vertices> vs;
   vs.at(0) = vertices.at(e[ind][0]);
   vs.at(1) = vertices.at(e[ind][1]);
   return   EdgeElementType(vs);
}

/*
int Tetra4<V>::GetEdgeIndex(Tetra4<V>::EdgeElementType edge)
Return the index for the edge as defined
*/
template <typename V>
int Tetra4<V>::GetEdgeIndex(Tetra4<V>::EdgeElementType edge) {
   for(int i=0; i<num_edges; i++) {
      if(edge == this->GetEdge(i)) {
         return i;
      }
   }
   return num_edges;
}

/*
std::vector<EdgeElementType> Tetra4<V>::GetEdges(int f)
Return the edges in an array which correspond to the given face, f
*/
template <typename V>
std::vector<typename Tetra4<V>::EdgeElementType> Tetra4<V>::GetEdges(int f) {
   int e[4][3][2] = { { {2,3}, {3,1}, {1,2} }, 
                       { {3,2}, {2,0}, {0,3} }, 
                       { {1,3}, {3,0}, {0,1} }, 
                       { {0,2}, {2,1}, {1,0} } };
   
   std::vector<EdgeElementType> edges;

   boost::array<VertexHandle,EdgeElementType::num_vertices> vertices_1, vertices_2, vertices_3;
   
   vertices_1.at(0) = this->vertices.at(e[f][0][0]);
   vertices_1.at(1) = this->vertices.at(e[f][0][1]);
   vertices_2.at(0) = this->vertices.at(e[f][1][0]);
   vertices_2.at(1) = this->vertices.at(e[f][1][1]);
   vertices_3.at(0) = this->vertices.at(e[f][2][0]);
   vertices_3.at(1) = this->vertices.at(e[f][2][1]);
   
   edges.push_back(EdgeElementType(vertices_1));
   edges.push_back(EdgeElementType(vertices_2));
   edges.push_back(EdgeElementType(vertices_3));
   
   return edges;
}

/*
std::vector<EdgeElementType> Tetra4<V>::GetAllEdges()
*/
template <typename V>
std::vector<typename Tetra4<V>::EdgeElementType> Tetra4<V>::GetAllEdges() {
   int e[4][3][2] = { { {2,3}, {3,1}, {1,2} }, 
                       { {3,2}, {2,0}, {0,3} }, 
                       { {1,3}, {3,0}, {0,1} }, 
                       { {0,2}, {2,1}, {1,0} } };
   
   std::vector<EdgeElementType> edges;
   
   boost::array<VertexHandle,EdgeElementType::num_vertices> vertices_1, vertices_2, vertices_3, vertices_4, vertices_5, vertices_6;
   
   vertices_1.at(0) = this->vertices.at(e[0][0][0]);
   vertices_1.at(1) = this->vertices.at(e[0][0][1]);
   vertices_2.at(0) = this->vertices.at(e[0][1][0]);
   vertices_2.at(1) = this->vertices.at(e[0][1][1]);
   vertices_3.at(0) = this->vertices.at(e[0][2][0]);
   vertices_3.at(1) = this->vertices.at(e[0][2][1]);
   vertices_4.at(0) = this->vertices.at(e[1][1][0]);
   vertices_4.at(1) = this->vertices.at(e[1][1][1]);
   vertices_5.at(0) = this->vertices.at(e[1][2][0]);
   vertices_5.at(1) = this->vertices.at(e[1][2][1]);
   vertices_6.at(0) = this->vertices.at(e[2][2][0]);
   vertices_6.at(1) = this->vertices.at(e[2][2][1]);
   
   edges.push_back(EdgeElementType(vertices_1));
   edges.push_back(EdgeElementType(vertices_2));
   edges.push_back(EdgeElementType(vertices_3));
   edges.push_back(EdgeElementType(vertices_4));
   edges.push_back(EdgeElementType(vertices_5));
   edges.push_back(EdgeElementType(vertices_6));
   
   return edges;
}

/*
Tetra4<V>::ElementHandle Tetra4<V>::NextOnEdge(int ind)
Return the next element on the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::ElementHandle Tetra4<V>::NextOnEdge(int ind) {
   int e1[6] = { 2, 3, 1, 0, 2, 0 };
   int e2[6] = { 3, 1, 2, 3, 0, 1 };
   return adjacents.at(e1[ind]);
}

/*
Tetra4<V>::ElementHandle Tetra4<V>::PrevOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::ElementHandle Tetra4<V>::PrevOnEdge(int ind) {
   int e1[6] = { 2, 3, 1, 0, 2, 0 };
   int e2[6] = { 3, 1, 2, 3, 0, 1 };
   return adjacents.at(e2[ind]);
}

/*
Tetra4<V>::ElementHandle Tetra4<V>::NextOnEdge(EdgeElementType edge)
Return the next element on the edge
*/
template <typename V>
typename Tetra4<V>::ElementHandle Tetra4<V>::NextOnEdge(EdgeElementType edge) {
   int e1[6] = { 2, 3, 1, 0, 2, 0 };
   int e2[6] = { 3, 1, 2, 3, 0, 1 };
   
   for(int i=0; i<num_edges; i++) {
      if(edge == this->GetEdge(i)) {
         return adjacents.at(e1[i]);
      }
      else if(edge == !this->GetEdge(i)) {
         return adjacents.at(e2[i]);
      }
   }
   
   return NULL;
}

/*
Tetra4<V>::ElementHandle Tetra4<V>::PrevOnEdge(EdgeElementType edge)
Return the previous element on the edge
*/
template <typename V>
typename Tetra4<V>::ElementHandle Tetra4<V>::PrevOnEdge(EdgeElementType edge) {
   int e1[6] = { 2, 3, 1, 0, 2, 0 };
   int e2[6] = { 3, 1, 2, 3, 0, 1 };
   
   for(int i=0; i<num_edges; i++) {
      if(edge == this->GetEdge(i)) {
         return adjacents.at(e2[i]);
      }
      else if(edge == !this->GetEdge(i)) {
         return adjacents.at(e1[i]);
      }
   }
   
   return NULL;}

/*
Tetra4<V>::CohElementHandle Tetra4<V>::NextCohOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::CohElementHandle Tetra4<V>::NextCohOnEdge(int ind) {
   int e[6] = { 2, 3, 1, 0, 2, 0 };
   return coh_adjacents.at(e[ind]);
}

/*
Tetra4<V>::CohElementHandle Tetra4<V>::PrevCohOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tetra4<V>::CohElementHandle Tetra4<V>::PrevCohOnEdge(int ind) {
   int e[6] = { 3, 1, 2, 3, 0, 1 };
   return coh_adjacents.at(e[ind]);
}

/*
typename Tetra4<V>::VertexHandle Tetra4<V>::GetCentre()
Return a the centre of the element (vertex will have it's global ID set to that of the element)
*/
template <typename V>
typename Tetra4<V>::VertexHandle Tetra4<V>::GetCentre() {
   typename V::CoordType new_x=0.0, new_y=0.0, new_z=0.0;
   for(int i=0; i<4; i++) {
      new_x += vertices.at(i)->coords.at(0);
      new_y += vertices.at(i)->coords.at(1);
      new_z += vertices.at(i)->coords.at(2);
   }
   new_x /= 4; new_y /= 4; new_z /= 4;
   
   return new V(new_x, new_y, new_z);
}

} // namespace DAPTA

---