From: Guido Kanschat Date: Tue, 7 Mar 2006 15:00:55 +0000 (+0000) Subject: compute row lengths for block patterns X-Git-Tag: v8.0.0~12140 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2440607fcf03893653c8912668b5019b040d7719;p=dealii.git compute row lengths for block patterns git-svn-id: https://svn.dealii.org/trunk@12550 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 1c7078430b..1ed08d628f 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -273,17 +273,32 @@ class DoFTools * couplings. Here, couplings due * to flux operators on faces are * considered. See #Coupling. - * - * @todo This function is only - * implemented for primitive - * elements. - */ + */ template static void compute_row_length_vector( const DH& dofs, std::vector& row_lengths, const Table<2,Coupling>& couplings, + const Table<2,Coupling>& flux_couplings); + + /** + * Same as the other function, + * but separates the blocks for a + * block system of equations. + * + * The result vector + * row_lengths contains + * FiniteElementData::n_blocks() + * vectors of length + * DH::n_dofs(). + */ + template + static + void compute_row_length_vector( + const DH& dofs, + std::vector >& row_lengths, + const Table<2,Coupling>& couplings, const Table<2,Coupling>& flux_couplings); /** diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 9d58f7c227..e96463df52 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -243,12 +243,12 @@ DoFTools::compute_row_length_vector( #else -template +template static void compute_cell_row_length_matrix( Table<2, unsigned int>& matrix, - const typename DH::cell_iterator& cell, + const CellIterator& cell, const FiniteElement& fe, const Table<2,DoFTools::Coupling>& couple_cell, const Table<2,DoFTools::Coupling>& couple_face) @@ -400,7 +400,8 @@ compute_cell_row_length_matrix( for (unsigned int mult=0;mult static void compute_face_row_length_matrix( Table<2, unsigned int>& matrix, Table<2, unsigned int>& nmatrix, - const typename DH::face_iterator& face, - const typename DH::cell_iterator& cell, - const typename DH::cell_iterator& neighbor, const FiniteElement& fe, const FiniteElement& nfe, const Table<2,DoFTools::Coupling>& couple_cell, @@ -715,6 +715,9 @@ DoFTools::compute_row_length_vector( std::vector > couple_face; convert_couplings_to_blocks(dofs, couplings, couple_cell); convert_couplings_to_blocks(dofs, flux_couplings, couple_face); + + Table<2, unsigned int> cell_couplings; + Table<2, unsigned int> neighbor_couplings; // We loop over cells and go from // cells to lower dimensional @@ -737,104 +740,127 @@ DoFTools::compute_row_length_vector( Assert (flux_couplings.n_cols()==fe.n_components(), ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components())); + cell_couplings.reinit(fe.dofs_per_cell, fe.n_blocks()); cell_indices.resize(fe.dofs_per_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; - while (i < fe.first_line_index) - { - for (unsigned int base=0;base( + cell_couplings, cell, fe, couple_cell[fe_index], couple_face[fe_index]); - // In all other cases we - // subtract adjacent faces to be - // added in the loop below. - while (i < fe.first_quad_index) + // 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) { - for (unsigned int base=0;base1) - ? (DH::dimension-1) - : GeometryInfo::faces_per_cell) - * fe.base_element(base).dofs_per_face; - row_lengths[cell_indices[i]] += increment; - } - ++i; + const typename DH::cell_iterator neighbor = cell->neighbor(iface); + const FiniteElement& nfe = neighbor->get_fe(); + typename DH::face_iterator face = cell->face(iface); + // Do this only once per + // face and not on the + // hanging faces. + if (face->user_flag_set() || neighbor->has_children()) + continue; + face->set_user_flag(); + + neighbor_couplings.reinit(nfe.dofs_per_cell, nfe.n_blocks()); + neighbor_indices.resize(nfe.dofs_per_cell); + neighbor->get_dof_indices(neighbor_indices); + + compute_face_row_length_matrix( + cell_couplings, neighbor_couplings,// face, cell, neighbor, + fe, nfe, couple_cell[fe_index], couple_face[fe_index]); + + for (unsigned int i=0;i +void +DoFTools::compute_row_length_vector( + const DH& dofs, + std::vector >& row_lengths, + const Table<2,Coupling>& couplings, + const Table<2,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 DH::cell_iterator end = dofs.end(); + typename DH::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; + convert_couplings_to_blocks(dofs, couplings, couple_cell); + convert_couplings_to_blocks(dofs, flux_couplings, couple_face); + + Table<2, unsigned int> cell_couplings; + Table<2, unsigned int> neighbor_couplings; + + // 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(); - // Now quads in 2D and 3D - while (i < fe.first_hex_index) - { - for (unsigned int base=0;base2) - ? (DH::dimension-2) - : GeometryInfo::faces_per_cell) - * fe.base_element(base).dofs_per_face; - row_lengths[cell_indices[i]] += increment; - } - ++i; - } + 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())); - // 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; - } + cell_couplings.reinit(fe.dofs_per_cell, fe.n_blocks()); + cell_indices.resize(fe.dofs_per_cell); + cell->get_dof_indices(cell_indices); + + compute_cell_row_length_matrix( + cell_couplings, cell, fe, couple_cell[fe_index], couple_face[fe_index]); // At this point, we have // counted all dofs @@ -849,79 +875,32 @@ DoFTools::compute_row_length_vector( // flux contributions. for (unsigned int iface=0;iface::faces_per_cell;++iface) { - if (cell->at_boundary(iface)) - { - for (unsigned int i=0;ineighbor(iface); const FiniteElement& nfe = neighbor->get_fe(); typename DH::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() || 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? + face->set_user_flag(); - // This will not work with - // different couplings on - // different cells. + neighbor_couplings.reinit(nfe.dofs_per_cell, nfe.n_blocks()); neighbor_indices.resize(nfe.dofs_per_cell); neighbor->get_dof_indices(neighbor_indices); - for (unsigned int base=0;base( + cell_couplings, neighbor_couplings,// face, cell, neighbor, + fe, nfe, couple_cell[fe_index], couple_face[fe_index]); + + for (unsigned int i=0;i