const Point<spacedim> &p,
Vector<double> &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 <int dim, int spacedim>
+ void create_point_source_vector(const Mapping<dim,spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &rhs_vector);
+
+ /**
+ * Calls the create_point_source_vector()
+ * function for vector-valued finite elements,
+ * see above, with
+ * <tt>mapping=MappingQ1@<dim@>()</tt>.
+ */
+ template <int dim, int spacedim>
+ void create_point_source_vector(const DoFHandler<dim,spacedim> &dof,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &rhs_vector);
+
+ /**
+ * Like the previous set of functions,
+ * but for hp objects.
+ */
+ template <int dim, int spacedim>
+ void create_point_source_vector(const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &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 <int dim, int spacedim>
+ void create_point_source_vector(const hp::DoFHandler<dim,spacedim> &dof,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &rhs_vector);
+
/**
* Create a right hand side
* vector from boundary
+
+ template <int dim, int spacedim>
+ void create_point_source_vector (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof_handler,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &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<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+ cell_point =
+ GridTools::find_active_cell_around_point (mapping, dof_handler, p);
+
+ Quadrature<dim> q(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+
+ const FEValuesExtractors::Vector vec (0);
+ FEValues<dim,spacedim> 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<unsigned int> local_dof_indices (dofs_per_cell);
+ cell_point.first->get_dof_indices (local_dof_indices);
+
+ for(unsigned int i=0; i<dofs_per_cell; i++)
+ rhs_vector(local_dof_indices[i]) = orientation * fe_values[vec].value(i,0);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void create_point_source_vector (const DoFHandler<dim,spacedim> &dof_handler,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &rhs_vector)
+ {
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_point_source_vector(StaticMappingQ1<dim,spacedim>::mapping, dof_handler,
+ p, orientation, rhs_vector);
+ }
+
+
+ template <int dim, int spacedim>
+ void create_point_source_vector (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof_handler,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &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<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+ cell_point =
+ GridTools::find_active_cell_around_point (mapping, dof_handler, p);
+
+ Quadrature<dim> q(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+
+ const FEValuesExtractors::Vector vec (0);
+ FEValues<dim> 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<unsigned int> local_dof_indices (dofs_per_cell);
+ cell_point.first->get_dof_indices (local_dof_indices);
+
+ for(unsigned int i=0; i<dofs_per_cell; i++)
+ rhs_vector(local_dof_indices[i]) = orientation * fe_values[vec].value(i,0);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void create_point_source_vector (const hp::DoFHandler<dim,spacedim> &dof_handler,
+ const Point<spacedim> &p,
+ const Point<dim> &orientation,
+ Vector<double> &rhs_vector)
+ {
+ Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+ create_point_source_vector(hp::StaticMappingQ1<dim>::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 <>
(const hp::DoFHandler<deal_II_dimension> &,
const Point<deal_II_dimension> &,
Vector<double> &);
+ template
+ void create_point_source_vector<deal_II_dimension>
+ (const Mapping<deal_II_dimension> &,
+ const DoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ Vector<double> &);
+ template
+ void create_point_source_vector<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ Vector<double> &);
+
+ template
+ void create_point_source_vector<deal_II_dimension>
+ (const hp::MappingCollection<deal_II_dimension> &,
+ const hp::DoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ Vector<double> &);
+ template
+ void create_point_source_vector<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ const Point<deal_II_dimension> &,
+ Vector<double> &);
#if deal_II_dimension > 1
template