From: wolf Date: Fri, 14 Sep 2001 10:49:13 +0000 (+0000) Subject: VectorTools::create_boundary_right_hand_side X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=14bf10b4c3becc178e98c8808be79d1ec46418b8;p=dealii-svn.git VectorTools::create_boundary_right_hand_side git-svn-id: https://svn.dealii.org/trunk@5002 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 49de131974..b3c9476462 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -14,6 +14,7 @@ #include #include #include +#include template class Function; template class FunctionMap; @@ -175,6 +176,15 @@ enum NormType { * @ref{MatrixCreator}@p{::create_*} functions which take a right hand side do, * but without assembling a matrix. * + * @item Creation of boundary right hand side vectors: The + * @p{create_boundary_right_hand_side} function computes the vector + * $f_i = \int_{\partial\Omega} g(x) \phi_i(x) dx$. This is the + * right hand side contribution of boundary forces when having + * inhomogeneous Neumann boundary values in Laplace's equation or + * other second order operators. This function also takes an + * optional argument denoting over which parts of the boundary the + * integration shall extend. + * * @item Interpolation of boundary values: * The @ref{MatrixTools}@p{::apply_boundary_values} function takes a list * of boundary nodes and their values. You can get such a list by interpolation @@ -438,7 +448,10 @@ class VectorTools const bool project_to_boundary_first = false); /** - * Create a right hand side vector. + * Create a right hand side + * vector. Prior content of the + * given @p{rhs_vector} vector is + * deleted. * * See the general documentation of this * class for further information. @@ -461,6 +474,51 @@ class VectorTools const Function &rhs, Vector &rhs_vector); + /** + * Create a right hand side + * vector from boundary + * forces. Prior content of the + * given @p{rhs_vector} vector is + * deleted. + * + * See the general documentation of this + * class for further information. + */ + template + static void create_boundary_right_hand_side (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &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 Mapping<1> &mapping, + const DoFHandler<1> &dof, + const Quadrature<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 + * @p{mapping=MappingQ1()}. + */ + template + static void create_boundary_right_hand_side (const DoFHandler &dof, + const Quadrature &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/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 6893354929..d13e051337 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -545,6 +545,139 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, } + +#if deal_II_dimension != 1 + +template +void +VectorTools::create_boundary_right_hand_side (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + const FiniteElement &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.clear(); + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + FEFaceValues fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points, + n_components = fe.n_components(); + + 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(); + + if (n_components==1) + { + std::vector rhs_values(n_q_points); + + 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())) + { + fe_values.reinit(cell, face); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector.clear(); + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); + + 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())) + { + fe_values.reinit(cell, face); + + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector.clear(); + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i &mapping, + const DoFHandler<1> &dof_handler, + const Quadrature<0> &quadrature, + const Function<1> &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) +{ + Assert (false, ExcImpossibleInDim(1)); +}; + +#endif + +template +void +VectorTools::create_boundary_right_hand_side (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + 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, + rhs_function, rhs_vector, + boundary_indicators); +} + + + #if deal_II_dimension == 1 void @@ -1361,6 +1494,25 @@ void VectorTools::create_right_hand_side (const DoFHandler &, const Quadrature &, const Function &, Vector &); + +#if deal_II_dimension != 1 +template +void +VectorTools::create_boundary_right_hand_side (const Mapping &, + const DoFHandler &, + const Quadrature &, + const Function &, + Vector &, + const std::set &); +#endif + +template +void +VectorTools::create_boundary_right_hand_side (const DoFHandler &, + const Quadrature &, + const Function &, + Vector &, + const std::set &); template void VectorTools::interpolate_boundary_values (const DoFHandler &, const unsigned char, diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 342dd22fc5..c4eea11daf 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -596,6 +596,13 @@ documentation, etc.

deal.II

    +

  1. + New: VectorTools::create_boundary_right_hand_side + integrates boundary forces for inhomogeneous Neumann boundary values. +
    + (WB 2001/09/13) +

    +

  2. New: DoFTools::make_flux_sparsity_pattern now exists also for 1d.