From: Wolfgang Bangerth Date: Tue, 29 Nov 2016 00:01:02 +0000 (-0700) Subject: Update GridTools::laplace_transform(). X-Git-Tag: v8.5.0-rc1~353^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F3645%2Fhead;p=dealii.git 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. --- 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 &S, - const std::map &m, - Vector &u) + void laplace_solve (const SparseMatrix &S, + const std::map &fixed_dofs, + Vector &u) { const unsigned int n_dofs=S.n(); FilteredMatrix > SF (S); @@ -681,7 +681,7 @@ namespace GridTools Vector 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::mapping, dof_handler, quadrature, S, coefficient); - // set up the boundary values for - // the laplace problem - std::vector > m(dim); + // set up the boundary values for the laplace problem + std::map fixed_dofs[dim]; typename std::map >::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::cell_iterator cell=dof_handler.begin_active(), endc=dof_handler.end(); for (; cell!=endc; ++cell) { - for (unsigned int face_no=0; face_no::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::vertices_per_cell; ++vertex_no) { - const typename DoFHandler::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::vertices_per_face; ++vertex_no) - { - const unsigned int vertex_index=face->vertex_index(vertex_no); - - const typename std::map >::const_iterator map_iter - = new_points.find(vertex_index); + const typename std::map >::const_iterator map_iter + = new_points.find(vertex_index); - if (map_iter!=map_end) - for (unsigned int i=0; i ( - face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); - } + if (map_iter!=map_end) + for (unsigned int i=0; i + (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 us[dim]; for (unsigned int i=0; i tasks; for (unsigned int i=0; i vertex_touched (triangulation.n_vertices(), false); for (cell=dof_handler.begin_active(); cell!=endc; ++cell) for (unsigned int vertex_no=0; vertex_no::vertices_per_cell; ++vertex_no) - { - Point &v=cell->vertex(vertex_no); - const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0); - for (unsigned int i=0; ivertex_index(vertex_no)] == false) + { + Point &v = cell->vertex(vertex_no); + + const types::global_dof_index dof_index = cell->vertex_dof_index(vertex_no, 0); + for (unsigned int i=0; ivertex_index(vertex_no)] = true; + } } template