From: guido Date: Sat, 16 Jul 2005 08:40:28 +0000 (+0000) Subject: function for estimating row lengths X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=535192262795ccc54e398da137a32c734faa6270;p=dealii-svn.git function for estimating row lengths git-svn-id: https://svn.dealii.org/trunk@11152 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 10a53886b7..65b8abc891 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -178,6 +178,38 @@ class DoFTools * @{ */ + /** + * On initializing a + * SparsityPattern, the number of + * entries required for each row + * must be known. This function + * estimates the maximum number + * of coupling degrees of freedom + * for each row of the sparsity + * pattern. + * + * @param dofs The DoFHandler + * @param row_lengths The vector + * containing the resulting row + * lengths. It must have the + * length DoFHandler::n_dofs(). + * + * @param flux_couplings + * Additional couplings between + * shape functions due to flux + * operators on faces are + * considered. If false, + * the default, only topological + * couplings between cells are + * counted. + */ + template + static + void compute_row_length_vector( + const DoFHandler& dofs, + std::vector& row_lengths, + const Coupling flux_couplings = none); + /** * On initializing a * SparsityPattern, the number of diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 05fcda862a..22e8db1d14 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -37,7 +37,213 @@ #include -//TODO:[GK] Traverse faces only once using flags + +#if deal_II_dimension == 1 + + +template +void +DoFTools::compute_row_length_vector( + const DoFHandler& dofs, + std::vector& row_lengths, + const Coupling) +{ + const FiniteElement& fe = dofs.get_fe(); + Assert (row_lengths.size() == dofs.n_dofs(), + ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); + + // Function starts here by + // resetting the counters. + std::fill(row_lengths.begin(), row_lengths.end(), 0); + + const typename DoFHandler::cell_iterator end = dofs.end(); + typename DoFHandler::active_cell_iterator cell; + std::vector cell_indices(fe.dofs_per_cell); + + for (cell = dofs.begin_active(); cell != end; ++cell) + { + cell->get_dof_indices(cell_indices); + unsigned int i = 0; + unsigned int increment = fe.dofs_per_cell; + while (i < fe.dofs_per_cell) + row_lengths[cell_indices[i++]] += increment; + } +} + + +#else + +template +void +DoFTools::compute_row_length_vector( + const DoFHandler& dofs, + std::vector& row_lengths, + const Coupling flux_coupling) +{ + const FiniteElement& fe = dofs.get_fe(); + Assert (row_lengths.size() == dofs.n_dofs(), + ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); + + // Function starts here by + // resetting the counters. + std::fill(row_lengths.begin(), row_lengths.end(), 0); + // We need the user flags, so we + // save them for later restoration + std::vector old_flags; + // We need a non-constant + // triangulation for the user + // flags. Since we restore them in + // the end, this cast is safe. + Triangulation& user_flags_triangulation = + const_cast&> (dofs.get_tria()); + user_flags_triangulation.save_user_flags(old_flags); + user_flags_triangulation.clear_user_flags(); + + const typename DoFHandler::cell_iterator end = dofs.end(); + typename DoFHandler::active_cell_iterator cell; + std::vector cell_indices(fe.dofs_per_cell); + std::vector neighbor_indices(fe.dofs_per_cell); + + // We loop over cells and go from + // cells to lower dimensional + // objects. This is the only way to + // cope withthe fact, that an + // unknown number of cells may + // share an object of dimension + // smaller than dim-1. + for (cell = dofs.begin_active(); cell != end; ++cell) + { + cell->get_dof_indices(cell_indices); + unsigned int i = 0; + // First, dofs on + // vertices. We assume that + // each vertex dof couples + // with all dofs on + // adjacent grid cells. + + // Adding all dofs of the cells + // will add dofs of the faces + // of the cell adjacent to the + // vertex twice. Therefore, we + // subtract these here and add + // them in a loop over the + // faces below. + + // in 1d, faces and vertices + // are identical. Nevertheless, + // this will only work if + // dofs_per_face is zero and + // dofs_per_vertex is + // arbitrary, not the other way + // round. + unsigned int increment = fe.dofs_per_cell - dim * fe.dofs_per_face; + while (i < fe.first_line_index) + row_lengths[cell_indices[i++]] += increment; + // From now on, if an object is + // a cell, its dofs only couple + // inside the cell. Since the + // faces are handled below, we + // have to subtract ALL faces + // in this case. + + // In all other cases we + // subtract adjacent faces to be + // added in the loop below. + increment = (dim>1) + ? fe.dofs_per_cell - (dim-1) * fe.dofs_per_face + : fe.dofs_per_cell - GeometryInfo::faces_per_cell * fe.dofs_per_face; + while (i < fe.first_quad_index) + row_lengths[cell_indices[i++]] += increment; + + // Now quads in 2D and 3D + increment = (dim>2) + ? fe.dofs_per_cell - (dim-2) * fe.dofs_per_face + : fe.dofs_per_cell - GeometryInfo::faces_per_cell * fe.dofs_per_face; + while (i < fe.first_hex_index) + row_lengths[cell_indices[i++]] += increment; + // Finally, cells in 3D + increment = fe.dofs_per_cell - GeometryInfo::faces_per_cell * fe.dofs_per_face; + while (i < fe.dofs_per_cell) + row_lengths[cell_indices[i++]] += increment; + + // At this point, we have + // counted all dofs + // contributiong from cells + // coupled topologically to the + // adjacent cells, but we + // subtracted some faces. + + // Now, let's go by the faces + // and add the missing + // contribution as well as the + // flux contributions. + for (unsigned int iface=0;iface::faces_per_cell;++iface) + { + if (cell->at_boundary(iface)) + { + for (unsigned int i=0;i::cell_iterator neighbor = cell->neighbor(iface); +// const bool neighbor_same_level = (neighbor->level() == cell->level()); +// const unsigned int +// nface = neighbor_same_level +// ? cell->neighbor_of_neighbor(iface) +// : cell->neighbor_of_coarser_neighbor(iface).first; + typename DoFHandler::face_iterator face = cell->face(iface); + + // Flux couplings are + // computed from both sides + // for simplicity. + + // The dofs on the common face + // will be handled below, + // therefore, we subtract them + // here. + if (flux_coupling != none) + { + unsigned int increment = fe.dofs_per_cell - fe.dofs_per_face; + for (unsigned int i=0;iuser_flag_set() || neighbor->has_children()) + continue; + face->set_user_flag(); + // At this point, we assume + // that each cell added its + // dofs minus the face to + // the couplings of the + // face dofs. Since we + // subtracted two faces, we + // have to re-add one. + + // If one side of the face + // is refined, all the fine + // face dofs couple with + // the coarse one. + + // Wolfgang, do they couple + // with each other by + // constraints? + neighbor->get_dof_indices(neighbor_indices); + for (unsigned int i=0;i @@ -48,6 +254,8 @@ DoFTools::compute_row_length_vector( Table<2,Coupling> couplings, Table<2,Coupling> flux_couplings) { + Assert(false, ExcNotImplemented()); + const FiniteElement& fe = dofs.get_fe(); const unsigned int ncomp = fe.n_components; @@ -69,56 +277,20 @@ DoFTools::compute_row_length_vector( // Function starts here by // resetting the counters. std::fill(row_lengths.begin(), row_lengths.end(), 0); - + // We need the user flags, so we + // save them for later restoration + std::vector old_flags; + user_flags_triangulation.save_user_flags(old_flags); + user_flags_triangulation.clear_user_flags(); + const typename DoFHandler::cell_iterator end = dofs.end(); typename DoFHandler::active_cell_iterator cell; - std::vector indices(fe.dofs_per_cell); + std::vector cell_indices(fe.dofs_per_cell); + std::vector face_indices(fe.dofs_per_face); for (cell = dofs.begin(); cell != end; ++cell) { - indices.resize(fe.dofs_per_cell); - cell->get_dof_indices(indices); - // First the simple case where - // all degrees of freedom on a - // cell couple. - if (couplings.n_rows() <= 1) - { - unsigned int i = 0; - // First, dofs on - // vertices. We assume that - // each vertex dof couples - // with all dofs on - // adjacent grid cells. - - // Since two cells share a - // facein 2d, we can subtract - // the dofs per face, - // unless we have hanging - // nodes. - - // In 3d, 8 cells share 12 - // faces. - unsigned int increment = fe.dofs_per_cell - fe.dofs_per_face; - while (i < fe.first_line_index) - row_lengths[indices[i++]] += increment; - // Lines are already cells - // in 1D - increment = (dim>1) - ? fe.dofs_per_cell - fe.dofs_per_face - : fe.dofs_per_cell; - while (i < fe.first_quad_index) - row_lengths[indices[i++]] += increment; - // Now quads in 2D and 3D - increment = (dim>2) - ? fe.dofs_per_cell - fe.dofs_per_face - : fe.dofs_per_cell; - while (i < fe.first_hex_index) - row_lengths[indices[i++]] += increment; - // Finally, cells in 3D - increment = fe.dofs_per_cell; - while (i < fe.ofs_per_cell) - row_lengths[indices[i++]] += increment; - } + cell->get_dof_indices(cell_indices); // At this point, we have // counted all dofs // contributiong from cells @@ -132,14 +304,32 @@ DoFTools::compute_row_length_vector( // flux contributions. for (unsigned int iface=0;iface::faces_per_cell;++iface) { + if (cell->at_boundary(iface)) + continue; + const bool neighbor_same_level = (cell->neighbor(iface)->level() + == cell->level()); + const unsigned int + nface = (neighbor_same_level) + ? cell->neighbor_of_neighbor(iface) + : cell->neighbor_of_coarser_neighbor(iface); typename DoFHandler::face_iterator face = cell->face(iface); - indices.resize(fe.dofs_per_face); - face->get_dof_indices(indices); + // Flux couplings are + // computed from both sides + // for simplicity. for (unsigned int i=0;i& dofs, std::vector& row_lengths, + const Coupling flux_coupling); + + template void DoFTools::make_sparsity_pattern (const DoFHandler &dof,