Commits

Anders Logg committed 8b273e1

Some improvements in bounding box code find while debugging

  • Participants
  • Parent commits 61b0565

Comments (0)

Files changed (8)

dolfin/function/Function.cpp

 // Modified by Andre Massing 2009
 //
 // First added:  2003-11-28
-// Last changed: 2013-09-25
+// Last changed: 2013-10-25
 
 #include <algorithm>
 #include <map>
   // Check that this is not a sub function.
   if (_vector->size() != _function_space->dofmap()->global_dimension())
   {
-    cout << "Size of vector: " << _vector->size() << endl;
-    cout << "Size of function space: "
-         << _function_space->dofmap()->global_dimension() << endl;
     dolfin_error("Function.cpp",
                  "access vector of degrees of freedom",
                  "Cannot access a non-const vector from a subfunction");
     }
     else
     {
+      cout << point << endl;
       dolfin_error("Function.cpp",
                    "evaluate function at point",
                    "The point is not inside the domain. Consider setting \"allow_extrapolation\" to allow extrapolation");
   const Mesh& mesh = *_function_space->mesh();
 
   const double* _x = x.data();
-  const std::size_t dim = mesh.geometry().dim();
-  const Point point(dim, _x);
+  const std::size_t gdim = mesh.geometry().dim();
+  const Point point(gdim, _x);
 
   // FIXME: Testing
   int ID = 0;
   // Check whether we are allowed to extrapolate to evaluate
   if (id == std::numeric_limits<unsigned int>::max() && !allow_extrapolation)
   {
+    cout << point << endl;
     dolfin_error("Function.cpp",
                  "evaluate function at point",
                  "The point is not inside the domain. Consider setting \"allow_extrapolation\" to allow extrapolation");
 
   // Alternative 2: Compute closest cell to point (x)
   if (id == std::numeric_limits<unsigned int>::max()
-      && allow_extrapolation && dim == 2)
+      && allow_extrapolation && gdim == 2)
   {
       id = mesh.bounding_box_tree()->compute_closest_entity(point).first;
   }
     const double * const * vertices = ufc_cell.coordinates;
 
     Point barycenter;
-    for (std::size_t i = 0; i <= dim; i++)
+    for (std::size_t i = 0; i <= gdim; i++)
     {
-      Point vertex(dim, vertices[i]);
+      Point vertex(gdim, vertices[i]);
       barycenter += vertex;
     }
-    barycenter /= (dim + 1);
+    barycenter /= (gdim + 1);
     id = mesh.bounding_box_tree()->compute_first_entity_collision(barycenter);
   }
 

dolfin/geometry/BoundingBoxTree1D.h

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-05-02
-// Last changed: 2013-09-12
+// Last changed: 2013-10-24
 
 #ifndef __BOUNDING_BOX_TREE_1D_H
 #define __BOUNDING_BOX_TREE_1D_H
       bbox[1] = b[1];
 
       // Compute min and max over remaining boxes
-      for (; it != end; ++it)
+      for (++it; it != end; ++it)
       {
         const double* b = leaf_bboxes.data() + 2*(*it);
         if (b[0] < bbox[0]) bbox[0] = b[0];

dolfin/geometry/BoundingBoxTree2D.h

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-05-02
-// Last changed: 2013-09-12
+// Last changed: 2013-10-25
 
 #ifndef __BOUNDING_BOX_TREE_2D_H
 #define __BOUNDING_BOX_TREE_2D_H
 
+
+#include <boost/range/algorithm/nth_element.hpp>
+
 #include <algorithm>
 #include <vector>
 #include <dolfin/common/constants.h>
       bbox[3] = b[3];
 
       // Compute min and max over remaining boxes
-      for (; it != end; ++it)
+      for (++it; it != end; ++it)
       {
         const double* b = leaf_bboxes.data() + 4*(*it);
         if (b[0] < bbox[0]) bbox[0] = b[0];

dolfin/geometry/BoundingBoxTree3D.h

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-04-09
-// Last changed: 2013-09-12
+// Last changed: 2013-10-24
 
 #ifndef __BOUNDING_BOX_TREE_3D_H
 #define __BOUNDING_BOX_TREE_3D_H
       bbox[5] = p[2];
 
       // Compute min and max over remaining points
-      for (; it != end; ++it)
+      for (++it; it != end; ++it)
       {
         const double* p = points[*it].coordinates();
         if (p[0] < bbox[0]) bbox[0] = p[0];

dolfin/geometry/GenericBoundingBoxTree.cpp

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-05-02
-// Last changed: 2013-09-12
+// Last changed: 2013-10-25
 
 // Define a maximum dimension used for a local array in the recursive
 // build function. Speeds things up compared to allocating it in each

dolfin/geometry/GenericBoundingBoxTree.h

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-04-23
-// Last changed: 2013-09-12
+// Last changed: 2013-10-25
 
 #ifndef __GENERIC_BOUNDING_BOX_TREE_H
 #define __GENERIC_BOUNDING_BOX_TREE_H

dolfin/mesh/TriangleCell.cpp

 // Modified by Jan Blechta 2013
 //
 // First added:  2006-06-05
-// Last changed: 2013-09-12
+// Last changed: 2013-10-25
 
 #include <algorithm>
 #include <dolfin/log/log.h>

test/unit/geometry/python/BoundingBoxTree.py

 # along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 #
 # First added:  2013-04-15
-# Last changed: 2013-10-22
+# Last changed: 2013-10-23
 
 import unittest
 import numpy
     print ""
     print "Testing BoundingBoxTree"
     print "------------------------------------------------"
-
-    print "FIXME: Temporary testing"
     unittest.main()