From: Wolfgang Bangerth Date: Sat, 23 Feb 2008 19:14:34 +0000 (+0000) Subject: Implement an entirely new scheme to deal with vector-valued problems. Also provide... X-Git-Tag: v8.0.0~9401 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=777a1fb2e21d6fea0b3016d07c5e25fc129bd242;p=dealii.git Implement an entirely new scheme to deal with vector-valued problems. Also provide an insane amount of documentation for it. git-svn-id: https://svn.dealii.org/trunk@15757 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 4c8006bdf3..6d93d2dab7 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -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 @@ -2684,7 +2684,12 @@ FiniteElement::is_primitive (const unsigned int i) const // std::vector is) is // probably more expensive than // just comparing against 1 - return (n_nonzero_components_table[i] == 1); + // + // for good measure, short circuit the test + // if the entire FE is primitive + return ((cached_primitivity == true) + || + (n_nonzero_components_table[i] == 1)); } diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index a8abf92fd7..78549296bb 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 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 @@ -25,11 +25,14 @@ DEAL_II_NAMESPACE_OPEN /** - * This class provides an interface to group several elements together - * into one. To the outside world, the resulting object looks just - * like a usual finite element object, which is composed of several - * other finite elements that are possibly of different type. The - * result is then a vector-valued finite element. + * This class provides an interface to group several elements together into + * one. To the outside world, the resulting object looks just like a usual + * finite element object, which is composed of several other finite elements + * that are possibly of different type. The result is then a vector-valued + * finite element. Vector valued elements are discussed in a number of + * tutorial programs, for example @ref step_8 "step-8", @ref step_20 + * "step-20", @ref step_21 "step-21", and in particular in the @ref vector_valued + * module. * * The overall numbering of degrees of freedom is as follows: for each * subobject (vertex, line, quad, or hex), the degrees of freedom are @@ -70,7 +73,7 @@ DEAL_II_NAMESPACE_OPEN * multiplicity. The number of blocks of the system is simply the sum * of all multiplicities. * - * @ingroup febase fe + * @ingroup febase fe vector_valued * * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002, 2003, 2006, Ralf Hartmann 2001. */ diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index f2e9e34210..2f81caf1f0 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -43,13 +44,380 @@ DEAL_II_NAMESPACE_OPEN template class Quadrature; +template class FEValuesBase; + + + +/** + * A namespace in which we declare "extractors", i.e. classes that when used + * as subscripts in operator[] expressions on FEValues, FEFaceValues, and + * FESubfaceValues objects extract certain components of a vector-valued + * element. There are extractors for single scalar components as well as + * vector components consisting of dim elements. + * + * See the description of the @ref vector_valued module for examples how to + * use the features of this namespace. + * + * @ingroup feaccess vector_valued + */ +namespace FEValuesExtractors +{ + /** + * Extractor for a single scalar component + * of a vector-valued element. The concept + * of extractors is defined in the + * documentation of the namespace + * FEValuesExtractors and in the @ref + * vector_valued module. + * + * @ingroup feaccess vector_valued + */ + struct Scalar + { + /** + * The selected scalar component of the + * vector. + */ + const unsigned int component; + + /** + * Constructor. Take the selected + * vector component as argument. + */ + Scalar (const unsigned int component); + }; + + + /** + * Extractor for a vector of + * dim components of a + * vector-valued element. The value of + * dim is defined by the + * FEValues object the extractor is applied + * to. + * + * The concept of + * extractors is defined in the + * documentation of the namespace + * FEValuesExtractors and in the @ref + * vector_valued module. + * + * @ingroup feaccess vector_valued + */ + struct Vector + { + /** + * The first component of the vector view. + */ + const unsigned int first_vector_component; + + /** + * Constructor. Take the first + * component of the selected vector + * inside the FEValues object as + * argument. + */ + Vector (const unsigned int first_vector_component); + }; +} + + + +/** + * A namespace for "views" on a FEValues, FEFaceValues, or FESubfaceValues + * object. A view represents only a certain part of the whole: whereas the + * FEValues object represents all values, gradients, or second + * derivatives of all components of a vector-valued element, views restrict + * the attention to only a single component or a subset of components. + * + * There are classes that present views for single scalar components as well + * as vector components consisting of dim elements. + * + * See the description of the @ref vector_valued module for examples how to + * use the features of this namespace. + * + * @ingroup feaccess vector_valued + */ +namespace FEValuesViews +{ + /** + * A class representing a view to a single + * scalar component of a possibly + * vector-valued finite element. Views are + * discussed in the @ref vector_valued module. + * + * @ingroup feaccess vector_valued + */ + template + class Scalar + { + public: + /** + * A typedef for the data type of + * values of the view this class + * represents. Since we deal with a + * single components, the value type is + * a scalar double. + */ + typedef double value_type; + + /** + * A typedef for the type of gradients + * of the view this class + * represents. Here, for a scalar + * component of the finite element, the + * gradient is a + * Tensor@<1,dim@>. + */ + typedef Tensor<1,dim> gradient_type; + + /** + * A typedef for the type of second + * derivatives of the view this class + * represents. Here, for a scalar + * component of the finite element, the + * Hessian is a + * Tensor@<2,dim@>. + */ + typedef Tensor<2,dim> hessian_type; + + /** + * Constructor for an object that + * represents a single scalar component + * of a FEValuesBase object (or of one + * of the classes derived from + * FEValuesBase). + */ + Scalar (const FEValuesBase &fe_values_base, + const unsigned int component); + + /** + * Return the value of the vector + * component selected by this view, for + * the shape function and quadrature + * point selected by the arguments. + */ + value_type + value (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the gradient (a tensor of + * rank 1) of the vector component + * selected by this view, for the shape + * function and quadrature point + * selected by the arguments. + */ + gradient_type + gradient (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the Hessian (the tensor of + * rank 2 of all second derivatives) of + * the vector component selected by + * this view, for the shape function + * and quadrature point selected by the + * arguments. + */ + hessian_type + hessian (const unsigned int shape_function, + const unsigned int q_point) const; + + + private: + /** + * A pointer to the FEValuesBase object + * we operator on. + */ + const SmartPointer > fe_values; + + /** + * The single scalar component this + * view represents of the FEValuesBase + * object. + */ + const unsigned int component; + }; + -//TODO: Add access to mapping values to FEValuesBase -//TODO: Several FEValuesBase of a system should share Mapping + /** + * A class representing a view to a set of + * dim components forming a + * vector part of a vector-valued finite + * element. Views are discussed in the + * @ref vector_valued module. + * + * @ingroup feaccess vector_valued + */ + template + class Vector + { + public: + /** + * A typedef for the data type of + * values of the view this class + * represents. Since we deal with a set + * of dim components, the + * value type is a Tensor<1,dim>. + */ + typedef Tensor<1,dim> value_type; + + /** + * A typedef for the type of gradients + * of the view this class + * represents. Here, for a set of + * dim components of the + * finite element, the gradient is a + * Tensor@<2,dim@>. + */ + typedef Tensor<2,dim> gradient_type; + + /** + * A typedef for the type of + * symmetrized gradients of the view + * this class represents. Here, for a + * set of dim components + * of the finite element, the + * symmetrized gradient is a + * SymmetricTensor@<2,dim@>. + */ + typedef SymmetricTensor<2,dim> symmetric_gradient_type; + + /** + * A typedef for the type of the + * divergence of the view this class + * represents. Here, for a set of + * dim components of the + * finite element, the divergence of + * course is a scalar. + */ + typedef double divergence_type; + + /** + * A typedef for the type of second + * derivatives of the view this class + * represents. Here, for a set of + * dim components of the + * finite element, the Hessian is a + * Tensor@<3,dim@>. + */ + typedef Tensor<3,dim> hessian_type; + + /** + * Constructor for an object that + * represents dim components of a + * FEValuesBase object (or of one of + * the classes derived from + * FEValuesBase), representing a + * vector-valued variable. + * + * The second argument denotes the + * index of the first component of the + * selected vector. + */ + Vector (const FEValuesBase &fe_values_base, + const unsigned int first_vector_component); + + /** + * Return the value of the vector + * components selected by this view, + * for the shape function and + * quadrature point selected by the + * arguments. Here, since the view + * represents a vector-valued part of + * the FEValues object with + * dim components, the + * return type is a tensor of rank 1 + * with dim components. + */ + value_type + value (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the gradient (a tensor of + * rank 2) of the vector component + * selected by this view, for the shape + * function and quadrature point + * selected by the arguments. + */ + gradient_type + gradient (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the symmetric gradient (a + * symmetric tensor of rank 2) of the + * vector component selected by this + * view, for the shape function and + * quadrature point selected by the + * arguments. + * + * The symmetric gradient is defined as + * $\frac 12 [(\nabla \phi_i(x_q)) + + * (\nabla \phi_i(x_q))^T]$, where + * $\phi_i$ represents the + * dim components selected + * from the FEValuesBase object, and + * $x_q$ is the location of the $q$th + * quadrature point. + */ + symmetric_gradient_type + symmetric_gradient (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the scalar divergence of + * the vector components selected by + * this view, for the shape function + * and quadrature point selected by the + * arguments. + */ + divergence_type + divergence (const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the Hessian (the tensor of + * rank 2 of all second derivatives) of + * the vector components selected by + * this view, for the shape function + * and quadrature point selected by the + * arguments. + */ + hessian_type + hessian (const unsigned int shape_function, + const unsigned int q_point) const; + + private: + /** + * A pointer to the FEValuesBase object + * we operator on. + */ + const SmartPointer > fe_values; + + /** + * The first component of the vector + * this view represents of the + * FEValuesBase object. + */ + const unsigned int first_vector_component; + }; +} + + + + /*!@addtogroup feaccess */ /*@{*/ + + +//TODO: Add access to mapping values to FEValuesBase +//TODO: Several FEValuesBase of a system should share Mapping + /** * Contains all data vectors for FEValues. * This class has been extracted from FEValuesBase to be handed @@ -455,6 +823,39 @@ class FEValuesBase : protected FEValuesData, * Destructor. */ ~FEValuesBase (); + + /// @name Extractors Methods to extract individual components + //@{ + + /** + * Create a view of the current FEValues + * object that represents a particular + * scalar component of the possibly + * vector-valued finite element. The + * concept of views is explained in the + * documentation of the namespace + * FEValuesViews and in particular + * in the @ref vector_valued module. + */ + FEValuesViews::Scalar + operator[] (const FEValuesExtractors::Scalar &scalar) const; + + /** + * Create a view of the current FEValues + * object that represents a set of + * dim scalar components + * (i.e. a vector) of the vector-valued + * finite element. The concept of views + * is explained in the documentation of + * the namespace FEValuesViews and in particular + * in the @ref vector_valued module. + */ + FEValuesViews::Vector + operator[] (const FEValuesExtractors::Vector &vector) const; + + //@} + + /// @name ShapeAccess Access to shape function values //@{ @@ -1793,7 +2194,15 @@ class FEValuesBase : protected FEValuesData, * copyable, we make it private, * and also do not implement it. */ - FEValuesBase & operator= (const FEValuesBase &); + FEValuesBase & operator= (const FEValuesBase &); + + /** + * Make the view classes friends of this + * class, since they access internal + * data. + */ + template friend class FEValuesViews::Scalar; + template friend class FEValuesViews::Vector; }; @@ -2487,15 +2896,724 @@ class FESubfaceValues : public FEFaceValuesBase /*@}*/ #ifndef DOXYGEN + + +/*------------------------ Inline functions: namespace FEValuesExtractors --------*/ + +namespace FEValuesExtractors +{ + inline + Scalar::Scalar (const unsigned int component) + : + component (component) + {} + + + + inline + Vector::Vector (const unsigned int first_vector_component) + : + first_vector_component (first_vector_component) + {} +} + + +/*------------------------ Inline functions: namespace FEValuesViews --------*/ + +namespace FEValuesViews +{ + template + inline + Scalar::Scalar (const FEValuesBase &fe_values, + const unsigned int component) + : + fe_values (&fe_values), + component (component) + { + Assert (component < this->fe_values->fe->n_components(), + ExcIndexRange(component, 0, this->fe_values->fe->n_components())); + } + + + + template + inline + typename Scalar::value_type + Scalar::value (const unsigned int shape_function, + const unsigned int q_point) const + { + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_values, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // an adaptation of the + // FEValuesBase::shape_value_component + // function except that here we know the + // component as fixed. see the + // comments there + // + // we can do away with some of the + // assertions since they are already + // taken care of in the constructor of + // this class + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + if (component == + fe_values->fe->system_to_component_index(shape_function).first) + return fe_values->shape_values(fe_values->shape_function_to_row_table[shape_function],q_point); + else + return 0; + } + else + { + if (fe_values->fe->get_nonzero_components(shape_function)[component] == false) + return 0; + + const unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin()+ + component, + true)); + return fe_values->shape_values(row, q_point); + } + } + + + + template + inline + typename Scalar::gradient_type + Scalar::gradient (const unsigned int shape_function, + const unsigned int q_point) const + { + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // an adaptation of the + // FEValuesBase::shape_grad_component + // function except that here we know the + // component as fixed. see the + // comments there + // + // we can do away with the assertions + // since they are already taken care of + // in the constructor of this class + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + if (component == + fe_values->fe->system_to_component_index(shape_function).first) + return fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + else + return gradient_type(); + } + else + { + if (fe_values->fe->get_nonzero_components(shape_function)[component] == false) + return gradient_type(); + + const unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin()+ + component, + true)); + return fe_values->shape_gradients[row][q_point]; + } + } + + + + template + inline + typename Scalar::hessian_type + Scalar::hessian (const unsigned int shape_function, + const unsigned int q_point) const + { + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_hessians, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // an adaptation of the + // FEValuesBase::shape_grad_component + // function except that here we know the + // component as fixed. see the + // comments there + // + // we can do away with the assertions + // since they are already taken care of + // in the constructor of this class + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + if (component == + fe_values->fe->system_to_component_index(shape_function).first) + return fe_values->shape_hessians[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + else + return hessian_type(); + } + else + { + if (fe_values->fe->get_nonzero_components(shape_function)[component] == false) + return hessian_type(); + + const unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin()+ + component, + true)); + return fe_values->shape_hessians[row][q_point]; + } + } + + + + template + inline + Vector::Vector (const FEValuesBase &fe_values, + const unsigned int first_vector_component) + : + fe_values (&fe_values), + first_vector_component (first_vector_component) + { + Assert (first_vector_component+dim-1 < this->fe_values->fe->n_components(), + ExcIndexRange(first_vector_component+dim-1, 0, + this->fe_values->fe->n_components())); + } + + + + template + inline + typename Vector::value_type + Vector::value (const unsigned int shape_function, + const unsigned int q_point) const + { + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_values, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // compared to the scalar case above, we + // can save some work here because we + // know that we are querying a contiguous + // range of components + value_type return_value; + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + // if this is a primitive shape + // function then at most one element + // of the output vector is + // nonzero. find out if indeed one is + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + dim)) + return_value[nonzero_component - first_vector_component] + = fe_values->shape_values(fe_values-> + shape_function_to_row_table[shape_function], + q_point); + } + else + { + unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin() + + first_vector_component, + true)); + for (unsigned int d=0; dfe->get_nonzero_components(shape_function)[first_vector_component+d] == + true) + { + return_value[d] = fe_values->shape_values(row, q_point); + + if ((d != dim-1) + && + (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] + == true)) + ++row; + } + } + + return return_value; + } + + + + template + inline + typename Vector::gradient_type + Vector::gradient (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + typename FEValuesBase::ExcAccessToUninitializedField()); + + gradient_type return_value; + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + dim)) + return_value[nonzero_component - first_vector_component] + = fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + } + else + { + unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin() + + first_vector_component, + true)); + for (unsigned int d=0; dfe->get_nonzero_components(shape_function)[first_vector_component+d] == + true) + { + return_value[d] = fe_values->shape_gradients[row][q_point]; + + if ((d != dim-1) + && + (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] + == true)) + ++row; + } + } + + return return_value; + } + + + + template + inline + typename Vector::divergence_type + Vector::divergence (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + typename FEValuesBase::ExcAccessToUninitializedField()); + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + dim)) + return + fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point] + [nonzero_component - first_vector_component]; + else + return 0; + } + else + { + unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin() + + first_vector_component, + true)); + + double div = 0; + for (unsigned int d=0; dfe->get_nonzero_components(shape_function)[first_vector_component+d] == + true) + { + div += fe_values->shape_gradients[row][q_point][d]; + + if ((d != dim-1) + && + (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] + == true)) + ++row; + } + return div; + } + } + + + + template + inline + typename Vector::hessian_type + Vector::hessian (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_hessians, + typename FEValuesBase::ExcAccessToUninitializedField()); + + hessian_type return_value; + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + dim)) + return_value[nonzero_component - first_vector_component] + = fe_values->shape_hessians[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + } + else + { + unsigned int + row = (fe_values->shape_function_to_row_table[shape_function] + + + std::count (fe_values->fe->get_nonzero_components(shape_function).begin(), + fe_values->fe->get_nonzero_components(shape_function).begin() + + first_vector_component, + true)); + for (unsigned int d=0; dfe->get_nonzero_components(shape_function)[first_vector_component+d] == + true) + { + return_value[d] = fe_values->shape_hessians[row][q_point]; + + if ((d != dim-1) + && + (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] + == true)) + ++row; + } + } + + return return_value; + } + + + + // we duplicate the following function for + // each dimension since we need to + // initialize a few arrays of dimension + // dependent size. even though this happens + // in a switch(dim) clause, we get warnings + // from the compiler if we don't separate + // things into different functions + template <> + inline + Vector<1>::symmetric_gradient_type + Vector<1>::symmetric_gradient (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + FEValuesBase<1>::ExcAccessToUninitializedField()); + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if (nonzero_component == first_vector_component) + { + // first get the one component of + // the nonsymmetrized gradient + // that is not zero + const Tensor<1,1> grad + = fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + + // then form a symmetric tensor + // out of it. note that access to + // individual elements of a + // SymmetricTensor object is + // fairly slow. if we implemented + // the following in the naive + // way, we would therefore incur + // a very significant penalty: a + // preliminary version of the + // Stokes tutorial program would + // slow down from 17 to 23 + // seconds because of this single + // function! + // + // consequently, we try to be a + // bit smarter by already laying + // out the data in the right + // format and creating a + // symmetric tensor of it + const double array[symmetric_gradient_type::n_independent_components] + = { grad[0] }; + return symmetric_gradient_type(array); + } + else + return symmetric_gradient_type(); + } + else + { + return symmetrize (gradient(shape_function, q_point)); + } + } + + + + template <> + inline + Vector<2>::symmetric_gradient_type + Vector<2>::symmetric_gradient (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + FEValuesBase<2>::ExcAccessToUninitializedField()); + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + 2)) + { + // first get the one component of + // the nonsymmetrized gradient + // that is not zero + const Tensor<1,2> grad + = fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + + // then form a symmetric tensor + // out of it. note that access to + // individual elements of a + // SymmetricTensor object is + // fairly slow. if we implemented + // the following in the naive + // way, we would therefore incur + // a very significant penalty: a + // preliminary version of the + // Stokes tutorial program would + // slow down from 17 to 23 + // seconds because of this single + // function! + // + // consequently, we try to be a + // bit smarter by already laying + // out the data in the right + // format and creating a + // symmetric tensor of it + switch (nonzero_component - first_vector_component) + { + case 0: + { + const double array[symmetric_gradient_type::n_independent_components] + = { grad[0], 0, grad[1]/2 }; + return symmetric_gradient_type(array); + } + + case 1: + { + const double array[symmetric_gradient_type::n_independent_components] + = { 0, grad[1], grad[0]/2 }; + return symmetric_gradient_type(array); + } + + default: + Assert (false, ExcInternalError()); + return symmetric_gradient_type(); + } + } + else + return symmetric_gradient_type(); + } + else + { + return symmetrize (gradient(shape_function, q_point)); + } + } + + + + template <> + inline + Vector<3>::symmetric_gradient_type + Vector<3>::symmetric_gradient (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + Assert (shape_function < fe_values->fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell)); + Assert (fe_values->update_flags & update_gradients, + FEValuesBase<3>::ExcAccessToUninitializedField()); + + if (fe_values->fe->is_primitive() || + fe_values->fe->is_primitive(shape_function)) + { + const unsigned int + nonzero_component + = fe_values->fe->system_to_component_index(shape_function).first; + + if ((nonzero_component >= first_vector_component) + && + (nonzero_component < first_vector_component + 3)) + { + // first get the one component of + // the nonsymmetrized gradient + // that is not zero + const Tensor<1,3> grad + = fe_values->shape_gradients[fe_values-> + shape_function_to_row_table[shape_function]][q_point]; + + // then form a symmetric tensor + // out of it. note that access to + // individual elements of a + // SymmetricTensor object is + // fairly slow. if we implemented + // the following in the naive + // way, we would therefore incur + // a very significant penalty: a + // preliminary version of the + // Stokes tutorial program would + // slow down from 17 to 23 + // seconds because of this single + // function! + // + // consequently, we try to be a + // bit smarter by already laying + // out the data in the right + // format and creating a + // symmetric tensor of it + switch (nonzero_component - first_vector_component) + { + case 0: + { + const double array[symmetric_gradient_type::n_independent_components] + = { grad[0], 0, 0, grad[1]/2 , grad[2]/2, 0}; + return symmetric_gradient_type(array); + } + + case 1: + { + const double array[symmetric_gradient_type::n_independent_components] + = { 0, grad[1], 0, grad[0]/2, 0, grad[2]/2 }; + return symmetric_gradient_type(array); + } + + case 2: + { + const double array[symmetric_gradient_type::n_independent_components] + = { 0, 0, grad[2], 0, grad[0]/2, grad[1]/2 }; + return symmetric_gradient_type(array); + } + + default: + Assert (false, ExcInternalError()); + return symmetric_gradient_type(); + } + } + else + return symmetric_gradient_type(); + } + else + { + return symmetrize (gradient(shape_function, q_point)); + } + } +} + + + + + + /*------------------------ Inline functions: FEValuesBase ------------------------*/ + + + +template +inline +FEValuesViews::Scalar +FEValuesBase:: +operator[] (const FEValuesExtractors::Scalar &scalar) const +{ + return FEValuesViews::Scalar (*this, scalar.component); +} + + + +template +inline +FEValuesViews::Vector +FEValuesBase:: +operator[] (const FEValuesExtractors::Vector &vector) const +{ + return + FEValuesViews::Vector (*this, vector.first_vector_component); +} + + + template inline double FEValuesBase::shape_value (const unsigned int i, const unsigned int j) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); Assert (fe->is_primitive (i), @@ -2527,6 +3645,8 @@ FEValuesBase::shape_value_component (const unsigned int i, const unsigned int j, const unsigned int component) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2560,7 +3680,7 @@ FEValuesBase::shape_value_component (const unsigned int i, // is non-zero at all within // this component: if (fe->get_nonzero_components(i)[component] == false) - return 0.; + return 0; // count how many non-zero // component the shape function @@ -2577,7 +3697,7 @@ FEValuesBase::shape_value_component (const unsigned int i, fe->get_nonzero_components(i).begin()+component, true)); return this->shape_values(row, j); - }; + } } @@ -2588,6 +3708,8 @@ const Tensor<1,dim> & FEValuesBase::shape_grad (const unsigned int i, const unsigned int j) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); Assert (fe->is_primitive (i), @@ -2623,6 +3745,8 @@ FEValuesBase::shape_grad_component (const unsigned int i, const unsigned int j, const unsigned int component) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2673,7 +3797,7 @@ FEValuesBase::shape_grad_component (const unsigned int i, fe->get_nonzero_components(i).begin()+component, true)); return this->shape_gradients[row][j]; - }; + } } @@ -2682,8 +3806,10 @@ template inline const Tensor<2,dim> & FEValuesBase::shape_hessian (const unsigned int i, - const unsigned int j) const + const unsigned int j) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField()); Assert (fe->is_primitive (i), @@ -2722,13 +3848,16 @@ FEValuesBase::shape_2nd_derivative (const unsigned int i, } + template inline Tensor<2,dim> FEValuesBase::shape_hessian_component (const unsigned int i, - const unsigned int j, - const unsigned int component) const + const unsigned int j, + const unsigned int component) const { + Assert (i < fe->dofs_per_cell, + ExcIndexRange (i, 0, fe->dofs_per_cell)); Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2779,10 +3908,11 @@ FEValuesBase::shape_hessian_component (const unsigned int i, fe->get_nonzero_components(i).begin()+component, true)); return this->shape_hessians[row][j]; - }; + } } + template inline Tensor<2,dim> @@ -2794,6 +3924,7 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i, } + template inline const FiniteElement & @@ -2870,19 +4001,24 @@ FEValuesBase::JxW (const unsigned int i) const } + template template +inline void FEValuesBase::get_function_grads (const InputVector &fe_function, - std::vector > &gradients) const + std::vector > &gradients) const { get_function_gradients(fe_function, gradients); } + template template -void FEValuesBase::get_function_grads ( +inline +void +FEValuesBase::get_function_grads ( const InputVector& fe_function, const VectorSlice >& indices, std::vector > &values) const @@ -2891,8 +4027,10 @@ void FEValuesBase::get_function_grads ( } + template template +inline void FEValuesBase:: get_function_grads (const InputVector &fe_function, @@ -2902,9 +4040,12 @@ get_function_grads (const InputVector &fe_function, } + template template -void FEValuesBase::get_function_grads ( +inline +void +FEValuesBase::get_function_grads ( const InputVector& fe_function, const VectorSlice >& indices, std::vector > >& values, @@ -2917,7 +4058,8 @@ void FEValuesBase::get_function_grads ( template template -inline void +inline +void FEValuesBase:: get_function_2nd_derivatives (const InputVector &fe_function, std::vector > &hessians) const @@ -2929,6 +4071,7 @@ get_function_2nd_derivatives (const InputVector &fe_function, template template +inline void FEValuesBase:: get_function_2nd_derivatives (const InputVector &fe_function, diff --git a/deal.II/doc/doxygen/headers/fe.h b/deal.II/doc/doxygen/headers/fe.h index 5fbb1c2024..7e429b4514 100644 --- a/deal.II/doc/doxygen/headers/fe.h +++ b/deal.II/doc/doxygen/headers/fe.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2003, 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 @@ -67,7 +67,7 @@ * * The classes in this module are used when one wants to assemble matrices or * vectors. They link finite elements, quadrature objects, and mappins: the - * finite element classes describes a finite element space on a unit cell + * finite element classes describe a finite element space on a unit cell * (i.e. the unit line segment, square, or cube [0,1]^d), the * quadrature classes describe where quadrature points are located and what * weight they have, and the mapping classes describe how to map a point from @@ -80,7 +80,10 @@ * functionality as the FEValues class does for cells. Finally, the * FESubfaceValues class offers the possibility to ingrate on parts of faces * if the neighboring cell is refined and the present cell shares only a part - * of its face with the neighboring cell. + * of its face with the neighboring cell. If vector-valued elements are used, + * the FEValues and related classes allow access to all vector components; if + * one wants to pick individual components, there are extractor classes that + * make this task simpler, as described in the @ref vector_valued module. * * The last member of this group, the UpdateFlags enumeration, is used as an * optimization: instead of letting the FEValues class compute every possible @@ -89,8 +92,8 @@ * UpdateFlags enumeration is used to offer symbolic names denoting what you * want the FEValues class to compute. * - * All these classes are used in all tutorial programs from step-3 onward, and - * are described there in significant detail. + * All these classes are used in all @ref Tutorial "tutorial programs" from + * step-3 onward, and are described there in significant detail. * * The actual workings of the FEValues class and friends is * complicated because it has to be general yet efficient. The page on diff --git a/deal.II/doc/doxygen/headers/glossary.h b/deal.II/doc/doxygen/headers/glossary.h index d925476530..d596548693 100644 --- a/deal.II/doc/doxygen/headers/glossary.h +++ b/deal.II/doc/doxygen/headers/glossary.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 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 @@ -29,8 +29,10 @@ *
Blocks were introduced in BlockVector, * BlockSparseMatrix and related classes. These are used to reflect the * structure of a PDE system in linear algebra, in particular allowing - * for modular solvers. In DoFHandler, this block structure is - * prepared by DoFRenumbering::block_wise(). + * for modular solvers for problems with multiple solution components. + * How to implement this is described in more detail in the + * @ref vector_valued report and the tutorial programs referenced + * therein. * * Originally, this concept was intermixed with the idea of the vector * @ref GlossComponent "component". Since the introduction of diff --git a/deal.II/doc/doxygen/headers/vector_valued.h b/deal.II/doc/doxygen/headers/vector_valued.h new file mode 100644 index 0000000000..b16aa273a4 --- /dev/null +++ b/deal.II/doc/doxygen/headers/vector_valued.h @@ -0,0 +1,740 @@ +//------------------------------------------------------------------------- +// $Id: fe.h 15458 2007-11-07 22:59:17Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//------------------------------------------------------------------------- + + +/** + * @defgroup vector_valued Handling vector valued problems + * + * + * Vector-valued problems are partial differential equations in which the + * solution variable is not a scalar function, but a vector-valued function or + * a set of functions. This includes, for example: + *
    + *
  • The elasticity equation discussed in @ref step_8 "step-8", + * @ref step_17 "step-17", and @ref step_18 "step-18" in which the + * solution is the vector-valued displacement at each point. + *
  • The mixed Laplace equation and its extensions discussed in + * @ref step_20 "step-20", and @ref step_21 "step-21" in which the + * solution is the scalar pressure and the vector-valued velocity + * at each point. + *
  • The Stokes equation and its extensions discussed in + * @ref step_22 "step-22", and @ref step_31 "step-31" in which again + * the solution is the scalar pressure and the vector-valued velocity + * at each point. + *
  • Complex-valued solutions consisting of real and imaginary parts, as + * discussed for example in @ref step_29 "step-29". + *
+ * + * This page gives an overview of how to implement such vector-valued problems + * efficiently in deal.II. + * + * + * + *
%Table of contents
+ *
    + *
  1. @ref VVPhilosophy "A philosophical view of vector-valued problems" + *
  2. @ref VVFEs "Describing finite element spaces" + *
  3. @ref VVAssembling "Assembling linear systems" + *
  4. @ref VVAlternative "An alternative approach" + *
  5. @ref VVBlockSolvers "Block solvers" + *
+ * + * + * + * @anchor VVPhilosophy + *

A philosophical view of vector-valued problems

+ * + * The way one deals systematically with vector-valued problems is not + * fundamentally different from scalar problems in that one needs to come up + * first with a weak (variational) formulation of the problem that takes into + * account all the solution variables. To understand how to do that, let us + * consider a simple example, the mixed Laplace equations discussed in + * @ref step_20 "step-20": +@f{eqnarray*} + \textbf{u} + \nabla p &=& 0, + \\ + -\textrm{div}\; \textbf{u} &=& f, +@f} + * + * Here, we have dim+1 solution components: the scalar pressure + * $p$ and the vector-valued velocity $\textbf u$ with dim vector + * components. + * + * A systematic way to get a weak or variational form for this and other + * vector problems is to first consider it as a problem where the operators + * and solution variables are written in vector and matrix form. For the + * example, this would read as follows: +@f{eqnarray*} + \left( + \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array} + \right) + \left( + \begin{array}{c} \mathbf u \\ p \end{array} + \right) + = + \left( + \begin{array}{c} \mathbf 0 \\ f \end{array} + \right) +@f} + * + * This makes it clear that the solution +@f{eqnarray*} + U = + \left( + \begin{array}{c} \mathbf u \\ p \end{array} + \right) +@f} + * indeed has dim+1 components. We note that we could change the + * ordering of the solution components $\textbf u$ and $p$ inside $U$ if we + * also change columns of the matrix operator. + * + * Next, we need to think about test functions $V$. We want to multiply both + * sides of the equation with them, then integrate over $\Omega$. The result + * should be a scalar equality. We can achieve this by choosing $V$ also + * vector valued as +@f{eqnarray*} + V = + \left( + \begin{array}{c} \mathbf v \\ q \end{array} + \right). +@f} + * + * It is convenient to multiply the matrix-vector equation by the test + * function from the left, since this way we automatically get the correct + * matrix later on (in the linear system, the matrix is also multiplied from + * the right with the solution variable, not from the left), whereas if we + * multiplied from the right then the matrix so assembled is the transpose of + * the one we really want. + * + * With this in mind, let us multiply by $V$ and integrate to get the + * following equation which has to hold for all test functions $V$: +@f{eqnarray*} + \int_\Omega + \left( + \begin{array}{c} \mathbf v \\ q \end{array} + \right)^T + \left( + \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array} + \right) + \left( + \begin{array}{c} \mathbf u \\ p \end{array} + \right) + \ dx + = + \int_\Omega + \left( + \begin{array}{c} \mathbf v \\ q \end{array} + \right)^T + \left( + \begin{array}{c} \mathbf 0 \\ f \end{array} + \right) + \ dx, +@f} + * or equivalently: +@f{eqnarray*} + (\mathbf v, \mathbf u) + + + (\mathbf v, \nabla p) + - + (q, \mathrm{div}\ \mathbf u) + = + (q,f). +@f} + * + * We get the final form by integrating by part the second term: +@f{eqnarray*} + (\mathbf v, \mathbf u) + - + (\mathrm{div}\ \mathbf v, p) + - + (q, \mathrm{div}\ \mathbf u) + = + (q,f) + (\mathbf n\cdot\mathbf v, p)_{\partial\Omega}. +@f} + * + * It is this form that we will later use in assembling the discrete weak form + * into a matrix and a right hand side vector: the form in which we have + * solution and test functions $U,V$ that each consist of a number of vector + * components that we can extract. + * + * + * @anchor VVFEs + *

Describing finite element spaces

+ * + * Once we have settled on this description, we need to find a way to describe + * the vector-valued finite element space from which we draw solution and test + * functions. This is where the FESystem class comes in: it composes + * vector-valued finite element spaces from simpler one. For example, if we + * were to attempt to use $Q_1$ elements for all dim components + * of $\mathbf u$ and the one pressure component $p$, we could use the + * following object: + * @code + * FESystem finite_element (FE_Q(1), dim, + * FE_Q(1), 1); + * @endcode + * + * This means that the final finite element will consist of dim + * components made up of FE_Q elements of degree 1, and another one also of + * degree 1. If course, a simpler (and more efficient) way to achieve the same + * is to use the following form instead: + * @code + * FESystem finite_element (FE_Q(1), dim+1); + * @endcode + * Another way (also not efficient, but making it clear which components + * belong to which element) would be + * @code + * FESystem finite_element (FESystem(FE_Q(1), dim), 1, + * FE_Q(1), 1); + * @endcode + * meaning that we first create a vector element with dim + * components, each consisting of FE_Q elements of order 1. We then couple + * one copy of this already vector-valued element to a single scalar + * copy of the FE_Q element of order 1 which will describe the pressure + * component. + * + * As it turns out these (equivalent) choices do not lead to a stable scheme + * for the mixed Laplace equation. In @ref step_20 "step-20", we therefore use + * a Raviart-Thomas element for the velocities. What exactly this means may be + * of less concern for now except that the FE_RaviartThomas class describes + * elements that already have dim components. For the pressure, + * we use piecewise bi-/trilinear elements that may be discontinuous between + * cells; this is done using the FE_DGQ class. The combined element will then + * be described by + * @code + * FESystem finite_element (FERaviartThomas(1), 1, + * FE_DGQ(1), 1); + * @endcode + * i.e. we combine a single copy of the Raviart-Thomas element with a single + * copy of the element used for the pressure $p$. + * + * In this manner, we can combine the whole vector-valued element from its + * individual components. On a sidenote it may also be worth mentioning why we + * want to also mirror our logical construction of the weak form in our + * program: by creating a finite element object as shown above and then + * handing it off to a DoFHandler object, we get an enumeration of all + * degrees of freedom at once, whether velocity or pressure. Initially, they + * are enumerated in a rather random order, but they can be arranged such that + * first all velocity components are enumerated and then all pressures by, for + * example, calling DoFRenumbering::component_wise . By and large, many things + * become much simpler if we can rely on the fact that there is only a single + * DoFHandler that describes everything. + * + * + * @anchor VVAssembling + *

Assembling linear systems

+ * + * The next step is to assemble the linear system. How to do this for the + * simple case of a scalar problem has been shown in many tutorial programs, + * starting with @ref step_3 "step-3". Here we will show how to do it for + * vector problems. + * + * How to do this is possibly best explained by showing an example + * illustrating how the local contribution of a cell to the weak form of above + * mixed Laplace equations could be assembled. This is essentially how @ref + * step_20 "step-20" does it: + * @code + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + ... + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; q + *
  • The first thing we do is to declare "extractors" (see the + * FEValuesExtractors namespace). These are objects that don't + * do much except store which components of a vector-valued finite + * element constitute a single scalar component, or a tensor of + * rank 1 (i.e. what we call a "physical vector", always consisting + * of dim components). Here, we declare + * an object that represents the velocities consisting of + * dim components starting at component zero, and the + * extractor for the pressure, which is a scalar component at + * position dim. + * + *
  • We then do our usual loop over all cells, shape functions, and + * quadrature points. In the innermost loop, we compute the local + * contribution of a pair of shape functions to the global matrix + * and right hand side vector. Recall that the cell contributions + * to the bilinear form (i.e. neglecting boundary terms) looked as + * follows, based on shape functions + * $V_i=\left(\begin{array}{c}\mathbf v_i q_i\end{array}\right), + * V_j=\left(\begin{array}{c}\mathbf v_j q_j\end{array}\right)$: + @f{eqnarray*} + (\mathbf v_i, \mathbf v_j) + - + (\mathrm{div}\ \mathbf v_i, q_j) + - + (q_i, \mathrm{div}\ \mathbf v_j) + @f} + * whereas the implementation looked like this: + * @code + local_matrix(i,j) += (fe_values[velocities].value (i, q) * + fe_values[velocities].value (j, q) + - + fe_values[velocities].divergence (i, q) * + fe_values[pressure].value (j, q) + - + fe_values[pressure].value (i, q) * + fe_values[velocities].divergence (j, q)) * + fe_values.JxW(q); + * @endcode + * The similarities are pretty obvious. + * + *
  • Essentially, what happens in above code is this: when you do + * fe_values[pressure], a so-called "view" is created, i.e. + * an object that unlike the full FEValues object represents not all + * components of a finite element, but only the one(s) represented by + * the extractor object pressure or + * velocities. + * + *
  • These views can then be asked for information about these individual + * components. For example, when you write + * fe_values[pressure].value(i,q) you get the + * value of the pressure component of the $i$th shape function $V_i$ at + * the $q$th quadrature point. Because the extractor + * pressure represents a scalar component, the results of + * the operator fe_values[pressure].value(i,q) is a scalar + * number. On the other hand, the call + * fe_values[velocities].value(i,q) would produce the + * value of a whole set of dim components, which would + * be of type Tensor@<1,dim@>. + * + *
  • Other things that can be done with views is to ask for the gradient + * of a particular shape function's components described by an + * extractor. For example, fe_values[pressure].gradient(i,q) + * would represent the gradient of the scalar pressure component, which + * is of type Tensor@<1,dim@>, whereas the gradient of the + * velocities components, fe_values[velocities].value(i,q) + * is a Tensor@<2,dim@>, i.e. a matrix $G_{ij}$ that consits + * of entries $G_{ij}=\frac{\partial\phi_i}{\partial x_j}$. Finally, + * both scalar and vector views can be asked for the second derivatives + * ("Hessians"); vector views can also be queried for the + * symmetric gradients as well as the divergence (as we do in the + * example), the latter former a tensor of rank 3, the latter a scalar. + * Finally, vector views can be asked for the symmetric gradient, + * defined as $S_{ij}=\frac 12 \left[\frac{\partial\phi_i}{\partial x_j} + * + \frac{\partial\phi_j}{\partial x_i}\right]$. + * + * Other examples of using extractors and views are shown in tutorial programs + * @ref step_21 "step-21", + * @ref step_31 "step-31", + * @ref step_22 "step-22" and a few other programs. + * + * + * @anchor VVAlternative + *

    An alternative approach

    + * + * There are situations in which one can optimize the assembly of a matrix or + * right hand side vector a bit, using knowledge of the finite element in + * use. Consider, for example, the bilinear form of the elasticity equations + * which we are concerned with first in @ref step_8 "step-8": + * +@f[ + a({\mathbf u}, {\mathbf v}) = + \left( + \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v} + \right)_\Omega + + + \sum_{i,j} + \left( + \mu \partial_i u_j, \partial_i v_j + \right)_\Omega, + + + \sum_{i,j} + \left( + \mu \partial_i u_j, \partial_j v_i + \right)_\Omega, +@f] + * + * Here, $\mathbf u$ is a vector function with dim components, + * $\mathbf v$ the corresponding test function, and $\lambda,\mu$ are material + * parameters. Given our discussions above, the obvious way to implement this + * bilinear form would be as follows, using an extractor object that interprets + * all dim components of the finite element as single vector, + * rather than disjoint components: + * + * @code + const FEValuesExtractors::Vector displacements (0); + + ... + + for (unsigned int q_point=0; q_point phi_i_grad + = fe_values[displacements].gradient (i,q_point); + const double phi_i_div + = fe_values[displacements].divergence (i,q_point); + + for (unsigned int j=0; j phi_j_grad + = fe_values[displacements].gradient (j,q_point); + const double phi_j_div + = fe_values[displacements].divergence (j,q_point); + + cell_matrix(i,j) + += (phi_i_div * phi_j_div * + lambda_values[q_point] + + + scalar_product(phi_i_grad, phi_j_grad) * + mu_values[q_point] + + + scalar_product(phi_i_grad, transpose(phi_j_grad)) * + mu_values[q_point]) + * + fe_values.JxW(q_point); + } + } + * @endcode + * + * The scalar product between two tensors used in this bilinear form is + * implemented as follows: + * + * @code +template +double +scalar_product (const Tensor<2,dim> &u, + const Tensor<2,dim> &v) +{ + double tmp = 0; + for (unsigned int i=0; i phi_i_symmgrad + = fe_values[displacements].symmetric_gradient (i,q_point); + const double phi_i_div + = fe_values[displacements].divergence (i,q_point); + + for (unsigned int j=0; j phi_j_symmgrad + = fe_values[displacements].symmetric_gradient (j,q_point); + const double phi_j_div + = fe_values[displacements].divergence (j,q_point); + + cell_matrix(i,j) + += (phi_i_div * phi_j_div * + lambda_values[q_point] + + + 2 * + (phi_i_symmgrad * phi_j_symmgrad) * + mu_values[q_point]) + * + fe_values.JxW(q_point)); + } + } + * @endcode + * + * So if, again, this is not the code we use in @ref step_8 "step-8", what do + * we do there? The answer rests on the finite element we use. There, we use the + * following element: + * @code + * FESystem finite_element (FE_Q(1), dim); + * @endcode + * In other words, the finite element we use consists of dim copies + * of the same scalar element. This is what we call a @ref GlossPrimitive + * "primitive" element: an element that may be vector-valued but where each + * shape function has exactly one non-zero component. In other words: if the + * $x$-component of a displacement shape function is nonzero, then the $y$- + * and $z$-components must be zero and so on. What this means is that also + * derived quantities based on shape functions inherit this sparsity property. + * For example: the divergence + * $\mathrm{div}\ \Phi(x,y,z)=\partial_x\varphi_x(x,y,z) + + * \partial_y\varphi_y(x,y,z) + \partial_z\varphi_z(x,y,z)$ + * of a vector-valued shape function + * $\Phi(x,y,z)=(\varphi_x(x,y,z), \varphi_y(x,y,z), \varphi_z(x,y,z))^T$ is, + * in the present case, either + * $\mathrm{div}\ \Phi(x,y,z)=\partial_x\varphi_x(x,y,z)$, + * $\mathrm{div}\ \Phi(x,y,z)=\partial_y\varphi_y(x,y,z)$, or + * $\mathrm{div}\ \Phi(x,y,z)=\partial_z\varphi_z(x,y,z)$, because exactly one + * of the $\varphi_\ast$ is nonzero. Knowing this means that we can save a + * number of computations that, if we were to do them, would only yield + * zeros to add up. + * + * In a similar vein, if only one component of a shape function is nonzero, + * then only one row of its gradient $\nabla\Phi$ is nonzero. What this means + * for terms like $(\mu \nabla\Phi_i,\nabla\Phi_j)$, where the scalar product + * between two tensors is defined as + * $(\tau, \gamma)_\Omega=\int_\Omega \sum_{i,j=1}^d \tau_{ij} \gamma_{ij}$, + * is that the term is only nonzero if both tensors have their nonzero + * entries in the same row, which means that the two shape functions have + * to have their single nonzero component in the same location. + * + * If we use this sort of knowledge, then we can in a first step avoid + * computing gradient tensors if we can determine up front that their + * scalar product will be nonzero, in a second step avoid + * building the entire tensors and only get its nonzero components, + * and in a final step simplify the scalar product by only considering + * that index $i$ for the one nonzero row, rather than multiplying and + * adding up zeros. + * + * The vehicle for all this is the ability to determine which vector + * component is going to be nonzero. This information is provided by the + * FiniteElement::system_to_component_index function. What can be done with + * it, using the example above, is explained in detail in + * @ref step_8 "step-8". + * + * + * @anchor VVBlockSolvers + *

    Block solvers

    + * + * Using techniques as shown above, it isn't particularly complicated to + * assemble the linear system, i.e. matrix and right hand side, for a + * vector-valued problem. However, then it also has to be solved. This is more + * complicated. Naively, one could just consider the matrix as a whole. For + * most problems, this matrix is not going to be definite (except for special + * cases like the elasticity equations covered in @ref step_8 "step-8" and + * @ref step_17 "step-17"). It will, often, also not be symmetric. This rather + * general class of matrices presents problems for iterative solvers: the lack + * of structural properties prevents the use of most efficient methods and + * preconditioners. While it can be done, the solution process will therefore + * most often be slower than necessary. + * + * The answer to this problem is to make use of the structure of the + * problem. For example, for the mixed Laplace equations discussed above, the + * operator has the form +@f{eqnarray*} + \left( + \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array} + \right) +@f} + * + * It would be nice if this structure could be recovered in the linear system + * as well. For example, after discretization, we would like to have a matrix + * with the following block structure: +@f{eqnarray*} + \left( + \begin{array}{cc} M & B \\ B^T & 0 \end{array} + \right), +@f} + * where $M$ represents the mass matrix that results from discretizing the + * identity operator $\mathbf 1$ and $B$ the equivalent of the gradient + * operator. + * + * By default, this is not what happens, however. Rather, deal.II assigns + * %numbers to degrees of freedom in a rather random manner. Consequently, if + * you form a vector out of the values of degrees of freedom will not be + * neatly ordered in a vector like +@f{eqnarray*} + \left( + \begin{array}{c} U \\ P \end{array} + \right). +@f} + * Rather, it will be a permutation of this, with %numbers of degrees of + * freedom corresponding to velocities and pressures intermixed. Consequently, + * the system matrix will also not have the nice structure mentioned above, + * but with the same permutation or rows and columns. + * + * What is needed is to re-enumerate degrees of freedom so that velocities + * come first and pressures last. This can be done using the + * DoFRenumbering::component_wise function, as explained in @ref step_20 + * "step-20", @ref step_21 "step-21", @ref step_22 "step-22", and @ref step_31 + * "step-31". After this, at least the degrees of freedom are partitioned + * properly. + * + * But then we still have to make use of it, i.e. we have to come up with a + * solver that uses the structure. For example, in @ref step_20 "step-20", we + * do a block elimination of the linear system +@f{eqnarray*} + \left( + \begin{array}{cc} M & B \\ B^T & 0 \end{array} + \right) + \left( + \begin{array}{c} U \\ P \end{array} + \right) + = + \left( + \begin{array}{c} F \\ G \end{array} + \right). +@f} + * What this system means, of course, is +@f{eqnarray*} + MU + BP &=& F,\\ + B^TU &=& G. +@f} + * + * So, if we multiply the first equation by $B^TM^{-1}$ and subtract the + * second from the result, we get +@f{eqnarray*} + B^TM^{-1}BP &=& B^TM^{-1}F-G. +@f} + * + * This is an equation that now only contains the pressure variables. If we + * can solve it, we can in a second step solve for the velocities using +@f{eqnarray*} + MU = F-BP. +@f} + * + * This has the advantage that the matrices $B^TM^{-1}B$ and $M$ that we have + * to solve with are both symmetric and positive definite, as opposed to the + * large whole matrix we had before. + * + * How a solver like this is implemented is explained in more detail in @ref + * step_20 "step-20", @ref step_31 "step-31", and a few other tutorial + * programs. What we would like to point out here is that we now need a way to + * extract certain parts of a matrix or vector: if we are to multiply, say, + * the $U$ part of the solution vector by the $M$ part of the global matrix, + * then we need to have a way to access these parts of the whole. + * + * This is where the BlockVector, BlockSparseMatrix, and similar classes come + * in. For all practical purposes, then can be used as regular vectors or + * sparse matrices, i.e. they offer element access, provide the usual vector + * operations and implement, for example, matrix-vector multiplications. In + * other words, assembling matrices and right hand sides works in exactly the + * same way as for the non-block versions. That said, internally they store + * the elements of vectors and matrices in "blocks"; for example, instead of + * using one large array, the BlockVector class stores it as a set of arrays + * each of which we call a block. The advantage is that, while the whole thing + * can be used as a vector, one can also access an individual block which + * then, again, is a vector with all the vector operations. + * + * To show how to do this, let us consider the second equation $MU=F-BP$ to be + * solved above. This can be achieved using the following sequence similar to + * what we have in @ref step_20 "step-20": + * @code + Vector tmp (solution.block(0).size()); + system_matrix.block(0,1).vmult (tmp, solution.block(1)); + tmp *= -1; + tmp += system_rhs.block(0); + + + SolverControl solver_control (solution.block(0).size(), + 1e-8*tmp.l2_norm()); + SolverCG<> cg (solver_control, vector_memory); + + cg.solve (system_matrix.block(0,0), + solution.block(0), + tmp, + PreconditionIdentity()); + * @endcode + * + * What's happening here is that we allocate a temporary vector with as many + * elements as the first block of the solution vector, i.e. the velocity + * component $U$, has. We then set this temporary vector equal to the $(0,1)$ + * block of the matrix, i.e. $B$, times component 1 of the solution which is + * the previously computed pressure $P$. The result is multiplied by $-1$, and + * component 0 of the right hand side, $F$ is added to it. The temporary + * vector now contains $F-BP$. The rest of the code snippet simply solves a + * linear system with $F-BP$ as right hand side and the $(0,0)$ block of the + * global matrix, i.e. $M$. Using block vectors and matrices in this way + * therefore allows us to quite easily write rather complicated solvers making + * use of the block structure of a linear system. + * + * @ingroup feall feaccess + */ + diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index a4674816d9..c2ea7c8cfe 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -109,6 +109,15 @@ inconvenience this causes. (Moritz Allmaras, 2007/10/31)

    +
  • New: A significantly simpler way to code the assembly of linear + systems for problems with more than one solution variable has been + implemented. This is explained in detail in the report on @ref vector_valued + and tutorial programs @ref step_20 "step-20" and @ref step_21 "step-21" + have been converted as well. +
    + (WB, 2008/02/23) +

    +
  • Improved: On Mac OS X, the operating system provides for "frameworks", which are essentially collections of shared libraries. We now link with the "Accelerate" framework (from Mac OS X 10.4 diff --git a/deal.II/examples/step-20/doc/intro.dox b/deal.II/examples/step-20/doc/intro.dox index 46f0756d03..97bcd8cda1 100644 --- a/deal.II/examples/step-20/doc/intro.dox +++ b/deal.II/examples/step-20/doc/intro.dox @@ -24,6 +24,9 @@ We are going to extend this tutorial program in @ref step_21 "step-21" to solve not only the mixed Laplace equation, but add another equation that describes the transport of a mixture of two fluids. +The equations covered here fall into the class of vector-valued problems. A +toplevel overview of this topic can be found in the @ref vector_valued module. +

    Formulation, weak form, and discrete problem

    @@ -202,90 +205,84 @@ This is, at best, tedious, error prone, and not dimension independent. There are obvious ways to make things dimension independent, but in the end, the code is simply not pretty. What would be much nicer is if we could simply extract the ${\mathbf u}$ and $p$ components of a shape function $x_h^i$. In the -program we do that, by writing functions like this one: +program we do that in the following way: @code -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - for (unsigned int d=0; dfe_values object, to extract +What this piece of code does is, given an fe_values object, to extract the values of the first $dim$ components of shape function i at quadrature points q, that is the velocity components of that shape function. Put differently, if we write shape functions $x_h^i$ as the tuple $\{{\mathbf u}_h^i,p_h^i\}$, then the function returns the velocity part of this tuple. Note that the velocity is of course a dim-dimensional tensor, and -that the function returns a corresponding object. - -Likewise, we have a function that extracts the pressure component of a shape -function: -@code -template -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - return fe_values.shape_value_component (i,q,dim); -} -@endcode -Finally, the bilinear form contains terms involving the gradients of the -velocity component of shape functions. To be more precise, we are not really -interested in the full gradient, but only the divergence of the velocity -components, i.e. ${\textrm{div}}\ {\mathbf u}_h^i = \sum_{d=0}^{dim-1} -\frac{\partial}{\partial x_d} ({\mathbf u}_h^i)_d$. Here's a function that returns -this quantity: -@code -template -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - double divergence = 0; - for (unsigned int d=0; d::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) { - const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); - const double phi_i_p = extract_p (fe_values, i, q); - - for (unsigned int j=0; j phi_j_u = extract_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); - - local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u - - div_phi_i_u * phi_j_p - - phi_i_p * div_phi_j_u) - * fe_values.JxW(q); - } - - local_rhs(i) += -(phi_i_p * - rhs_values[q] * - fe_values.JxW(q)); - } + fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); + k_inverse.value_list (fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; q phi_i_u = fe_values[velocities].value (i, q); + const double div_phi_i_u = fe_values[velocities].divergence (i, q); + const double phi_i_p = fe_values[pressure].value (i, q); + + for (unsigned int j=0; j phi_j_u = fe_values[velocities].value (j, q); + const double div_phi_j_u = fe_values[velocities].divergence (j, q); + const double phi_j_p = fe_values[pressure].value (j, q); + + local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u + - div_phi_i_u * phi_j_p + - phi_i_p * div_phi_j_u) * + fe_values.JxW(q); + } + + local_rhs(i) += -phi_i_p * + rhs_values[q] * + fe_values.JxW(q); + } @endcode This very closely resembles the form in which we have originally written down the bilinear form and right hand side. @@ -294,40 +291,29 @@ There is one final term that we have to take care of: the right hand side contained the term $(g,{\mathbf v}\cdot {\mathbf n})_{\partial\Omega}$, constituting the weak enforcement of pressure boundary conditions. We have already seen in @ref step_7 "step-7" how to deal with face integrals: essentially exactly the same as with -domain integrals, except that we have to use the FEFaceValues class +domain integrals, except that we have to use the FEFaceValues class instead of FEValues. To compute the boundary term we then simply have -to loop over all boundary faces and integrate there. If you look closely at -the definitions of the extract_* functions above, you will realize -that it isn't even necessary to write new functions that extract the velocity -and pressure components of shape functions from FEFaceValues objects: -both FEValues and FEFaceValues are derived from a common -base class, FEValuesBase, and the extraction functions above can -therefore deal with both in exactly the same way. Assembling the missing -boundary term then takes on the following form: +to loop over all boundary faces and integrate there. The mechanism works in +the same way as above, i.e. the extractor classes also work on FEFaceValues objects: @code -for (unsigned int face_no=0; - face_no::faces_per_cell; - ++face_no) - if (cell->at_boundary(face_no)) - { - fe_face_values.reinit (cell, face_no); - - pressure_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - boundary_values); - - for (unsigned int q=0; q - phi_i_u = extract_u (fe_face_values, i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } - } + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q::value_list (const std::vector > &points, } - // @sect3{extract_u and friends} - - // The next three functions are - // needed for matrix and right hand - // side assembly. They are described - // in detail in the introduction to - // this program, so that we do not - // need to discuss them here again: -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - for (unsigned int d=0; d -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - double divergence = 0; - for (unsigned int d=0; d -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - return fe_values.shape_value_component (i,q,dim); -} - - // @sect3{MixedLaplaceProblem class implementation} @@ -646,6 +598,23 @@ void MixedLaplaceProblem::assemble_system () std::vector boundary_values (n_face_q_points); std::vector > k_inverse_values (n_q_points); + // Finally, we need a couple of extractors + // that we will use to get at the velocity + // and pressure components of vector-valued + // shape functions. Their function and use + // is described in detail in the @ref + // vector_valued report. Essentially, we + // will use them as subscripts on the + // FEValues objects below: the FEValues + // object describes all vector components + // of shape functions, while after + // subscription, it will only refer to the + // velocities (a set of dim + // components starting at component zero) + // or the pressure (a scalar component + // located at position dim): + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); // With all this in place, we can // go on with the loop over all @@ -670,15 +639,15 @@ void MixedLaplaceProblem::assemble_system () for (unsigned int q=0; q phi_i_u = extract_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); - const double phi_i_p = extract_p (fe_values, i, q); + const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); + const double div_phi_i_u = fe_values[velocities].divergence (i, q); + const double phi_i_p = fe_values[pressure].value (i, q); for (unsigned int j=0; j phi_j_u = extract_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); + const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q); + const double div_phi_j_u = fe_values[velocities].divergence (j, q); + const double phi_j_p = fe_values[pressure].value (j, q); local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u - div_phi_i_u * phi_j_p @@ -704,15 +673,10 @@ void MixedLaplaceProblem::assemble_system () for (unsigned int q=0; q - phi_i_u = extract_u (fe_face_values, i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } + local_rhs(i) += -(fe_face_values[velocities].value (i, q) * + fe_face_values.normal_vector(q) * + boundary_values[q] * + fe_face_values.JxW(q)); } // The final step in the loop diff --git a/deal.II/examples/step-21/doc/intro.dox b/deal.II/examples/step-21/doc/intro.dox index 66761658ac..d08c107a41 100644 --- a/deal.II/examples/step-21/doc/intro.dox +++ b/deal.II/examples/step-21/doc/intro.dox @@ -10,6 +10,11 @@ equation. This is therefore also the first time-dependent tutorial program (besides the somewhat strange time-dependence of @ref step_18 "step-18"). +The equations covered here are an extension of the material already covered in +@ref step_20 "step-20". In particular, they fall into the class of +vector-valued problems. A toplevel overview of this topic can be found in the +@ref vector_valued module. +

    The two phase flow problem

    diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index b81f519fca..a5d093599b 100644 --- a/deal.II/examples/step-21/step-21.cc +++ b/deal.II/examples/step-21/step-21.cc @@ -495,80 +495,6 @@ double f_saturation (const double S, - // @sect3{extract_u and friends} - - // More tools: We need methods to extract the - // velocity, pressure, and saturation - // components of finite element shape - // functions. These functions here are - // completely analogous to the ones we have - // already used in step-20: -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - for (unsigned int d=0; d -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - double divergence = 0; - for (unsigned int d=0; d -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - return fe_values.shape_value_component (i,q,dim); -} - - - -template -double extract_s (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - return fe_values.shape_value_component (i,q,dim+1); -} - - - -template -Tensor<1,dim> -extract_grad_s (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - for (unsigned int d=0; d::assemble_system () std::vector > old_solution_values(n_q_points, Vector(dim+2)); std::vector > > old_solution_grads(n_q_points, std::vector > (dim+2)); - + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Scalar saturation (dim+1); + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -908,18 +838,18 @@ void TwoPhaseFlowProblem::assemble_system () { const double old_s = old_solution_values[q](dim+1); - const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); - const double phi_i_p = extract_p (fe_values, i, q); - const double phi_i_s = extract_s (fe_values, i, q); - const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q); + const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); + const double div_phi_i_u = fe_values[velocities].divergence (i, q); + const double phi_i_p = fe_values[pressure].value (i, q); + const double phi_i_s = fe_values[saturation].value (i, q); + const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient(i, q); for (unsigned int j=0; j phi_j_u = extract_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); - const double phi_j_s = extract_s (fe_values, j, q); + const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q); + const double div_phi_j_u = fe_values[velocities].divergence (j, q); + const double phi_j_p = fe_values[pressure].value (j, q); + const double phi_j_s = fe_values[saturation].value (j, q); local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * mobility_inverse(old_s,viscosity) * phi_j_u @@ -952,7 +882,7 @@ void TwoPhaseFlowProblem::assemble_system () for (unsigned int i=0; i - phi_i_u = extract_u (fe_face_values, i, q); + phi_i_u = fe_face_values[velocities].value (i, q); local_rhs(i) += -(phi_i_u * fe_face_values.normal_vector(q) * @@ -1023,6 +953,8 @@ void TwoPhaseFlowProblem::assemble_rhs_S () std::vector local_dof_indices (dofs_per_cell); SaturationBoundaryValues saturation_boundary_values; + + const FEValuesExtractors::Scalar saturation (dim+1); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1049,8 +981,8 @@ void TwoPhaseFlowProblem::assemble_rhs_S () for (unsigned int d=0; d grad_phi_i_s = extract_grad_s(fe_values, i, q); + const double phi_i_s = fe_values[saturation].value (i, q); + const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient (i, q); local_rhs(i) += (time_step * f_saturation(old_s,viscosity) * @@ -1128,7 +1060,7 @@ void TwoPhaseFlowProblem::assemble_rhs_S () : neighbor_saturation[q]), viscosity) * - extract_s(fe_face_values,i,q) * + fe_face_values[saturation].value (i,q) * fe_face_values.JxW(q); } } diff --git a/deal.II/examples/step-22/doc/intro.dox b/deal.II/examples/step-22/doc/intro.dox index 358f177418..dbb770c0de 100644 --- a/deal.II/examples/step-22/doc/intro.dox +++ b/deal.II/examples/step-22/doc/intro.dox @@ -92,6 +92,9 @@ mixing over time scales of many million years, a time scale much shorter than for the same amount of heat to be distributed by thermal conductivity. +The equations covered here fall into the class of vector-valued problems. A +toplevel overview of this topic can be found in the @ref vector_valued module. +

    %Boundary and initial conditions

    @@ -264,119 +267,6 @@ temperature field using the just computed velocity field. -@code -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - tmp[component] = fe_values.shape_value (i,q); - - return tmp; -} - - - -template -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - return fe_values.shape_grad (i,q)[component]; - else - return 0; -} - - - -template -Tensor<2,dim> -extract_grad_s_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<2,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - { - const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2; - - for (unsigned int e=0; e -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim) - return fe_values.shape_value (i,q); - else - return 0; -} - - - -template -double extract_T (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim+1) - return fe_values.shape_value (i,q); - else - return 0; -} - - - -template -Tensor<1,dim> -extract_grad_T (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim+1) - tmp = fe_values.shape_grad (i,q); - - return tmp; -} -@endcode - @code { Vector xx(dof_handler.n_dofs()); diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index fd58cac8ea..c19973e269 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -44,7 +44,6 @@ #include #include #include -#include #include #include @@ -232,122 +231,6 @@ RightHandSide::vector_value (const Point &p, - - -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - tmp[component] = fe_values.shape_value (i,q); - - return tmp; -} - - - -template -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - return fe_values.shape_grad (i,q)[component]; - else - return 0; -} - - - -template -Tensor<2,dim> -extract_grad_s_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<2,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - { - const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2; - - for (unsigned int e=0; e -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim) - return fe_values.shape_value (i,q); - else - return 0; -} - - - -template -double extract_T (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim+1) - return fe_values.shape_value (i,q); - else - return 0; -} - - - -template -Tensor<1,dim> -extract_grad_T (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim+1) - tmp = fe_values.shape_grad (i,q); - - return tmp; -} - - - - template class InverseMatrix : public Subscriptor { @@ -534,19 +417,6 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) } -template -double -scalar_product (const Tensor<2,dim> &a, - const Tensor<2,dim> &b) -{ - double tmp = 0; - for (unsigned int i=0; i void BoussinesqFlowProblem::assemble_system () @@ -592,6 +462,10 @@ void BoussinesqFlowProblem::assemble_system () const double Raleigh_number = 10; + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Scalar temperature (dim+1); + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -610,24 +484,27 @@ void BoussinesqFlowProblem::assemble_system () for (unsigned int i=0; i phi_i_u = extract_u (fe_values, i, q); + const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); if (rebuild_matrices) { - const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); - const double phi_i_p = extract_p (fe_values, i, q); - const double phi_i_T = extract_T (fe_values, i, q); - const Tensor<1,dim> grad_phi_i_T = extract_grad_T(fe_values, i, q); + const SymmetricTensor<2,dim> + phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q); + const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); + const double div_phi_i_u = fe_values[velocities].divergence (i, q); + const double phi_i_p = fe_values[pressure].value (i, q); + const double phi_i_T = fe_values[temperature].value (i, q); + const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient(i, q); for (unsigned int j=0; j phi_j_grads_u = extract_grad_s_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); - const double phi_j_T = extract_T (fe_values, j, q); + const SymmetricTensor<2,dim> + phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q); + const double div_phi_j_u = fe_values[velocities].divergence (j, q); + const double phi_j_p = fe_values[pressure].value (j, q); + const double phi_j_T = fe_values[temperature].value (j, q); - local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u) + local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u + phi_i_p * phi_j_p @@ -660,7 +537,7 @@ void BoussinesqFlowProblem::assemble_system () for (unsigned int i=0; i - phi_i_u = extract_u (fe_face_values, i, q); + phi_i_u = fe_face_values[velocities].value (i, q); local_rhs(i) += -(phi_i_u * fe_face_values.normal_vector(q) * @@ -791,6 +668,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () std::vector local_dof_indices (dofs_per_cell); TemperatureBoundaryValues temperature_boundary_values; + const FEValuesExtractors::Scalar temperature (dim+1); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -817,8 +695,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () for (unsigned int d=0; d grad_phi_i_T = extract_grad_T(fe_values, i, q); + const double phi_i_T = fe_values[temperature].value (i, q); + const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient (i, q); const Point p = fe_values.quadrature_point(q); @@ -901,7 +779,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () old_solution_values_face[q](dim+1) : neighbor_temperature[q]) * - extract_T(fe_face_values,i,q) * + fe_face_values[temperature].value (i,q) * fe_face_values.JxW(q); } } @@ -954,7 +832,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () old_solution_values_face[q](dim+1) : neighbor_temperature[q]) * - extract_T(fe_face_values,i,q) * + fe_face_values[temperature].value (i,q) * fe_face_values.JxW(q); } } @@ -1014,7 +892,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () old_solution_values_face[q](dim+1) : neighbor_temperature[q]) * - extract_T(fe_face_values,i,q) * + fe_face_values[temperature].value (i,q) * fe_face_values.JxW(q); } } @@ -1348,7 +1226,7 @@ void BoussinesqFlowProblem::run () if (timestep_number % 10 == 0) refine_mesh (); } - while (time <= 500); + while (time <= 5); } diff --git a/deal.II/examples/step-29/doc/intro.dox b/deal.II/examples/step-29/doc/intro.dox index 2486a6c5f0..fedf82c278 100644 --- a/deal.II/examples/step-29/doc/intro.dox +++ b/deal.II/examples/step-29/doc/intro.dox @@ -29,6 +29,10 @@ ParameterHandler class first used in @ref step_19 "step-19", which provides a convenient way for reading parameters from a configuration file at runtime without the need to recompile the program code. +The equations covered here fall into the class of vector-valued problems. A +toplevel overview of this topic can be found in the @ref vector_valued module. + +

    Problem setting

    The original purpose of this program is to simulate the focussing properties diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index 34a85e260b..54f260ad8b 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -29,6 +29,9 @@ To be well-posed, we will have to add boundary conditions to the equations. What boundary conditions are readily possible here will become clear once we discuss the weak form of the equations. +The equations covered here fall into the class of vector-valued problems. A +toplevel overview of this topic can be found in the @ref vector_valued module. +

    Weak form

    diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index b486221bc7..a726cb4ed3 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -249,8 +249,8 @@ class RightHandSide : public Function template double -RightHandSide::value (const Point &p, - const unsigned int component) const +RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const { return 0; } @@ -266,120 +266,6 @@ RightHandSide::vector_value (const Point &p, } - - // @sect3{extract_u and friends} - - // The next four functions are needed for - // the assembly of the system matrix and - // the right hand side. They are very similar - // to the ones used in step-20, except - // that we are going to use Q(p+1)Qp elements - // instead of divergence-free Raviart-Thomas - // elements, which simplifies this procedure. - // The only function that is new is - // extract_grad_s_u, which - // gets the symmetric gradient of u. - // As discussed in the introduction, this - // is a second-rank tensor, formed by - // contributions from the gradient and its - // transpose. -template -Tensor<1,dim> -extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<1,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - tmp[component] = fe_values.shape_value (i,q); - - return tmp; -} - - - -template -double -extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - return fe_values.shape_grad (i,q)[component]; - else - return 0; -} - - -template -Tensor<2,dim> -extract_grad_s_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - Tensor<2,dim> tmp; - - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component < dim) - { - const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2; - - for (unsigned int e=0; e -double extract_p (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) -{ - const unsigned int component - = fe_values.get_fe().system_to_component_index(i).first; - - if (component == dim) - return fe_values.shape_value (i,q); - else - return 0; -} - - // @sect4{Inner product of second-rank tensors} - - // In the assembly process, we will need - // to form inner products of second-rank - // tensors. The way how to do this was - // discussed in the introduction - just - // take the sum of the product of the - // individual entries. -template -double -scalar_product (const Tensor<2,dim> &a, - const Tensor<2,dim> &b) -{ - double tmp = 0; - for (unsigned int i=0; i::setup_dofs () // hand side, and global // numbers of the degrees of freedom // for the present cell. - template void StokesProblem::assemble_system () { @@ -787,11 +672,18 @@ void StokesProblem::assemble_system () const PressureBoundaryValues pressure_boundary_values; std::vector boundary_values (n_face_q_points); - + const RightHandSide right_hand_side; std::vector > rhs_values (n_q_points, Vector(dim+1)); + // Next, we need two objects that work as + // extractors for the FEValues + // object. Their use is explained in detail + // in the report on @ref vector_valued : + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + // This starts the loop over all // cells. With the abbreviations // extract_u etc. @@ -810,20 +702,21 @@ void StokesProblem::assemble_system () rhs_values); for (unsigned int q=0; q phi_i_grads_u= extract_grad_s_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); - const double phi_i_p = extract_p (fe_values, i, q); + { + for (unsigned int i=0; i + phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q); + const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); + const double div_phi_i_u = fe_values[velocities].divergence (i, q); + const double phi_i_p = fe_values[pressure].value (i, q); for (unsigned int j=0; j phi_j_grads_u= extract_grad_s_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); - const double phi_j_p = extract_p (fe_values, j, q); - - + const SymmetricTensor<2,dim> + phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q); + const double div_phi_j_u = fe_values[velocities].divergence (j, q); + const double phi_j_p = fe_values[pressure].value (j, q); // Note the way we write the // contributions // phi_i_p * phi_j_p , @@ -837,11 +730,20 @@ void StokesProblem::assemble_system () // is only non-zero when all the // other terms vanish (and the other // way around). - local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u) + // + // Note also that operator* + // is overloaded for + // symmetric tensors, + // yielding the scalar + // product between the two + // tensors in the first + // line: + local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u + phi_i_p * phi_j_p) * fe_values.JxW(q); + } const unsigned int component_i = fe.system_to_component_index(i).first; @@ -871,15 +773,15 @@ void StokesProblem::assemble_system () for (unsigned int q=0; q - phi_i_u = extract_u (fe_face_values, i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } + { + const Tensor<1,dim> + phi_i_u = fe_face_values[velocities].value (i, q); + + local_rhs(i) += -(phi_i_u * + fe_face_values.normal_vector(q) * + boundary_values[q] * + fe_face_values.JxW(q)); + } } // The final step is, as usual, diff --git a/deal.II/examples/step-8/doc/intro.dox b/deal.II/examples/step-8/doc/intro.dox index c7e916b8df..7822a7bcc5 100644 --- a/deal.II/examples/step-8/doc/intro.dox +++ b/deal.II/examples/step-8/doc/intro.dox @@ -4,7 +4,8 @@ In real life, most partial differential equations are really systems of equations. Accordingly, the solutions are usually -vector-valued. The deal.II library supports such problems, and we will show +vector-valued. The deal.II library supports such problems (see the +extensive documentation in the @ref vector_valued module), and we will show that that is mostly rather simple. The only more complicated problems are in assembling matrix and right hand side, but these are easily understood as well. @@ -56,7 +57,7 @@ and the respective bilinear form is then \sum_{i,j} \left( \mu \partial_i u_j, \partial_i v_j - \right)_\Omega, + \right)_\Omega + \sum_{i,j} \left( @@ -74,7 +75,7 @@ or also writing the first term a sum over components: \sum_{k,l} \left( \mu \partial_i u_j, \partial_i v_j - \right)_\Omega, + \right)_\Omega + \sum_{i,j} \left( @@ -83,17 +84,27 @@ or also writing the first term a sum over components: @f] -How do we now assemble the matrix for such an equation? The first thing we -need is some knowledge about how the shape functions work in the case of -vector-valued finite elements. Basically, this comes down to the following: +How do we now assemble the matrix for such an equation? A very long answer +with a number of different alternatives is given in the documentation of the +@ref vector_valued module. Historically, the solution shown below was the only +one available in the early years of the library. It turns out to also be the +fastest. On the other hand, if a few per cent of compute time do not matter, +there are simpler and probably more intuitive ways to assemble the linear +system than the one discussed below but that weren't available until several +years after this tutorial program was first written; if you are interested in +them, take a look at the @ref vector_valued module. + +Let us go back to the question of how to assemble the linear system. The first +thing we need is some knowledge about how the shape functions work in the case +of vector-valued finite elements. Basically, this comes down to the following: let $n$ be the number of shape functions for the scalar finite element of which we build the vector element (for example, we will use bilinear functions for each component of the vector-valued finite element, so the scalar finite -element is the FE_Q(1) element which we have used in previous examples -already, and $n=4$ in two space dimensions). Further, let $N$ be the number of -shape functions for the vector element; in two space dimensions, we need $n$ -shape functions for each component of the vector, so $N=2n$. Then, the $i$th -shape function of the vector element has the form +element is the FE_Q(1) element which we have used in previous +examples already, and $n=4$ in two space dimensions). Further, let $N$ be the +number of shape functions for the vector element; in two space dimensions, we +need $n$ shape functions for each component of the vector, so $N=2n$. Then, +the $i$th shape function of the vector element has the form @f[ \Phi_i({\mathbf x}) = \varphi_{base(i)}({\mathbf x})\ {\mathbf e}_{comp(i)}, @f] @@ -311,13 +322,13 @@ Likewise, the contribution of cell $K$ to the right hand side vector is This is the form in which we will implement the local stiffness matrix and right hand side vectors. -As a final note: in the @ref step_17 "step-17" example program, we will revisit the elastic -problem laid out here, and will show how to solve it in parallel on a cluster -of computers. The resulting program will thus be able to solve this problem to -significantly higher accuracy, and more efficiently if this is -required. In addition, in @ref step_20 "step-20", @ref step_21 "step-21", and -@ref step_22 "step-22", we will revisit some -vector-valued problems and show a few techniques that may make it +As a final note: in the @ref step_17 "step-17" example program, we will +revisit the elastic problem laid out here, and will show how to solve it in +parallel on a cluster of computers. The resulting program will thus be able to +solve this problem to significantly higher accuracy, and more efficiently if +this is required. In addition, in @ref step_20 "step-20", @ref step_21 +"step-21", as well as a few other of the later tutorial programs, we will +revisit some vector-valued problems and show a few techniques that may make it simpler to actually go through all the stuff shown above, with -FiniteElement::system_to_component_index etc. +FiniteElement::system_to_component_index etc. diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile index 849142261c..6eea11aaf9 100644 --- a/tests/deal.II/Makefile +++ b/tests/deal.II/Makefile @@ -74,6 +74,7 @@ tests_x = block_matrices \ create_* \ line_coarsening_3d \ no_flux_* \ + fe_values_view_* \ rt_hessian \ rt_covariant diff --git a/tests/deal.II/fe_values_view_01.cc b/tests/deal.II/fe_values_view_01.cc new file mode 100644 index 0000000000..9a6b189685 --- /dev/null +++ b/tests/deal.II/fe_values_view_01.cc @@ -0,0 +1,115 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 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 +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + + +// test the FEValues views and extractor classes. these tests use a primitive +// finite element and scalar extractors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + deallog << "FE=" << fe.get_name() + << std::endl; + + const QGauss quadrature(2); + FEValues fe_values (fe, quadrature, + update_values | update_gradients | update_hessians); + for (typename DoFHandler::active_cell_iterator + cell = dof.begin_active(); cell != dof.end(); ++cell) + { + fe_values.reinit (cell); + + for (unsigned int c=0; c +void test_hyper_sphere() +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + FESystem fe (FE_Q(1), 1, + FE_Q(2), 2, + FE_DGQ(3), dim); + test(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("fe_values_view_01/output"); + deallog << std::setprecision (2); + + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-12); + + test_hyper_sphere<2>(); + test_hyper_sphere<3>(); +} diff --git a/tests/deal.II/fe_values_view_02.cc b/tests/deal.II/fe_values_view_02.cc new file mode 100644 index 0000000000..23ca64d614 --- /dev/null +++ b/tests/deal.II/fe_values_view_02.cc @@ -0,0 +1,137 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 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 +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + + +// test the FEValues views and extractor classes. these tests use a primitive +// finite element and vector extractors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + deallog << "FE=" << fe.get_name() + << std::endl; + + const QGauss quadrature(2); + FEValues fe_values (fe, quadrature, + update_values | update_gradients | update_hessians); + for (typename DoFHandler::active_cell_iterator + cell = dof.begin_active(); cell != dof.end(); ++cell) + { + fe_values.reinit (cell); + + for (unsigned int c=0; c +void test_hyper_sphere() +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + FESystem fe (FE_Q(1), 1, + FE_Q(2), 2, + FE_DGQ(3), dim); + test(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("fe_values_view_02/output"); + deallog << std::setprecision (2); + + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-12); + + test_hyper_sphere<2>(); + test_hyper_sphere<3>(); +} diff --git a/tests/deal.II/fe_values_view_03.cc b/tests/deal.II/fe_values_view_03.cc new file mode 100644 index 0000000000..cbb9c620ba --- /dev/null +++ b/tests/deal.II/fe_values_view_03.cc @@ -0,0 +1,117 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 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 +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + + +// test the FEValues views and extractor classes. these tests use a non-primitive +// finite element and scalar extractors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + deallog << "FE=" << fe.get_name() + << std::endl; + + const QGauss quadrature(2); + FEValues fe_values (fe, quadrature, + update_values | update_gradients | update_hessians); + for (typename DoFHandler::active_cell_iterator + cell = dof.begin_active(); cell != dof.end(); ++cell) + { + fe_values.reinit (cell); + + for (unsigned int c=0; c +void test_hyper_sphere() +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + FESystem fe (FE_Q(1), 1, + FE_RaviartThomas(1), 1, + FE_Nedelec(1), 1); + test(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("fe_values_view_03/output"); + deallog << std::setprecision (2); + + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-12); + + test_hyper_sphere<2>(); + test_hyper_sphere<3>(); +} diff --git a/tests/deal.II/fe_values_view_04.cc b/tests/deal.II/fe_values_view_04.cc new file mode 100644 index 0000000000..9368301ba2 --- /dev/null +++ b/tests/deal.II/fe_values_view_04.cc @@ -0,0 +1,133 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 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 +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + + +// test the FEValues views and extractor classes. these tests use a non-primitive +// finite element and vector extractors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + deallog << "FE=" << fe.get_name() + << std::endl; + + const QGauss quadrature(2); + FEValues fe_values (fe, quadrature, + update_values | update_gradients | update_hessians); + for (typename DoFHandler::active_cell_iterator + cell = dof.begin_active(); cell != dof.end(); ++cell) + { + fe_values.reinit (cell); + + for (unsigned int c=0; c +void test_hyper_sphere() +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + FESystem fe (FE_Q(1), 1, + FE_RaviartThomas(1), 1, + FE_Nedelec(1), 1); + test(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("fe_values_view_04/output"); + deallog << std::setprecision (2); + + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-12); + + test_hyper_sphere<2>(); + test_hyper_sphere<3>(); +}