From f2552d7f12e7ecd443f3bb9e669ede121a1d2244 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 28 Nov 2016 17:01:02 -0700
Subject: [PATCH] Update GridTools::laplace_transform().

Incomprehensively, the function looped over all faces of a cell, and then
all vertices of the face. This visited each vertex exactly 'dim' times.
This can be done easier.

While there, also fix a number of issues where we confuse 'unsigned int'
with 'types::global_dof_index'. This doesn't matter here because the
function only works in sequential settings anyway, but we should be
consistent.

Finally, instead of setting vertex locations from each adjacent cell, only
do it once by keeping track when we touch the vertex the first time.
---
 source/grid/grid_tools.cc | 69 +++++++++++++++++++--------------------
 1 file changed, 33 insertions(+), 36 deletions(-)

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 5d667876a6..7d09959df7 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -665,9 +665,9 @@ namespace GridTools
      * of the @p dim space dimensions. Factorized into a function of its own
      * in order to allow parallel execution.
      */
-    void laplace_solve (const SparseMatrix<double> &S,
-                        const std::map<unsigned int,double> &m,
-                        Vector<double> &u)
+    void laplace_solve (const SparseMatrix<double>                     &S,
+                        const std::map<types::global_dof_index,double> &fixed_dofs,
+                        Vector<double>                                 &u)
     {
       const unsigned int n_dofs=S.n();
       FilteredMatrix<Vector<double> > SF (S);
@@ -681,7 +681,7 @@ namespace GridTools
 
       Vector<double> f(n_dofs);
 
-      SF.add_constraints(m);
+      SF.add_constraints(fixed_dofs);
       SF.apply_constraints (f, true);
       solver.solve(SF, u, f, PF);
     }
@@ -728,41 +728,34 @@ namespace GridTools
 
     MatrixCreator::create_laplace_matrix(StaticMappingQ1<dim>::mapping, dof_handler, quadrature, S, coefficient);
 
-    // set up the boundary values for
-    // the laplace problem
-    std::vector<std::map<unsigned int,double> > m(dim);
+    // set up the boundary values for the laplace problem
+    std::map<types::global_dof_index,double> fixed_dofs[dim];
     typename std::map<unsigned int,Point<dim> >::const_iterator map_end=new_points.end();
 
-    // fill these maps using the data
-    // given by new_points
+    // fill these maps using the data given by new_points
     typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin_active(),
                                             endc=dof_handler.end();
     for (; cell!=endc; ++cell)
       {
-        for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+        // loop over all vertices of the cell and see if it is listed in the map
+        // given as first argument of the function
+        for (unsigned int vertex_no=0;
+             vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
           {
-            const typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+            const unsigned int vertex_index=cell->vertex_index(vertex_no);
 
-            // loop over all vertices of the cell and see if it is listed in the map
-            // given as first argument of the function
-            for (unsigned int vertex_no=0;
-                 vertex_no<GeometryInfo<dim>::vertices_per_face; ++vertex_no)
-              {
-                const unsigned int vertex_index=face->vertex_index(vertex_no);
-
-                const typename std::map<unsigned int,Point<dim> >::const_iterator map_iter
-                  = new_points.find(vertex_index);
+            const typename std::map<unsigned int,Point<dim> >::const_iterator map_iter
+              = new_points.find(vertex_index);
 
-                if (map_iter!=map_end)
-                  for (unsigned int i=0; i<dim; ++i)
-                    m[i].insert(std::pair<unsigned int,double> (
-                                  face->vertex_dof_index(vertex_no, 0), map_iter->second(i)));
-              }
+            if (map_iter!=map_end)
+              for (unsigned int i=0; i<dim; ++i)
+                fixed_dofs[i].insert(std::pair<types::global_dof_index,double>
+                                     (cell->vertex_dof_index(vertex_no, 0),
+                                      map_iter->second(i)));
           }
       }
 
-    // solve the dim problems with
-    // different right hand sides.
+    // solve the dim problems with different right hand sides.
     Vector<double> us[dim];
     for (unsigned int i=0; i<dim; ++i)
       us[i].reinit (dof_handler.n_dofs());
@@ -771,21 +764,25 @@ namespace GridTools
     Threads::TaskGroup<> tasks;
     for (unsigned int i=0; i<dim; ++i)
       tasks += Threads::new_task (&laplace_solve,
-                                  S, m[i], us[i]);
+                                  S, fixed_dofs[i], us[i]);
     tasks.join_all ();
 
-    // change the coordinates of the
-    // points of the triangulation
+    // change the coordinates of the points of the triangulation
     // according to the computed values
+    std::vector<bool> vertex_touched (triangulation.n_vertices(), false);
     for (cell=dof_handler.begin_active(); cell!=endc; ++cell)
       for (unsigned int vertex_no=0;
            vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
-        {
-          Point<dim> &v=cell->vertex(vertex_no);
-          const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0);
-          for (unsigned int i=0; i<dim; ++i)
-            v(i)=us[i](dof_index);
-        }
+        if (vertex_touched[cell->vertex_index(vertex_no)] == false)
+          {
+            Point<dim> &v = cell->vertex(vertex_no);
+
+            const types::global_dof_index dof_index = cell->vertex_dof_index(vertex_no, 0);
+            for (unsigned int i=0; i<dim; ++i)
+              v(i) = us[i](dof_index);
+
+            vertex_touched[cell->vertex_index(vertex_no)] = true;
+          }
   }
 
   template <int dim, int spacedim>
-- 
2.39.5