Garth Wells avatar Garth Wells committed 07e4dcf Merge

Merge branch 'garth/scotch-graph-reordering'

Comments (0)

Files changed (8)

cmake/scripts/download-demo-data

             output=`wget http://fenicsproject.org/pub/data/meshes/$MESH 2>&1`
         fi
 
-	# If download fails
-	if [ $? -ne 0 ]; then
+	    # If download fails
+	    if [ $? -ne 0 ]; then
             echo "$output" 1>&2
-	    exit 1
-	fi
+	        exit 1
+	    fi
+
         mv $MESH $DEMO
     fi
 }

demo/documented/mesh-generation/python/demo_mesh_generation.py

-# Copyright (C) 2012 Garth N. Wells
-#
-# This file is part of DOLFIN.
-#
-# DOLFIN is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# DOLFIN is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
-#
-# First added:  2012-02-02
-# Last changed:
-#
-# This program generates a mesh for a polygonal domain that is
-# represented by a list of its vertices.
-# Modified 02-08-2013 Solveig Masvie
-# Begin demo
-
-
-from dolfin import *
-
-if not has_cgal():
-    print "DOLFIN must be compiled with CGAL to run this demo."
-    exit(0)
-
-# Create empty Mesh
-mesh = Mesh()
-
-# Create list of polygonal domain vertices
-domain_vertices = [Point(0.0, 0.0),
-                   Point(10.0, 0.0),
-                   Point(10.0, 2.0),
-                   Point(8.0, 2.0),
-                   Point(7.5, 1.0),
-                   Point(2.5, 1.0),
-                   Point(2.0, 4.0),
-                   Point(0.0, 4.0),
-                   Point(0.0, 0.0)]
-
-# Generate mesh and plot
-PolygonalMeshGenerator.generate(mesh, domain_vertices, 0.25);
-plot(mesh, interactive=True)
-
-# Generate 3D mesh from OFF file input (tetrahedron)
-PolyhedralMeshGenerator.generate(mesh, "../tetrahedron.off", 0.05)
-plot(mesh, interactive=True)
-
-# Generate 3D mesh from OFF file input (cube)
-PolyhedralMeshGenerator.generate(mesh, "../cube.off", 0.05)
-plot(mesh, interactive=True)

dolfin/fem/DofMapBuilder.cpp

 #include <dolfin/common/Timer.h>
 #include <dolfin/graph/BoostGraphOrdering.h>
 #include <dolfin/graph/GraphBuilder.h>
+#include <dolfin/graph/SCOTCH.h>
 #include <dolfin/log/log.h>
 #include <dolfin/mesh/BoundaryMesh.h>
 #include <dolfin/mesh/DistributedMeshTools.h>
   // Build dofmap based on UFC-provided map. This function does not
   // set local_range
   map restricted_dofs_inverse;
-  build_ufc_dofmap(dofmap, restricted_dofs_inverse, mesh, slave_master_entities,
-            restriction);
+  build_ufc_dofmap(dofmap, restricted_dofs_inverse, mesh,
+                   slave_master_entities, restriction);
 
   // Check if dofmap is distributed
   const bool distributed = MPI::num_processes() > 1;
 
   // Extract ufc sub-dofmap from parent and get offset
   dolfin_assert(parent_dofmap._ufc_dofmap);
-  sub_dofmap._ufc_dofmap = extract_ufc_sub_dofmap(*(parent_dofmap._ufc_dofmap),
-                                                  offset, component,
-                                                  parent_dofmap.num_global_mesh_entities);
+  sub_dofmap._ufc_dofmap
+    = extract_ufc_sub_dofmap(*(parent_dofmap._ufc_dofmap),
+                             offset, component,
+                             parent_dofmap.num_global_mesh_entities);
   dolfin_assert(sub_dofmap._ufc_dofmap);
 
   // Check for dimensional consistency between the dofmap and mesh
       {
         // Get dof index
         ufc_to_current_dof = parent_dofmap.ufc_map_to_dofmap.find(*dof);
-        dolfin_assert(ufc_to_current_dof != parent_dofmap.ufc_map_to_dofmap.end());
+        dolfin_assert(ufc_to_current_dof
+                      != parent_dofmap.ufc_map_to_dofmap.end());
 
         // Add to ufc-to-current dof map
         sub_dofmap.ufc_map_to_dofmap.insert(*ufc_to_current_dof);
         if (parent_shared != parent_dofmap._shared_dofs.end())
         {
           sub_dofmap._shared_dofs.insert(*parent_shared);
-          sub_dofmap._neighbours.insert(parent_shared->second.begin(), parent_shared->second.end());
+          sub_dofmap._neighbours.insert(parent_shared->second.begin(),
+                                        parent_shared->second.end());
         }
       }
     }
   std::vector<bool> vertex_shared(mesh.num_vertices(), false);
   boost::unordered_map<unsigned int, std::vector<std::pair<unsigned int,
       unsigned int> > >::const_iterator shared_vertex;
-  for (shared_vertex = shared_vertices.begin(); shared_vertex != shared_vertices.end(); ++shared_vertex)
+  for (shared_vertex = shared_vertices.begin();
+       shared_vertex != shared_vertices.end(); ++shared_vertex)
   {
     dolfin_assert(shared_vertex->first < vertex_shared.size());
     vertex_shared[shared_vertex->first] = true;
 
   // Mark slave vertices
   std::vector<bool> slave_vertex(mesh.num_vertices(), false);
-  std::map<unsigned int, std::pair<unsigned int, unsigned int> >::const_iterator slave;
-  for (slave = slave_to_master_vertices.begin(); slave != slave_to_master_vertices.end(); ++slave)
+  std::map<unsigned int, std::pair<unsigned int,
+                                   unsigned int> >::const_iterator slave;
+  for (slave = slave_to_master_vertices.begin();
+       slave != slave_to_master_vertices.end(); ++slave)
   {
     dolfin_assert(slave->first < slave_vertex.size());
     slave_vertex[slave->first] = true;
 
   // Compute modified global vertex indices
   std::size_t new_index = 0;
-  modified_global_indices = std::vector<std::size_t>(mesh.num_vertices(), std::numeric_limits<std::size_t>::max());
+  modified_global_indices
+    = std::vector<std::size_t>(mesh.num_vertices(),
+                               std::numeric_limits<std::size_t>::max());
   for (VertexIterator vertex(mesh); !vertex.end(); ++vertex)
   {
     const std::size_t local_index = vertex->index();
       boost::unordered_map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::const_iterator
         it = shared_vertices.find(local_index);
       dolfin_assert(it != shared_vertices.end());
-      const std::vector<std::pair<unsigned int, unsigned int> >& sharing_procs = it->second;
+      const std::vector<std::pair<unsigned int, unsigned int> >& sharing_procs
+        = it->second;
 
       // Figure out if this is the lowest rank process sharing the vertex
       std::vector<std::pair<unsigned int, unsigned int> >::const_iterator
-       min_sharing_rank = std::min_element(sharing_procs.begin(), sharing_procs.end());
+       min_sharing_rank = std::min_element(sharing_procs.begin(),
+                                           sharing_procs.end());
       std::size_t _min_sharing_rank = proc_num + 1;
       if (min_sharing_rank != sharing_procs.end())
         _min_sharing_rank = min_sharing_rank->first;
   // Set index for shared vertices that have been numbered by another process
   for (std::size_t p = 0; p < received_vertex_data.size(); ++p)
   {
-    const std::vector<std::size_t>& received_vertex_data_p = received_vertex_data[p];
+    const std::vector<std::size_t>& received_vertex_data_p
+      = received_vertex_data[p];
     for (std::size_t i = 0; i < received_vertex_data_p.size(); i += 2)
     {
       const unsigned int local_index = received_vertex_data_p[i];
   }
 
   // Request master vertex index from master owner
-  std::vector<std::vector<std::size_t> > master_send_buffer(MPI::num_processes());
-  std::vector<std::vector<std::size_t> > local_slave_index(MPI::num_processes());
-  std::map<unsigned int, std::pair<unsigned int, unsigned int> >::const_iterator master;
-  for (master = slave_to_master_vertices.begin(); master != slave_to_master_vertices.end(); ++master)
+  std::vector<std::vector<std::size_t> >
+    master_send_buffer(MPI::num_processes());
+  std::vector<std::vector<std::size_t> >
+    local_slave_index(MPI::num_processes());
+  std::map<unsigned int,
+           std::pair<unsigned int, unsigned int> >::const_iterator master;
+  for (master = slave_to_master_vertices.begin();
+       master != slave_to_master_vertices.end(); ++master)
   {
     const unsigned int local_index = master->first;
     const unsigned int master_proc = master->second.first;
   MPI::all_to_all(master_send_buffer, received_slave_vertex_indices);
 
   // Send back new master vertex index
-  std::vector<std::vector<std::size_t> > master_vertex_indices(MPI::num_processes());
+  std::vector<std::vector<std::size_t> >
+    master_vertex_indices(MPI::num_processes());
   for (std::size_t p = 0; p < received_slave_vertex_indices.size(); ++p)
   {
-    const std::vector<std::size_t>& local_master_indices = received_slave_vertex_indices[p];
+    const std::vector<std::size_t>& local_master_indices
+      = received_slave_vertex_indices[p];
     for (std::size_t i = 0; i < local_master_indices.size(); ++i)
     {
       std::size_t master_local_index = local_master_indices[i];
   // Set index for slave vertices
   for (std::size_t p = 0; p < received_new_slave_vertex_indices.size(); ++p)
   {
-    const std::vector<std::size_t>& new_indices = received_new_slave_vertex_indices[p];
+    const std::vector<std::size_t>& new_indices
+      = received_new_slave_vertex_indices[p];
     const std::vector<std::size_t>& local_indices = local_slave_index[p];
     for (std::size_t i = 0; i < new_indices.size(); ++i)
     {
       const std::size_t new_global_index   = new_indices[i];
 
       dolfin_assert(local_index < modified_global_indices.size());
-      //dolfin_assert(modified_global_indices[local_index] == std::numeric_limits<std::size_t>::max());
       modified_global_indices[local_index] = new_global_index;
     }
   }
 
-  // Sanity check
-  //for (std::size_t i = 0; i < modified_global_indices.size(); ++i)
-  //{
-  //  dolfin_assert(modified_global_indices[i] != std::numeric_limits<std::size_t>::max());
-  //}
-
-  // Serial hack
-  //for (slave = slave_to_master_vertices.begin(); slave != slave_to_master_vertices.end(); ++slave)
-  //{
-  //  dolfin_assert(slave->first < modified_global_indices.size());
-  //  modified_global_indices[slave->first] = modified_global_indices[slave->second.second];
-  //}
-
   // Send new indices to process that share a vertex but were not
   // responsible for re-numbering
   return MPI::sum(new_index);
   // Build local graph for blocks
   for (CellIterator cell(mesh); !cell.end(); ++cell)
   {
-    const std::vector<dolfin::la_index>& dofs0 = dofmap.cell_dofs(cell->index());
-    const std::vector<dolfin::la_index>& dofs1 = dofmap.cell_dofs(cell->index());
+    const std::vector<dolfin::la_index>& dofs0
+      = dofmap.cell_dofs(cell->index());
+    const std::vector<dolfin::la_index>& dofs1
+      = dofmap.cell_dofs(cell->index());
 
     dolfin_assert(dofs0.size() % block_size == 0);
     const std::size_t nodes_per_cell = dofs0.size()/block_size;
           graph[dofs0[i] % num_nodes].insert(dofs1[j] % num_nodes);
   }
 
-  // Reorder block graph (reverse Cuthill-McKee)
-  const std::vector<std::size_t> block_remap
-    = BoostGraphOrdering::compute_cuthill_mckee(graph, true);
+  // Reorder block graph
+  const std::string ordering_library
+    = dolfin::parameters["dof_ordering_library"];
+  std::vector<std::size_t> block_remap;
+  if (ordering_library == "Boost")
+    block_remap = BoostGraphOrdering::compute_cuthill_mckee(graph, true);
+  else if (ordering_library == "SCOTCH")
+    block_remap = SCOTCH::compute_gps(graph);
+  else if (ordering_library == "random")
+  {
+    // NOTE: Randomised dof ordering should only be used for testing/benchmarking
+    block_remap.resize(graph.size());
+    for (std::size_t i = 0; i < block_remap.size(); ++i)
+      block_remap[i] = i;
+    std::random_shuffle(block_remap.begin(), block_remap.end());
+  }
+  else
+  {
+    dolfin_error("DofMapBuilder.cpp",
+                 "reorder degrees of freedom",
+                 "The requested ordering library '%s' is unknown", ordering_library.c_str());
+  }
 
   // Re-number dofs for each cell
   std::vector<std::vector<dolfin::la_index> >::iterator cell_map;
   std::vector<dolfin::la_index>::iterator dof;
-  for (cell_map = dofmap._dofmap.begin(); cell_map != dofmap._dofmap.end(); ++cell_map)
+  for (cell_map = dofmap._dofmap.begin(); cell_map != dofmap._dofmap.end();
+       ++cell_map)
   {
     for (dof = cell_map->begin(); dof != cell_map->end(); ++dof)
     {
 
   // Sanity checks on UFC dofmap
   dolfin_assert(dofmap._ufc_dofmap);
-  dolfin_assert(dofmap._ufc_dofmap->geometric_dimension() == mesh.geometry().dim());
-  dolfin_assert(dofmap._ufc_dofmap->topological_dimension() == mesh.topology().dim());
+  dolfin_assert(dofmap._ufc_dofmap->geometric_dimension()
+                == mesh.geometry().dim());
+  dolfin_assert(dofmap._ufc_dofmap->topological_dimension()
+                == mesh.topology().dim());
 
   // Clear ufc-dofs-to-actual-dofs
   dofmap.ufc_map_to_dofmap.clear();
 
   // Global enity indices
-  std::vector<std::vector<std::size_t> > global_entity_indices(mesh.topology().dim() + 1);
+  std::vector<std::vector<std::size_t> >
+    global_entity_indices(mesh.topology().dim() + 1);
 
   // Generate and number required mesh entities. Mesh indices are modified
   // for periodic bcs
   const std::size_t D = mesh.topology().dim();
-  dofmap.num_global_mesh_entities = std::vector<std::size_t>(mesh.topology().dim() + 1, 0);
+  dofmap.num_global_mesh_entities
+    = std::vector<std::size_t>(mesh.topology().dim() + 1, 0);
   if (!slave_master_entities)
   {
     // Compute number of mesh entities
       if (dofmap._ufc_dofmap->needs_mesh_entities(d))
       {
         // Get master-slave map
-        dolfin_assert(slave_master_entities->find(d) != slave_master_entities->end());
+        dolfin_assert(slave_master_entities->find(d)
+                      != slave_master_entities->end());
         const std::map<unsigned int, std::pair<unsigned int, unsigned int> >&
           slave_to_master_entities = slave_master_entities->find(d)->second;
 
         if (d == 0)
         {
           // Compute modified global vertex indices
-          const std::size_t num_vertices = build_constrained_vertex_indices(mesh,
+          const std::size_t num_vertices
+            = build_constrained_vertex_indices(mesh,
                 slave_to_master_entities, global_entity_indices[0]);
           dofmap.num_global_mesh_entities[0] = num_vertices;
         }
           // Get number of entities
           std::map<unsigned int, std::set<unsigned int> > shared_entities;
           const std::size_t num_entities
-            = DistributedMeshTools::number_entities(mesh, slave_to_master_entities,
+            = DistributedMeshTools::number_entities(mesh,
+                                                    slave_to_master_entities,
                                                     global_entity_indices[d],
                                                     shared_entities, d);
           dofmap.num_global_mesh_entities[d] = num_entities;
       if (!global_entity_indices[d].empty())
       {
         for (std::size_t i = 0; i < ufc_cell.num_cell_entities[d]; ++i)
-          ufc_cell.entity_indices[d][i] = global_entity_indices[d][cell->entities(d)[i]];
+        {
+          ufc_cell.entity_indices[d][i]
+            = global_entity_indices[d][cell->entities(d)[i]];
+        }
       }
     }
 
     // Tabulate standard UFC dof map
     ufc_dofs.resize(local_dim);
     dofmap._ufc_dofmap->tabulate_dofs(ufc_dofs.data(),
-                                      dofmap.num_global_mesh_entities, ufc_cell);
+                                      dofmap.num_global_mesh_entities,
+                                      ufc_cell);
     std::copy(ufc_dofs.begin(), ufc_dofs.end(), cell_dofs.begin());
 
     // Renumber dofs if mesh is restricted
         continue;
 
       // Tabulate dofs on cell
-      const std::vector<dolfin::la_index>& cell_dofs = dofmap.cell_dofs(c.index());
+      const std::vector<dolfin::la_index>& cell_dofs
+        = dofmap.cell_dofs(c.index());
 
       // Tabulate which dofs are on the facet
       dofmap.tabulate_facet_dofs(facet_dofs, c.index(f));
   std::vector<std::size_t> recv_buffer;
   for (std::size_t k = 1; k < num_prococesses; ++k)
   {
-    const std::size_t src  = (process_number - k + num_prococesses) % num_prococesses;
+    const std::size_t src
+      = (process_number - k + num_prococesses) % num_prococesses;
     const std::size_t dest = (process_number + k) % num_prococesses;
     MPI::send_recv(send_buffer, dest, recv_buffer, src);
 
           shared_unowned_nodes.insert(received_node);
           shared_owned_nodes.erase(received_node);
         }
-        else if (received_vote == node_vote[received_node] && process_number > src)
+        else if (received_vote == node_vote[received_node]
+                 && process_number > src)
         {
           // If votes are equal, let lower rank process take ownership
           shared_unowned_nodes.insert(received_node);
         // Remember the sharing of the node
         shared_node_processes[received_node].push_back(src);
       }
-      else if (shared_unowned_nodes.find(received_node) != shared_unowned_nodes.end())
+      else if (shared_unowned_nodes.find(received_node)
+               != shared_unowned_nodes.end())
       {
         // Remember the sharing of the node
         shared_node_processes[received_node].push_back(src);
   if (process_number == 0)
   {
     shared_owned_nodes.insert(global_dofs.begin(), global_dofs.end());
-    for (set::const_iterator dof = global_dofs.begin(); dof != global_dofs.end(); ++dof)
+    for (set::const_iterator dof = global_dofs.begin();
+         dof != global_dofs.end(); ++dof)
     {
       set::const_iterator _dof = shared_unowned_nodes.find(*dof);
       if (_dof != shared_unowned_nodes.end())
   else
   {
     shared_unowned_nodes.insert(global_dofs.begin(), global_dofs.end());
-    for (set::const_iterator dof = global_dofs.begin(); dof != global_dofs.end(); ++dof)
+    for (set::const_iterator dof = global_dofs.begin();
+         dof != global_dofs.end(); ++dof)
     {
       set::const_iterator _dof = shared_owned_nodes.find(*dof);
       if (_dof != shared_owned_nodes.end())
   // Mark all shared-and-owned dofs as owned by the processes
   for (CellIterator cell(mesh); !cell.end(); ++cell)
   {
-    const std::vector<dolfin::la_index>& cell_dofs = dofmap.cell_dofs(cell->index());
-    //const std::size_t cell_dimension = dofmap.cell_dimension(cell->index());
+    const std::vector<dolfin::la_index>& cell_dofs
+      = dofmap.cell_dofs(cell->index());
     for (std::size_t i = 0; i < cell_dofs.size(); ++i)
     {
       // Get cell node
   log(TRACE, "Finished determining dof ownership for parallel dof map");
 }
 //-----------------------------------------------------------------------------
-void DofMapBuilder::parallel_renumber(const boost::array<set, 3>& node_ownership,
-                                      const vec_map& shared_node_processes,
-                                      DofMap& dofmap,
-                                      const Mesh& mesh,
-                                      boost::shared_ptr<const Restriction> restriction,
-                                      const map& restricted_nodes_inverse,
-                                      std::size_t block_size)
+void
+DofMapBuilder::parallel_renumber(const boost::array<set, 3>& node_ownership,
+                                 const vec_map& shared_node_processes,
+                                 DofMap& dofmap,
+                                 const Mesh& mesh,
+                                 boost::shared_ptr<const Restriction> restriction,
+                                 const map& restricted_nodes_inverse,
+                                 std::size_t block_size)
 {
   log(TRACE, "Renumber dofs for parallel dof map");
 
                  "The degree of freedom mapping cannot be renumbered twice");
   }
 
-  const std::vector<std::vector<dolfin::la_index> >& old_dofmap = dofmap._dofmap;
+  const std::vector<std::vector<dolfin::la_index> >& old_dofmap
+    = dofmap._dofmap;
   std::vector<std::vector<dolfin::la_index> > new_dofmap(old_dofmap.size());
   dolfin_assert(old_dofmap.size() == mesh.num_cells());
 
   // Compute offset for owned and non-shared nodes
-  const std::size_t process_offset = MPI::global_offset(owned_nodes.size(), true);
+  const std::size_t process_offset = MPI::global_offset(owned_nodes.size(),
+                                                        true);
 
   // Clear some data
   dofmap._off_process_owner.clear();
   }
 
   // Reorder nodes locally
-  const std::vector<std::size_t> node_remap
-      = BoostGraphOrdering::compute_cuthill_mckee(graph, true);
+  // Reorder block graph
+  const std::string ordering_library
+    = dolfin::parameters["dof_ordering_library"];
+  std::vector<std::size_t> node_remap;
+  if (ordering_library == "Boost")
+    node_remap = BoostGraphOrdering::compute_cuthill_mckee(graph, true);
+  else if (ordering_library == "SCOTCH")
+   node_remap = SCOTCH::compute_gps(graph);
+  else if (ordering_library == "random")
+  {
+    // NOTE: Randomised dof ordering should only be used for testing/benchmarking
+    node_remap.resize(graph.size());
+    for (std::size_t i = 0; i < node_remap.size(); ++i)
+      node_remap[i] = i;
+    std::random_shuffle(node_remap.begin(), node_remap.end());
+  }
+  else
+  {
+    dolfin_error("DofMapBuilder.cpp",
+                 "reorder degrees of freedom",
+                 "The requested ordering library '%s' is unknown", ordering_library.c_str());
+  }
 
   // Map from old to new index for dofs
   boost::unordered_map<std::size_t, std::size_t> old_to_new_node_index;
     for (std::size_t i = 0; i < block_size; ++i)
     {
       std::size_t ufc_dof_index = *owned_node + i*num_nodes;
-      std::size_t new_dof_index = (process_offset + node_remap[counter])*block_size + i;
+      std::size_t new_dof_index = (process_offset
+                                   + node_remap[counter])*block_size + i;
 
       dofmap.ufc_map_to_dofmap[ufc_dof_index] = new_dof_index;
     }
   std::vector<std::size_t> recv_buffer;
   for (std::size_t k = 1; k < num_processes; ++k)
   {
-    const std::size_t src  = (process_number - k + num_processes) % num_processes;
+    const std::size_t src
+      = (process_number - k + num_processes) % num_processes;
     const std::size_t dest = (process_number + k) % num_processes;
     MPI::send_recv(send_buffer, dest, recv_buffer, src);
 
       const std::size_t received_new_node_index = recv_buffer[i + 1];
 
       // Check if this process has shared dof (and is not the owner)
-      if (shared_unowned_nodes.find(received_old_node_index) != shared_unowned_nodes.end())
+      if (shared_unowned_nodes.find(received_old_node_index)
+          != shared_unowned_nodes.end())
       {
         // Add to old-to-new node map
-        old_to_new_node_index[received_old_node_index] = received_new_node_index;
+        old_to_new_node_index[received_old_node_index]
+          = received_new_node_index;
 
         // Store map from off-process dof to owner and update
         // UFC-to-renumbered map
           dofmap._off_process_owner[new_dof_index] = src;
 
           // Update UFC-to-renumbered map
-          //const std::size_t old_dof_index = received_new_node_index*block_size + i;
           dofmap.ufc_map_to_dofmap[ufc_dof_index] = new_dof_index;
         }
-        //dofmap._off_process_owner[received_new_node_index] = src;
-
-        // Update UFC-to-renumbered map
-        //dofmap.ufc_map_to_dofmap[received_old_node_index] = received_new_node_index;
       }
     }
   }
         const std::size_t dof = it->first*block_size + i;
         dofmap._shared_dofs.insert(std::make_pair(dof, it->second));
       }
-      //dofmap._shared_dofs.insert(*it);
     }
     else
     {
         const std::size_t dof = (new_index->second)*block_size + i;
         dofmap._shared_dofs.insert(std::make_pair(dof, it->second));
       }
-      //dofmap._shared_dofs.insert(std::make_pair(new_index->second, it->second));
     }
     dofmap._neighbours.insert(it->second.begin(), it->second.end());
   }
     new_dofmap[cell_index].resize(cell_dimension);
     for (std::size_t i = 0; i < cell_dimension; ++i)
     {
-      // Get old dof
-      //std::size_t old_index = old_dofmap[cell_index][i];
-      //std::size_t old_node = old_dofmap[cell_index][i] % num_nodes;
-
-      // Map back to original (and common) numbering for restricted space
-      //if (restriction)
-      //{
-      //  const map_iterator it = restricted_nodes_inverse.find(old_index);
-      //  dolfin_assert(it != restricted_nodes_inverse.end());
-      //  old_index = it->second;
-      //}
-
-      // Insert dof
-      //new_dofmap[cell_index][i] = old_to_new_node_index[old_index];
-      //new_dofmap[cell_index][i] = old_to_new_node_index[old_index];
-
       const std::size_t old_index = old_dofmap[cell_index][i];
       const std::size_t old_node  = old_index % num_nodes;
       const std::size_t new_node  = old_to_new_node_index[old_node];
       // Create dummy cell argument to tabulate single global dof
       boost::scoped_ptr<ufc::cell> ufc_cell(new ufc::cell);
       std::size_t dof = 0;
-      ufc_dofmap->tabulate_dofs(&dof, dofmap.num_global_mesh_entities, *ufc_cell);
+      ufc_dofmap->tabulate_dofs(&dof, dofmap.num_global_mesh_entities,
+                                *ufc_cell);
 
       // Insert global dof index
       std::pair<DofMapBuilder::set::iterator, bool> ret
     for (std::size_t i = 0; i < ufc_dofmap->num_sub_dofmaps(); ++i)
     {
       // Extract sub-dofmap and intialise
-      boost::shared_ptr<ufc::dofmap> sub_dofmap(ufc_dofmap->create_sub_dofmap(i));
+      boost::shared_ptr<ufc::dofmap>
+        sub_dofmap(ufc_dofmap->create_sub_dofmap(i));
 
       compute_global_dofs(global_dofs, offset, sub_dofmap, dofmap);
 
   for (std::size_t i = 0; i < component[0]; i++)
   {
     // Extract sub dofmap
-    boost::scoped_ptr<ufc::dofmap> ufc_tmp_dofmap(ufc_dofmap.create_sub_dofmap(i));
+    boost::scoped_ptr<ufc::dofmap>
+      ufc_tmp_dofmap(ufc_dofmap.create_sub_dofmap(i));
     dolfin_assert(ufc_tmp_dofmap);
 
     // Check dimensional consistency between UFC dofmap and the mesh
   }
 
   // Create UFC sub-system
-  boost::shared_ptr<ufc::dofmap> sub_dofmap(ufc_dofmap.create_sub_dofmap(component[0]));
+  boost::shared_ptr<ufc::dofmap>
+    sub_dofmap(ufc_dofmap.create_sub_dofmap(component[0]));
   dolfin_assert(sub_dofmap);
 
   // Return sub-system if sub-sub-system should not be extracted,
   if (ufc_dofmap.num_sub_dofmaps() > 1)
   {
     // Create UFC first sub-dofmap
-    boost::scoped_ptr<ufc::dofmap> ufc_sub_dofmap0(ufc_dofmap.create_sub_dofmap(0));
+    boost::scoped_ptr<ufc::dofmap>
+      ufc_sub_dofmap0(ufc_dofmap.create_sub_dofmap(0));
     dolfin_assert(ufc_sub_dofmap0);
 
     // Create UFC sub-dofmaps and check that all sub dofmaps have the
       // same number of dofs per entity
       for (std::size_t i = 1; i < ufc_dofmap.num_sub_dofmaps(); ++i)
       {
-        boost::scoped_ptr<ufc::dofmap> ufc_sub_dofmap(ufc_dofmap.create_sub_dofmap(i));
+        boost::scoped_ptr<ufc::dofmap>
+          ufc_sub_dofmap(ufc_dofmap.create_sub_dofmap(i));
         dolfin_assert(ufc_sub_dofmap);
         for (std::size_t d = 0; d <= ufc_dofmap.topological_dimension(); ++d)
         {
-          if (ufc_sub_dofmap->num_entity_dofs(d) != ufc_sub_dofmap0->num_entity_dofs(d))
+          if (ufc_sub_dofmap->num_entity_dofs(d)
+              != ufc_sub_dofmap0->num_entity_dofs(d))
           {
             has_block_structure = false;
             break;

dolfin/graph/BoostGraphOrdering.cpp

 #include <boost/graph/cuthill_mckee_ordering.hpp>
 #include <boost/graph/properties.hpp>
 
+#include <dolfin/common/Timer.h>
 #include "Graph.h"
 #include "BoostGraphOrdering.h"
 
 std::vector<std::size_t>
   BoostGraphOrdering::compute_cuthill_mckee(const Graph& graph, bool reverse)
 {
+  Timer timer("Boost Cuthill-McKee graph ordering");
+
   // Number of vertices
   const std::size_t n = graph.size();
 

dolfin/graph/SCOTCH.cpp

 #include <map>
 #include <numeric>
 #include <set>
+#include <boost/lexical_cast.hpp>
 
 #include <dolfin/common/Set.h>
 #include <dolfin/common/Timer.h>
 #include "GraphBuilder.h"
 #include "SCOTCH.h"
 
+#include <dolfin/common/dolfin_common.h>
+
 #ifdef HAS_SCOTCH
 extern "C"
 {
 
   // Compute partitions
   const std::size_t num_global_vertices = mesh_data.num_global_cells;
-  const std::vector<std::size_t>& global_cell_indices = mesh_data.global_cell_indices;
+  const std::vector<std::size_t>& global_cell_indices
+    = mesh_data.global_cell_indices;
   partition(local_graph, ghost_vertices, global_cell_indices,
             num_global_vertices, cell_partition);
 }
 //-----------------------------------------------------------------------------
-std::vector<std::size_t> SCOTCH::compute_reordering(const Graph& graph)
+std::vector<std::size_t> SCOTCH::compute_gps(const Graph& graph,
+                                             std::size_t num_passes)
+{
+  // Create strategy string for Gibbs-Poole-Stockmayer ordering
+  std::string strategy = "g{pass= "
+    + boost::lexical_cast<std::string>(num_passes) + "}";
+
+  return compute_reordering(graph, strategy);
+}
+//-----------------------------------------------------------------------------
+std::vector<std::size_t> SCOTCH::compute_reordering(const Graph& graph,
+                                                    std::string scotch_strategy)
 {
   std::vector<std::size_t> permutation, inverse_permutation;
-  compute_reordering(graph, permutation, inverse_permutation);
+  compute_reordering(graph, permutation, inverse_permutation, scotch_strategy);
   return permutation;
 }
 //-----------------------------------------------------------------------------
 void SCOTCH::compute_reordering(const Graph& graph,
-                               std::vector<std::size_t>& permutation,
-                               std::vector<std::size_t>& inverse_permutation)
+                                std::vector<std::size_t>& permutation,
+                                std::vector<std::size_t>& inverse_permutation,
+                                std::string scotch_strategy)
 {
-  // Remove graph loops
-  Graph _graph = graph;
-  for(std::size_t i = 0; i < _graph.size(); ++i)
-  {
-    _graph[i].set().erase(std::remove(_graph[i].set().begin(), _graph[i].set().end(), i), _graph[i].set().end());
-  }
+  Timer timer("SCOTCH graph ordering");
 
   // Number of local graph vertices (cells)
-  const SCOTCH_Num vertnbr = _graph.size();
+  const SCOTCH_Num vertnbr = graph.size();
 
-  // Data structures for graph input to SCOTCH (add 1 for case that graph size is zero)
+  // Data structures for graph input to SCOTCH (add 1 for case that
+  // graph size is zero)
   std::vector<SCOTCH_Num> verttab;
   verttab.reserve(vertnbr + 1);
   std::vector<SCOTCH_Num> edgetab;
-  edgetab.reserve(10*vertnbr);
+  edgetab.reserve(20*vertnbr);
 
   // Build local graph input for SCOTCH
   // (number of local + ghost graph vertices (cells),
   SCOTCH_Num edgenbr = 0;
   verttab.push_back(0);
   Graph::const_iterator vertex;
-  for(vertex = _graph.begin(); vertex != _graph.end(); ++vertex)
+  for(vertex = graph.begin(); vertex != graph.end(); ++vertex)
   {
     edgenbr += vertex->size();
     verttab.push_back(verttab.back() + vertex->size());
   }
 
   // Check graph data for consistency
+  /*
   #ifdef DEBUG
   if (SCOTCH_graphCheck(&scotch_graph))
   {
                  "Consistency error in SCOTCH graph");
   }
   #endif
+  */
 
   // Re-ordering strategy
   SCOTCH_Strat strat;
   SCOTCH_stratInit(&strat);
 
+  // Set SCOTCH strategy (if provided)
+  //SCOTCH_stratGraphOrderBuild(&strat, SCOTCH_STRATQUALITY, 0, 0);
+  //SCOTCH_stratGraphOrderBuild(&strat, SCOTCH_STRATSPEED, 0, 0);
+  if (!scotch_strategy.empty())
+    SCOTCH_stratGraphOrder(&strat, scotch_strategy.c_str());
+
   // Vector to hold permutation vectors
   std::vector<SCOTCH_Num> permutation_indices(vertnbr);
   std::vector<SCOTCH_Num> inverse_permutation_indices(vertnbr);
 
-  // Reset SCOTCH random number generator to produce deterministic partitions
+  // Reset SCOTCH random number generator to produce deterministic
+  // partitions
   SCOTCH_randomReset();
 
   // Compute re-ordering
   // Number of local graph vertices (cells)
   const SCOTCH_Num vertlocnbr = local_graph.size();
 
-  // Data structures for graph input to SCOTCH (add 1 for case that local graph size is zero)
+  // Data structures for graph input to SCOTCH (add 1 for case that
+  // local graph size is zero)
   std::vector<SCOTCH_Num> vertloctab;
   vertloctab.reserve(local_graph.size() + 1);
   std::vector<SCOTCH_Num> edgeloctab;
   // Array containing . . . . (some sanity checks)
   std::vector<std::size_t> procvrttab(num_processes + 1);
   for (std::size_t i = 0; i < num_processes; ++i)
-    procvrttab[i] = std::accumulate(proccnttab.begin(), proccnttab.begin() + i, 0);
-  procvrttab[num_processes] = procvrttab[num_processes - 1] + proccnttab[num_processes - 1];
+  {
+    procvrttab[i] = std::accumulate(proccnttab.begin(),
+                                    proccnttab.begin() + i, 0);
+  }
+  procvrttab[num_processes] = procvrttab[num_processes - 1]
+    + proccnttab[num_processes - 1];
 
   // Sanity check
   for (std::size_t i = 1; i <= proc_num; ++i)
 
   // Copy partition
   cell_partition.resize(_cell_partition.size());
-  std::copy(_cell_partition.begin(), _cell_partition.end(), cell_partition.begin());
+  std::copy(_cell_partition.begin(), _cell_partition.end(),
+            cell_partition.begin());
 }
 //-----------------------------------------------------------------------------
 #else
                "DOLFIN has been configured without support for SCOTCH");
 }
 //-----------------------------------------------------------------------------
-void SCOTCH::partition(const std::vector<std::set<std::size_t> >& local_graph,
-                       const std::set<std::size_t>& ghost_vertices,
-                       const std::vector<std::size_t>& global_cell_indices,
-                       std::size_t num_global_vertices,
-                       std::vector<std::size_t>& cell_partition)
+std::vector<std::size_t> SCOTCH::compute_gps(const Graph& graph,
+                                             std::size_t num_passes)
 {
   dolfin_error("SCOTCH.cpp",
-               "partition mesh using SCOTCH",
+               "re-order graph using SCOTCH",
                "DOLFIN has been configured without support for SCOTCH");
+  return std::vector<std::size_t>();
 }
 //-----------------------------------------------------------------------------
-std::vector<std::size_t> SCOTCH::compute_reordering(const Graph& graph)
+std::vector<std::size_t>
+SCOTCH::compute_reordering(const Graph& graph,
+                           std::string scotch_strategy)
 {
   dolfin_error("SCOTCH.cpp",
                "re-order graph using SCOTCH",
                "DOLFIN has been configured without support for SCOTCH");
-  std::vector<std::size_t> x;
-  return x;
+  return std::vector<std::size_t>();
 }
 //-----------------------------------------------------------------------------
 void SCOTCH::compute_reordering(const Graph& graph,
-                                 std::vector<std::size_t>& permutation,
-                                 std::vector<std::size_t>& inverse_permutation)
+                                std::vector<std::size_t>& permutation,
+                                std::vector<std::size_t>& inverse_permutation,
+                                std::string scotch_strategy)
+
 {
   dolfin_error("SCOTCH.cpp",
                "re-order graph using SCOTCH",
                "DOLFIN has been configured without support for SCOTCH");
 }
 //-----------------------------------------------------------------------------
+void SCOTCH::partition(const std::vector<std::set<std::size_t> >& local_graph,
+                       const std::set<std::size_t>& ghost_vertices,
+                       const std::vector<std::size_t>& global_cell_indices,
+                       std::size_t num_global_vertices,
+                       std::vector<std::size_t>& cell_partition)
+{
+  dolfin_error("SCOTCH.cpp",
+               "partition mesh using SCOTCH",
+               "DOLFIN has been configured without support for SCOTCH");
+}
+//-----------------------------------------------------------------------------
 
 #endif

dolfin/graph/SCOTCH.h

 
 #include <cstddef>
 #include <set>
+#include <string>
 #include <vector>
 
 #include "Graph.h"
     static void compute_partition(std::vector<std::size_t>& cell_partition,
                                   const LocalMeshData& mesh_data);
 
+    /// Compute reordering (map[old] -> new) using
+    /// Gibbs-Poole-Stockmeyer re-ordering
+    static std::vector<std::size_t> compute_gps(const Graph& graph,
+                                                std::size_t num_passes=5);
+
     // Compute graph re-ordering
-    static std::vector<std::size_t> compute_reordering(const Graph& graph);
+    static std::vector<std::size_t> compute_reordering(const Graph& graph,
+                                                       std::string scotch_strategy="");
 
     // Compute graph re-ordering
-    static void compute_reordering(const Graph& graph,
-                                   std::vector<std::size_t>& permutation,
-                                   std::vector<std::size_t>& inverse_permutation);
+    static
+      void compute_reordering(const Graph& graph,
+                              std::vector<std::size_t>& permutation,
+                              std::vector<std::size_t>& inverse_permutation,
+                              std::string scotch_strategy="");
 
   private:
 
     // Compute cell partitions from distribted dual graph
-    static void partition(const std::vector<std::set<std::size_t> >& local_graph,
-                          const std::set<std::size_t>& ghost_vertices,
-                          const std::vector<std::size_t>& global_cell_indices,
-                          const std::size_t num_global_vertices,
-                          std::vector<std::size_t>& cell_partition);
+    static
+      void partition(const std::vector<std::set<std::size_t> >& local_graph,
+                     const std::set<std::size_t>& ghost_vertices,
+                     const std::vector<std::size_t>& global_cell_indices,
+                     const std::size_t num_global_vertices,
+                     std::vector<std::size_t>& cell_partition);
 
   };
 

dolfin/la/PETScKrylovSolver.cpp

 #ifdef HAS_PETSC
 
 #include <petsclog.h>
-
 #include <boost/assign/list_of.hpp>
-#include <dolfin/common/NoDeleter.h>
-#include <dolfin/log/dolfin_log.h>
+
 #include <dolfin/common/MPI.h>
+#include <dolfin/common/NoDeleter.h>
+#include <dolfin/common/Timer.h>
 #include "GenericMatrix.h"
 #include "GenericVector.h"
 #include "KrylovSolver.h"
 //-----------------------------------------------------------------------------
 std::size_t PETScKrylovSolver::solve(PETScVector& x, const PETScVector& b)
 {
+  Timer timer("PETSc Krylov solver");
+
   dolfin_assert(_A);
   dolfin_assert(_ksp);
 

dolfin/parameter/GlobalParameters.h

       // Line width relative to edge length in SVG output
       p.add("relative_line_width", 0.025);
 
-
-      // Threaded computation
-
       // Number of threads to run, 0 = run serial version
       p.add("num_threads", 0);
 
       // DOF reordering when running in serial
       p.add("reorder_dofs_serial", true);
 
+      // Allowed dofs ordering libraries and set default
+      std::set<std::string> allowed_dof_ordering_libraries;
+      allowed_dof_ordering_libraries.insert("Boost");
+      allowed_dof_ordering_libraries.insert("random");
+      allowed_dof_ordering_libraries.insert("SCOTCH");
+      std::string default_dof_ordering_library = "Boost";
+      #ifdef HAS_SCOTCH
+      default_dof_ordering_library = "SCOTCH";
+      #endif
+
+      // Add dof ordering library
+      p.add("dof_ordering_library", default_dof_ordering_library,
+            allowed_dof_ordering_libraries);
+
       // Print the level of thread support provided by the MPI library
       p.add("print_mpi_thread_support_level", false);
 
       #endif
 
       // Add mesh/graph partitioner
-      p.add("mesh_partitioner",
-            default_mesh_partitioner,
+      p.add("mesh_partitioner", default_mesh_partitioner,
             allowed_mesh_partitioners);
 
       // Approaches to partitioning (following Zoltan syntax)
             allowed_partitioning_approaches);
 
       #ifdef HAS_PARMETIS
-      // Repartitioning parameter, determines how strongly to hold on to cells
-      // when shifting between processes
+      // Repartitioning parameter, determines how strongly to hold on
+      // to cells when shifting between processes
       p.add("ParMETIS_repartitioning_weight", 1000.0);
       #endif
 
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.