Mesh.cpp

Go to the documentation of this file.

/* Generated by Together */

#include <string.h>   /* added by Zaczek,M 6/7/2001 ... for memcpy() */
#include "Mesh.hpp"
#include "Element.hpp"
#include "Face.hpp"
#include "Node.hpp"
#include "IntList.hpp"
#include "NodeList.hpp"
#include "ElementList.hpp"
#include "FaceList.hpp"
#include "FaceListList.hpp"

Mesh::Mesh() :
  d_nodes(0),
  d_numNodes(0),
  d_nodeArraySize(0),
  d_faces(0),
  d_numFaces(0),
  d_faceArraySize(0),
  d_elements(0),
  d_numElements(0),
  d_elementArraySize(0)
{

  Element::setMesh( this );
  Face::setMesh( this );
  minEdge = 1E+25;

}

Mesh::~Mesh(){

  if( d_elements ){
    int i;
    for( i = 0; i < d_numElements; i++ ){
      delete d_elements[i];
    }
    delete [] d_elements;
  }

  if( d_faces ){
    int i;
    for( i = 0; i < d_numFaces; i++ ){
      delete d_faces[i];
    }
    delete [] d_faces;
  }

  if( d_nodes ){
    int i;
    for( i = 0; i < d_numNodes; i++ ){
      delete d_nodes[i];
    }
    delete [] d_nodes;
  }
}


void Mesh::addNode(Node* node){

  if( d_nodeArraySize == d_numNodes ){
    d_nodeArraySize =  (int)( d_nodeArraySize * 1.25 + 3);
    Node **new_arr = new Node*[d_nodeArraySize];
    if( d_numNodes > 0 ){
      memcpy(new_arr, d_nodes, d_numNodes * sizeof(Node*) );
      delete [] d_nodes;
    }
    d_nodes = new_arr;
  }
  d_nodes[d_numNodes] = node;
  node->setID(d_numNodes + 1);
  d_numNodes++;

}

void Mesh::removeNode(Node * node){

  int i;
  int dif = 0;
  for( i = 0; i < d_numNodes; i++ ){
    if( d_nodes[i] == node ){
      dif++;
      continue;
    }
    if( dif > 0 ){
      d_nodes[i - dif] = d_nodes[i];
      d_nodes[i - dif]->setID( i - dif + 1 ); 
    }
  }
  d_numNodes -= dif;

}


void Mesh::addFace(Face* face){

  if( d_faceArraySize == d_numFaces ){
    d_faceArraySize =  (int)( d_faceArraySize * 1.25 + 3);
    Face **new_arr = new Face*[d_faceArraySize];
    if( d_numFaces > 0 ){
      memcpy(new_arr, d_faces, d_numFaces * sizeof(Face*) );
      delete [] d_faces;
    }
    d_faces = new_arr;
  }
  d_faces[d_numFaces] = face;
  face->setID(d_numFaces + 1);
  d_numFaces++;
}

void Mesh::removeFace(Face * face){

  int i;
  int dif = 0;
  for( i = 0; i < d_numFaces; i++ ){
    if( d_faces[i] == face ){
      dif++;
      continue;
    }
    if( dif > 0 ){
      d_faces[i - dif] = d_faces[i];
      d_faces[i - dif]->setID( i - dif + 1 ); 
    }
  }
  d_numFaces -= dif;


}

void Mesh::addElement(Element* element){

  if( d_elementArraySize == d_numElements ){
    d_elementArraySize =  (int)( d_elementArraySize * 1.25 + 3);
    Element **new_arr = new Element*[d_elementArraySize];
    if( d_numElements > 0 ){
      memcpy(new_arr, d_elements, d_numElements * sizeof(Element*) );
      delete [] d_elements;
    }
    d_elements = new_arr;
  }
  d_elements[d_numElements] = element;
  element->setID(d_numElements + 1);
  d_numElements++;

}

void Mesh::removeElement(Element * element){

  int i;
  int dif = 0;
  for( i = 0; i < d_numElements; i++ ){
    if( d_elements[i] == element ){
      dif++;
      continue;
    }
    if( dif > 0 ){
      d_elements[i - dif] = d_elements[i];
      d_elements[i - dif]->setID( i - dif + 1 ); 
    }
  }
  d_numElements -= dif;

}


boolean Mesh::addCohesive( int material1, int material2, int new_material ){

  // go over all faces & for each one which has two elems
  // with the mat1 & mat2 - generate a cohesive (if neither of
  // the two  elems is cohesive) - have a list of those (to split
  // at nodes later)
  
  FaceList faces;

  // collect the faces
  int i;
  for( i = 0; i < d_numFaces; i++ ){
    Element *e1 = d_faces[i]->getElement1();
    Element *e2 = d_faces[i]->getElement2();
    if( e1 && e2 && !e1->isCohesive()
        && !e2->isCohesive()
        && ( ( e1->getMaterialType() == material1 
               && e2->getMaterialType() == material2 )
             || ( e1->getMaterialType() == material2 
                  && e2->getMaterialType() == material1 ) ) ){
      // fits
      faces.insert_first( d_faces[i] );
    }
  }

  // Collect nodes  & their faces

  NodeList sep_nodes;
  FaceListList node_faces;
  int size = faces.size();
  faces.reset();
  for( i = 0; i < size; i++ ){
    Face* face = faces.get();
    faces.next();
    int num = face->getNumNodes();
    Node** fnodes = face->getNodes();
    int j;
    for( j = 0; j < num; j++ ){
      if( sep_nodes.move_to( fnodes[j] ) ){
                int ind = sep_nodes.index();
                node_faces.index( ind );
                node_faces.get()->insert_first( face );
      }
      else {
                sep_nodes.insert_first( fnodes[j] );
                FaceList* newlist = new FaceList;
                newlist->insert_first( face );
                node_faces.insert_first( newlist );
      }
    }
  }

  // now go over the nodes & check if separable...

  // check at each node 
  // if have two element groups separated by cohesives (or our faces)
  // if yes separate the two. 
  
  int sep_count = 0;
  int node_size = sep_nodes.size();
  sep_nodes.reset();
  node_faces.reset();
  for( i = 0; i < node_size; i++ ) {
        int j = i;
    Node* node = sep_nodes.get();
    FaceList* flist = node_faces.get();
    sep_nodes.next();
    node_faces.next();
    boolean separated = node->separate( this, flist, new_material );
    if( separated ) sep_count++;
    delete flist;
  }
  cerr << " separated " << sep_count << " nodes, out of " << node_size 
       << " involved" << endl;
  
  return ( sep_count >0 ? TRUE : FALSE );
}

void  Mesh::replaceNode( Node* node, Node* new_node, 
                         ElementList *sep_elements, int new_material ){

  // add any cohesives which have one side inside the sep_elements

  int size = sep_elements->size();
  //out_faces will ultimately contain the list of faces that form the outside
  //of the sep_elements group of elements. 
  FaceList out_faces;
  int i;
  //for each element in the sep_elements list...
  for( i = 0; i < size; i++ ) {
    sep_elements->index(i);
    Element* elem = sep_elements->get();
    int numf = elem->getNumFaces();
    Face** faces = elem->getFaces();
    int j;
    //for each face in each element...
    for( j = 0; j < numf; j++ ){
      Face* face = faces[j];
      
      //see of this face contains the node we're working on...
      if( !face->containsNode( node ) ){
                continue;  //try next face
      }
      
      //if yes, get the element on the other side of this face...
      Element* oelem= (face->getElement1() == elem 
                    ? face->getElement2()
                    : face->getElement1() );
                    
          //if its an outside face (no second element), then go to the next face...
      if( !oelem ){
                continue;
      }
      
      //if the second element is already on our element list, go to next face...
      if( sep_elements->move_to( oelem ) ){
                continue;
      }
      
      //otherwise, if it's a cohesive element on the other side of the face...
      if( oelem->isCohesive() ) {
        //get the second face of the cohesive element (first is in "face")
                Face** cofaces = oelem->getFaces();
                Face* oface = ( cofaces[0] == face ? cofaces[1] : cofaces[0] );
                //now get the element that is on the other side of the cohesive element
                Element* other = (oface->getElement1() == oelem 
                          ? oface->getElement2()
                          : oface->getElement1() );
                //if that other element is already on our list of elements, then append 
                //the first one to the list (it wasn't on our list - we checked above). 
                if( sep_elements->move_to( other ) ) {
                        sep_elements->append( oelem );
                }
                else { 
                        //otherwise add the face we're working with to the out_faces list. 
                        //(if the opposite element does not belong in our element list, then it is
                        //an outside "separator" face)
                        out_faces.insert_first( face );
                }
      } // if( oelem->isCohesive() )
      else { //if the other element is not cohesive, and is not in the sep_elements list, 
                //then it is a boundary.
                out_faces.insert_first( face );
      }
      
    } //for( j = 0 ...
  } //  for( i = 0 ...
  
  int num_out_faces = out_faces.size();
   
  //Task:  if the face does not belong to a cohesive element, then we have to make one. 
  //Otherwise, we need to replace the node in the current face.  Somewhere here there
  //is a bug right now (8/13/01)
  
  //for each face in the out_faces array, get the elements on either side.
  for( i = 0; i < num_out_faces; i++ ){
    Face *face = out_faces.get();
    out_faces.next();
    Element* e1 = face->getElement1();
    Element* e2 = face->getElement2();
    if( !e2 ) { // the sep replacement will handle
      continue;
    }
    //if neither of them are cohesive, add a cohesive element between them.
    if( !e1->isCohesive() && !e2->isCohesive() ){
    //elem is used to ensure correct direction of the normals. 
      Element* elem = ( sep_elements->move_to( e1 ) ? e1 : e2 );
      Element* new_elem = face->buildCohesive( elem, node, new_node );
      new_elem->setMaterialType( new_material );
      addElement( new_elem );
   }
    else if( e1->isCohesive() ){ // e2 is in sep
      e1->replaceFaceNode( node, new_node, face );
    }
    else {
      e2->replaceFaceNode( node, new_node, face );
    }
  }

  // go and replace in all the elements
  // including inside the faces 
  //replace nodes in the elements in the sep_element list.  Bug may be here as well. 
  sep_elements->reset();
  size = sep_elements->size();
  for( i = 0; i < size; i++ ){

    sep_elements->get()->replaceNode( node, new_node );
    sep_elements->next();
  }
}


ostream & operator<<(ostream & stream, Mesh & mesh){

  int j;
  int numcohesive = 0;
  for (j=0; j<mesh.d_numElements; j++) {
        if ((*mesh.d_elements[j]).isCohesive()) numcohesive += 1;
  }

  stream << (mesh.d_numElements - numcohesive) << '\t' <<  mesh.d_numNodes 
         << '\t' << numcohesive << '\t' << mesh.getMinEdge() <<endl;

  stream.setf(ios::uppercase);
  stream.setf(ios::scientific,ios::floatfield);
  stream.precision(9);
  int numbc_nodes = 0;
  int i;
  for( i = 0; i < mesh.d_numNodes;i++ ){
    stream << (*mesh.d_nodes[i]); 
    if( mesh.d_nodes[i]->getFlag() != Node::e_unset_flag ){
      numbc_nodes++;
    }
  }
  stream << numbc_nodes << endl;
  for( i = 0; i < mesh.d_numNodes;i++ ){
    if( mesh.d_nodes[i]->getFlag() != Node::e_unset_flag ){
      stream <<  mesh.d_nodes[i]->getID() << '\t' << mesh.d_nodes[i]->getFlag()<< endl;
    }
  }
  for( i = 0; i < mesh.d_numElements; i++ ){
    stream << (*mesh.d_elements[i]);
  }
  
  return stream;
}


//Returns the length of the shortest edge in an element.
double Mesh::calcEdgeLength(int elementID){

   return(d_elements[elementID]->getMinEdgeLength());

}

istream & operator>>(istream & stream, Mesh & mesh){

  int numn;
  int nume;
  stream >> nume >> numn;
  
  mesh.d_nodeArraySize = numn;
  mesh.d_nodes = new Node*[numn];
  // read nodes
  int i;
  for( i = 0; i < numn;i++ ){
    mesh.d_nodes[i] = new Node();
    stream >> (*mesh.d_nodes[i]);
  }
  mesh.d_numNodes = numn;

  mesh.d_elementArraySize = nume;
  mesh.d_elements = new Element*[nume];

  // assume ~ twice as elements

  mesh.d_faceArraySize = nume * 2;
  mesh.d_faces = new Face*[mesh.d_faceArraySize];


  for( i = 0; i < nume;i++ ){
    int id, type;
    double currentEdge;
    stream >> id >> type;
    mesh.d_elements[i] = Element::create( id, (Element::Type)type );
    stream >> (*mesh.d_elements[i]);
    currentEdge = (mesh.calcEdgeLength(i));
    if (currentEdge < mesh.getMinEdge()) mesh.setMinEdge(currentEdge);
  }
  mesh.d_numElements = nume;
  // read groups

  int numbcn;
  stream >> numbcn;
  // read nodes with flag
  for( i = 0; i < numbcn; i++ ){
    int num;
    int flag;
    stream >> num >> flag;
    mesh.d_nodes[num-1]->setFlag( flag );
  }
  
  //read faces with flag (three nodes and the flag)
  int numFaceFlags, numFaceNodes, elementtype, theflag, totalFaces;
  
  stream >> numFaceFlags;  //number of faces in input file with flags
  totalFaces = mesh.getNumFaces(); //total number of existing faces in mesh
  
  for (i=0 ; i < numFaceFlags ; i++) {
        stream >> elementtype;  //type of element the face is associated with
        //make a temporary working face to hold the node info
        Face* theFace = Face::create((Face::eType)elementtype);
        stream >> *theFace;
        stream >> theflag;
        numFaceNodes = theFace->getNumNodes(); //number of nodes in this face
        
        //take "theFace" and compare it's nodes against the nodes in every face in the 
        //Mesh list of faces until we find a match.  When we find a match, set it's flag,
        //and then go on to the next face in the list.  This is very inefficient, but
        //without knowing what element the supplied face is in, I don't know how else to
        //go about it. 
        int j;
        for (j=0; j < totalFaces; j++) {
                Face* trialFace = mesh.getFace(j+1);  //getFace assumes ID's starting with 1
                Node** theNodes = theFace->getNodes();
                Node** trialNodes = trialFace->getNodes();
                boolean nodeFound;
                int numFound;
                numFound=0;
                
                //for each node in the face to be found...
                int k;
                for (k=0; k<numFaceNodes; k++) {
                        nodeFound=FALSE;
                        int numExistingFaceNodes;
                        numExistingFaceNodes = trialFace->getNumNodes();
                        if (numFaceNodes != numExistingFaceNodes) break; //faces aren't compatible
                        
                        //...see if it's in the current face in the mesh that we're looking at.
                        int m;
                        for (m=0; m<numExistingFaceNodes; m++) {
                                if(theNodes[k]==trialNodes[m]) {
                                        nodeFound=TRUE;
                                        numFound++;
                                        break;
                                }
                        }
                        if(!nodeFound) break;  //if we didn't find the node, try the next face
                        //if we found all the nodes in any face, then set the flag and read the next face in.
                        if(numFound==numFaceNodes) {
                                trialFace->setFlag(theflag);
                                break;
                        }
                }
        }
        delete theFace;
  }
  return stream;
}

void  Mesh::write_boundary(ostream & stream){

  int count = 0;
  int i;
  for( i = 0; i < d_numFaces; i++ ){
    if( d_faces[i]->isExterior() ){
      count++;
    }
  }
  stream << count << endl;

  for( i = 0; i < d_numFaces; i++ ){
    if( d_faces[i]->isExterior() ){
      Element* e1 = d_faces[i]->getElement1();
      stream << e1->getID() << '\t';
      int numf = e1->getNumFaces();
      Face** efaces = e1->getFaces();
      int j;
      for( j = 0; j < numf; j++ ){
                if( d_faces[i] == efaces[j] ){
                        stream << j + 1 << '\t' << d_faces[i]->getFlag() << endl;
                        break;
                }
      }
    }
  }
}