/**
* 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
* 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.
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<dim> * 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<double> &matrix,
- const FunctionMap<1>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<1> * const a = 0);
-
-
+ const Function<dim> * const weight = 0,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
+
/**
* Calls the
* create_boundary_mass_matrix()
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<dim> * const a = 0);
+ const Function<dim> * const a = 0,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
* Same function as above, but for hp
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<dim> * 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<double> &matrix,
- const FunctionMap<1>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<1> * const a = 0);
+ const Function<dim> * const a = 0,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
* Same function as above, but for hp
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
- const Function<dim> * const a = 0);
+ const Function<dim> * const a = 0,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
* Assemble the Laplace
*/
template <int dim>
static
- void create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
+ void create_boundary_mass_matrix_1 (boost::tuple<const Mapping<dim>&,
+ const DoFHandler<dim>&,
+ const Quadrature<dim-1>&> commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
const Function<dim> * const a,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
*/
template <int dim>
static
- void create_boundary_mass_matrix_1 (const hp::MappingCollection<dim> &mapping,
- const hp::DoFHandler<dim> &dof,
- const hp::QCollection<dim-1> &q,
+ void create_boundary_mass_matrix_1 (boost::tuple<const hp::MappingCollection<dim>&,
+ const hp::DoFHandler<dim>&,
+ const hp::QCollection<dim-1>&> commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
const Function<dim> * const a,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
};
* given function may be, taking into account that a virtual function has
* to be called.
*
- * <li> 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.
+ * <li> Projection: compute the <i>L</i><sup>2</sup>-projection of the
+ * given function onto the finite element space, i.e. if <i>f</i> is
+ * the function to be projected, compute <i>f<sub>h</sub></i> in
+ * <i>V<sub>h</sub></i> such that
+ * (<i>f<sub>h</sub></i>,<i>v<sub>h</sub></i>)=(<i>f</i>,<i>v<sub>h</sub></i>)$
+ * for all discrete test functions <i>v<sub>h</sub></i>. This is done
+ * through the solution of the linear system of equations <i> M v =
+ * f</i> where <i>M</i> 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 <i>f<sub>h</sub></i>. 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 <i>L</i><sup>2</sup>-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
*
* 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
+ * <i>L</i><sup>2</sup>-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
+ * <i>L</i><sup>2</sup>-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.
*
* detail may change in the future.
*
* <li> 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
* <tt>MatrixCreator::create_*</tt> functions which take a right hand side do,
* but without assembling a matrix.
*
* <li> 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$.
*
* <li> 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
* integration shall extend.
*
* <li> 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 <tt>unsigned char</tt>;
* see the section in the documentation of the Triangulation class for more
* 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
* the place of the respective boundary point.
*
* <li> 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 <i>L</i><sup>2</sup>-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.
*
* <li> 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
*
* 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.
*
* 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 <i>H</i><sup>1</sup> seminorm is the <i>L</i><sup>2</sup>
* 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.
+ * <i>H</i><sup>1</sup> norm is the sum of the square of seminorm
+ * and the square of the <i>L</i><sup>2</sup> norm.
*
* To get the global <i>L<sup>1</sup></i> error, you have to sum up the
* entries in @p difference, e.g. using
- * <tt>Vector<double>::l1_norm</tt> function. For the global $L^2$
+ * Vector::l1_norm() function. For the global <i>L</i><sup>2</sup>
* difference, you have to sum up the squares of the entries and
* take the root of the sum, e.g. using
- * <tt>Vector<double>::l2_norm</tt>. 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 <i>l</i><sub>1</sub> and <i>l</i><sub>2</sub> 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 <tt>Vector<double>::linfty_norm</tt> 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 <i>H</i><sup>1</sup> norm and seminorm, the same rule applies as for the
+ * <i>L</i><sup>2</sup> norm: compute the <i>l</i><sub>2</sub> norm of the cell error vector.
* </ul>
*
* All functions use the finite element given to the DoFHandler object the last
VECTOR &vec);
/**
- * Calls the @p interpolate
+ * Calls the @p interpolate()
* function above with
* <tt>mapping=MappingQ1@<dim>@()</tt>.
*/
/**
* Calls the other
- * @p interpolate_boundary_values
+ * interpolate_boundary_values()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
/**
* Calls the other
- * @p interpolate_boundary_values
+ * interpolate_boundary_values()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
/**
- * 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
* 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 <int dim>
static void project_boundary_values (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &boundary_functions,
const Quadrature<dim-1> &q,
- std::map<unsigned int,double> &boundary_values);
+ std::map<unsigned int,double> &boundary_values,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
- * Calls the @p project_boundary_values
+ * Calls the project_boundary_values()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
static void project_boundary_values (const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &boundary_function,
const Quadrature<dim-1> &q,
- std::map<unsigned int,double> &boundary_values);
+ std::map<unsigned int,double> &boundary_values,
+ std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
Vector<double> &rhs_vector);
/**
- * Calls the @p create_right_hand_side
+ * Calls the create_right_hand_side()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
Vector<double> &rhs_vector);
/**
- * Calls the @p create_point_source_vector
+ * Calls the create_point_source_vector()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
/**
* Calls the
- * @p create_boundary_right_hand_side
+ * create_boundary_right_hand_side()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
const std::set<unsigned char> &boundary_indicators = std::set<unsigned char>());
/**
- * 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 <int dim>
static void create_boundary_right_hand_side (const hp::DoFHandler<dim> &dof,
const double exponent = 2.);
/**
- * Calls the @p integrate_difference
+ * Calls the integrate_difference()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
const double exponent = 2.);
/**
- * Calls the @p integrate_difference
+ * Calls the integrate_difference()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
const unsigned int component);
/**
- * Calls the @p compute_mean_value
+ * Calls the compute_mean_value()
* function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
{
#if deal_II_dimension == 1
+ template <int dim>
void
- interpolate_zero_boundary_values (const dealii::DoFHandler<1> &dof_handler,
+ interpolate_zero_boundary_values (const dealii::DoFHandler<dim> &dof_handler,
std::map<unsigned int,double> &boundary_values)
{
// we only need to find the
// dof indices. that's easy :-)
for (unsigned direction=0; direction<2; ++direction)
{
- dealii::DoFHandler<1>::cell_iterator
+ typename dealii::DoFHandler<dim>::cell_iterator
cell = dof_handler.begin(0);
while (!cell->at_boundary(direction))
cell = cell->neighbor(direction);
const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &boundary_functions,
const Quadrature<dim-1> &,
- std::map<unsigned int,double> &boundary_values)
+ std::map<unsigned int,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
{
+ Assert (component_mapping.size() == 0, ExcNotImplemented());
// projection in 1d is equivalent
// to interpolation
interpolate_boundary_values (mapping, dof, boundary_functions,
const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &boundary_functions,
const Quadrature<dim-1> &q,
- std::map<unsigned int,double> &boundary_values)
+ std::map<unsigned int,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
{
//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity
// structures, matrices and right hand sides or distribution of
// 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<unsigned int> dof_to_boundary_mapping;
std::set<unsigned char> selected_boundary_components;
MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
mass_matrix, boundary_functions,
- rhs, dof_to_boundary_mapping);
+ rhs, dof_to_boundary_mapping, (const Function<dim>*) 0,
+ component_mapping);
// For certain weird elements,
// there might be degrees of
// 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<Vector<double> > filtered_mass_matrix(mass_matrix, true);
FilteredMatrix<Vector<double> > filtered_precondition;
VectorTools::project_boundary_values (const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &boundary_functions,
const Quadrature<dim-1> &q,
- std::map<unsigned int,double> &boundary_values)
+ std::map<unsigned int,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- project_boundary_values(StaticMappingQ1<dim>::mapping, dof, boundary_functions, q, boundary_values);
+ project_boundary_values(StaticMappingQ1<dim>::mapping, dof, boundary_functions, q,
+ boundary_values, component_mapping);
}
// $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
#if deal_II_dimension == 1
-void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const)
+template <int dim>
+void MatrixCreator::create_boundary_mass_matrix (const Mapping<dim>&,
+ const DoFHandler<dim>&,
+ const Quadrature<dim-1>&,
+ SparseMatrix<double>&,
+ const typename FunctionMap<dim>::type&,
+ Vector<double>&,
+ std::vector<unsigned int>&,
+ const Function<dim>* const,
+ std::vector<unsigned int>)
{
// what would that be in 1d? the
// identity matrix on the boundary
}
-#endif
+#else
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const coefficient)
+ const Function<dim> * const coefficient,
+ std::vector<unsigned int> component_mapping)
{
+ const FiniteElement<dim> &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;i<n_components;++i)
+ component_mapping.push_back(i);
+ }
+ else
+ AssertDimension (n_components, component_mapping.size());
+
+ Assert (coefficient ==0 ||
+ coefficient->n_components==1 ||
+ coefficient->n_components==n_components, ExcComponentMismatch());
+
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadGroup<> threads;
// mutex to synchronise access to
// the matrix
Threads::ThreadMutex mutex;
+
+ typedef boost::tuple<const Mapping<dim>&,
+ const DoFHandler<dim>&,
+ const Quadrature<dim-1>&> Commons;
// then assemble in parallel
typedef void (*create_boundary_mass_matrix_1_t)
- (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
+ (Commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<dim> * const coefficient,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1<dim>;
for (unsigned int thread=0; thread<n_threads; ++thread)
- threads += Threads::spawn (p)(mapping, dof, q, matrix,
+ threads += Threads::spawn (p)(Commons(mapping, dof, q), matrix,
boundary_functions, rhs_vector,
dof_to_boundary_mapping, coefficient,
+ component_mapping,
thread_ranges[thread], mutex);
threads.join_all ();
}
template <int dim>
void
MatrixCreator::
-create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
+create_boundary_mass_matrix_1 (boost::tuple<const Mapping<dim>&,
+ const DoFHandler<dim>&,
+ const Quadrature<dim-1>&> commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<dim> * const coefficient,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<DoFHandler<dim> > range,
Threads::ThreadMutex &mutex)
{
+ // All assertions for this function
+ // are in the calling function
+ // before creating threads.
+ const Mapping<dim>& mapping = boost::get<0>(commons);
+ const DoFHandler<dim>& dof = boost::get<1>(commons);
+ const Quadrature<dim-1>& q = boost::get<2>(commons);
+
const FiniteElement<dim> &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;
UpdateFlags update_flags = UpdateFlags (update_values |
update_JxW_values |
+ update_normal_vectors |
update_quadrature_points);
FEFaceValues<dim> fe_values (mapping, fe, q, update_flags);
std::vector<double> rhs_values_scalar (fe_values.n_quadrature_points);
std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
- Vector<double>(n_components));
+ Vector<double>(n_function_components));
std::vector<unsigned int> dofs (dofs_per_cell);
std::vector<unsigned int> dofs_on_face_vector (dofs_per_face);
for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
coefficient_vector_values[point] = coefficient_values[point];
}
+
+ // Special treatment
+ // for Hdiv and Hcurl
+ // elements, where only
+ // the normal or
+ // tangential component
+ // should be projected.
+ std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points,
+ std::vector<double>(n_components, 1.));
+
+ for (unsigned int comp = 0;comp<n_components;++comp)
+ {
+ const FiniteElement<dim>& 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<dim>::H1) &&
+ base.conforms(FiniteElementData<dim>::Hdiv))
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp)
+ * fe_values.normal_vector(point)(bcomp);
+ }
for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
{
{
for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
{
- if (fe.system_to_component_index(j).first ==
- fe.system_to_component_index(i).first)
+ if (fe.system_to_component_index(j).first
+ == fe.system_to_component_index(i).first)
{
cell_matrix(i,j)
+= weight
}
}
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
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;
}
}
}
}
+#endif
template <int dim>
void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &rhs,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const a)
+ const Function<dim> * const a,
+ std::vector<unsigned int> component_mapping)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
create_boundary_mass_matrix(StaticMappingQ1<dim>::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<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const)
+template <int dim>
+void
+MatrixCreator::create_boundary_mass_matrix (const hp::MappingCollection<dim>&,
+ const hp::DoFHandler<dim>&,
+ const hp::QCollection<dim-1>&,
+ SparseMatrix<double>&,
+ const typename FunctionMap<dim>::type&,
+ Vector<double>&,
+ std::vector<unsigned int>&,
+ const Function<dim>* const,
+ std::vector<unsigned int>)
{
// what would that be in 1d? the
// identity matrix on the boundary
}
-#endif
-
+#else
template <int dim>
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const coefficient)
+ const Function<dim> * const coefficient,
+ std::vector<unsigned int> component_mapping)
{
+ const hp::FECollection<dim> &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<n_components;++i)
+ component_mapping.push_back(i);
+ }
+ else
+ AssertDimension (n_components, component_mapping.size());
+
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadGroup<> threads;
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
+ typedef boost::tuple<const hp::MappingCollection<dim>&,
+ const hp::DoFHandler<dim>&,
+ const hp::QCollection<dim-1>&> 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<dim> &mapping,
- const hp::DoFHandler<dim> &dof,
- const hp::QCollection<dim-1> &q,
+ (Commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<dim> * const coefficient,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex);
create_boundary_mass_matrix_1_t p = &MatrixCreator::template create_boundary_mass_matrix_1<dim>;
for (unsigned int thread=0; thread<n_threads; ++thread)
- threads += Threads::spawn (p)(mapping, dof, q, matrix,
+ threads += Threads::spawn (p)(Commons(mapping, dof, q), matrix,
boundary_functions, rhs_vector,
dof_to_boundary_mapping, coefficient,
+ component_mapping,
thread_ranges[thread], mutex);
threads.join_all ();
}
template <int dim>
void
MatrixCreator::
-create_boundary_mass_matrix_1 (const hp::MappingCollection<dim> &mapping,
- const hp::DoFHandler<dim> &dof,
- const hp::QCollection<dim-1> &q,
+create_boundary_mass_matrix_1 (boost::tuple<const hp::MappingCollection<dim>&,
+ const hp::DoFHandler<dim>&,
+ const hp::QCollection<dim-1>&> commons,
SparseMatrix<double> &matrix,
const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<dim> * const coefficient,
+ const std::vector<unsigned int>& component_mapping,
const IteratorRange<hp::DoFHandler<dim> > range,
Threads::ThreadMutex &mutex)
-{
+{
+ const hp::MappingCollection<dim>& mapping = boost::get<0>(commons);
+ const hp::DoFHandler<dim>& dof = boost::get<1>(commons);
+ const hp::QCollection<dim-1>& q = boost::get<2>(commons);
const hp::FECollection<dim> &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)
{
// FE has several components
{
rhs_values_system.resize (fe_values.n_quadrature_points,
- Vector<double>(n_components));
+ Vector<double>(n_function_components));
boundary_functions.find(cell->face(face)->boundary_indicator())
->second->vector_value_list (fe_values.get_quadrature_points(),
rhs_values_system);
}
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;
}
}
}
(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;
}
}
}
}
}
+#endif
template <int dim>
void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler<dim> &dof,
const typename FunctionMap<dim>::type &rhs,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const a)
+ const Function<dim> * const a,
+ std::vector<unsigned int> component_mapping)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
create_boundary_mass_matrix(hp::StaticMappingQ1<dim>::mapping_collection, 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
template
void
MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
const FunctionMap<deal_II_dimension>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<deal_II_dimension> * const a);
-#endif
+ const Function<deal_II_dimension> * const a,
+ std::vector<unsigned int>);
+
template
void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
const FunctionMap<deal_II_dimension>::type &rhs,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<deal_II_dimension> * const a);
+ const Function<deal_II_dimension> * const a,
+ std::vector<unsigned int>);
-#if deal_II_dimension != 1
template
void
MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
-(const hp::MappingCollection<deal_II_dimension> &mapping,
- const hp::DoFHandler<deal_II_dimension> &dof,
- const hp::QCollection<deal_II_dimension-1> &q,
- SparseMatrix<double> &matrix,
- const FunctionMap<deal_II_dimension>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<deal_II_dimension> * const a);
-#endif
+(const hp::MappingCollection<deal_II_dimension>&,
+ const hp::DoFHandler<deal_II_dimension>&,
+ const hp::QCollection<deal_II_dimension-1>&,
+ SparseMatrix<double>&,
+ const FunctionMap<deal_II_dimension>::type&,
+ Vector<double>&,
+ std::vector<unsigned int>&,
+ const Function<deal_II_dimension> * const,
+ std::vector<unsigned int>);
+
template
void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension>
-(const hp::DoFHandler<deal_II_dimension> &dof,
- const hp::QCollection<deal_II_dimension-1> &q,
- SparseMatrix<double> &matrix,
- const FunctionMap<deal_II_dimension>::type &rhs,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<deal_II_dimension> * const a);
+(const hp::DoFHandler<deal_II_dimension>&,
+ const hp::QCollection<deal_II_dimension-1>&,
+ SparseMatrix<double>&,
+ const FunctionMap<deal_II_dimension>::type&,
+ Vector<double>&,
+ std::vector<unsigned int>&,
+ const Function<deal_II_dimension> * const,
+ std::vector<unsigned int>);
const DoFHandler<deal_II_dimension> &,
const FunctionMap<deal_II_dimension>::type &,
const Quadrature<deal_II_dimension-1>&,
- std::map<unsigned int,double> &);
+ std::map<unsigned int,double>&, std::vector<unsigned int>);
template
void VectorTools::project_boundary_values<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &,
const FunctionMap<deal_II_dimension>::type &,
const Quadrature<deal_II_dimension-1>&,
- std::map<unsigned int,double> &);
+ std::map<unsigned int,double>&, std::vector<unsigned int>);
#if deal_II_dimension != 1