From 78a562f74f60312f9828bdbdf4358d8802479c3f Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 6 Aug 2008 17:15:08 +0000 Subject: [PATCH] Implement a second function for partitioning a triangulation. git-svn-id: https://svn.dealii.org/trunk@16502 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_tools.h | 86 ++++++++++++++- deal.II/deal.II/source/grid/grid_tools.cc | 127 ++++++++++++++-------- deal.II/doc/news/changes.h | 12 +- 3 files changed, 179 insertions(+), 46 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h index a4bce05c5d..9e6a88e70e 100644 --- a/deal.II/deal.II/include/grid/grid_tools.h +++ b/deal.II/deal.II/include/grid/grid_tools.h @@ -33,6 +33,8 @@ namespace hp template class MappingCollection; } +class SparsityPattern; + /** * This class is a collection of algorithms working on triangulations, @@ -446,7 +448,7 @@ class GridTools * values between zero and * @p n_partitions-1. You can access the * subdomain id of a cell by using - * cell->subdomain_id(). + * cell-@>subdomain_id(). * * This function will generate an error * if METIS is not installed unless @@ -461,6 +463,88 @@ class GridTools static void partition_triangulation (const unsigned int n_partitions, Triangulation &triangulation); + + /** + * This function does the same as the + * previous one, i.e. it partitions a + * triangulation using METIS into a + * number of subdomains identified by the + * cell-@>subdomain_id() + * flag. + * + * The difference to the previous + * function is the second argument, a + * sparsity pattern that represents the + * connectivity pattern between cells. + * + * While the function above builds it + * directly from the triangulation by + * considering which cells neighbor each + * other, this function can take a more + * refined connectivity graph. The + * sparsity pattern needs to be of size + * $N\times N$, where $N$ is the number + * of active cells in the + * triangulation. If the sparsity pattern + * contains an entry at position $(i,j)$, + * then this means that cells $i$ and $j$ + * (in the order in which they are + * traversed by active cell iterators) + * are to be considered connected; METIS + * will then try to partition the domain + * in such a way that (i) the subdomains + * are of roughly equal size, and (ii) a + * minimal number of connections are + * broken. + * + * This function is mainly useful in + * cases where connections between cells + * exist that are not present in the + * triangulation alone (otherwise the + * previous function would be the simpler + * one to use). Such connections may + * include that certain parts of the + * boundary of a domain are coupled + * through symmetric boundary conditions + * or integrals (e.g. friction contact + * between the two sides of a crack in + * the domain), or if a numerical scheme + * is used that not only connects + * immediate neighbors but a larger + * neighborhood of cells (e.g. when + * solving integral equations). + * + * In addition, this function may be + * useful in cases where the default + * sparsity pattern is not entirely + * sufficient. This can happen because + * the default is to just consider face + * neighbors, not neighboring cells that + * are connected by edges or + * vertices. While the latter couple when + * using continuous finite elements, they + * are typically still closely connected + * in the neighborship graph, and METIS + * will not usually cut important + * connections in this case. However, if + * there are vertices in the mesh where + * many cells (many more than the common + * 4 or 6 in 2d and 3d, respectively) + * come together, then there will be a + * significant number of cells that are + * connected across a vertex, but several + * degrees removed in the connectivity + * graph built only using face + * neighbors. In a case like this, METIS + * may sometimes make bad decisions and + * you may want to build your own + * connectivity graph. + */ + template + static + void partition_triangulation (const unsigned int n_partitions, + const SparsityPattern &cell_connection_graph, + Triangulation &triangulation); /** * For each active cell, return in the diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index c552eaf894..868ea6c964 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -839,61 +839,94 @@ partition_triangulation (const unsigned int n_partitions, // we decompose the domain by first // generating the connection graph of all // cells with their neighbors, and then - // passing this graph off to METIS. To make - // things a little simpler and more - // general, we let the function - // DoFTools:make_flux_sparsity_pattern - // function generate the connection graph - // for us and reuse the SparsityPattern - // data structure for the connection - // graph. The connection structure of the - // mesh is obtained by using a fake - // piecewise constant finite element - // - // Since in 3d the generation of a - // sparsity pattern can be expensive, we - // take the detour of the compressed - // sparsity pattern, which is a little - // slower but more efficient in terms of - // memory - FE_DGQ fake_q0(0); - DoFHandler dof_handler (triangulation); - dof_handler.distribute_dofs (fake_q0); - Assert (dof_handler.n_dofs() == triangulation.n_active_cells(), - ExcInternalError()); + // passing this graph off to METIS. + // + // as built in this function, we only + // consider face neighbors, which leads to + // a fixed number of entries per row (don't + // forget that each cell couples with + // itself) + SparsityPattern cell_connectivity (triangulation.n_active_cells(), + triangulation.n_active_cells(), + GeometryInfo::faces_per_cell+1); + + // next we have to build a mapping from the + // list of cells to their indices. for + // this, use the user_index field + std::vector saved_user_indices; + triangulation.save_user_indices (saved_user_indices); + unsigned int index = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell, ++index) + cell->set_user_index (index); + + // next loop over all cells and their + // neighbors to build the sparsity pattern + index = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell, ++index) + { + cell_connectivity.add (index, index); + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f) == false) + cell_connectivity.add (index, + cell->neighbor(f)->user_index()); + } - CompressedSparsityPattern csp (dof_handler.n_dofs(), - dof_handler.n_dofs()); - DoFTools::make_flux_sparsity_pattern (dof_handler, csp); + // now compress the so-built connectivity + // pattern and restore user indices + cell_connectivity.compress (); + triangulation.save_user_indices (saved_user_indices); - SparsityPattern sparsity_pattern; - sparsity_pattern.copy_from (csp); + // finally defer to the other function for + // partitioning and assigning subdomain ids + partition_triangulation (n_partitions, + cell_connectivity, + triangulation); +} + + + +template +void +GridTools:: +partition_triangulation (const unsigned int n_partitions, + const SparsityPattern &cell_connection_graph, + Triangulation &triangulation) +{ + Assert (n_partitions > 0, ExcInvalidNumberOfPartitions(n_partitions)); + Assert (cell_connection_graph.n_rows() == triangulation.n_active_cells(), + ExcMessage ("Connectivity graph has wrong size")); + Assert (cell_connection_graph.n_cols() == triangulation.n_active_cells(), + ExcMessage ("Connectivity graph has wrong size")); + + // check for an easy return + if (n_partitions == 1) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + cell->set_subdomain_id (0); + return; + } // partition this connection graph and get // back a vector of indices, one per degree // of freedom (which is associated with a // cell) std::vector partition_indices (triangulation.n_active_cells()); - sparsity_pattern.partition (n_partitions, partition_indices); + cell_connection_graph.partition (n_partitions, partition_indices); // finally loop over all cells and set the - // subdomain ids. for this, get the DoF - // index of each cell and extract the - // subdomain id from the vector obtained - // above + // subdomain ids std::vector dof_indices(1); - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - { - cell->get_dof_indices(dof_indices); - Assert (dof_indices[0] < triangulation.n_active_cells(), - ExcInternalError()); - Assert (partition_indices[dof_indices[0]] < n_partitions, - ExcInternalError()); - - cell->set_subdomain_id (partition_indices[dof_indices[0]]); - } + unsigned int index = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell, ++index) + cell->set_subdomain_id (partition_indices[index]); } @@ -1278,6 +1311,12 @@ void GridTools::partition_triangulation (const unsigned int, Triangulation &); +template +void +GridTools::partition_triangulation (const unsigned int, + const SparsityPattern &, + Triangulation &); + template void GridTools:: diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0394b0109d..eda6995d4f 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -234,18 +234,28 @@ inconvenience this causes.

deal.II

    +
  1. +

    + New: There is a second GridTools::partition_triangulation + function that takes a cell connectivity pattern as argument, rather + than computing it itself as the existing function. Use cases are + discussed in the documentation of the new function. + (WB 2008/08/06) +

  2. +
  3. Fixed: GridTools::find_cells_adjacent_to_vertex had a bug that prevented its correct functioning in three dimensions. Some cases were left out due to uncorrect assumptions on the various refinement possibilities. +
    (Luca Heltai 2008/07/17)

  4. New: There is now a new - Triangulation::prevent_distorted_boundary_cells function which + Triangulation::prevent_distorted_boundary_cells function which is only useful in case of anisotropic refinement. At the boundary of the domain, the new point on the face may be far inside the current cell, if the boundary has a strong curvature. If we -- 2.39.5