From: guido Date: Wed, 24 Aug 2005 07:47:39 +0000 (+0000) Subject: compute_row_length_vector for systems X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1b48765a057772776bd5cf399d55e9f950df9323;p=dealii-svn.git compute_row_length_vector for systems git-svn-id: https://svn.dealii.org/trunk@11312 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 65b8abc891..7b5a2bc4e0 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -225,8 +225,15 @@ class DoFTools * not coupling in the * differential equation can be * eliminated by the two optional - * tables. + * tables, which may reduce the + * amount of pre-allocated + * memory dramatically. * + * @todo This function is only + * implemented for primitive + * elements. Implementation in 1D + * is missing completely. + * * @param dofs The DoFHandler * @param row_lengths The vector * containing the resulting row @@ -254,8 +261,8 @@ class DoFTools void compute_row_length_vector( const DoFHandler& dofs, std::vector& row_lengths, - Table<2,Coupling> couplings /*= typename Table<2,Coupling>()*/, - Table<2,Coupling> flux_couplings /*= Table<2,Coupling>()*/); + const Table<2,Coupling>& couplings /*= typename Table<2,Coupling>()*/, + const Table<2,Coupling>& flux_couplings /*= Table<2,Coupling>()*/); /** * Locate non-zero entries of the @@ -330,7 +337,7 @@ class DoFTools /** * Locate non-zero entries for - * ector valued finite elements. + * vector valued finite elements. * This function does mostly the * same as the previous * @p make_sparsity_pattern, but @@ -372,16 +379,16 @@ class DoFTools * contained in there. * * This function is designed to - * accept a mask, like the one + * accept a coupling pattern, like the one * shown above, through the - * @p mask parameter, which - * contains boolean values. It + * @p couplings parameter, which + * contains values of type #Coupling. It * builds the matrix structure * just like the previous * function, but does not create * matrix elements if not - * specified by the mask. If the - * mask is symmetric, then so + * specified by the coupling pattern. If the + * couplings are symmetric, then so * will be the resulting sparsity * pattern. * @@ -402,11 +409,24 @@ class DoFTools * more than one component (in * deal.II speak: they are * non-primitive). In this case, - * the mask element correspoding - * to the first non-zero - * component is taken, and the - * mask elements of all following - * components are ignored. + * the coupling element + * correspoding to the first + * non-zero component is taken + * and additional ones for this + * component are ignored. + */ + template + static + void + make_sparsity_pattern (const DoFHandler& dof, + const Table<2, Coupling>& coupling, + SparsityPattern& sparsity_pattern); + /** + * @deprecated This is the old + * form of the previous + * function. It generates a table + * of #Coupling values and calls + * it. */ template static @@ -1554,12 +1574,39 @@ class DoFTools // ---------------------- inline and template functions -------------------- +template +inline +void +DoFTools::make_sparsity_pattern ( + const DoFHandler &dof, + const std::vector > &mask, + SparsityPattern &sparsity_pattern) +{ + const unsigned int ncomp = dof.get_fe().n_components(); + + Assert (mask.size() == ncomp, + ExcDimensionMismatch(mask.size(), ncomp)); + for (unsigned int i=0; i couplings(ncomp, ncomp); + for (unsigned int i=0;i +inline void -DoFTools:: -map_support_points_to_dofs (const Mapping &mapping, - const DoFHandler &dof_handler, - std::map, unsigned int, Comp> &point_to_index_map) +DoFTools::map_support_points_to_dofs ( + const Mapping &mapping, + const DoFHandler &dof_handler, + std::map, unsigned int, Comp> &point_to_index_map) { // let the checking of arguments be // done by the function first diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 74a3a82a58..92fb934be2 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -40,13 +40,13 @@ #if deal_II_dimension == 1 - +// Specialization for 1D template void DoFTools:: compute_row_length_vector(const DoFHandler &dofs, std::vector &row_lengths, - const Coupling) + const Coupling flux_coupling) { const FiniteElement& fe = dofs.get_fe(); Assert (row_lengths.size() == dofs.n_dofs(), @@ -68,12 +68,87 @@ compute_row_length_vector(const DoFHandler &dofs, // dof on this cell for (unsigned int i=0; i::faces_per_cell;++face) + { + if (cell->at_boundary(face)) continue; + for (unsigned int i=0; i +void +DoFTools::compute_row_length_vector( + const DoFHandler& dofs, + std::vector& row_lengths, + const Table<2,Coupling>& couplings, + const Table<2,Coupling>& flux_couplings) +{ + const FiniteElement& fe = dofs.get_fe(); + Assert (fe.is_primitive(), FiniteElementBase::ExcFENotPrimitive()); + Assert (row_lengths.size() == dofs.n_dofs(), + ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); + Assert (couplings.n_rows()==fe.n_components(), + ExcDimensionMismatch(couplings.n_rows(), fe.n_components())); + Assert (couplings.n_cols()==fe.n_components(), + ExcDimensionMismatch(couplings.n_cols(), fe.n_components())); + Assert (flux_couplings.n_rows()==fe.n_components(), + ExcDimensionMismatch(flux_couplings.n_rows(), fe.n_components())); + Assert (flux_couplings.n_cols()==fe.n_components(), + ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components())); + + // Function starts here by + // resetting the counters. + std::fill(row_lengths.begin(), row_lengths.end(), 0); + + // Two vectors containing the + // number of dofs of a cell, split + // by components. + std::vector dofs_cell(fe.n_components()); + + for (unsigned int comp=0;comp::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); + + // each dof can couple with each other + // dof on this cell + for (unsigned int i=0; i::faces_per_cell;++face) + { + if (cell->at_boundary(face)) continue; + for (unsigned int i=0; i void DoFTools::compute_row_length_vector( @@ -243,21 +318,236 @@ DoFTools::compute_row_length_vector( user_flags_triangulation.load_user_flags(old_flags); } -#endif - -//TODO:[GK] This function is not finished yet!!! +//TODO:[GK] Extend this to non-primitive elements +// This is the template for 2D and 3D. See version for 1D above template void DoFTools::compute_row_length_vector( const DoFHandler& dofs, std::vector& row_lengths, - Table<2,Coupling> couplings, - Table<2,Coupling> flux_couplings) + const Table<2,Coupling>& couplings, + const Table<2,Coupling>& flux_couplings) { - Assert(false, ExcNotImplemented()); + const FiniteElement& fe = dofs.get_fe(); + Assert (fe.is_primitive(), FiniteElementBase::ExcFENotPrimitive()); + Assert (row_lengths.size() == dofs.n_dofs(), + ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); + Assert (couplings.n_rows()==fe.n_components(), + ExcDimensionMismatch(couplings.n_rows(), fe.n_components())); + Assert (couplings.n_cols()==fe.n_components(), + ExcDimensionMismatch(couplings.n_cols(), fe.n_components())); + Assert (flux_couplings.n_rows()==fe.n_components(), + ExcDimensionMismatch(flux_couplings.n_rows(), fe.n_components())); + Assert (flux_couplings.n_cols()==fe.n_components(), + ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components())); + + // 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); + + // Two vectors containing the + // number of dofs of a cell and of + // a face, split by components. + std::vector dofs_cell(fe.n_components()); + std::vector dofs_face(fe.n_components()); + + for (unsigned int comp=0;compget_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; + while (i < fe.first_line_index) + { + for (unsigned int comp=0;comp1) + ? dofs_cell[comp] - (dim-1) * dofs_face[comp] + : dofs_cell[comp] - GeometryInfo::faces_per_cell * dofs_face[comp]; + row_lengths[cell_indices[i]] += increment; + } + ++i; + } + + // Now quads in 2D and 3D + while (i < fe.first_hex_index) + { + for (unsigned int comp=0;comp2) + ? dofs_cell[comp] - (dim-2) * dofs_face[comp] + : dofs_cell[comp] - GeometryInfo::faces_per_cell * dofs_face[comp]; + row_lengths[cell_indices[i]] += increment; + } + ++i; + } + + // Finally, cells in 3D + while (i < fe.dofs_per_cell) + { + for (unsigned int comp=0;comp::faces_per_cell * dofs_face[comp]; + row_lengths[cell_indices[i]] += increment; + } + ++i; + } + + // 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. + for (unsigned int comp=0;compuser_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 comp=0;comp void @@ -292,9 +582,9 @@ DoFTools::make_sparsity_pattern ( template void DoFTools::make_sparsity_pattern ( - const DoFHandler &dof, - const std::vector > &mask, - SparsityPattern &sparsity) + const DoFHandler& dof, + const Table<2,Coupling>& couplings, + SparsityPattern& sparsity) { const unsigned int n_dofs = dof.n_dofs(); const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; @@ -303,11 +593,10 @@ DoFTools::make_sparsity_pattern ( ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); - Assert (mask.size() == dof.get_fe().n_components(), - ExcDimensionMismatch(mask.size(), dof.get_fe().n_components())); - for (unsigned int i=0; i& dofs, std::vector& row_lengths, const Coupling flux_coupling); +template void +DoFTools::compute_row_length_vector( + const DoFHandler& dofs, std::vector& row_lengths, + const Table<2,Coupling>& couplings, const Table<2,Coupling>& flux_couplings); template void DoFTools::make_sparsity_pattern @@ -3128,27 +3420,19 @@ DoFTools::make_sparsity_pattern -(const DoFHandler &dof, - const std::vector > &mask, - SparsityPattern &sparsity); +(const DoFHandler&, const Table<2,Coupling>&, SparsityPattern&); template void DoFTools::make_sparsity_pattern -(const DoFHandler &dof, - const std::vector > &mask, - CompressedSparsityPattern &sparsity); +(const DoFHandler&, const Table<2,Coupling>&, CompressedSparsityPattern&); template void DoFTools::make_sparsity_pattern -(const DoFHandler &dof, - const std::vector > &mask, - BlockSparsityPattern &sparsity); +(const DoFHandler&, const Table<2,Coupling>&, BlockSparsityPattern&); template void DoFTools::make_sparsity_pattern -(const DoFHandler &dof, - const std::vector > &mask, - CompressedBlockSparsityPattern &sparsity); +(const DoFHandler&, const Table<2,Coupling>&, CompressedBlockSparsityPattern&); #if deal_II_dimension > 1 template void