From eedbb402bc621abaeafba02ee57f9fc556e2c83f Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 11 Jan 2007 03:41:18 +0000 Subject: [PATCH] Implement a whole host of new functions for hp as well. git-svn-id: https://svn.dealii.org/trunk@14296 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 89 +++- .../include/numerics/vectors.templates.h | 447 ++++++++++++++++-- deal.II/deal.II/source/numerics/vectors.cc | 50 +- 3 files changed, 540 insertions(+), 46 deletions(-) diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index f823aa6f36..ee1d1cb6d2 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -567,6 +567,27 @@ class VectorTools const Function &rhs, Vector &rhs_vector); + /** + * Like the previous set of functions, + * but for hp objects. + */ + template + static void create_right_hand_side (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + const Function &rhs, + Vector &rhs_vector); + + /** + * Like the previous set of functions, + * but for hp objects. + */ + template + static void create_right_hand_side (const hp::DoFHandler &dof, + const hp::QCollection &q, + const Function &rhs, + Vector &rhs_vector); + /** * Create a right hand side * vector for a point source at point @p p. @@ -593,6 +614,31 @@ class VectorTools const Point &p, Vector &rhs_vector); + /** + * Like the previous set of functions, + * but for hp objects. + */ + template + static void create_point_source_vector(const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const Point &p, + Vector &rhs_vector); + + /** + * Like the previous set of functions, + * but for hp objects. The function uses + * the default Q1 mapping object. Note + * that if your hp::DoFHandler uses any + * active fe index other than zero, then + * you need to call the function above + * that provides a mapping object for + * each active fe index. + */ + template + static void create_point_source_vector(const hp::DoFHandler &dof, + const Point &p, + Vector &rhs_vector); + /** * Create a right hand side * vector from boundary @@ -638,6 +684,47 @@ class VectorTools Vector &rhs_vector, const std::set &boundary_indicators = std::set()); + /** + * Same as the set of functions above, + * but for hp objects. + */ + template + static void create_boundary_right_hand_side (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, + const Function &rhs, + Vector &rhs_vector, + const std::set &boundary_indicators = std::set()); + + /** + * Specialization of above + * function for 1d. Since the + * computation is not useful in + * 1d, this function simply + * throws an exception. + */ + static void create_boundary_right_hand_side (const hp::MappingCollection<1> &mapping, + const hp::DoFHandler<1> &dof, + const hp::QCollection<0> &q, + const Function<1> &rhs, + Vector &rhs_vector, + const std::set &boundary_indicators = std::set()); + + /** + * Calls the @p + * create_boundary_right_hand_side + * function, see above, with a single Q1 + * mapping as collection. This function + * therefore will only work if the only + * active fe index in use is zero. + */ + template + static void create_boundary_right_hand_side (const hp::DoFHandler &dof, + const hp::QCollection &q, + const Function &rhs, + Vector &rhs_vector, + const std::set &boundary_indicators = std::set()); + /** * Prepare Dirichlet boundary * conditions. Make up the list diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 9b9004629d..e1540eb8ed 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name: $ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -43,7 +43,7 @@ DEAL_II_NAMESPACE_OPEN -//TODO[RH]: Use StaticQ1Mapping where appropriate +//TODO[RH]: Use StaticMappingQ1 where appropriate template void VectorTools::interpolate (const Mapping &mapping, @@ -524,8 +524,9 @@ void VectorTools::create_right_hand_side (const Mapping &mapping, std::vector dofs (dofs_per_cell); Vector cell_vector (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); if (n_components==1) { @@ -536,7 +537,8 @@ void VectorTools::create_right_hand_side (const Mapping &mapping, fe_values.reinit(cell); const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); cell_vector = 0; for (unsigned int point=0; point &mapping, } else { - std::vector > rhs_values(n_q_points, Vector(n_components)); + std::vector > rhs_values(n_q_points, + Vector(n_components)); - // Use the faster code if the FiniteElement is primitive - if (fe.is_primitive ()) + for (; cell!=endc; ++cell) { - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), + rhs_values); - cell_vector = 0; + cell_vector = 0; + // Use the faster code if the + // FiniteElement is primitive + if (fe.is_primitive ()) + { for (unsigned int point=0; point &mapping, fe_values.shape_value(i,point) * weights[point]; } - - cell->get_dof_indices (dofs); - - for (unsigned int i=0; i &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; + // Otherwise do it the way + // proposed for vector valued + // elements for (unsigned int point=0; point &mapping, fe_values.shape_value_component(i,point,comp_i) * weights[point]; } - - cell->get_dof_indices (dofs); - - for (unsigned int i=0; iget_dof_indices (dofs); + + for (unsigned int i=0; i void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, const Quadrature &quadrature, @@ -627,6 +621,148 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, } + + +template +void VectorTools::create_right_hand_side (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector) +{ + const hp::FECollection &fe = dof_handler.get_fe(); + Assert (fe.n_components() == rhs_function.n_components, + ExcComponentMismatch()); + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + rhs_vector = 0; + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + hp::FEValues x_fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int n_components = fe.n_components(); + + std::vector dofs (fe.max_dofs_per_cell()); + Vector cell_vector (fe.max_dofs_per_cell()); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + if (n_components==1) + { + std::vector rhs_values; + + for (; cell!=endc; ++cell) + { + x_fe_values.reinit(cell); + + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points); + dofs.resize (dofs_per_cell); + cell_vector.reinit (dofs_per_cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values; + + for (; cell!=endc; ++cell) + { + x_fe_values.reinit(cell); + + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points, + Vector(n_components)); + dofs.resize (dofs_per_cell); + cell_vector.reinit (dofs_per_cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), + rhs_values); + + cell_vector = 0; + + // Use the faster code if the + // FiniteElement is primitive + if (cell->get_fe().is_primitive ()) + { + for (unsigned int point=0; pointget_fe().system_to_component_index(i).first; + + cell_vector(i) += rhs_values[point](component) * + fe_values.shape_value(i,point) * + weights[point]; + } + } + else + { + // Otherwise do it the way proposed + // for vector valued elements + for (unsigned int point=0; pointget_fe().get_nonzero_components(i)[comp_i]) + { + cell_vector(i) += rhs_values[point](comp_i) * + fe_values.shape_value_component(i,point,comp_i) * + weights[point]; + } + } + + cell->get_dof_indices (dofs); + + for (unsigned int i=0; i +void VectorTools::create_right_hand_side (const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_right_hand_side(hp::StaticMappingQ1::mapping_collection, + dof_handler, quadrature, + rhs_function, rhs_vector); +} + + + + template void VectorTools::create_point_source_vector (const Mapping &mapping, const DoFHandler &dof_handler, @@ -635,6 +771,9 @@ void VectorTools::create_point_source_vector (const Mapping &mapping, { Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + Assert (dof_handler.get_fe().n_components() == 1, + ExcMessage ("This function only works for scalar finite elements")); + rhs_vector = 0; std::pair::active_cell_iterator, Point > @@ -643,7 +782,8 @@ void VectorTools::create_point_source_vector (const Mapping &mapping, Quadrature q(GeometryInfo::project_to_unit_cell(cell_point.second)); - FEValues fe_values(dof_handler.get_fe(), q, UpdateFlags(update_values)); + FEValues fe_values(mapping, dof_handler.get_fe(), + q, UpdateFlags(update_values)); fe_values.reinit(cell_point.first); const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; @@ -655,6 +795,8 @@ void VectorTools::create_point_source_vector (const Mapping &mapping, rhs_vector(local_dof_indices[i]) = fe_values.shape_value(i,0); } + + template void VectorTools::create_point_source_vector (const DoFHandler &dof_handler, const Point &p, @@ -665,6 +807,53 @@ void VectorTools::create_point_source_vector (const DoFHandler &dof_hand p, rhs_vector); } + +template +void VectorTools::create_point_source_vector (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) +{ + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + Assert (dof_handler.get_fe().n_components() == 1, + ExcMessage ("This function only works for scalar finite elements")); + + rhs_vector = 0; + + std::pair::active_cell_iterator, Point > + cell_point = + GridTools::find_active_cell_around_point (mapping, dof_handler, p); + + Quadrature q(GeometryInfo::project_to_unit_cell(cell_point.second)); + + FEValues fe_values(mapping[cell_point.first->active_fe_index()], + cell_point.first->get_fe(), q, UpdateFlags(update_values)); + fe_values.reinit(cell_point.first); + + const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell; + + std::vector local_dof_indices (dofs_per_cell); + cell_point.first->get_dof_indices (local_dof_indices); + + for(unsigned int i=0; i +void VectorTools::create_point_source_vector (const hp::DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_point_source_vector(hp::StaticMappingQ1::mapping_collection, + dof_handler, + p, rhs_vector); +} + + #if deal_II_dimension != 1 template @@ -746,7 +935,8 @@ VectorTools::create_boundary_right_hand_side (const Mapping &mapping, cell_vector = 0; - // Use the faster code if the FiniteElement is primitive + // Use the faster code if the + // FiniteElement is primitive if (fe.is_primitive ()) { for (unsigned int point=0; point &mapping, } else { - // And the full featured code, if vector valued FEs are used + // And the full featured + // code, if vector valued + // FEs are used for (unsigned int point=0; point &dof_handl const std::set &boundary_indicators) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - create_boundary_right_hand_side(mapping, dof_handler, quadrature, + + create_boundary_right_hand_side(StaticMappingQ1::mapping, dof_handler, + quadrature, rhs_function, rhs_vector, boundary_indicators); } +#if deal_II_dimension != 1 + +template +void +VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + const hp::FECollection &fe = dof_handler.get_fe(); + Assert (fe.n_components() == rhs_function.n_components, + ExcComponentMismatch()); + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + + rhs_vector = 0; + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + hp::FEFaceValues x_fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int n_components = fe.n_components(); + + std::vector dofs (fe.max_dofs_per_cell()); + Vector cell_vector (fe.max_dofs_per_cell()); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + if (n_components==1) + { + std::vector rhs_values; + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary () && + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != + boundary_indicators.end())) + { + x_fe_values.reinit(cell, face); + + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values; + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary () && + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != + boundary_indicators.end())) + { + x_fe_values.reinit(cell, face); + + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + rhs_values.resize (n_q_points, Vector(n_components)); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + + // Use the faster code if the + // FiniteElement is primitive + if (cell->get_fe().is_primitive ()) + { + for (unsigned int point=0; pointget_fe().system_to_component_index(i).first; + + cell_vector(i) += rhs_values[point](component) * + fe_values.shape_value(i,point) * + weights[point]; + } + } + else + { + // And the full featured + // code, if vector valued + // FEs are used + for (unsigned int point=0; pointget_fe().get_nonzero_components(i)[comp_i]) + { + cell_vector(i) + += rhs_values[point](comp_i) * + fe_values.shape_value_component(i,point,comp_i) * + weights[point]; + } + } + + cell->get_dof_indices (dofs); + + for (unsigned int i=0; i &, + const hp::DoFHandler<1> &, + const hp::QCollection<0> &, + const Function<1> &, + Vector &, + const std::set &) +{ + Assert (false, ExcImpossibleInDim(1)); +} + +#endif + +template +void +VectorTools::create_boundary_right_hand_side (const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_boundary_right_hand_side(hp::StaticMappingQ1::mapping_collection, + dof_handler, quadrature, + rhs_function, rhs_vector, + boundary_indicators); +} + + + + #if deal_II_dimension == 1 //TODO[?] Actually the Mapping object should be a MappingCollection object for the diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index f096705170..bc888cdd07 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -50,6 +50,21 @@ void VectorTools::create_right_hand_side const Quadrature &, const Function &, Vector &); + +template +void VectorTools::create_right_hand_side +(const hp::MappingCollection &, + const hp::DoFHandler &, + const hp::QCollection &, + const Function &, + Vector &); +template +void VectorTools::create_right_hand_side +(const hp::DoFHandler &, + const hp::QCollection &, + const Function &, + Vector &); + template void VectorTools::create_point_source_vector (const Mapping &, @@ -62,6 +77,18 @@ void VectorTools::create_point_source_vector const Point &, Vector &); +template +void VectorTools::create_point_source_vector +(const hp::MappingCollection &, + const hp::DoFHandler &, + const Point &, + Vector &); +template +void VectorTools::create_point_source_vector +(const hp::DoFHandler &, + const Point &, + Vector &); + #if deal_II_dimension != 1 template void @@ -83,6 +110,27 @@ VectorTools::create_boundary_right_hand_side Vector &, const std::set &); +#if deal_II_dimension != 1 +template +void +VectorTools::create_boundary_right_hand_side +(const hp::MappingCollection &, + const hp::DoFHandler &, + const hp::QCollection &, + const Function &, + Vector &, + const std::set &); +#endif + +template +void +VectorTools::create_boundary_right_hand_side +(const hp::DoFHandler &, + const hp::QCollection &, + const Function &, + Vector &, + const std::set &); + template void VectorTools::interpolate_boundary_values ( const DoFHandler &, -- 2.39.5