--- Begin Message ---
Merge authors:
Joachim Haga (jobh)
Related merge proposals:
https://code.launchpad.net/~jobh/dolfin/misc/+merge/96079
proposed by: Joachim Haga (jobh)
------------------------------------------------------------
revno: 6626 [merge]
committer: Garth N. Wells <gnw20@xxxxxxxxx>
branch nick: staging
timestamp: Tue 2012-03-06 13:55:20 +0000
message:
merge bc changes from Joachim.
added:
dolfin/common/RangedIndexSet.h
modified:
demo/undocumented/sym-dirichlet-bc/python/speed-up-test.py
dolfin/fem/DirichletBC.cpp
dolfin/mesh/SubDomain.cpp
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/trunk
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/trunk/+edit-subscription
=== 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-05 15:18:32 +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-05 15:18:32 +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 13:55: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 13:55: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;
--- End Message ---
