From 9b30af4b37cf67e8abd51069cf7ee44f0654be79 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 9 Aug 2006 03:28:13 +0000 Subject: [PATCH] Implement VectorTools::project also for 1d. git-svn-id: https://svn.dealii.org/trunk@13620 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 64 ++------ .../include/numerics/vectors.templates.h | 153 +++++++++--------- deal.II/doc/news/changes.html | 7 + 3 files changed, 96 insertions(+), 128 deletions(-) diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index b50bf9c89b..ca7561820a 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -496,6 +496,13 @@ class VectorTools * * See the general documentation of this * class for further information. + * + * In 1d, the default value of + * the boundary quadrature + * formula is an invalid object + * since integration on the + * boundary doesn't happen in + * 1d. */ template static void project (const Mapping &mapping, @@ -505,28 +512,9 @@ class VectorTools const Function &function, VECTOR &vec, const bool enforce_zero_boundary = false, - const Quadrature &q_boundary = QGauss(2), - const bool project_to_boundary_first = false); - - /** - * Declaration of specialization - * of the previous function for - * 1d. At present, it is not - * implemented. - * - * The default value of the boundary - * quadrature formula is an invalid - * object since it makes no sense in 1d. - */ - template - static void project (const Mapping<1> &mapping, - const DoFHandler<1> &dof, - const ConstraintMatrix &constraints, - const Quadrature<1> &quadrature, - const Function<1> &function, - VECTOR &vec, - const bool enforce_zero_boundary = false, - const Quadrature<0> &q_boundary = *invalid_face_quadrature, + const Quadrature &q_boundary = (dim > 1 ? + QGauss(2) : + Quadrature(0)), const bool project_to_boundary_first = false); /** @@ -541,29 +529,11 @@ class VectorTools const Function &function, VECTOR &vec, const bool enforce_zero_boundary = false, - const Quadrature &q_boundary = QGauss(2), + const Quadrature &q_boundary = (dim > 1 ? + QGauss(2) : + Quadrature(0)), const bool project_to_boundary_first = false); - /** - * Declaration of specialization - * of the previous function for - * 1d. At present, it is not - * implemented. - * - * The default value of the boundary - * quadrature formula is an invalid - * object since it makes no sense in 1d. - */ - template - static void project (const DoFHandler<1> &dof, - const ConstraintMatrix &constraints, - const Quadrature<1> &quadrature, - const Function<1> &function, - VECTOR &vec, - const bool enforce_zero_boundary = false, - const Quadrature<0> &q_boundary = *invalid_face_quadrature, - const bool project_to_boundary_first = false); - /** * Create a right hand side * vector. Prior content of the @@ -1243,14 +1213,6 @@ class VectorTools * Exception */ DeclException0 (ExcNoComponentSelected); - - private: - /** - * Null pointer used to - * denote invalid face - * quadrature formulas in 1d. - */ - static const Quadrature<0> * const invalid_face_quadrature; }; diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 98760a90c9..e956d7248e 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -276,79 +276,41 @@ VectorTools::interpolate (const DoFHandler &dof_1, } -#if deal_II_dimension == 1 - -template -void VectorTools::project (const Mapping<1> &, - const DoFHandler<1> &, - const ConstraintMatrix &, - const Quadrature<1> &, - const Function<1> &, - VECTOR &, - const bool , - const Quadrature<0> &, - const bool ) -{ - // this function should easily be - // implemented using the template - // below. However some changes have - // to be made since faces don't - // exist in 1D. Maybe integrate the - // creation of zero boundary values - // into the project_boundary_values - // function? - Assert (false, ExcNotImplemented()); -} - - -template -void VectorTools::project (const DoFHandler<1> &dof_handler, - const ConstraintMatrix &constraints, - const Quadrature<1> &quadrature, - const Function<1> &function, - VECTOR &vec_result, - const bool enforce_zero_boundary, - const Quadrature<0> &q_boundary, - const bool project_to_boundary_first) -{ - Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1<1> mapping; - project (mapping, dof_handler, constraints, quadrature, function, vec_result, - enforce_zero_boundary, q_boundary, project_to_boundary_first); -} - - -#endif - - -template -void VectorTools::project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VECTOR &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) +namespace internal { - Assert (dof.get_fe().n_components() == function.n_components, - ExcInvalidFE()); + namespace VectorTools + { +#if deal_II_dimension == 1 - Assert (vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); - - const FiniteElement &fe = dof.get_fe(); + void + interpolate_zero_boundary_values (const ::DoFHandler<1> &dof_handler, + std::map &boundary_values) + { + // we only need to find the + // left-most and right-most + // vertex and query its vertex + // dof indices. that's easy :-) + for (unsigned direction=0; direction<2; ++direction) + { + ::DoFHandler<1>::cell_iterator + cell = dof_handler.begin(0); + while (!cell->at_boundary(direction)) + cell = cell->neighbor(direction); - // make up boundary values - std::map boundary_values; + for (unsigned int i=0; ivertex_dof_index (direction, i)] = 0.; + } + } - if (enforce_zero_boundary == true) - // no need to project boundary - // values, but enforce - // homogeneous boundary values - // anyway +#else + + template + void + interpolate_zero_boundary_values (const ::DoFHandler &dof_handler, + std::map &boundary_values) { + const FiniteElement &fe = dof_handler.get_fe(); + // loop over all boundary faces // to get all dof indices of // dofs on the boundary. note @@ -370,8 +332,9 @@ void VectorTools::project (const Mapping &mapping, // that is actually wholly on // the boundary, not only by // one line or one vertex - typename DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); + typename ::DoFHandler::active_face_iterator + face = dof_handler.begin_active_face(), + endf = dof_handler.end_face(); std::vector face_dof_indices (fe.dofs_per_face); for (; face!=endf; ++face) if (face->at_boundary()) @@ -385,24 +348,60 @@ void VectorTools::project (const Mapping &mapping, // vector valued elements here, // since we set all components boundary_values[face_dof_indices[i]] = 0.; - }; + } } + +#endif + } +} + + + +template +void VectorTools::project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VECTOR &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) +{ + Assert (dof.get_fe().n_components() == function.n_components, + ExcInvalidFE()); + + Assert (vec_result.size() == dof.n_dofs(), + ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); + + // make up boundary values + std::map boundary_values; + + if (enforce_zero_boundary == true) + // no need to project boundary + // values, but enforce + // homogeneous boundary values + // anyway + internal::VectorTools:: + interpolate_zero_boundary_values (dof, boundary_values); + else // no homogeneous boundary values if (project_to_boundary_first == true) // boundary projection required { - // set up a list of boundary functions for - // the different boundary parts. We want the - // @p{function} to hold on all parts of the - // boundary + // set up a list of boundary + // functions for the + // different boundary + // parts. We want the + // @p{function} to hold on + // all parts of the boundary typename FunctionMap::type boundary_functions; for (unsigned char c=0; c<255; ++c) boundary_functions[c] = &function; project_boundary_values (dof, boundary_functions, q_boundary, boundary_values); - }; - + } // set up mass matrix and right hand side Vector vec (dof.n_dofs()); diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index a4426aee0d..4b25b8524b 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -786,6 +786,13 @@ inconvenience this causes.

deal.II

    +
  1. + Extended: The VectorTools::project functions + are now also implemented for 1d. +
    + (WB 2006/08/08) +

    +
  2. Extended: DerivativeApproximation now offers access to the -- 2.39.5