Merge lp:~jobh/dolfin/misc into lp:~fenics-core/dolfin/trunk

[Joachim Haga]

This branch finishes the DirichletBC speed-up effort for topological and geometric search (and, incidentally, SubDomain::apply_markers, since that is used when initialising a DirichletBC from a sub domain). These search methods (from a sub domain) are now roughly as fast as pointwise search again.

It also fixes the problem with setting BCs on a subdomain which has been reported (it checks if dofmap is a view before calling ownership_range()).

I seem to have created some problems earlier which may be due to criss-cross merges, so I haven't merged with trunk this time; however, I have performed a test-merge which was clean so there should be no problems.

It passes all tests, and the test-merge has passed the quicktests.

[Garth Wells (garth-wells)]

review: Approve

=== modified file 'demo/undocumented/sym-dirichlet-bc/python/speed-up-test.py'
--- demo/undocumented/sym-dirichlet-bc/python/speed-up-test.py	2011-10-10 00:06:19 +0000
+++ demo/undocumented/sym-dirichlet-bc/python/speed-up-test.py	2012-03-06 09:08:20 +0000
@@ -22,8 +22,8 @@
 # First added:  2008-08-13
 # Last changed: 2008-08-13
 
+from dolfin import *
 import time
-from dolfin import *
 
 
 # Sub domain for Dirichlet boundary condition
@@ -31,14 +31,14 @@
     def inside(self, x, on_boundary):
         return bool(on_boundary and x[0] < DOLFIN_EPS)
 
-mesh = UnitSquare(500,500)
+mesh = UnitCube(32,32,32)
 V = FunctionSpace(mesh, "CG", 1)
 
 # Define variational problem
 v = TestFunction(V)
 u = TrialFunction(V)
 
-f = Function(V, "500.0*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
+f = Expression("500.0*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
 
 
 a = dot(grad(v), grad(u))*dx
@@ -71,5 +71,10 @@
     t1 = time.time()
     print "time for new assembly      ", t1-t0, " using ", backend
 
+    t0 = time.time()
+    A, Aa = symmetric_assemble(a, bcs=bc)
+    b = assemble(L, bcs=bc)
+    t1 = time.time()
+    print "time for symm assembly     ", t1-t0, " using ", backend
 
 #summary()
 
=== added file 'dolfin/common/RangedIndexSet.h'
--- dolfin/common/RangedIndexSet.h	1970-01-01 00:00:00 +0000
+++ dolfin/common/RangedIndexSet.h	2012-03-06 09:08:20 +0000
@@ -0,0 +1,108 @@
+// Copyright (C) 2012 Joachim B Haga
+//
+// 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-03-02
+// Last changed: 2012-03-02
+
+#ifndef __RANGED_INDEX_SET_H
+#define __RANGED_INDEX_SET_H
+
+#include <vector>
+#include <dolfin/log/log.h>
+#include "types.h"
+
+namespace dolfin
+{
+
+  /// This class provides an special-purpose data structure for testing if a given
+  /// index within a range is set. It is very fast.
+  ///
+  /// The memory requirements are one bit per item in range, since it uses a
+  /// (packed) std::vector<bool> for storage.
+
+  class RangedIndexSet
+  {
+
+  public:
+
+    /// Create a ranged set with range given as a (lower, upper) pair.
+    RangedIndexSet(std::pair<uint, uint> range)
+      : _range(range), _is_set(range.second - range.first)
+    {
+      clear();
+    }
+
+    /// Create a ranged set with 0 as lower range
+    RangedIndexSet(uint upper_range)
+      : _range(std::pair<uint, uint>(0, upper_range)), _is_set(upper_range)
+    {
+      clear();
+    }
+
+    /// Return true if a given index is within range, i.e., if it can be stored in the set.
+    bool in_range(uint i) const
+    {
+      return (i >= _range.first && i < _range.second);
+    }
+
+    /// Check is the set contains the given index.
+    bool has_index(uint i) const
+    {
+      dolfin_assert(in_range(i));
+      return _is_set[i - _range.first];
+    }
+
+    /// Insert a given index into the set. Returns true if the index was
+    /// inserted (i.e., the index was not already in the set).
+    bool insert(uint i)
+    {
+      dolfin_assert(in_range(i));
+      std::vector<bool>::reference entry = _is_set[i - _range.first];
+      if (entry)
+      {
+        return false;
+      }
+      else
+      {
+        entry = true;
+        return true;
+      }
+    }
+
+    /// Erase an index from the set.
+    void erase(uint i)
+    {
+      dolfin_assert(in_range(i));
+      _is_set[i - _range.first] = false;
+    }
+
+    /// Erase all indices from the set.
+    void clear()
+    {
+      std::fill(_is_set.begin(), _is_set.end(), false);
+    }
+
+  private:
+
+    const std::pair<uint, uint> _range;
+    std::vector<bool> _is_set;
+
+  };
+
+}
+
+#endif
 
=== modified file 'dolfin/fem/DirichletBC.cpp'
--- dolfin/fem/DirichletBC.cpp	2012-03-02 15:40:22 +0000
+++ dolfin/fem/DirichletBC.cpp	2012-03-06 09:08:20 +0000
@@ -32,6 +32,7 @@
 #include <dolfin/common/constants.h>
 #include <dolfin/common/Array.h>
 #include <dolfin/common/NoDeleter.h>
+#include <dolfin/common/RangedIndexSet.h>
 #include <dolfin/function/GenericFunction.h>
 #include <dolfin/function/FunctionSpace.h>
 #include <dolfin/function/Constant.h>
@@ -818,6 +819,11 @@
   log(TRACE, "Computing facets, needed for geometric application of boundary conditions.");
   mesh.init(mesh.topology().dim() - 1);
 
+  // Speed up the computations by only visiting (most) dofs once
+  RangedIndexSet already_visited(dofmap.is_view()
+                                 ? std::pair<uint, uint>(0,0)
+                                 : dofmap.ownership_range());
+
   // Iterate over facets
   Progress p("Computing Dirichlet boundary values, geometric search", facets.size());
   for (uint f = 0; f < facets.size(); ++f)
@@ -841,17 +847,28 @@
       {
         ufc_cell.update(*c, facet_number);
 
+        bool tabulated = false;
         bool interpolated = false;
 
-        // Tabulate coordinates of dofs on cell
-        dofmap.tabulate_coordinates(data.coordinates, ufc_cell);
-
         // Tabulate dofs on cell
         const std::vector<uint>& cell_dofs = dofmap.cell_dofs(c->index());
 
         // Loop over all dofs on cell
-        for (uint i = 0; i < dofmap.cell_dimension(c->index()); ++i)
+        for (uint i = 0; i < cell_dofs.size(); ++i)
         {
+          const uint global_dof = cell_dofs[i];
+
+          // Skip already checked dofs
+          if (already_visited.in_range(global_dof) && !already_visited.insert(global_dof))
+            continue;
+
+          if (!tabulated)
+          {
+            // Tabulate coordinates of dofs on cell
+            dofmap.tabulate_coordinates(data.coordinates, ufc_cell);
+            tabulated = true;
+          }
+
           // Check if the coordinates are on current facet and thus on boundary
           if (!on_facet(&(data.coordinates[i][0]), facet))
             continue;
@@ -860,10 +877,10 @@
           {
             // Restrict coefficient to cell
             g->restrict(&data.w[0], *_function_space->element(), cell, ufc_cell);
+            interpolated = true;
           }
 
           // Set boundary value
-          const uint global_dof = cell_dofs[i];
           const double value = data.w[i];
           boundary_values[global_dof] = value;
         }
@@ -896,10 +913,10 @@
   // Create UFC cell object
   UFCCell ufc_cell(mesh);
 
-  // Speeder-upper
-  std::pair<uint,uint> local_range = dofmap.ownership_range();
-  std::vector<bool> already_visited(local_range.second - local_range.first);
-  std::fill(already_visited.begin(), already_visited.end(), false);
+  // Speed up the computations by only visiting (most) dofs once
+  RangedIndexSet already_visited(dofmap.is_view()
+                                 ? std::pair<uint, uint>(0,0)
+                                 : dofmap.ownership_range());
 
   // Iterate over cells
   Progress p("Computing Dirichlet boundary values, pointwise search", mesh.num_cells());
@@ -921,13 +938,10 @@
     for (uint i = 0; i < dofmap.cell_dimension(cell->index()); ++i)
     {
       const uint global_dof = cell_dofs[i];
-      if (global_dof >= local_range.first && global_dof < local_range.second)
-      {
-        const uint dof_index = global_dof - local_range.first;
-        if (already_visited[dof_index])
-          continue;
-        already_visited[dof_index] = true;
-      }
+
+      // Skip already checked dofs
+      if (already_visited.in_range(global_dof) && !already_visited.insert(global_dof))
+        continue;
 
       // Check if the coordinates are part of the sub domain (calls user-defined 'inside' function)
       Array<double> x(gdim, &data.coordinates[i][0]);
 
=== modified file 'dolfin/mesh/SubDomain.cpp'
--- dolfin/mesh/SubDomain.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/mesh/SubDomain.cpp	2012-03-06 09:08:20 +0000
@@ -21,6 +21,7 @@
 // Last changed: 2011-08-31
 
 #include <dolfin/common/Array.h>
+#include <dolfin/common/RangedIndexSet.h>
 #include <dolfin/log/log.h>
 #include "Mesh.h"
 #include "MeshData.h"
@@ -164,6 +165,13 @@
   // Array for vertex coordinate
   Array<double> x;
 
+  // Speed up the computation by only checking each vertex once (or twice if it
+  // is on the boundary for some but not all facets).
+  RangedIndexSet boundary_visited(mesh.num_vertices());
+  RangedIndexSet interior_visited(mesh.num_vertices());
+  std::vector<bool> boundary_inside(mesh.num_vertices());
+  std::vector<bool> interior_inside(mesh.num_vertices());
+
   // Compute sub domain markers
   Progress p("Computing sub domain markers", mesh.num_entities(dim));
   for (MeshEntityIterator entity(mesh, dim); !entity.end(); ++entity)
@@ -175,6 +183,10 @@
         (exterior.empty() || exterior[*entity]);
     }
 
+    // Select the visited-cache to use for this facet (or entity)
+    RangedIndexSet&    is_visited = (on_boundary ? boundary_visited : interior_visited);
+    std::vector<bool>& is_inside  = (on_boundary ? boundary_inside  : interior_inside);
+
     // Start by assuming all points are inside
     bool all_points_inside = true;
 
@@ -183,8 +195,13 @@
     {
       for (VertexIterator vertex(*entity); !vertex.end(); ++vertex)
       {
-        x.update(_geometric_dimension, const_cast<double*>(vertex->x()));
-        if (!inside(x, on_boundary))
+        if (is_visited.insert(vertex->index()))
+        {
+          x.update(_geometric_dimension, const_cast<double*>(vertex->x()));
+          is_inside[vertex->index()] = inside(x, on_boundary);
+        }
+
+        if (!is_inside[vertex->index()])
         {
           all_points_inside = false;
           break;