elem/Tri3.cpp

#include "DAPTA.hpp"

namespace DAPTA { // Define namespace DAPTA

/*
Tri3<V>::Tri3(VertexHandle i, VertexHandle j, VertexHandle k)
*/
template <typename V>
Tri3<V>::Tri3(VertexHandle i, VertexHandle j, VertexHandle k) {
   assert(i != NULL);
   assert(j != NULL);
   assert(k != NULL);

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

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

/*
Tri3<V>::Tri3(boost::array<VertexHandle, num_vertices> nodes)
*/
template <typename V>
Tri3<V>::Tri3(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 Tri3<V>::FindVertex(VertexHandle v)
Returns the index in the vertices vector where v is located.
*/
template <typename V>
int Tri3<V>::FindVertex(VertexHandle v) {
   assert(v != NULL);
   
   int j=0;
   while(j<vertices.size() && vertices.at(j)!=v)
      j++;
   return j;
}

/*
bool Tri3<V>::IsBoundaryNode(VertexHandle v1)
*/
template <typename V>
bool Tri3<V>::IsBoundaryNode(VertexHandle v1) {
   return this->IsBoundaryNode(FindVertex(v1));
}

/*
bool Tri3<V>::IsBoundaryNode(int v1)
*/
template <typename V>
bool Tri3<V>::IsBoundaryNode(int v1) {
   assert(v1 < vertices.size());

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

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

   do {
      if((cur_tri->adjacents.at(next_face) == NULL) || (cur_tri->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_tri = cur_tri->NextTriOnNode(v1, next_face);
   } while(cur_tri != this);
   
   // We didn't find a boundary face, so this isn't a boundary edge
   return false;
}

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

/*
bool Tri3<V>::IsBoundaryEdge(int v1, int v2)
*/
template <typename V>
bool Tri3<V>::IsBoundaryEdge(int v1, int v2) {
   if((this->adjacents.at(this->FaceToVertex(v1, v2)) == NULL) || (this->coh_adjacents.at(this->FaceToVertex(v1, v2)) != NULL))
      return true;
   else
      return false;
}

/*
typename Tri3<V>::ElementHandle Tri3<V>::NextTriOnNode(int &v1, int &f)
Update v1 and f indexes and return pointer to the next tri along when iterating
around a node defined by v1 going in the direction of f.
*/
template <typename V>
typename Tri3<V>::ElementHandle Tri3<V>::NextTriOnNode(int &v1, int &f) {
   ElementHandle new_tri = adjacents.at(f);

   assert(new_tri != NULL);
   
   int new_f=this->FaceToVertex(v1, f);

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

   v1 = new_tri->FindVertex(v1_vertex);
   f = new_tri->FindVertex(new_f_vertex);
   
   assert(v1 < new_tri->vertices.size());
   assert(f < new_tri->vertices.size());
   
   return new_tri;
}

/*
int Tri3<V>::FaceToVertex(VertexHandle v1, VertexHandle v2)
*/
template <typename V>
int Tri3<V>::FaceToVertex(VertexHandle v1, VertexHandle v2) {
   return this->FaceToVertex(FindVertex(v1), FindVertex(v2));
}

/*
int Tri3<V>::FaceToVertex(int v1, int v2)
*/
template <typename V>
int Tri3<V>::FaceToVertex(int v1, int v2) {
   int ret = 0;
   while((v1 == ret) || (v2 == ret))
      ret++;
   return ret;
}

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

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

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

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

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

/*
std::vector<EdgeElementType> Tri3<V>::GetEdges(int f)
Return the edges in an array which correspond to the given face, f
*/
template <typename V>
std::vector<typename Tri3<V>::EdgeElementType> Tri3<V>::GetEdges(int f) {
   std::vector<EdgeElementType> edges;

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

/*
std::vector<EdgeElementType> Tri3<V>::GetAllEdges()
*/
template <typename V>
std::vector<typename Tri3<V>::EdgeElementType> Tri3<V>::GetAllEdges() {
   std::vector<EdgeElementType> edges;
   
   boost::array<VertexHandle,EdgeElementType::num_vertices> vertices_1, vertices_2, vertices_3;
   
   vertices_1.at(0) = this->vertices.at(0);
   vertices_2.at(0) = this->vertices.at(1);
   vertices_3.at(0) = this->vertices.at(2);
   edges.push_back(EdgeElementType(vertices_1));
   edges.push_back(EdgeElementType(vertices_2));
   edges.push_back(EdgeElementType(vertices_3));
   
   return edges;
}

/*
Tri3<V>::ElementHandle Tri3<V>::NextOnEdge(int ind)
Return the next element on the edge defined by index, ind
*/
template <typename V>
typename Tri3<V>::ElementHandle Tri3<V>::NextOnEdge(int ind) {
   return adjacents.at((ind+1)%3);
}

/*
Tri3<V>::ElementHandle Tri3<V>::PrevOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tri3<V>::ElementHandle Tri3<V>::PrevOnEdge(int ind) {
   return adjacents.at((ind+2)%3);
}

/*
Tri3<V>::CohElementHandle Tri3<V>::NextCohOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tri3<V>::CohElementHandle Tri3<V>::NextCohOnEdge(int ind) {
   return coh_adjacents.at((ind+1)%3);
}

/*
Tri3<V>::CohElementHandle Tri3<V>::PrevCohOnEdge(int ind)
Return the previous element on the edge defined by index, ind
*/
template <typename V>
typename Tri3<V>::CohElementHandle Tri3<V>::PrevCohOnEdge(int ind) {
   return coh_adjacents.at((ind+2)%3);
}

/*
typename Tri3<V>::VertexHandle Tri3<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 Tri3<V>::VertexHandle Tri3<V>::GetCentre() {
   typename V::CoordType new_x=0.0, new_y=0.0, new_z=0.0;
   for(int i=0; i<3; 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 /= 3; new_y /= 3; new_z /= 3;
   
   return new V(new_x, new_y, new_z);
}

} // namespace DAPTA

---