mesh/Mesh.cpp

#include "DAPTA.hpp"

namespace DAPTA { // Define namespace DAPTA

/*
void Mesh<Traits>::FractureFace(ElementHandle e, int f)
Fracture the f'th face on e and updates the vertices and cohesive elements lists as appropriate.
*/
template <typename Traits>
void Mesh<Traits>::FractureFace(ElementHandle e, int f) {
   // Setup some variables for later
   int f_2=ElementType::num_adjacents;
   ElementHandle e_2 = e->adjacents.at(f);
   
   // If no adjacent then we can't fracture, so return NULL
   if(e_2 == NULL)
      return;
   
   // If there's already a cohesive element here then return
   if(e->coh_adjacents.at(f) != NULL)
      return;
   
   FaceElementType face_1 = e->GetFace(f);
   FaceElementType face_2 = face_1;
   bool found_face_2 = false;
   
   for(int i=0; i<ElementType::num_adjacents; i++) {
      face_2 = e_2->GetFace(i);
//    if((face_2 == face_1) || (!face_2 == face_1)) {
      if((!face_2 == face_1)) {
         found_face_2 = true;
         f_2 = i;
         break;
      }
   }
   assert(found_face_2);
   
/*    if(e->globalID == 255) {
         std::cout << "   FACE = " << f_2 << "(" << face_2.vertices[0]->globalID << ":" << face_2.vertices[1]->globalID << ")" << std::endl;
      }
   while(!(face_1 == face_2) && (f_2 != ElementType::num_adjacents)) {
      f_2++;
      face_2 = e_2->GetFace(f_2);
      if(e->globalID == 255) {
         std::cout << "   " << f_2 << std::endl;
      }
      if(e->globalID == 255) {
         std::cout << "   FACE = " << f_2 << "(" << face_2.vertices[0]->globalID << ":" << face_2.vertices[1]->globalID << ")" << std::endl;
      }
   }
   if(e->globalID == 255) {
      std::cout << "   " << f_2 << std::endl;
   }
   assert(f_2 != ElementType::num_adjacents); // Stops the whole thing looping in the while loop above */
   
   // Get list of edges for this face
   std::vector<EdgeElementType> frac_edges = e->GetEdges(f);
   
/* for(typename std::vector<EdgeElementType>::iterator iter=frac_edges.begin(); iter!=frac_edges.end(); ++iter) {
      std::sort(iter->vertices.begin(), iter->vertices.end(), IDComparison<VertexType>());
   }*/
   
   // Loop through each edge
   for(int i=0; i<frac_edges.size(); i++) {
      std::sort(frac_edges.at(i).vertices.begin(), frac_edges.at(i).vertices.end(), IDComparison<VertexType>());
      edges[frac_edges.at(i)]->break_chain(e->globalID, e_2->globalID);
   }
   
   CohElementHandle coh_elem = new CohElementType(e, f, e_2, f_2);
   e->coh_adjacents.at(f) = coh_elem;
   e_2->coh_adjacents.at(f_2) = coh_elem;
   
   AddCohElement(coh_elem);
   
   /*return coh_elem;*/
}

template <typename Traits>
void Mesh<Traits>::UpdateEdges() {
   for(EdgeMapIterator iter=EdgeMapBegin(); iter!=EdgeMapEnd(); ++iter) {
      // Check if this edge needs splitting and new nodes created
      if(iter->second->needSplit()) {
         std::vector<typename ElementConListType::iterator> del_iters;
         for(typename ElementConListType::iterator edge_iter=iter->second->begin(); edge_iter!=iter->second->end(); ++edge_iter) {
            // Check if broken or not linked
            if(edge_iter->link_state != ElementConListType::con) {
               typename ElementConListType::iterator save=edge_iter; save++;
               del_iters.push_back(save);
               edge_iter->link_state = ElementConListType::discon;
            }
         }
         
         boost::array<VertexHandle, EdgeElementType::num_vertices> old_vs;
         for(int j=0; j<EdgeElementType::num_vertices; j++) {
            old_vs.at(j) = iter->first.vertices.at(j);
         }
         
         // Now we have a split list of where the new nodes need to be - so now update them!
         for(int i=0; i<del_iters.size()-1; i++) {
            typename ElementConListType::iterator start = del_iters.at(i), end = del_iters.at(i+1);
            
            // Create the new element connectivity list
            ElementConListHandle new_con_list = new ElementConListType();
            
            // Create the new vertices and edge data structure
            boost::array<VertexHandle, EdgeElementType::num_vertices> new_vs;
            for(int j=0; j<EdgeElementType::num_vertices; j++) {
               new_vs.at(j) = new VertexType(old_vs.at(j));
               AddVertex(new_vs.at(j));
            }
            EdgeElementType edge_2(new_vs);
            edges[edge_2] = new_con_list;
            
            for(typename ElementConListType::iterator edge_iter=start; edge_iter!=end; ++edge_iter) {
               // Update the element
               /* For parallel, the element might not be on this processor, so we don't need to update it, 
                * so we try to find the element, but if it's NULL then we know it's not on this processor.
                */
               ElementHandle elem = FindElementFromIndex((*edge_iter).elem_globalID);
               if(elem != NULL) {
                  for(int j=0; j<EdgeElementType::num_vertices; j++) {
                     for(int k=0; k<ElementType::num_vertices; k++) {
                        if(elem->vertices.at(k) == old_vs.at(j)) {
                           elem->vertices.at(k) = new_vs.at(j);
                        }
                     }
                  }
               }
            }
            
            new_con_list->insert(new_con_list->begin(), start, end);
            iter->second->erase(start, end);
         }
      }
   }
}

/*
bool Mesh<Traits>::ReadFile (std::istream in)
Setup the mesh by reading from an input file
*/
template <typename Traits>
bool Mesh<Traits>::ReadFile (std::istream &in) {
   std::string elem_keyword = ElementType::deck_keyword;
   std::string node_keyword = VertexType::deck_keyword;
   
   std::string cur_keyword;
   bool startdata=false;
   bool enddata=false;
   
   std::map<int, VertexHandle> temp_vertices;

   while(!in.eof()) {
      in >> cur_keyword;
      if(startdata && !enddata) {
         if(cur_keyword == node_keyword) {
            int node_num;
            in >> node_num;
         
            double x=0.0, y=0.0, z=0.0;
            in >> x >> y >> z;
         
            VertexHandle v = new VertexType(x,y,z,node_num);
            this->AddVertex(v);
            
            // Note - not checking to see if node is already defined!
            temp_vertices[node_num] = v;
         }
         else if(cur_keyword == elem_keyword) {
            int elem_num;
            in >> elem_num;
         
            boost::array<VertexHandle, ElementType::num_vertices> nodes;

            int node_num=0;
            for(int i=0; i<ElementType::num_vertices; i++) {
               in >> node_num; nodes.at(i) = temp_vertices[node_num];
            }
            
            ElementHandle e = new ElementType(nodes);
            e->globalID = elem_num;
            this->AddElement(e);
         }
         else {
            //return false;
         }
      }
      if(cur_keyword == "STARTDATA") {
         startdata = true;
      }
      if(cur_keyword == "ENDDATA") {
         enddata = true;
      }
   }
   
   // Connect all the elements
   std::for_each(ElementMapBegin(), ElementMapEnd(), bind1st(std::mem_fun(&Mesh<Traits>::ConnectElement),this));
   
   // Form the element connectivity lists
   std::for_each(ElementMapBegin(), ElementMapEnd(), bind1st(std::mem_fun(&Mesh<Traits>::AddElementToConLists),this));
   
   return true;
}

/*
std::string Mesh<Traits>::ExportFile()
Export the mesh to a file
*/
template <typename Traits>
std::string Mesh<Traits>::ExportFile() {
   std::stringstream export_file;
   
   export_file << "STARTDATA\n";
   
   for(VertexMapIterator v_iter=VertexMapBegin(); v_iter!=VertexMapEnd(); ++v_iter) {
      export_file << VertexType::deck_keyword << "\t";
      export_file << " " << (*v_iter)->globalID;
      for(int i=0; i<VertexType::dimension; i++) {
         export_file << "\t" << (*v_iter)->coords.at(i);
      }
      export_file << "\n";
   }
   for(ElementMapIterator e_iter=ElementMapBegin(); e_iter!=ElementMapEnd(); ++e_iter) {
      export_file << ElementType::deck_keyword << "\t";
      export_file << " " << (*e_iter)->globalID;
      for(int i=0; i<ElementType::num_vertices; i++) {
         export_file << "\t" << (*e_iter)->vertices.at(i)->globalID;
      }
      export_file << "\n";
   }
   
   export_file << "ENDDATA\n";
   
   return export_file.str();
}

/*
std::string Mesh<Traits>::ExportDomainFile (std::set<Mesh<Traits>::VertexHandle> part_vert, std::set<Mesh<Traits>::ElementHandle> part_elem)
Export part of the mesh to a file
*/
template <typename Traits>
std::string Mesh<Traits>::ExportDomainFile (std::set<Mesh<Traits>::VertexHandle> part_vert, std::set<Mesh<Traits>::ElementHandle> part_elem) {
   std::stringstream export_file;
   
   export_file << "STARTDATA\n";
   
   for(typename std::set<VertexHandle>::iterator v_iter=part_vert.begin(); v_iter!=part_vert.end(); ++v_iter) {
      export_file << VertexType::deck_keyword << "\t";
      export_file << " " << (*v_iter)->globalID;
      for(int i=0; i<VertexType::dimension; i++) {
         export_file << "\t" << (*v_iter)->coords.at(i);
      }
      export_file << "\n";
   }
   for(typename std::set<ElementHandle>::iterator e_iter=part_elem.begin(); e_iter!=part_elem.end(); ++e_iter) {
      export_file << ElementType::deck_keyword << "\t";
      export_file << " " << (*e_iter)->globalID;
      for(int i=0; i<ElementType::num_vertices; i++) {
         export_file << "\t" << (*e_iter)->vertices.at(i)->globalID;
      }
      export_file << "\n";
   }
   
   export_file << "ENDDATA\n";
   
   return export_file.str();
}

/*
bool Mesh<Traits>::AddVertex(VertexHandle v)
Add a vertex to the Mesh.
*/
template <typename Traits>
bool Mesh<Traits>::AddVertex(VertexHandle v) {
   vertices.push_back(v);
   if(v->globalID == 0)
      v->globalID = this->maxVertex+1;
   this->maxVertex = std::max(v->globalID, this->maxVertex);
   std::stringstream s; s << newPrefix << v->globalID;
   int new_id; s >> new_id;
   v->globalID = new_id;
   return true;
}

/*
bool Mesh<Traits>::AddElement(ElementHandle e)
Add an element to the Mesh.
*/
template <typename Traits>
bool Mesh<Traits>::AddElement(ElementHandle e) {
   //ConnectElement(e);
   elements.push_back(e);
   if(e->globalID == 0)
      e->globalID = this->maxElement+1;
   this->maxElement = std::max(e->globalID, this->maxElement);
   return true;
}

/*
bool Mesh<Traits>::AddCohElement(CohElementHandle e)
Add a cohesive element to the Mesh.
*/
template <typename Traits>
bool Mesh<Traits>::AddCohElement(CohElementHandle e) {
   cohelements.push_back(e);
   if(e->globalID == 0)
      e->globalID = this->maxCohElement+1;
   this->maxCohElement = std::max(e->globalID, this->maxCohElement);
   return true;
}

/*
bool Mesh<Traits>::ClearMesh()
Clear all vertices and elements from the mesh
*/
template <typename Traits>
bool Mesh<Traits>::ClearMesh() {
   std::for_each(VertexMapBegin(),     VertexMapEnd(),     zap<VertexHandle>);
   std::for_each(ElementMapBegin(),    ElementMapEnd(),    zap<ElementHandle>);
   std::for_each(CohElementMapBegin(), CohElementMapEnd(), zap<CohElementHandle>);
   std::for_each(EdgeMapBegin(),       EdgeMapEnd(),       TakeSecond<typename EdgeMapType::value_type>(zap<ElementConListHandle>));

   vertices.clear();
   elements.clear();
   cohelements.clear();
   faces.clear();
   edges.clear();
   return true;
}

/*
bool Mesh<Traits>::RemoveElement(ElementHandle e)
*/
template <typename Traits>
bool Mesh<Traits>::RemoveElement(ElementHandle e) {
   ElementMapIterator e_iter = std::find(ElementMapBegin(), ElementMapEnd(), e);
   assert(e_iter != ElementMapEnd());
   elements.erase(e_iter);
   elements_inactive.push_back(e);
   //delete e;
   return true;
}

/*
bool Mesh<Traits>::RemoveElements(std::vector<int> elems)
*/
template <typename Traits>
bool Mesh<Traits>::RemoveElements(std::vector<int> elems) {
   for(std::vector<int>::iterator iter=elems.begin(); iter!=elems.end(); ++iter) {
      RemoveElement(FindElementFromIndex(*iter));
   }
   return true;
}

/*
bool Mesh<Traits>::RemoveVertex(VertexHandle v)
*/
template <typename Traits>
bool Mesh<Traits>::RemoveVertex(VertexHandle v) {
   VertexMapIterator v_iter = std::find(VertexMapBegin(), VertexMapEnd(), v);
   assert(v_iter != VertexMapEnd());
   vertices.erase(v_iter);
   vertices_inactive.push_back(v);
   //delete v;
   return true;
}

/*
bool Mesh<Traits>::RemoveVertices(std::vector<int> verts)
*/
template <typename Traits>
bool Mesh<Traits>::RemoveVertices(std::vector<int> verts) {
   for(std::vector<int>::iterator iter=verts.begin(); iter!=verts.end(); ++iter) {
      RemoveVertex(FindVertexFromIndex(*iter));
   }
   return true;
}

/*
bool SetActive(std::vector<int> elems)
*/
template <typename Traits>
bool Mesh<Traits>::SetActive(std::vector<int> act_elems_vec) {
   std::copy(elements.begin(), elements.end(), std::back_inserter(elements_inactive));
   elements.clear();
   
   for(std::vector<int>::iterator act_iter=act_elems_vec.begin(); act_iter!=act_elems_vec.end(); ++act_iter) {
      ElementMapIterator found = std::find_if(elements_inactive.begin(), 
                                    elements_inactive.end(), 
                                    std::bind1st(IDEqual<ElementType>(), *act_iter));
      if(found != elements_inactive.end()) {
         elements.push_back(*found);
         elements_inactive.erase(found);
      }
   }
   
   std::sort(elements.begin(), elements.end(), IDComparison<ElementType>());
   
   std::copy(vertices.begin(), vertices.end(), std::back_inserter(vertices_inactive));
   vertices.clear();
   
   for(ElementMapIterator e_iter=elements.begin(); e_iter!=elements.end(); ++e_iter) {
      for(int i=0; i<ElementType::num_vertices; i++) {
         VertexMapIterator found = std::find(vertices_inactive.begin(), 
                                    vertices_inactive.end(), 
                                    (*e_iter)->vertices.at(i));
         if(found != vertices_inactive.end()) {
            vertices.push_back(*found);
            vertices_inactive.erase(found);
         }
      }
   }
   
   std::sort(vertices.begin(), vertices.end(), IDComparison<VertexType>());
   
   /*std::vector<int> all_elems_vec;
   for(ElementMapIterator e_iter=ElementMapBegin(); e_iter!=ElementMapEnd(); ++e_iter)
      all_elems_vec.push_back((*e_iter)->globalID);
   
   std::set<int> all_elems(all_elems_vec.begin(), all_elems_vec.end());
   std::set<int> act_elems(act_elems_vec.begin(), act_elems_vec.end());
   
   std::set<int> del_elems;
   std::set_difference(all_elems.begin(), all_elems.end(), 
                       act_elems.begin(), act_elems.end(), 
                  std::insert_iterator<std::set<int> >(del_elems, del_elems.begin()));
   std::vector<int> del_elems_vec(del_elems.begin(), del_elems.end());
   
   std::vector<int> all_verts_vec;
   for(VertexMapIterator v_iter=VertexMapBegin(); v_iter!=VertexMapEnd(); ++v_iter)
      all_verts_vec.push_back((*v_iter)->globalID);
   
   std::set<int> all_verts(all_verts_vec.begin(), all_verts_vec.end());
   
   std::set<int> act_verts;
   for(std::vector<int>::iterator de_iter=del_elems_vec.begin(); de_iter!=del_elems_vec.end(); ++de_iter) {
      ElementHandle e = FindElementFromIndex(*de_iter);
      for(int i=0; i<ElementType::num_vertices; i++) {
         act_verts.insert(e->vertices.at(i)->globalID);
      }
   }
   
   std::set<int> del_verts;
   std::set_difference(all_verts.begin(), all_verts.end(), 
                       act_verts.begin(), act_verts.end(), 
                  std::insert_iterator<std::set<int> >(del_verts, del_verts.begin()););
   std::vector<int> del_verts_vec(del_verts.begin(), del_verts.end());
   
   RemoveElements(del_elems_vec);
   RemoveVertices(del_verts_vec);*/
   
   return true;
}

/*
typename Mesh<Traits>::VertexHandle Mesh<Traits>::FindVertexFromIndex (int index)
*/
template <typename Traits>
typename Mesh<Traits>::VertexHandle Mesh<Traits>::FindVertexFromIndex (int index) {
   VertexMapIterator found = std::find_if(VertexMapBegin(), VertexMapEnd(), std::bind1st(IDEqual<VertexType>(), index));
   /*for(VertexMapIterator v_iter=VertexMapBegin(); v_iter!=VertexMapEnd(); ++v_iter) {
      if((*v_iter)->globalID == index)
         return (*v_iter);
   }
   return NULL*/
   if(found != VertexMapEnd())
      return (*found);
   else return NULL;
}

/*
typename Mesh<Traits>::ElementHandle Mesh<Traits>::FindElementFromIndex (int index)
*/
template <typename Traits>
typename Mesh<Traits>::ElementHandle Mesh<Traits>::FindElementFromIndex (int index) {
   ElementMapIterator found = std::find_if(ElementMapBegin(), ElementMapEnd(), std::bind1st(IDEqual<ElementType>(), index));
   /*for(ElementMapIterator e_iter=ElementMapBegin(); e_iter!=ElementMapEnd(); ++e_iter) {
      if((*e_iter)->globalID == index)
         return (*e_iter);
   }
   return NULL;*/
   if(found != ElementMapEnd()) return (*found);
   else return NULL;
}

/*
void Mesh<Traits>::ConnectElement(ElementHandle e)
*/
template <typename Traits>
void Mesh<Traits>::ConnectElement(ElementHandle e) {
   for(int i=0; i<ElementType::num_adjacents; i++) {
      FaceElementType f = e->GetFace(i);
      
      std::sort(f.vertices.begin(), f.vertices.end(), IDComparison<VertexType>());
      FaceMapIterator find_face = faces.find(f);
      
      if(find_face != FaceMapEnd()) {
         e->MakeConnection(i, faces[f]);
         faces[f]->MakeConnection(faces[f]->GetFaceIndex(f), e);
      }
      else {
         faces[f] = e;
      }
   }
}

/*
void Mesh<Traits>::AddElementToConLists(ElementHandle e)
*/
template <typename Traits>
void Mesh<Traits>::AddElementToConLists(ElementHandle e) {
   for(int i=0; i<ElementType::num_edges; i++) {
      EdgeElementType edge = e->GetEdge(i);
      int next_ind = 0, prev_ind = 0;
      
      std::sort(edge.vertices.begin(), edge.vertices.end(), IDComparison<VertexType>());
      
      //if(e->NextOnEdge(i) != NULL) next_ind = e->NextOnEdge(i)->globalID;
      //if(e->PrevOnEdge(i) != NULL) prev_ind = e->PrevOnEdge(i)->globalID;
      
      if(e->NextOnEdge(edge) != NULL) next_ind = e->NextOnEdge(edge)->globalID;
      if(e->PrevOnEdge(edge) != NULL) prev_ind = e->PrevOnEdge(edge)->globalID;
      
      EdgeMapIterator find_edge = edges.find(edge);
      if(find_edge == EdgeMapEnd()) {
         edges[edge] = new ElementConListType();
      }
      edges[edge]->add_elem(e->globalID, prev_ind, next_ind);
      e->con_lists.at(i) = edges[edge];
      /*if(EdgeElementType::num_vertices > 1) {
         EdgeMapIterator find_not_edge = edges.find(!edge);
         if(find_not_edge == EdgeMapEnd()) {
            edges[!edge] = new ElementConListType();
         }
         edges[!edge]->add_elem(e->globalID, next_ind, prev_ind);
      }*/
      
      /* if(edge.vertices[0]->globalID == 950) {
         std::cout << "ELEM: " << e->globalID << ", NEXT: " << next_ind << ", PREV: " << prev_ind << std::endl;
         std::cout << " " << edge.vertices.at(0)->globalID << "\t";
         std::cout << edges[edge]->isClosed() << "\t";
         std::cout << edges[edge]->needSplit() << "\t";
         for(typename ElementConListType::iterator list_iter=edges[edge]->begin(); list_iter!=edges[edge]->end(); ++list_iter) {
            std::cout << (*list_iter).elem_globalID << "(" << (*list_iter).link_state << ")" << "->";
         }
         std::cout << std::endl;
      } */
   }
}

/*
std::vector<typename Mesh<Traits>::VertexHandle> Mesh<Traits>::GetCentres()
Return a list of vertices denoting the centres of all the elements in the mesh.
Each vertex will have its global ID set to the global ID of the element
*/
template <typename Traits>
std::vector<typename Mesh<Traits>::VertexHandle> Mesh<Traits>::GetCentres() {
   std::vector<VertexHandle> centres;
   
   for(ElementMapIterator e_iter=ElementMapBegin(); e_iter!=ElementMapEnd(); ++e_iter) {
      VertexHandle centre_vertex = (*e_iter)->GetCentre();
      centre_vertex->globalID = (*e_iter)->globalID;
      centres.push_back(centre_vertex);
   }
   
   return centres;
}

} // namespace DAPTA

---