From 0994b9289c3e9cacd1bf6e668ab679bbaae5a911 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Thu, 22 Dec 2005 09:47:32 +0000 Subject: [PATCH] compute_row_length_vector for hpDoFHandler git-svn-id: https://svn.dealii.org/trunk@11916 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_tools.h | 11 +- deal.II/deal.II/source/dofs/dof_tools.cc | 250 ++++++++++++----------- 2 files changed, 133 insertions(+), 128 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index e38c2615b6..3a97a1f3df 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -223,10 +223,10 @@ class DoFTools * couplings between cells are * counted. */ - template + template static void compute_row_length_vector( - const DoFHandler& dofs, + const DH& dofs, std::vector& row_lengths, const Coupling flux_couplings = none); @@ -273,13 +273,12 @@ class DoFTools * * @todo This function is only * implemented for primitive - * elements. Implementation in 1D - * is missing completely. + * elements. */ - template + template static void compute_row_length_vector( - const DoFHandler& dofs, + 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>()*/); diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 4a0fa4f1da..e11ec399d4 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -135,14 +135,14 @@ namespace #if deal_II_dimension == 1 // Specialization for 1D -template +template void DoFTools:: -compute_row_length_vector(const DoFHandler &dofs, - std::vector &row_lengths, - const Coupling flux_coupling) +compute_row_length_vector( + const DH& 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())); @@ -150,12 +150,14 @@ compute_row_length_vector(const DoFHandler &dofs, // 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); + const typename DH::cell_iterator end = dofs.end(); + typename DH::active_cell_iterator cell; + std::vector cell_indices; for (cell = dofs.begin_active(); cell != end; ++cell) { + const FiniteElement& fe = cell->get_fe(); + cell_indices.resize(fe.dofs_per_cell); cell->get_dof_indices(cell_indices); // each dof can couple with each other @@ -166,58 +168,52 @@ compute_row_length_vector(const DoFHandler &dofs, // If fluxes couple, add // coupling to neighbor cells if (flux_coupling != none) - for (unsigned int face=0;face::faces_per_cell;++face) + for (unsigned int face=0;face::faces_per_cell;++face) { if (cell->at_boundary(face)) continue; - + const FiniteElement& nfe = cell->get_fe(); for (unsigned int i=0; i +template void DoFTools::compute_row_length_vector( - const DoFHandler& dofs, + const DH& 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(), typename FiniteElement::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()); + std::fill(row_lengths.begin(), row_lengths.end(), 0); - 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); + const typename DH::cell_iterator end = dofs.end(); + typename DH::active_cell_iterator cell; + std::vector cell_indices; for (cell = dofs.begin_active(); cell != end; ++cell) { + const FiniteElement& fe = cell->get_fe(); + 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(), + 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_dof_indices(cell_indices); // each dof can couple with each other @@ -225,17 +221,20 @@ DoFTools::compute_row_length_vector( for (unsigned int i=0; i::faces_per_cell;++face) + for (unsigned int face=0;face::faces_per_cell;++face) { if (cell->at_boundary(face)) continue; + const FiniteElement& nfe = cell->neighbor(face)->get_fe(); for (unsigned int i=0; i +template void -DoFTools::compute_row_length_vector(const DoFHandler &dofs, - std::vector &row_lengths, - const Coupling flux_coupling) +DoFTools::compute_row_length_vector( + const DH& 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())); @@ -264,15 +263,15 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // 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()); + 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); + const typename DH::cell_iterator end = dofs.end(); + typename DH::active_cell_iterator cell; + std::vector cell_indices; + std::vector neighbor_indices; // We loop over cells and go from // cells to lower dimensional @@ -283,6 +282,8 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // smaller than dim-1. for (cell = dofs.begin_active(); cell != end; ++cell) { + const FiniteElement& fe = cell->get_fe(); + cell_indices.resize(fe.dofs_per_cell); cell->get_dof_indices(cell_indices); unsigned int i = 0; // First, dofs on @@ -306,7 +307,8 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // dofs_per_vertex is // arbitrary, not the other way // round. - unsigned int increment = fe.dofs_per_cell - dim * fe.dofs_per_face; +//TODO: This assumes that even in hp context, the dofs per face coincide! + unsigned int increment = fe.dofs_per_cell - DH::dimension * fe.dofs_per_face; while (i < fe.first_line_index) row_lengths[cell_indices[i++]] += increment; // From now on, if an object is @@ -319,20 +321,20 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // 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; + increment = (DH::dimension>1) + ? fe.dofs_per_cell - (DH::dimension-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; + increment = (DH::dimension>2) + ? fe.dofs_per_cell - (DH::dimension-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; + 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; @@ -347,7 +349,7 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // and add the missing // contribution as well as the // flux contributions. - for (unsigned int iface=0;iface::faces_per_cell;++iface) + for (unsigned int iface=0;iface::faces_per_cell;++iface) { if (cell->at_boundary(iface)) { @@ -356,13 +358,9 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, continue; } - const typename DoFHandler::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); + const typename DH::cell_iterator neighbor = cell->neighbor(iface); + const FiniteElement& nfe = neighbor->get_fe(); + typename DH::face_iterator face = cell->face(iface); // Flux couplings are // computed from both sides @@ -374,7 +372,7 @@ DoFTools::compute_row_length_vector(const DoFHandler &dofs, // here. if (flux_coupling != none) { - unsigned int increment = fe.dofs_per_cell - fe.dofs_per_face; + unsigned int increment = nfe.dofs_per_cell - nfe.dofs_per_face; for (unsigned int i=0;i &dofs, // is refined, all the fine // face dofs couple with // the coarse one. + neighbor_indices.resize(nfe.dofs_per_cell); neighbor->get_dof_indices(neighbor_indices); for (unsigned int i=0;i &dofs, //TODO:[GK] Extend this to non-primitive elements // This is the template for 2D and 3D. See version for 1D above -template +template void DoFTools::compute_row_length_vector( - const DoFHandler& dofs, + const DH& 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(), typename FiniteElement::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. @@ -442,28 +430,16 @@ DoFTools::compute_row_length_vector( // 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()); + 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); + const typename DH::cell_iterator end = dofs.end(); + typename DH::active_cell_iterator cell; + std::vector cell_indices; + std::vector neighbor_indices; - // 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;comp& fe = cell->get_fe(); + 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(), + 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->get_dof_indices(cell_indices); + cell_indices.resize(fe.dofs_per_cell); unsigned int i = 0; // First, dofs on // vertices. We assume that @@ -502,7 +490,8 @@ DoFTools::compute_row_length_vector( for (unsigned int comp=0;comp1) - ? dofs_cell[comp] - (dim-1) * dofs_face[comp] - : dofs_cell[comp] - GeometryInfo::faces_per_cell * dofs_face[comp]; + increment = fe.base_element(fe.component_to_base(comp).first).dofs_per_cell + - ((DH::dimension>1) + ? (DH::dimension-1) + : GeometryInfo::faces_per_cell) + * fe.base_element(fe.component_to_base(comp).first).dofs_per_face; row_lengths[cell_indices[i]] += increment; } ++i; @@ -536,9 +527,11 @@ DoFTools::compute_row_length_vector( for (unsigned int comp=0;comp2) - ? dofs_cell[comp] - (dim-2) * dofs_face[comp] - : dofs_cell[comp] - GeometryInfo::faces_per_cell * dofs_face[comp]; + increment = fe.base_element(fe.component_to_base(comp).first).dofs_per_cell + - ((DH::dimension>2) + ? (DH::dimension-2) + : GeometryInfo::faces_per_cell) + * fe.base_element(fe.component_to_base(comp).first).dofs_per_face; row_lengths[cell_indices[i]] += increment; } ++i; @@ -550,7 +543,9 @@ DoFTools::compute_row_length_vector( for (unsigned int comp=0;comp::faces_per_cell * dofs_face[comp]; + increment = fe.base_element(fe.component_to_base(comp).first).dofs_per_cell + - GeometryInfo::faces_per_cell + * fe.base_element(fe.component_to_base(comp).first).dofs_per_face; row_lengths[cell_indices[i]] += increment; } ++i; @@ -567,7 +562,7 @@ DoFTools::compute_row_length_vector( // and add the missing // contribution as well as the // flux contributions. - for (unsigned int iface=0;iface::faces_per_cell;++iface) + for (unsigned int iface=0;iface::faces_per_cell;++iface) { if (cell->at_boundary(iface)) { @@ -576,13 +571,9 @@ DoFTools::compute_row_length_vector( continue; } - const typename DoFHandler::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); + const typename DH::cell_iterator neighbor = cell->neighbor(iface); + const FiniteElement& nfe = neighbor->get_fe(); + typename DH::face_iterator face = cell->face(iface); // Flux couplings are // computed from both sides @@ -596,7 +587,8 @@ DoFTools::compute_row_length_vector( for (unsigned int i=0;iget_dof_indices(neighbor_indices); for (unsigned int comp=0;comp& dofs, std::vector& row_lengths, const Coupling flux_coupling); +template void +DoFTools::compute_row_length_vector( + const hpDoFHandler& 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::compute_row_length_vector( + const hpDoFHandler& dofs, std::vector& row_lengths, + const Table<2,Coupling>& couplings, const Table<2,Coupling>& flux_couplings); + template void DoFTools::make_sparsity_pattern (const DoFHandler &dof, -- 2.39.5