Mesh structure code for elements implemented.

AMDiS/src/Element.h

@@ -370,6 +370,13 @@ namespace AMDiS {
     /// Returns whether Element has sideElem as one of its sides.
     virtual bool hasSide(Element *sideElem) const = 0;
 
+    /** \brief
+     * Returns for a given element type number the element type number of the children.
+     * For 1d and 2d this is always 0, because element type number are used in the 
+     * 3d case only.
+     */
+    virtual int getChildType(int elType) const = 0;
+
     /** \brief
      * Traverses an edge/face of a given element (this includes also all children of 
      * the element having the same edge/face). All vertex dofs alonge this edge/face

AMDiS/src/Line.h

@@ -141,6 +141,12 @@ namespace AMDiS {
     { 
       return false; 
     }
+
+    /// Element type number is not used in 1d, so return 0.
+    inline int getChildType(int) const
+    {
+      return 0;
+    }
   
     std::string getTypeName() const 
     { 

AMDiS/src/MeshStructure.cc

@@ -54,6 +54,32 @@ namespace AMDiS {
     commit();
   }
 
+  void MeshStructure::init(Element *el, int ithSide, int elType)
+  {
+    FUNCNAME("MeshStructure::init()");
+
+    clear();
+
+    addAlongSide(el, ithSide, elType);
+
+    commit();    
+  }
+
+  void MeshStructure::addAlongSide(Element *el, int ithSide, int elType)
+  {
+    insertElement(el->isLeaf());
+    
+    if (!el->isLeaf()) {
+      int s1 = el->getSideOfChild(0, ithSide, elType);
+      int s2 = el->getSideOfChild(1, ithSide, elType);
+
+      if (s1 != -1) 
+	addAlongSide(el->getFirstChild(), s1, el->getChildType(elType));
+      if (s2 != -1)
+	addAlongSide(el->getSecondChild(), s2, el->getChildType(elType));
+    }
+  }
+
   void MeshStructure::reset() 
   {
     currentIndex = 0;

AMDiS/src/MeshStructure.h

@@ -41,9 +41,11 @@ namespace AMDiS {
   
     void clear();
 
-    /// Creates a mesh structure code from a Mesh object by traversing it in preorder.
+    /// Creates a mesh structure code from a mesh object by traversing it in preorder.
     void init(Mesh *mesh);
 
+    void init(Element *el, int ithSide, int elType);
+
     void init(const std::vector<unsigned long int>& initCode, int n) 
     {
       code = initCode;
@@ -102,13 +104,13 @@ namespace AMDiS {
       bool cont = true;
       while (cont) {
 	if (isLeafElement())
-	  MSG("0");
+	  std::cout << "0";
 	else
-	  MSG("1");
+	  std::cout << "1";
 	
 	cont = nextElement();
       }
-      MSG("\n");
+      std::cout << std::endl;
     }
 
     /// Returns the mesh structure code.
@@ -131,6 +133,8 @@ namespace AMDiS {
     /// Insert a new element to the structure code. Is used by the init function.
     void insertElement(bool isLeaf);
 
+    void addAlongSide(Element *el, int ithSide, int elType);
+
     /// Merges two mesh structure codes to one structure code.     
     void merge(MeshStructure *structure1,
 	       MeshStructure *structure2,

AMDiS/src/ParallelDomainBase.cc

@@ -18,6 +18,7 @@
 #include "ElementFileWriter.h"
 #include "VertexVector.h"
 #include "StdMpi.h"
+#include "MeshStructure.h"
 
 #include "petscksp.h"
 
@@ -694,7 +695,7 @@ namespace AMDiS {
       stdMpi.send(sendIt->first, dofs);
     }
 
-    stdMpi.startCommunication();
+    stdMpi.startCommunication<int>();
 #endif
 
 #if 1
@@ -758,9 +759,56 @@ namespace AMDiS {
 
   void ParallelDomainBase::checkMeshChange()
   {
+    // === If mesh has not been changed, return. ===
+
     if (mesh->getChangeIndex() == lastMeshChangeIndex)
       return;
 
+
+    // === Create mesh structure codes for all ranks boundary elements. ===
+
+    typedef std::vector<MeshStructure> MeshCodeVec;
+    std::map<int, MeshCodeVec > sendCodes;
+
+    for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
+	 it != myIntBoundary.boundary.end(); ++it) {    
+
+      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+	   boundIt != it->second.end(); ++boundIt) {
+
+	MeshStructure elCode;
+	elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, 
+		    boundIt->rankObj.elType);
+	sendCodes[it->first].push_back(elCode);
+      }
+    }
+
+    StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm);
+    stdMpi.prepareCommunication(true);
+    stdMpi.send(sendCodes);
+
+    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
+	 it != otherIntBoundary.boundary.end(); ++it)
+      stdMpi.recv(it->first);
+
+    stdMpi.startCommunication<unsigned long int>();
+
+      
+#if 0
+    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
+	 it != otherIntBoundary.boundary.end(); ++it) {
+
+      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+	   boundIt != it->second.end(); ++boundIt) {
+
+	MeshStructure elCode;
+	elCode.init(boundIt->rankObj.el, 
+		    boundIt->rankObj.ithObj, 
+		    boundIt->rankObj.elType);
+      }
+    }
+#endif
+
     std::cout << "MESH HAS BEEN CHANGED!" << std::endl;
     exit(0);    
 
@@ -1262,7 +1310,7 @@ namespace AMDiS {
     for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
 	 recvIt != recvNewDofs.end(); ++recvIt, i++)
       stdMpi.recv(recvIt->first, recvIt->second * 2); 
-    stdMpi.startCommunication();
+    stdMpi.startCommunication<int>();
     std::map<int, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > > &dofMap = stdMpi.getRecvData();
 
     // === Change dof indices at boundary from other ranks. ===

AMDiS/src/StdMpi.h

@@ -24,6 +24,7 @@
 
 #include <map>
 #include "mpi.h"
+#include "MeshStructure.h"
 
 namespace AMDiS {
 
@@ -32,7 +33,16 @@ namespace AMDiS {
     return data.size() * 2;
   }
 
-  void makeIntBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf)
+  int intSizeOf(std::vector<MeshStructure> &data)
+  {
+    int s = 0;
+    for (int i = 0; i < data.size(); i++)
+      s += data[i].getCode().size() + 2;
+
+    return s;
+  }
+
+  void makeBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf)
   {
     int i = 0;
     for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator it = data.begin();
@@ -42,8 +52,8 @@ namespace AMDiS {
     }
   }
   
-  void makeFromIntBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data, 
-		      int *buf, int bufSize)
+  void makeFromBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data, 
+		   int *buf, int bufSize)
   {    
     if (bufSize == 0)
       return;
@@ -58,7 +68,35 @@ namespace AMDiS {
     } while (i < bufSize);
   }
 
+  void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf)
+  {
+    int pos = 0;
+    for (int i = 0; i < data.size(); i++) {
+      buf[pos++] = data[i].getCode().size();
+      buf[pos++] = data[i].getNumElements();
+      for (int j = 0; j < data[i].getCode().size(); j++)
+	buf[pos++] = data[i].getCode()[j];
+    }
+  }
 
+  void makeFromBuf(std::vector<MeshStructure> &data, unsigned long int *buf, int bufSize)
+  {
+    int pos = 0;
+
+    while (pos < bufSize) {       
+      int codeSize = buf[pos++];
+      int nElements = buf[pos++];
+      std::vector<unsigned long int> code;
+      code.reserve(codeSize);
+      for (int i = 0; i < codeSize; i++)
+	code[i] = buf[pos++];
+
+      MeshStructure meshCode;
+      meshCode.init(code, nElements);
+
+      data.push_back(meshCode);
+    }	
+  }
 
   template<typename SendT, typename RecvT>
   class StdMpi 
@@ -116,21 +154,48 @@ namespace AMDiS {
       return recvData;
     }
 
+    void commDataSize()
+    {
+      MPI::Request request[sendData.size() + recvDataSize.size()];
+      std::vector<int> sendBuffers, recvBuffers;
+      int requestCounter = 0;
+
+      int i = 0;
+      for (typename std::map<int, int>::iterator sendIt = sendDataSize.begin();
+	   sendIt != sendDataSize.end(); ++sendIt) {
+	sendBuffers.push_back(sendIt->second);
+	request[requestCounter++] =
+	  mpiComm.Isend(&(sendBuffers[i]), 1, MPI_INT, sendIt->first, 0);
+	i++;
+      }
+
+      for (std::map<int, int>::iterator recvIt = recvDataSize.begin();
+	   recvIt != recvDataSize.end(); ++recvIt)
+	request[requestCounter++] =
+	  mpiComm.Irecv(&(recvIt->second), 1, MPI_INT, recvIt->first, 0);
+      
+      MPI::Request::Waitall(requestCounter, request);      
+    }
+
+    template<class T>
     void startCommunication()
     {
       FUNCNAME("StdMpi::startCommunication()");
 
       TEST_EXIT_DBG(commPrepared)("Communication is not prepared!\n");
 
+      if (exchangeDataSize)
+	commDataSize();
+
       MPI::Request request[sendData.size() + recvDataSize.size()];
       int requestCounter = 0;
-      std::vector<int*> sendBuffers, recvBuffers;
+      std::vector<T*> sendBuffers, recvBuffers;
 
       for (typename std::map<int, SendT>::iterator sendIt = sendData.begin();
 	   sendIt != sendData.end(); ++sendIt) {
-	int bufferSize = intSizeOf(sendIt->second);
-	int* buf = new int[bufferSize];
-	makeIntBuf(sendIt->second, buf);
+	int bufferSize = sendDataSize[sendIt->first];
+	T* buf = new T[bufferSize];
+	makeBuf(sendIt->second, buf);
 
   	request[requestCounter++] = 
   	  mpiComm.Isend(buf, bufferSize, MPI_INT, sendIt->first, 0);
@@ -140,7 +205,7 @@ namespace AMDiS {
 
       for (std::map<int, int>::iterator recvIt = recvDataSize.begin();
 	   recvIt != recvDataSize.end(); ++recvIt) {
-	int* buf = new int[recvIt->second];
+	T* buf = new T[recvIt->second];
 
   	request[requestCounter++] =
   	  mpiComm.Irecv(buf, recvIt->second, MPI_INT, recvIt->first, 0);
@@ -157,7 +222,7 @@ namespace AMDiS {
       int i = 0;
       for (std::map<int, int>::iterator recvIt = recvDataSize.begin();
 	   recvIt != recvDataSize.end(); ++recvIt) {
-	makeFromIntBuf(recvData[recvIt->first], recvBuffers[i], recvIt->second);
+	makeFromBuf(recvData[recvIt->first], recvBuffers[i], recvIt->second);
 	delete [] recvBuffers[i];
 	i++;
       }      

AMDiS/src/Tetrahedron.h

@@ -123,6 +123,12 @@ namespace AMDiS {
     { 
       return true; 
     }
+
+    /// Return the new element type of the children.
+    inline int getChildType(int elType) const
+    {
+      return (elType + 1) % 3;
+    }
  
     std::string getTypeName() const 
     { 

AMDiS/src/Triangle.h

@@ -114,6 +114,12 @@ namespace AMDiS {
       return false; 
     }
 
+    /// Element type number is not used in 2d, so return 0.
+    inline int getChildType(int) const
+    {
+      return 0;
+    }
+
     /// implements Element::getSideOfChild()
     virtual int getSideOfChild(int child, int side, int) const 
     {