Commits

Dmitry Karpeyev committed d440d83 Merge

Comments (0)

Files changed (7)

src/dm/mesh/mesh.c

 #define __FUNCT__ "preallocateMatrix"
 PetscErrorCode preallocateMatrix(ALE::Mesh *mesh, const ALE::Obj<ALE::Mesh::section_type>& field, const ALE::Obj<ALE::Mesh::numbering_type>& globalOrder, Mat A)
 {
-  const ALE::Obj<ALE::Mesh::sieve_type>     serialGraph   = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
-  const ALE::Obj<ALE::Mesh::sieve_type>     parallelGraph = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
+        ALE::Obj<ALE::Mesh::topology_type>  parallelTopology = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
+  const ALE::Obj<ALE::Mesh::sieve_type>     serialGraph      = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
+  const ALE::Obj<ALE::Mesh::topology_type>  serialTopology   = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
   const ALE::Mesh::section_type::patch_type patch = 0;
   const ALE::Obj<ALE::Mesh::sieve_type>&    sieve = mesh->getTopologyNew()->getPatch(patch);
   PetscInt       numLocalRows, firstRow, lastRow;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
+  serialTopology->setPatch(patch, serialGraph);
   ierr = MatGetLocalSize(A, &numLocalRows, PETSC_NULL);CHKERRQ(ierr);
   ierr = MatGetOwnershipRange(A, &firstRow, &lastRow);CHKERRQ(ierr);
   ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr);
     serialGraph->addCone(sieve->cone(sieve->support(point)), point);
   }
   /* Distribute adjacency graph */
-  ALE::New::Distribution<serialGraph::topology_type>::scatterSieve(serialGraph, parallelGraph);
+  const ALE::Obj<ALE::Mesh::send_overlap_type>&              vertexSendOverlap = mesh->getVertexSendOverlap();
+  const ALE::Obj<ALE::Mesh::recv_overlap_type>&              vertexRecvOverlap = mesh->getVertexRecvOverlap();
+  const ALE::Obj<ALE::Mesh::numbering_type>&                vNumbering  = mesh->getLocalNumbering(0);
+  const ALE::Obj<ALE::Mesh::topology_type::label_sequence>& vertices    = field->getAtlas()->getTopology()->depthStratum(patch, 0);
+  int                                                       numVertices = vertices->size();
+  short                                                    *assignment  = new short[numVertices];
+  int                                                       rank        = mesh->commRank();
+
+  for(ALE::Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
+    // FIX if (vertexOverlap->capContains(*v_iter)) {
+    if ((rank && vertexRecvOverlap->baseContains(*v_iter)) || (!rank && vertexSendOverlap->capContains(*v_iter))) {
+      int minRank = field->commSize();
+
+      if (!rank) {
+        const Obj<ALE::Mesh::send_overlap_type::traits::supportSequence>& sendPatches = vertexSendOverlap->support(*v_iter);
+
+        for(ALE::Mesh::send_overlap_type::traits::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
+          if (*p_iter < minRank) minRank = *p_iter;
+        }
+      } else {
+        const Obj<ALE::Mesh::recv_overlap_type::traits::coneSequence>& recvPatches = vertexRecvOverlap->cone(*v_iter);
+
+        for(ALE::Mesh::recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != recvPatches->end(); ++p_iter) {
+          if (*p_iter < minRank) minRank = *p_iter;
+        }
+      }
+      if (minRank < rank) {
+        assignment[vNumbering->getIndex(*v_iter)] = minRank;
+      } else {
+        assignment[vNumbering->getIndex(*v_iter)] = rank;
+      }
+    } else {
+      assignment[vNumbering->getIndex(*v_iter)] = rank;
+    }
+  }
+  ALE::New::Distribution<ALE::Mesh::topology_type>::scatterTopology(serialTopology, numVertices, assignment, parallelTopology);
   /* Read out adjacency graph */
+  const ALE::Obj<ALE::Mesh::sieve_type> parallelGraph = parallelTopology->getPatch(patch);
+
   ierr = PetscMemzero(dnz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
   ierr = PetscMemzero(onz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
   for(ALE::Mesh::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
     const ALE::Mesh::atlas_type::point_type& point = c_iter->first;
 
     if (globalOrder->isLocal(point)) {
-      const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = graph->cone(point);
+      const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = parallelGraph->cone(point);
       const int                                               row   = globalOrder->getIndex(point);
       const int                                               rSize = c_iter->second.index;
 
       }
     }
   }
+  for(int r = 0; r < numLocalRows; r++) {
+    std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
+  }
   ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
   ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
   ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);

src/dm/mesh/sieve/CoSieve.hh

       const Obj<capSequence>& cap() {return this->_points;};
       const Obj<baseSequence>& leaves() {return this->_points;};
       const Obj<baseSequence>& base() {return this->_points;};
+      template<typename Color>
+      void addArrow(const point_type& p, const point_type& q, const Color& color) {
+        throw ALE::Exception("Cannot add an arrow to a DiscreteSieve");
+      };
+      void stratify() {};
       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
         ostringstream txt;
         int rank;

src/dm/mesh/sieve/Completion.hh

       };
       const value_type *restrict(const patch_type& patch, const point_type& p) {return this->restrictPoint(patch, p);};
       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {
-        if (patch != p) {
-          throw ALE::Exception("Point must be identical to patch in a PartitionSizeSection");
-        }
-        return &this->_sizes[patch];
+        return &this->_sizes[p];
       };
       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
         throw ALE::Exception("Cannot update a PartitionSizeSection");
       };
       const value_type *restrict(const patch_type& patch, const point_type& p) {return this->restrictPoint(patch, p);};
       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {
-        if (patch != p) {
-          throw ALE::Exception("Point must be identical to patch in a PartitionSection");
-        }
-        return this->_points[patch];
+        return this->_points[p];
       };
       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
         throw ALE::Exception("Cannot update a PartitionSection");
       };
       template<typename Filler, typename Section>
       static void fillSend(const Obj<Filler>& sendFiller, const Obj<Section>& sendSection) {
-        const topology_type::sheaf_type& patches = sendSection->getAtlas()->getTopology()->getPatches();
-        const topology_type::patch_type  patch   = 0; // FIX: patch should come from overlap
+        const topology_type::sheaf_type& ranks = sendSection->getAtlas()->getTopology()->getPatches();
+        const topology_type::patch_type  patch = 0; // FIX: patch should come from overlap
 
-        for(topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
+        for(topology_type::sheaf_type::const_iterator p_iter = ranks.begin(); p_iter != ranks.end(); ++p_iter) {
+          const int&                                          rank = p_iter->first;
           const Obj<topology_type::sieve_type::baseSequence>& base = p_iter->second->base();
 
           for(topology_type::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
-            sendSection->update(p_iter->first, *b_iter, sendFiller->restrict(patch, *b_iter));
+            sendSection->update(rank, *b_iter, sendFiller->restrict(patch, *b_iter));
           }
         }
       };

src/dm/mesh/sieve/Distribution.hh

         Obj<Mesh::sieve_type>  sieve         = topology->getPatch(0);
         const Obj<typename send_overlap_type::traits::capSequence> cap = cellOverlap->cap();
 
+        parallelMesh->setVertexSendOverlap(vertexOverlap);
         for(typename send_overlap_type::traits::baseSequence::iterator p_iter = cap->begin(); p_iter != cap->end(); ++p_iter) {
           const Obj<typename send_overlap_type::traits::supportSequence>& ranks = cellOverlap->support(*p_iter);
 
         Obj<Mesh::sieve_type>  parallelSieve = parallelTopology->getPatch(0);
         const Obj<typename send_overlap_type::traits::baseSequence> base = cellOverlap->base();
 
+        parallelMesh->setVertexRecvOverlap(vertexOverlap);
         for(typename send_overlap_type::traits::baseSequence::iterator p_iter = base->begin(); p_iter != base->end(); ++p_iter) {
           const Obj<typename send_overlap_type::traits::coneSequence>& ranks = cellOverlap->cone(*p_iter);
 
       #undef __FUNCT__
       #define __FUNCT__ "createScatterOverlap"
       static void createScatterOverlap(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
-        // There are arrows to each rank whose color is the partition point (also the rank)
-        for(int p = 1; p < sendOverlap->commSize(); p++) {
-          sendOverlap->addCone(p, p, p);
+        const int rank = sendOverlap->commRank();
+
+        std::cout << "Rank: " << rank << " Size: " << sendOverlap->commSize() << std::endl;
+        for(int p = 0; p < sendOverlap->commSize(); p++) {
+          if (p != rank) {
+            // There are arrows to each rank whose color is the partition point (also the rank)
+            sendOverlap->addCone(p, p, p);
+            //  The arrow is from rank 0 with partition point 0
+            recvOverlap->addCone(p, rank, rank);
+          }
         }
-        if (sendOverlap->debug()) {sendOverlap->view(std::cout, "Send overlap for sieve scatter");}
-        //  The arrow is from rank 0 with partition point 0
-        recvOverlap->addCone(0, recvOverlap->commRank(), recvOverlap->commRank());
-        if (recvOverlap->debug()) {recvOverlap->view(std::cout, "Receive overlap for sieve scatter");}
+        if (sendOverlap->debug) {sendOverlap->view(std::cout, "Send overlap for sieve scatter");}
+        if (recvOverlap->debug) {recvOverlap->view(std::cout, "Receive overlap for sieve scatter");}
       };
-      static void localSieve(const Obj<typename Mesh::sieve_type>& oldSieve, int numElements, short assignment[], const Obj<Mesh::sieve_type>& newSieve) {
+      #undef __FUNCT__
+      #define __FUNCT__ "localSieve"
+      static void localSieve(const Obj<typename Mesh::sieve_type>& oldSieve, const int numElements, const short assignment[], const Obj<Mesh::sieve_type>& newSieve) {
         int rank = oldSieve->commRank();
 
         for(int e = 0; e < numElements; e++) {
             }
           }
         }
-        if (overlap->debug()) {overlap->view(std::cout, "Overlap for points");}
+        if (overlap->debug) {overlap->view(std::cout, "Overlap for points");}
       };
       #undef __FUNCT__
       #define __FUNCT__ "updateSieve"
       template<typename RecvSection>
       static void updateSieve(const Obj<RecvSection>& recvSection, const Obj<topology_type>& topology) {
-        const typename RecvSection::patch_type                      patch = 0;
-        const typename RecvSection::topology_type::sheaf_type&      ranks = recvSection->getAtlas()->getTopology()->getPatches();
-        const Obj<typename RecvSection::topology_type::sieve_type>& sieve = topology->getPatch(patch);
+        const typename RecvSection::patch_type                 patch = 0;
+        const typename RecvSection::topology_type::sheaf_type& ranks = recvSection->getAtlas()->getTopology()->getPatches();
+        const Obj<typename topology_type::sieve_type>&         sieve = topology->getPatch(patch);
 
         for(typename RecvSection::topology_type::sheaf_type::const_iterator p_iter = ranks.begin(); p_iter != ranks.end(); ++p_iter) {
           int                                                              rank = p_iter->first;
       }
       #undef __FUNCT__
       #define __FUNCT__ "scatterSieve"
-      template<typename Mesh>
-      static void scatterSieve(const Obj<Mesh>& oldMesh, const Obj<Mesh>& newMesh) {
-        if (oldMesh->commSize() == 1) {
-          newMesh = oldMesh;
+      static void scatterTopology(const Obj<topology_type>& oldTopology, const int numElements, const short assignment[], Obj<topology_type>& newTopology) {
+        if (oldTopology->commSize() == 1) {
+          newTopology = oldTopology;
           return;
         }
         typedef typename Mesh::point_type value_type;
         typedef typename ALE::New::Completion<typename Mesh::topology_type,value_type> completion_type;
         typedef typename completion_type::atlas_type atlas_type;
-        Obj<send_overlap_type> sendOverlap = new send_overlap_type(oldMesh->comm(), oldMesh->debug());
-        Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(oldMesh->comm(), oldMesh->debug());
-        const typename Mesh::patch_type patch       = 0;
+        Obj<send_overlap_type> sendOverlap = new send_overlap_type(oldTopology->comm(), oldTopology->debug());
+        Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(oldTopology->comm(), oldTopology->debug());
+        const Obj<typename topology_type::sieve_type> newSieve = new typename topology_type::sieve_type(newTopology->comm(), newTopology->debug());
+        const typename topology_type::patch_type patch    = 0;
 
         typedef typename ALE::New::OverlapValues<send_overlap_type, atlas_type, value_type> send_section_type;
         typedef typename ALE::New::OverlapValues<recv_overlap_type, atlas_type, value_type> recv_section_type;
-        const Obj<send_section_type> sendSection = new send_section_type(oldMesh->comm(), oldMesh->debug());
-        const Obj<recv_section_type> recvSection = new recv_section_type(oldMesh->comm(), oldMesh->debug());
+        const Obj<send_section_type> sendSection = new send_section_type(oldTopology->comm(), oldTopology->debug());
+        const Obj<recv_section_type> recvSection = new recv_section_type(oldTopology->comm(), oldTopology->debug());
 
-        // Create Partition
-        short *assignment = ALE::New::Partitioner<typename Mesh::topology_type>::partitionSieve_Chaco(oldMesh->getTopology(), oldMesh->getDimension());
-        int numElements = oldMesh->getTopology()->heightStratum(patch, 0)->size();
         // Create Overlap - Could be calculated from partition, but would need communication
         createScatterOverlap(sendOverlap, recvOverlap);
         // Create local sieve
-        localSieve(oldMesh->getTopology()->getPatch(patch), numElements, assignment, newMesh->getTopology()->getPatch(patch));
+        newTopology->setPatch(patch, newSieve);
+        localSieve(oldTopology->getPatch(patch), numElements, assignment, newTopology->getPatch(patch));
         // Distribute points
-        Obj<topology_type>          secTopology          = completion_type::createSendTopology(sendOverlap);
+        Obj<completion_type::topology_type>                   secTopology          = completion_type::createSendTopology(sendOverlap);
         Obj<typename completion_type::partition_size_section> partitionSizeSection = new typename completion_type::partition_size_section(secTopology, numElements, assignment);
         Obj<typename completion_type::partition_section>      partitionSection     = new typename completion_type::partition_section(secTopology, numElements, assignment);
-        completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
+        sectionCompletion::completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
         // Create point overlap
-        createPointOverlap(sendSection, sendOverlap);
-        createPointOverlap(recvSection, recvOverlap);
+        createPointOverlap(sendOverlap, sendSection);
+        createPointOverlap(recvOverlap, recvSection);
         // Distribute cones
         secTopology = completion_type::createSendTopology(sendOverlap);
-        Obj<typename completion_type::cone_size_section> coneSizeSection = new typename completion_type::cone_size_section(secTopology, oldMesh->getTopology()->getPatch(patch));
-        Obj<typename completion_type::cone_section>      coneSection     = new typename completion_type::cone_section(secTopology, oldMesh->getTopology()->getPatch(patch));
-        completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
+        Obj<typename completion_type::cone_size_section> coneSizeSection = new typename completion_type::cone_size_section(secTopology, oldTopology->getPatch(patch));
+        Obj<typename completion_type::cone_section>      coneSection     = new typename completion_type::cone_section(secTopology, oldTopology->getPatch(patch));
+        sectionCompletion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
         // Update sieve
-        updateSieve(recvSection, newMesh->getTopology());
+        updateSieve(recvSection, newTopology);
       };
     };
   }

src/dm/mesh/sieve/Mesh.hh

     typedef ALE::New::Section<atlas_type, ALE::pair<int,double> > foliated_section_type;
     typedef struct {double x, y, z;}                                           split_value;
     typedef ALE::New::Section<atlas_type, ALE::pair<point_type, split_value> > split_section_type;
+    typedef ALE::New::Completion<topology_type, point_type>::send_overlap_type send_overlap_type;
+    typedef ALE::New::Completion<topology_type, point_type>::recv_overlap_type recv_overlap_type;
     int debug;
   private:
     Obj<sieve_type>            topology;
     Obj<topology_type>         _topology;
     Obj<foliated_section_type> _boundaries;
     Obj<split_section_type>    _splitField;
+    Obj<send_overlap_type>     _vertexSendOverlap;
+    Obj<recv_overlap_type>     _vertexRecvOverlap;
     MPI_Comm        _comm;
     int             _commRank;
     int             _commSize;
     void setTopologyNew(const Obj<topology_type>& topology) {this->_topology = topology;};
     const Obj<split_section_type>& getSplitSection() const {return this->_splitField;};
     void                           setSplitSection(const Obj<split_section_type>& splitField) {this->_splitField = splitField;};
+    const Obj<send_overlap_type>&  getVertexSendOverlap() const {return this->_vertexSendOverlap;};
+    void                           setVertexSendOverlap(const Obj<send_overlap_type>& vertexOverlap) {this->_vertexSendOverlap = vertexOverlap;};
+    const Obj<recv_overlap_type>&  getVertexRecvOverlap() const {return this->_vertexRecvOverlap;};
+    void                           setVertexRecvOverlap(const Obj<recv_overlap_type>& vertexOverlap) {this->_vertexRecvOverlap = vertexOverlap;};
     // Printing
     template <typename Stream_>
     friend Stream_& operator<<(Stream_& os, const split_value& v) {

src/dm/mesh/sieve/Sifter.hh

 
     template<typename ostream_type>
     void view(ostream_type& os, const char* label = NULL, bool rawData = false){
+      int rank = this->commRank();
+
       if(label != NULL) {
-        os << "Viewing Sifter '" << label << "':" << std::endl;
+        os << "["<<rank<<"]Viewing Sifter '" << label << "':" << std::endl;
       } 
       else {
-        os << "Viewing a Sifter:" << std::endl;
+        os << "["<<rank<<"]Viewing a Sifter:" << std::endl;
       }
       if(!rawData) {
         os << "cap --> base:" << std::endl;

src/docs/website/documentation/installation.html

       </p>
       <font color="#ff0000">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
 ./config/configure.py --with-cc='win32fe cl' --with-fc='win32fe f90'
---download-f-blas-lapack=1 LIBS="-L/cygdrive/c/Program Files/Microsoft
-Visual Studio/DF98/LIB"</font><br>
+--download-f-blas-lapack=1 LIBS="-L'/cygdrive/c/Program Files/Microsoft
+Visual Studio/DF98/LIB'"</font><br><br>
+Note the usage of both single quotes and double quotes in the above line. <br>
       <p><font color="#551a8b"><b>Configure</b></font><font
  color="#551a8b"><b>:</b> </font>--download-package
 option does not
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.