From: buerg Date: Sun, 6 May 2012 10:03:52 +0000 (+0000) Subject: Add create_point_source_vector for vector-valued FE. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0bcf2d26a2cf6155577cc4d07629415c6b7f11f8;p=dealii-svn.git Add create_point_source_vector for vector-valued FE. git-svn-id: https://svn.dealii.org/trunk@25497 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h index f2823f9300..f1fc9bf688 100644 --- a/deal.II/include/deal.II/numerics/vectors.h +++ b/deal.II/include/deal.II/numerics/vectors.h @@ -1714,6 +1714,63 @@ namespace VectorTools const Point &p, Vector &rhs_vector); + /** + * Create a right hand side + * vector for a point source at point @p p + * for vector-valued finite elements. + * Prior content of the + * given @p rhs_vector vector is + * deleted. + * + * See the general documentation of this + * class for further information. + */ + template + void create_point_source_vector(const Mapping &mapping, + const DoFHandler &dof, + const Point &p, + const Point &orientation, + Vector &rhs_vector); + + /** + * Calls the create_point_source_vector() + * function for vector-valued finite elements, + * see above, with + * mapping=MappingQ1@(). + */ + template + void create_point_source_vector(const DoFHandler &dof, + const Point &p, + const Point &orientation, + Vector &rhs_vector); + + /** + * Like the previous set of functions, + * but for hp objects. + */ + template + void create_point_source_vector(const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const Point &p, + const Point &orientation, + 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 + void create_point_source_vector(const hp::DoFHandler &dof, + const Point &p, + const Point &orientation, + Vector &rhs_vector); + /** * Create a right hand side * vector from boundary diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 2fb9e657db..72b02f2b7b 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -1008,6 +1008,105 @@ namespace VectorTools + + template + void create_point_source_vector (const Mapping &mapping, + const DoFHandler &dof_handler, + const Point &p, + const Point &orientation, + 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() == dim, + ExcMessage ("This function only works for vector-valued 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)); + + const FEValuesExtractors::Vector vec (0); + 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; + + std::vector local_dof_indices (dofs_per_cell); + cell_point.first->get_dof_indices (local_dof_indices); + + for(unsigned int i=0; i + void create_point_source_vector (const DoFHandler &dof_handler, + const Point &p, + const Point &orientation, + Vector &rhs_vector) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_point_source_vector(StaticMappingQ1::mapping, dof_handler, + p, orientation, rhs_vector); + } + + + template + void create_point_source_vector (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const Point &p, + const Point &orientation, + 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() == dim, + ExcMessage ("This function only works for vector-valued 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)); + + const FEValuesExtractors::Vector vec (0); + 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 create_point_source_vector (const hp::DoFHandler &dof_handler, + const Point &p, + const Point &orientation, + Vector &rhs_vector) + { + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + create_point_source_vector(hp::StaticMappingQ1::mapping_collection, + dof_handler, + p, orientation, rhs_vector); + } + + + // separate implementation for 1D because otherwise we get linker errors since // FEFaceValues<1> is not compiled template <> diff --git a/deal.II/source/numerics/vectors.inst.in b/deal.II/source/numerics/vectors.inst.in index 1062bbc524..c7fd3f27a1 100644 --- a/deal.II/source/numerics/vectors.inst.in +++ b/deal.II/source/numerics/vectors.inst.in @@ -462,6 +462,33 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (const hp::DoFHandler &, const Point &, Vector &); + template + void create_point_source_vector + (const Mapping &, + const DoFHandler &, + const Point &, + const Point &, + Vector &); + template + void create_point_source_vector + (const DoFHandler &, + const Point &, + const Point &, + Vector &); + + template + void create_point_source_vector + (const hp::MappingCollection &, + const hp::DoFHandler &, + const Point &, + const Point &, + Vector &); + template + void create_point_source_vector + (const hp::DoFHandler &, + const Point &, + const Point &, + Vector &); #if deal_II_dimension > 1 template