From: kanschat Date: Fri, 30 May 2008 16:19:48 +0000 (+0000) Subject: allow component mapping in VectorTools::project_boundary_values() X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a62fc062ba9595e756e5ba7e08b43ef5738ff3a1;p=dealii-svn.git allow component mapping in VectorTools::project_boundary_values() remove some 1-d specializations in favor of more readable documentation git-svn-id: https://svn.dealii.org/trunk@16274 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 1f1f5f6735..b4c05d17c7 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -309,9 +309,7 @@ class MatrixCreator /** * Assemble the mass matrix and a * right hand side vector along - * the boundary. If no - * coefficient is given, it is - * assumed to be constant one. + * the boundary. * * The matrix is assumed to * already be initialized with a @@ -323,6 +321,26 @@ class MatrixCreator * to use multithreading, this * function works in parallel. * + * @arg @p weight: an optional + * weight for the computation of + * the mass matrix. If no weight + * is given, it is set to one. + * + * @arg @p component_mapping: if + * the components in @p + * boundary_functions and @p dof + * do not coincide, this vector + * allows them to be remapped. If + * the vector is not empty, it + * has to have one entry for each + * component in @p dof. This + * entry is the component number + * in @p boundary_functions that + * should be used for this + * component in @p dof. By + * default, no remapping is + * applied. + * * @todo This function does not * work for finite elements with * cell-dependent shape functions. @@ -336,22 +354,9 @@ class MatrixCreator const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, - const Function * const a = 0); - - /** - * Same function, but for 1d. - */ - static - void create_boundary_mass_matrix (const Mapping<1> &mapping, - const DoFHandler<1> &dof, - const Quadrature<0> &q, - SparseMatrix &matrix, - const FunctionMap<1>::type &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function<1> * const a = 0); - - + const Function * const weight = 0, + std::vector component_mapping = std::vector()); + /** * Calls the * create_boundary_mass_matrix() @@ -366,7 +371,8 @@ class MatrixCreator const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, - const Function * const a = 0); + const Function * const a = 0, + std::vector component_mapping = std::vector()); /** * Same function as above, but for hp @@ -381,21 +387,8 @@ class MatrixCreator const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, - const Function * const a = 0); - - /** - * Same function as above, but for hp - * objects. - */ - static - void create_boundary_mass_matrix (const hp::MappingCollection<1> &mapping, - const hp::DoFHandler<1> &dof, - const hp::QCollection<0> &q, - SparseMatrix &matrix, - const FunctionMap<1>::type &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function<1> * const a = 0); + const Function * const a = 0, + std::vector component_mapping = std::vector()); /** * Same function as above, but for hp @@ -409,7 +402,8 @@ class MatrixCreator const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, - const Function * const a = 0); + const Function * const a = 0, + std::vector component_mapping = std::vector()); /** * Assemble the Laplace @@ -746,14 +740,15 @@ class MatrixCreator */ template static - void create_boundary_mass_matrix_1 (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, + void create_boundary_mass_matrix_1 (boost::tuple&, + const DoFHandler&, + const Quadrature&> commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, const Function * const a, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex); @@ -762,14 +757,15 @@ class MatrixCreator */ template static - void create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, + void create_boundary_mass_matrix_1 (boost::tuple&, + const hp::DoFHandler&, + const hp::QCollection&> commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector&dof_to_boundary_mapping, const Function * const a, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex); }; diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index f19fa4e131..49285ed8d2 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -96,19 +96,22 @@ class ConstraintMatrix; * given function may be, taking into account that a virtual function has * to be called. * - *
  • Projection: compute the $L^2$-projection of the given function - * onto the finite element space, i.e. if $f$ is the function to be - * projected, compute $f_h\in V_h$ such that $(f_h,v_h)=(f,v_h)$ for - * all discrete test functions $v_h$. This is done through the - * solution of the linear system of equations $M v = f$ where $M$ is - * the mass matrix $m_{ij} = \int_\Omega \phi_i(x) \phi_j(x) dx$ and - * $f_i = \int_\Omega f(x) \phi_i(x) dx$. The solution vector $v$ then - * is the nodal representation of the projection $f_h$. The project() - * functions are used in the @ref step_23 "step-23" tutorial program. + *
  • Projection: compute the L2-projection of the + * given function onto the finite element space, i.e. if f is + * the function to be projected, compute fh in + * Vh such that + * (fh,vh)=(f,vh)$ + * for all discrete test functions vh. This is done + * through the solution of the linear system of equations M v = + * f where M is the mass matrix $m_{ij} = \int_\Omega + * \phi_i(x) \phi_j(x) dx$ and $f_i = \int_\Omega f(x) \phi_i(x) + * dx$. The solution vector $v$ then is the nodal representation of + * the projection fh. The project() functions are + * used in the @ref step_23 "step-23" tutorial program. * * In order to get proper results, it be may necessary to treat * boundary conditions right. Below are listed some cases where this - * may be needed. If needed, this is done by $L^2$-projection of + * may be needed. If needed, this is done by L2-projection of * the trace of the given function onto the finite element space * restricted to the boundary of the domain, then taking this * information and using it to eliminate the boundary nodes from the @@ -138,10 +141,10 @@ class ConstraintMatrix; * * Obviously, the results of the two schemes for projection are * different. Usually, when projecting to the boundary first, the - * $L^2$-norm of the difference between original + * L2-norm of the difference between original * function and projection over the whole domain will be larger * (factors of five have been observed) while the - * $L^2$-norm of the error integrated over the + * L2-norm of the error integrated over the * boundary should of course be less. The reverse should also hold * if no projection to the boundary is performed. * @@ -165,17 +168,17 @@ class ConstraintMatrix; * detail may change in the future. * *
  • Creation of right hand side vectors: - * The @p create_right_hand_side function computes the vector + * The create_right_hand_side() function computes the vector * $f_i = \int_\Omega f(x) \phi_i(x) dx$. This is the same as what the * MatrixCreator::create_* functions which take a right hand side do, * but without assembling a matrix. * *
  • Creation of right hand side vectors for point sources: - * The @p create_point_source_vector function computes the vector + * The create_point_source_vector() function computes the vector * $f_i = \int_\Omega \delta_0(x-x_0) \phi_i(x) dx$. * *
  • Creation of boundary right hand side vectors: The - * @p create_boundary_right_hand_side function computes the vector + * 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 @@ -184,9 +187,9 @@ class ConstraintMatrix; * integration shall extend. * *
  • Interpolation of boundary values: - * The MatrixTools@p ::apply_boundary_values function takes a list + * The MatrixTools::apply_boundary_values() function takes a list * of boundary nodes and their values. You can get such a list by interpolation - * of a boundary function using the @p interpolate_boundary_values function. + * of a boundary function using the interpolate_boundary_values() function. * To use it, you have to * specify a list of pairs of boundary indicators (of type unsigned char; * see the section in the documentation of the Triangulation class for more @@ -200,7 +203,7 @@ class ConstraintMatrix; * Within this function, boundary values are interpolated, i.e. a node is given * the point value of the boundary function. In some cases, it may be necessary * to use the L2-projection of the boundary function or any other method. For - * this purpose we refer to the VectorTools@p ::project_boundary_values + * this purpose we refer to the project_boundary_values() * function below. * * You should be aware that the boundary function may be evaluated at nodes @@ -219,24 +222,24 @@ class ConstraintMatrix; * the place of the respective boundary point. * *
  • Projection of boundary values: - * The @p project_boundary_values function acts similar to the - * @p interpolate_boundary_values function, apart from the fact that it does + * The project_boundary_values() function acts similar to the + * interpolate_boundary_values() function, apart from the fact that it does * not get the nodal values of boundary nodes by interpolation but rather - * through the $L^2$-projection of the trace of the function to the boundary. + * through the L2-projection of the trace of the function to the boundary. * * The projection takes place on all boundary parts with boundary - * indicators listed in the map (FunctioMap@p ::FunctionMap) + * indicators listed in the map (FunctioMap::FunctionMap) * of boundary functions. These boundary parts may or may not be * continuous. For these boundary parts, the mass matrix is * assembled using the - * MatrixTools@p ::create_boundary_mass_matrix function, as + * MatrixTools::create_boundary_mass_matrix() function, as * well as the appropriate right hand side. Then the resulting * system of equations is solved using a simple CG method (without * preconditioning), which is in most cases sufficient for the * present purpose. * *
  • Computing errors: - * The function @p integrate_difference performs the calculation of + * The function integrate_difference() performs the calculation of * the error between a given (continuous) reference function and the * finite element solution in different norms. The integration is * performed using a given quadrature formula and assumes that the @@ -256,7 +259,7 @@ class ConstraintMatrix; * * Presently, there is the possibility to compute the following values from the * difference, on each cell: @p mean, @p L1_norm, @p L2_norm, @p Linfty_norm, - * @p H1_seminorm and @p H1_norm, see @p VectorTools::NormType. + * @p H1_seminorm and @p H1_norm, see VectorTools::NormType. * For the mean difference value, the reference function minus the numerical * solution is computed, not the other way round. * @@ -279,27 +282,27 @@ class ConstraintMatrix; * use of the wrong quadrature formula may show a significantly wrong result * and care should be taken to chose the right formula. * - * The $H^1$ seminorm is the $L^2$ + * The H1 seminorm is the L2 * norm of the gradient of the difference. The square of the full - * $H^1$ norm is the sum of the square of seminorm - * and the square of the $L^2$ norm. + * H1 norm is the sum of the square of seminorm + * and the square of the L2 norm. * * To get the global L1 error, you have to sum up the * entries in @p difference, e.g. using - * Vector::l1_norm function. For the global $L^2$ + * Vector::l1_norm() function. For the global L2 * difference, you have to sum up the squares of the entries and * take the root of the sum, e.g. using - * Vector::l2_norm. These two operations - * represent the $l_1$ and $l_2$ norms of the vectors, but you need + * Vector::l2_norm(). These two operations + * represent the l1 and l2 norms of the vectors, but you need * not take the absolute value of each entry, since the cellwise * norms are already positive. * * To get the global mean difference, simply sum up the elements as above. * To get the $L_\infty$ norm, take the maximum of the vector elements, e.g. - * using the Vector::linfty_norm function. + * using the Vector::linfty_norm() function. * - * For the global $H^1$ norm and seminorm, the same rule applies as for the - * $L^2$ norm: compute the $l_2$ norm of the cell error vector. + * For the global H1 norm and seminorm, the same rule applies as for the + * L2 norm: compute the l2 norm of the cell error vector. * * * All functions use the finite element given to the DoFHandler object the last @@ -438,7 +441,7 @@ class VectorTools VECTOR &vec); /** - * Calls the @p interpolate + * Calls the @p interpolate() * function above with * mapping=MappingQ1@@(). */ @@ -652,7 +655,7 @@ class VectorTools /** * Calls the other - * @p interpolate_boundary_values + * interpolate_boundary_values() * function, see above, with * mapping=MappingQ1@(). */ @@ -668,7 +671,7 @@ class VectorTools /** * Calls the other - * @p interpolate_boundary_values + * interpolate_boundary_values() * function, see above, with * mapping=MappingQ1@(). */ @@ -682,7 +685,7 @@ class VectorTools /** - * Project @p function to the boundary + * Project a function to the boundary * of the domain, using the given quadrature * formula for the faces. If the * @p boundary_values contained values @@ -691,29 +694,53 @@ class VectorTools * of the boundary part to be projected * on already was in the variable. * - * It is assumed that the number - * of components of - * @p boundary_function - * matches that of the finite - * element used by @p dof. + * If @p component_mapping is + * empty, it is assumed that the + * number of components of @p + * boundary_function matches that + * of the finite element used by + * @p dof. * * In 1d, projection equals * interpolation. Therefore, * interpolate_boundary_values is * called. * - * See the general documentation of this - * class for further information. + * @arg @p boundary_values: the + * result of this function, a map + * containing all indices of + * degrees of freedom at the + * boundary (as covered by the + * boundary parts in @p + * boundary_functions) and the + * computed dof value for this + * degree of freedom. + * + * @arg @p component_mapping: if + * the components in @p + * boundary_functions and @p dof + * do not coincide, this vector + * allows them to be remapped. If + * the vector is not empty, it + * has to have one entry for each + * component in @p dof. This + * entry is the component number + * in @p boundary_functions that + * should be used for this + * component in @p dof. By + * default, no remapping is + * applied. */ template static void project_boundary_values (const Mapping &mapping, const DoFHandler &dof, const typename FunctionMap::type &boundary_functions, const Quadrature &q, - std::map &boundary_values); + std::map &boundary_values, + std::vector component_mapping = std::vector()); /** - * Calls the @p project_boundary_values + * Calls the project_boundary_values() * function, see above, with * mapping=MappingQ1@(). */ @@ -721,7 +748,8 @@ class VectorTools static void project_boundary_values (const DoFHandler &dof, const typename FunctionMap::type &boundary_function, const Quadrature &q, - std::map &boundary_values); + std::map &boundary_values, + std::vector component_mapping = std::vector()); /** @@ -1053,7 +1081,7 @@ class VectorTools Vector &rhs_vector); /** - * Calls the @p create_right_hand_side + * Calls the create_right_hand_side() * function, see above, with * mapping=MappingQ1@(). */ @@ -1101,7 +1129,7 @@ class VectorTools Vector &rhs_vector); /** - * Calls the @p create_point_source_vector + * Calls the create_point_source_vector() * function, see above, with * mapping=MappingQ1@(). */ @@ -1155,7 +1183,7 @@ class VectorTools /** * Calls the - * @p create_boundary_right_hand_side + * create_boundary_right_hand_side() * function, see above, with * mapping=MappingQ1@(). */ @@ -1179,12 +1207,14 @@ class VectorTools 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. + * Calls the + * 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, @@ -1278,7 +1308,7 @@ class VectorTools const double exponent = 2.); /** - * Calls the @p integrate_difference + * Calls the integrate_difference() * function, see above, with * mapping=MappingQ1@(). */ @@ -1304,7 +1334,7 @@ class VectorTools const double exponent = 2.); /** - * Calls the @p integrate_difference + * Calls the integrate_difference() * function, see above, with * mapping=MappingQ1@(). */ @@ -1530,7 +1560,7 @@ class VectorTools const unsigned int component); /** - * Calls the @p compute_mean_value + * Calls the compute_mean_value() * function, see above, with * mapping=MappingQ1@(). */ diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index afabd9eb88..d330023d8b 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -311,8 +311,9 @@ namespace internal { #if deal_II_dimension == 1 + template void - interpolate_zero_boundary_values (const dealii::DoFHandler<1> &dof_handler, + interpolate_zero_boundary_values (const dealii::DoFHandler &dof_handler, std::map &boundary_values) { // we only need to find the @@ -321,7 +322,7 @@ namespace internal // dof indices. that's easy :-) for (unsigned direction=0; direction<2; ++direction) { - dealii::DoFHandler<1>::cell_iterator + typename dealii::DoFHandler::cell_iterator cell = dof_handler.begin(0); while (!cell->at_boundary(direction)) cell = cell->neighbor(direction); @@ -1660,8 +1661,10 @@ VectorTools::project_boundary_values (const Mapping &mapping, const DoFHandler &dof, const typename FunctionMap::type &boundary_functions, const Quadrature &, - std::map &boundary_values) + std::map &boundary_values, + std::vector component_mapping) { + Assert (component_mapping.size() == 0, ExcNotImplemented()); // projection in 1d is equivalent // to interpolation interpolate_boundary_values (mapping, dof, boundary_functions, @@ -1677,7 +1680,8 @@ VectorTools::project_boundary_values (const Mapping &mapping, const DoFHandler &dof, const typename FunctionMap::type &boundary_functions, const Quadrature &q, - std::map &boundary_values) + std::map &boundary_values, + std::vector component_mapping) { //TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity // structures, matrices and right hand sides or distribution of @@ -1685,9 +1689,18 @@ VectorTools::project_boundary_values (const Mapping &mapping, // there are no constrained nodes on the boundary, but is not // acceptable for higher dimensions. Fix this. - Assert (dof.get_fe().n_components() == boundary_functions.begin()->second->n_components, - ExcDimensionMismatch(dof.get_fe().n_components(), - boundary_functions.begin()->second->n_components)); + if (component_mapping.size() == 0) + { + AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components); + // I still do not see why i + // should create another copy + // here + component_mapping.resize(dof.get_fe().n_components()); + for (unsigned int i= 0 ;i < component_mapping.size() ; ++i) + component_mapping[i] = i; + } + else + AssertDimension (dof.get_fe().n_components(), component_mapping.size()); std::vector dof_to_boundary_mapping; std::set selected_boundary_components; @@ -1755,7 +1768,8 @@ VectorTools::project_boundary_values (const Mapping &mapping, MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, mass_matrix, boundary_functions, - rhs, dof_to_boundary_mapping); + rhs, dof_to_boundary_mapping, (const Function*) 0, + component_mapping); // For certain weird elements, // there might be degrees of @@ -1764,7 +1778,7 @@ VectorTools::project_boundary_values (const Mapping &mapping, // have support there. Let's // eliminate them here. -//TODO: Maybe we should figure out ith the element really needs this +//TODO: Maybe we should figure out if the element really needs this FilteredMatrix > filtered_mass_matrix(mass_matrix, true); FilteredMatrix > filtered_precondition; @@ -1824,10 +1838,12 @@ void VectorTools::project_boundary_values (const DoFHandler &dof, const typename FunctionMap::type &boundary_functions, const Quadrature &q, - std::map &boundary_values) + std::map &boundary_values, + std::vector component_mapping) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - project_boundary_values(StaticMappingQ1::mapping, dof, boundary_functions, q, boundary_values); + project_boundary_values(StaticMappingQ1::mapping, dof, boundary_functions, q, + boundary_values, component_mapping); } diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 6e6c9264ae..185f813235 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -909,14 +909,16 @@ void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, #if deal_II_dimension == 1 -void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &, - const DoFHandler<1> &, - const Quadrature<0> &, - SparseMatrix &, - const FunctionMap<1>::type&, - Vector &, - std::vector &, - const Function<1> * const) +template +void MatrixCreator::create_boundary_mass_matrix (const Mapping&, + const DoFHandler&, + const Quadrature&, + SparseMatrix&, + const typename FunctionMap::type&, + Vector&, + std::vector&, + const Function* const, + std::vector) { // what would that be in 1d? the // identity matrix on the boundary @@ -925,7 +927,7 @@ void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &, } -#endif +#else @@ -938,8 +940,33 @@ MatrixCreator::create_boundary_mass_matrix (const Mapping &mapping, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const coefficient) + const Function * const coefficient, + std::vector component_mapping) { + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0;in_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + const unsigned int n_threads = multithread_info.n_default_threads; Threads::ThreadGroup<> threads; @@ -953,24 +980,28 @@ MatrixCreator::create_boundary_mass_matrix (const Mapping &mapping, // mutex to synchronise access to // the matrix Threads::ThreadMutex mutex; + + typedef boost::tuple&, + const DoFHandler&, + const Quadrature&> Commons; // then assemble in parallel typedef void (*create_boundary_mass_matrix_1_t) - (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, + (Commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, const Function * const coefficient, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex); create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1; for (unsigned int thread=0; thread &mapping, template void MatrixCreator:: -create_boundary_mass_matrix_1 (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, +create_boundary_mass_matrix_1 (boost::tuple&, + const DoFHandler&, + const Quadrature&> commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, const Function * const coefficient, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex) { + // All assertions for this function + // are in the calling function + // before creating threads. + const Mapping& mapping = boost::get<0>(commons); + const DoFHandler& dof = boost::get<1>(commons); + const Quadrature& q = boost::get<2>(commons); + const FiniteElement &fe = dof.get_fe(); const unsigned int n_components = fe.n_components(); + const unsigned int n_function_components = boundary_functions.begin()->second->n_components; const bool fe_is_system = (n_components != 1); - const bool fe_is_primitive = fe.is_primitive(); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (n_components == boundary_functions.begin()->second->n_components, - ExcComponentMismatch()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); + const bool fe_is_primitive = fe.is_primitive(); const unsigned int dofs_per_cell = fe.dofs_per_cell, dofs_per_face = fe.dofs_per_face; @@ -1018,6 +1045,7 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, UpdateFlags update_flags = UpdateFlags (update_values | update_JxW_values | + update_normal_vectors | update_quadrature_points); FEFaceValues fe_values (mapping, fe, q, update_flags); @@ -1032,7 +1060,7 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, std::vector rhs_values_scalar (fe_values.n_quadrature_points); std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_components)); + Vector(n_function_components)); std::vector dofs (dofs_per_cell); std::vector dofs_on_face_vector (dofs_per_face); @@ -1085,6 +1113,27 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, for (unsigned int point=0; point > normal_adjustment(fe_values.n_quadrature_points, + std::vector(n_components, 1.)); + + for (unsigned int comp = 0;comp& base = fe.base_element(fe.component_to_base_index(comp).first); + const unsigned int bcomp = fe.component_to_base_index(comp).second; + + if (!base.conforms(FiniteElementData::H1) && + base.conforms(FiniteElementData::Hdiv)) + for (unsigned int point=0; point &mapping, { for (unsigned int j=0; j &mapping, } } cell_vector(i) += fe_values.shape_value(i,point) - * rhs_values_system[point](fe.system_to_component_index(i).first) + * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first]) * weight; } else @@ -1116,9 +1165,11 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, cell_matrix(i,j) += fe_values.shape_value_component(j,point,comp) * fe_values.shape_value_component(i,point,comp) + * normal_adjustment[point][comp] * weight * coefficient_vector_values[point](comp); cell_vector(i) += fe_values.shape_value_component(i,point,comp) * - rhs_values_system[point](comp) + rhs_values_system[point](component_mapping[comp]) + * normal_adjustment[point][comp] * weight; } } @@ -1222,6 +1273,7 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, } } +#endif template void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, @@ -1230,25 +1282,29 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, const typename FunctionMap::type &rhs, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a) + const Function * const a, + std::vector component_mapping) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); create_boundary_mass_matrix(StaticMappingQ1::mapping, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a); + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); } #if deal_II_dimension == 1 -void MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection<1> &, - const hp::DoFHandler<1> &, - const hp::QCollection<0> &, - SparseMatrix &, - const FunctionMap<1>::type&, - Vector &, - std::vector &, - const Function<1> * const) +template +void +MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection&, + const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const typename FunctionMap::type&, + Vector&, + std::vector&, + const Function* const, + std::vector) { // what would that be in 1d? the // identity matrix on the boundary @@ -1257,8 +1313,7 @@ void MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection<1> } -#endif - +#else template @@ -1270,8 +1325,32 @@ MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const coefficient) + const Function * const coefficient, + std::vector component_mapping) { + const hp::FECollection &fe_collection = dof.get_fe(); + const unsigned int n_components = fe_collection.n_components(); + + Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), + ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (boundary_functions.size() != 0, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), + ExcInternalError()); + Assert (coefficient ==0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcComponentMismatch()); + + if (component_mapping.size() == 0) + { + AssertDimension (n_components, boundary_functions.begin()->second->n_components); + for (unsigned int i=0;i threads; @@ -1282,27 +1361,31 @@ MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection = Threads::split_range (dof.begin_active(), dof.end(), n_threads); + typedef boost::tuple&, + const hp::DoFHandler&, + const hp::QCollection&> Commons; + // mutex to synchronise access to // the matrix Threads::ThreadMutex mutex; // then assemble in parallel typedef void (*create_boundary_mass_matrix_1_t) - (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, + (Commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, const Function * const coefficient, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex); create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1; for (unsigned int thread=0; thread template void MatrixCreator:: -create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, +create_boundary_mass_matrix_1 (boost::tuple&, + const hp::DoFHandler&, + const hp::QCollection&> commons, SparseMatrix &matrix, const typename FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, const Function * const coefficient, + const std::vector& component_mapping, const IteratorRange > range, Threads::ThreadMutex &mutex) -{ +{ + const hp::MappingCollection& mapping = boost::get<0>(commons); + const hp::DoFHandler& dof = boost::get<1>(commons); + const hp::QCollection& q = boost::get<2>(commons); const hp::FECollection &fe_collection = dof.get_fe(); const unsigned int n_components = fe_collection.n_components(); + const unsigned int n_function_components = boundary_functions.begin()->second->n_components; const bool fe_is_system = (n_components != 1); - - Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), - ExcInternalError()); - Assert (matrix.n() == matrix.m(), ExcInternalError()); - Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); - Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), - ExcInternalError()); - Assert (n_components == boundary_functions.begin()->second->n_components, - ExcComponentMismatch()); - Assert (coefficient ==0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcComponentMismatch()); #ifdef DEBUG if (true) { @@ -1406,7 +1481,7 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, // FE has several components { rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(n_components)); + Vector(n_function_components)); boundary_functions.find(cell->face(face)->boundary_indicator()) ->second->vector_value_list (fe_values.get_quadrature_points(), rhs_values_system); @@ -1436,7 +1511,7 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, } cell_vector(i) += v * rhs_values_system[point]( - fe.system_to_component_index(i).first) * weight; + component_mapping[fe.system_to_component_index(i).first]) * weight; } } } @@ -1464,7 +1539,7 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, (u * v * weight * coefficient_vector_values[point](component_i)); } } - cell_vector(i) += v * rhs_values_system[point](component_i) * weight; + cell_vector(i) += v * rhs_values_system[point](component_mapping[component_i]) * weight; } } } @@ -1687,6 +1762,7 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, } } +#endif template void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler &dof, @@ -1695,11 +1771,12 @@ void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler & const typename FunctionMap::type &rhs, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a) + const Function * const a, + std::vector component_mapping) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); create_boundary_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix,rhs, rhs_vector, dof_to_boundary_mapping, a); + matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping); } @@ -2619,7 +2696,6 @@ void MatrixCreator::create_mass_matrix -#if deal_II_dimension != 1 template void MatrixCreator::create_boundary_mass_matrix @@ -2630,8 +2706,9 @@ MatrixCreator::create_boundary_mass_matrix const FunctionMap::type &boundary_functions, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a); -#endif + const Function * const a, + std::vector); + template void MatrixCreator::create_boundary_mass_matrix @@ -2641,32 +2718,34 @@ void MatrixCreator::create_boundary_mass_matrix const FunctionMap::type &rhs, Vector &rhs_vector, std::vector &dof_to_boundary_mapping, - const Function * const a); + const Function * const a, + std::vector); -#if deal_II_dimension != 1 template void MatrixCreator::create_boundary_mass_matrix -(const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a); -#endif +(const hp::MappingCollection&, + const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const FunctionMap::type&, + Vector&, + std::vector&, + const Function * const, + std::vector); + template void MatrixCreator::create_boundary_mass_matrix -(const hp::DoFHandler &dof, - const hp::QCollection &q, - SparseMatrix &matrix, - const FunctionMap::type &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function * const a); +(const hp::DoFHandler&, + const hp::QCollection&, + SparseMatrix&, + const FunctionMap::type&, + Vector&, + std::vector&, + const Function * const, + std::vector); diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 1fec1f5211..7c79573098 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -143,14 +143,14 @@ void VectorTools::project_boundary_values const DoFHandler &, const FunctionMap::type &, const Quadrature&, - std::map &); + std::map&, std::vector); template void VectorTools::project_boundary_values (const DoFHandler &, const FunctionMap::type &, const Quadrature&, - std::map &); + std::map&, std::vector); #if deal_II_dimension != 1