From: Guido Kanschat Date: Fri, 27 Jan 2006 06:27:45 +0000 (+0000) Subject: compute_row_length_vector for MGDoFHandler X-Git-Tag: v8.0.0~12487 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ab22b0cdfdc2b7bae233894c2a0aad629d04f448;p=dealii.git compute_row_length_vector for MGDoFHandler git-svn-id: https://svn.dealii.org/trunk@12188 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 7adef23e6b..b4a3c7f86a 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -283,8 +283,8 @@ class DoFTools void compute_row_length_vector( const DH& dofs, std::vector& row_lengths, - const Table<2,Coupling>& couplings /*= typename Table<2,Coupling>()*/, - const Table<2,Coupling>& flux_couplings /*= Table<2,Coupling>()*/); + const Table<2,Coupling>& couplings, + const Table<2,Coupling>& flux_couplings); /** * Locate non-zero entries of the @@ -1085,7 +1085,6 @@ class DoFTools std::vector target_component = std::vector()); - /** * @deprecated See the previous * function with the same name diff --git a/deal.II/deal.II/include/multigrid/mg_tools.h b/deal.II/deal.II/include/multigrid/mg_tools.h index 43484488c5..9514c1bb2f 100644 --- a/deal.II/deal.II/include/multigrid/mg_tools.h +++ b/deal.II/deal.II/include/multigrid/mg_tools.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005 by the deal.II authors +// Copyright (C) 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -46,6 +46,33 @@ template class FullMatrix; class MGTools { public: + /** + * Compute row length vector for + * multilevel methods. + */ + template + static + void compute_row_length_vector( + const MGDoFHandler& dofs, + const unsigned int level, + std::vector& row_lengths, + const DoFTools::Coupling flux_couplings = none); + + /** + * Compute row length vector for + * multilevel methods with + * optimization for block + * couplings. + */ + template + static + void compute_row_length_vector( + const MGDoFHandler& dofs, + const unsigned int level, + std::vector& row_lengths, + const Table<2,DoFTools::Coupling>& couplings, + const Table<2,DoFTools::Coupling>& flux_couplings); + /** * Write the sparsity structure * of the matrix belonging to the diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 1c1478b2e0..0c345d6bd0 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -408,7 +408,6 @@ DoFTools::compute_row_length_vector( } -//TODO:[GK] Extend this to non-primitive elements // This is the template for 2D and 3D. See version for 1D above template void @@ -462,7 +461,6 @@ DoFTools::compute_row_length_vector( const FiniteElement& fe = cell->get_fe(); const unsigned int fe_index = cell->active_fe_index(); - Assert (fe.is_primitive(), typename FiniteElement::ExcFENotPrimitive()); Assert (couplings.n_rows()==fe.n_components(), ExcDimensionMismatch(couplings.n_rows(), fe.n_components())); Assert (couplings.n_cols()==fe.n_components(), @@ -472,8 +470,8 @@ DoFTools::compute_row_length_vector( Assert (flux_couplings.n_cols()==fe.n_components(), ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components())); - cell->get_dof_indices(cell_indices); cell_indices.resize(fe.dofs_per_cell); + cell->get_dof_indices(cell_indices); unsigned int i = 0; // First, dofs on // vertices. We assume that @@ -640,7 +638,7 @@ DoFTools::compute_row_length_vector( // This will not work with // different couplings on // different cells. - neighbor_indices.resize(fe.dofs_per_cell); + neighbor_indices.resize(nfe.dofs_per_cell); neighbor->get_dof_indices(neighbor_indices); for (unsigned int base=0;base +#if deal_II_dimension == 1 + +template +void +MGTools::compute_row_length_vector( + const MGDoFHandler&, + const unsigned int, + std::vector&, + const DoFTools::Coupling) +{ + Assert(false, ExcNotImplemented()); +} + + +template +void +MGTools::compute_row_length_vector( + const MGDoFHandler&, + const unsigned int, + std::vector&, + const Table<2,DoFTools::Coupling>&, + const Table<2,DoFTools::Coupling>&) +{ + Assert(false, ExcNotImplemented()); +} + + +#else + +// Template for 2D and 3D. For 1D see specialization above +template +void +MGTools::compute_row_length_vector( + const MGDoFHandler& dofs, + const unsigned int level, + std::vector& row_lengths, + const DoFTools::Coupling flux_coupling) +{ + 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 MGDoFHandler::cell_iterator end = dofs.end(level); + typename MGDoFHandler::active_cell_iterator cell; + std::vector cell_indices; + std::vector neighbor_indices; + + // We loop over cells and go from + // cells to lower dimensional + // objects. This is the only way to + // cope with the fact, that an + // unknown number of cells may + // share an object of dimension + // smaller than dim-1. + for (cell = dofs.begin(level); cell != end; ++cell) + { + const FiniteElement& fe = cell->get_fe(); + cell_indices.resize(fe.dofs_per_cell); + cell->get_mg_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. +//TODO: This assumes that even in hp context, the dofs per face coincide! + 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) + { + bool level_boundary = cell->at_boundary(iface); + typename MGDoFHandler::cell_iterator neighbor; + if (!level_boundary) + { + neighbor = cell->neighbor(iface); + if (static_cast(neighbor->level()) != level) + level_boundary = true; + } + + if (level_boundary) + { + for (unsigned int i=0;i& nfe = neighbor->get_fe(); + typename MGDoFHandler::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 != DoFTools::none) + { + unsigned int increment = nfe.dofs_per_cell - nfe.dofs_per_face; + for (unsigned int i=0;iuser_flag_set()) + 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. + neighbor_indices.resize(nfe.dofs_per_cell); + neighbor->get_mg_dof_indices(neighbor_indices); + for (unsigned int i=0;i +void +MGTools::compute_row_length_vector( + const MGDoFHandler& dofs, + const unsigned int level, + std::vector& row_lengths, + const Table<2,DoFTools::Coupling>& couplings, + const Table<2,DoFTools::Coupling>& flux_couplings) +{ + 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 MGDoFHandler::cell_iterator end = dofs.end(level); + typename MGDoFHandler::active_cell_iterator cell; + std::vector cell_indices; + std::vector neighbor_indices; + + // We have to translate the + // couplings from components to + // blocks, so this works for + // nonprimitive elements as well. + std::vector > couple_cell; + std::vector > couple_face; + DoFTools::convert_couplings_to_blocks(dofs, couplings, couple_cell); + DoFTools::convert_couplings_to_blocks(dofs, flux_couplings, couple_face); + + // 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) + { + const FiniteElement& fe = cell->get_fe(); + const unsigned int fe_index = cell->active_fe_index(); + + 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())); + + cell_indices.resize(fe.dofs_per_cell); + cell->get_mg_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 base=0;base1) + ? (dim-1) + : GeometryInfo::faces_per_cell) + * fe.base_element(base).dofs_per_face; + row_lengths[cell_indices[i]] += increment; + } + ++i; + } + + // Now quads in 2D and 3D + while (i < fe.first_hex_index) + { + for (unsigned int base=0;base2) + ? (dim-2) + : GeometryInfo::faces_per_cell) + * fe.base_element(base).dofs_per_face; + row_lengths[cell_indices[i]] += increment; + } + ++i; + } + + // Finally, cells in 3D + while (i < fe.dofs_per_cell) + { + for (unsigned int base=0;base::faces_per_cell + * fe.base_element(base).dofs_per_face; + 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) + { + bool level_boundary = cell->at_boundary(iface); + typename MGDoFHandler::cell_iterator neighbor; + if (!level_boundary) + { + neighbor = cell->neighbor(iface); + if (static_cast(neighbor->level()) != level) + level_boundary = true; + } + + if (level_boundary) + { + for (unsigned int i=0;i& nfe = neighbor->get_fe(); + typename MGDoFHandler::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 base=0;baseuser_flag_set()) + 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? + + // This will not work with + // different couplings on + // different cells. + neighbor_indices.resize(nfe.dofs_per_cell); + neighbor->get_mg_dof_indices(neighbor_indices); + for (unsigned int base=0;base void MGTools::make_sparsity_pattern ( const MGDoFHandler &dof, @@ -190,8 +654,8 @@ MGTools::make_flux_sparsity_pattern ( const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level, - const Table<2,DoFTools::Coupling> &int_mask, - const Table<2,DoFTools::Coupling> &flux_mask) + const Table<2,DoFTools::DoFTools::Coupling> &int_mask, + const Table<2,DoFTools::DoFTools::Coupling> &flux_mask) { const FiniteElement& fe = dof.get_fe(); const unsigned int n_dofs = dof.n_dofs(level); @@ -218,8 +682,8 @@ MGTools::make_flux_sparsity_pattern ( typename MGDoFHandler::cell_iterator cell = dof.begin(level), endc = dof.end(level); - Table<2,DoFTools::Coupling> int_dof_mask(total_dofs, total_dofs); - Table<2,DoFTools::Coupling> flux_dof_mask(total_dofs, total_dofs); + Table<2,DoFTools::DoFTools::Coupling> int_dof_mask(total_dofs, total_dofs); + Table<2,DoFTools::DoFTools::Coupling> flux_dof_mask(total_dofs, total_dofs); DoFTools::compute_dof_couplings(int_dof_mask, int_mask, fe); DoFTools::compute_dof_couplings(flux_dof_mask, flux_mask, fe); @@ -369,7 +833,7 @@ MGTools::make_flux_sparsity_pattern_edge ( const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level, - const Table<2,DoFTools::Coupling> &flux_mask) + const Table<2,DoFTools::DoFTools::Coupling> &flux_mask) { const FiniteElement& fe = dof.get_fe(); const unsigned int n_comp = fe.n_components(); @@ -399,7 +863,7 @@ MGTools::make_flux_sparsity_pattern_edge ( typename MGDoFHandler::cell_iterator cell = dof.begin(level), endc = dof.end(level); - Table<2,DoFTools::Coupling> flux_dof_mask(dofs_per_cell, dofs_per_cell); + Table<2,DoFTools::DoFTools::Coupling> flux_dof_mask(dofs_per_cell, dofs_per_cell); DoFTools::compute_dof_couplings(flux_dof_mask, flux_mask, fe); @@ -702,6 +1166,16 @@ MGTools::reinit_vector (const MGDoFHandler &mg_dof, #undef PATTERN +template void +MGTools::compute_row_length_vector( + const MGDoFHandler&, unsigned int, + std::vector&, const DoFTools::Coupling); +template void +MGTools::compute_row_length_vector( + const MGDoFHandler&, unsigned int, + std::vector&, + const Table<2,DoFTools::Coupling>&, const Table<2,DoFTools::Coupling>&); + template void MGTools::reinit_vector ( const MGDoFHandler&, MGLevelObject >&);