]> https://gitweb.dealii.org/ - dealii.git/commitdiff
First part of huge changes to allow projection of the trace of a function to the...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 20 May 1998 17:09:45 +0000 (17:09 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 20 May 1998 17:09:45 +0000 (17:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@324 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/function.h
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/include/numerics/error_estimator.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc

index 8c530205eb7d22815966e89ada100b756756c952..4cdfe88e4db153d0f6e7a3c33c2edd1adc4992bc 100644 (file)
@@ -75,6 +75,7 @@
  *  Once again the warning: do not use the #time# variable for any other
  *  purpose than the intended one! This will inevitably lead to confusion.
  *
+ *
  *  @author Wolfgang Bangerth, 1998
  */
 template <int dim>
index 7d3cd74743eaab742336671f9ce89baf4cd8f6cb..bac5a62d5609e1380f779df16cbb5d07c0cd5cbc 100644 (file)
@@ -5,6 +5,7 @@
 /*----------------------------   dof.h     ---------------------------*/
 
 #include <vector>
+#include <map>
 #include <base/exceptions.h>
 
 
@@ -25,8 +26,8 @@ template <int dim, class Accessor> class TriaIterator;
 template <int dim, class Accessor> class TriaActiveIterator;
 
 template <int dim> class Triangulation;
-
 template <int dim> class FiniteElementBase;
+template <int dim> class Function;
 
 class dVector;
 class dSMatrix;
@@ -447,6 +448,17 @@ class DoFHandler : public DoFDimensionInfo<dim> {
     typedef typename DoFDimensionInfo<dim>::face_iterator face_iterator;
     typedef typename DoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
 
+                                    /**
+                                     * Declare a data type which denotes a
+                                     * mapping between a boundary indicator
+                                     * and the function denoting the boundary
+                                     * values on this part of the boundary.
+                                     * Only one boundary function may be given
+                                     * for each boundary indicator, which is
+                                     * guaranteed by the #map# data type.
+                                     */
+    typedef map<unsigned char,const Function<dim>*> FunctionMap;
+
     
                                     /**
                                      * Constructor. Take #tria# as the
@@ -531,6 +543,57 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      */
     void make_sparsity_pattern (dSMatrixStruct &) const; 
 
+                                    /**
+                                     * Write the sparsity structure of the
+                                     * matrix composed of the basis functions
+                                     * on the boundary into the
+                                     * matrix structure. The sparsity
+                                     * pattern is not compressed, since if
+                                     * you want to call
+                                     * #ConstraintMatrix::condense(1)#
+                                     * afterwards, new entries have to be
+                                     * added. However, if you want to call
+                                     * #ConstraintMatrix::condense(1)#, you
+                                     * have to compress the matrix yourself,
+                                     * using #dSMatrixStruct::compress()#.
+                                     *
+                                     * Since this function is obviously useless
+                                     * in one spatial dimension, it is not
+                                     * implemented.
+                                     */
+    void make_boundary_sparsity_pattern (const vector<int> &dof_to_boundary_mapping,
+                                        dSMatrixStruct &) const; 
+
+                                    /**
+                                     * Write the sparsity structure of the
+                                     * matrix composed of the basis functions
+                                     * on the boundary into the
+                                     * matrix structure. In contrast to the
+                                     * previous function, only those parts
+                                     * of the boundary are considered of which
+                                     * the boundary indicator is listed in the
+                                     * set of numbers passed to this function.
+                                     *
+                                     * In fact, rather than a #set# of boundary
+                                     * indicators, a #map# needs to be passed,
+                                     * since most of the functions handling with
+                                     * boundary indicators take a mapping of
+                                     * boundary indicators and the respective
+                                     * boundary functions. The boundary function,
+                                     * however, is ignored in this function.
+                                     * If you have no functions at hand, but only
+                                     * the boundary indicators, set the function
+                                     * pointers to null pointers.
+                                     *
+                                     * Since this function is obviously useless
+                                     * in one spatial dimension, it is not
+                                     * implemented.
+                                     */
+    void make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators,
+                                        const vector<int> &dof_to_boundary_mapping,
+                                        dSMatrixStruct    &sparsity) const; 
+
+    
                                     /**
                                      * Make up the transfer matrix which
                                      * transforms the data vectors from one
@@ -596,6 +659,17 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      */
     unsigned int max_couplings_between_dofs () const;
 
+                                    /**
+                                     * Return the number of degrees of freedom
+                                     * located on the boundary another dof on
+                                     * the boundary can couple with.
+                                     *
+                                     * The number is the same as for
+                                     * #max_coupling_between_dofs# in one
+                                     * dimension less.
+                                     */
+    unsigned int max_couplings_between_boundary_dofs () const;
+    
                                     /**
                                      * Return the maximum number of entries
                                      * a row in a transfer matrix may contain
@@ -626,6 +700,39 @@ class DoFHandler : public DoFDimensionInfo<dim> {
     void distribute_cell_to_dof_vector (const dVector &cell_data,
                                        dVector       &dof_data) const;
 
+                                    /**
+                                     * Create a mapping from degree of freedom
+                                     * indices to the index of that degree
+                                     * of freedom on the boundary. After this
+                                     * operation, #mapping[dof]# gives the
+                                     * index of the the degree of freedom with
+                                     * global number #dof# in the list of
+                                     * degrees of freedom on the boundary.
+                                     * If the degree of freedom requested is
+                                     * not on the boundary, the value of
+                                     * #mapping[dof]# is #-1#. This function is
+                                     * mainly used when setting up matrices and
+                                     * vectors on the boundary from the ansatz
+                                     * functions, which have global numbers,
+                                     * while the matrices and vectors use
+                                     * numbers of the ansatz functions local
+                                     * to the boundary.
+                                     *
+                                     * Prior content of #mappin# is deleted.
+                                     *
+                                     * This function is not implemented for
+                                     * one dimension.
+                                     *
+                                     * DOC FIXME: algorithm
+                                     */
+    void map_dof_to_boundary_indices (vector<int> &mapping) const;
+
+                                    /**
+                                     * DOC FIXME: algorithm
+                                     */
+    void map_dof_to_boundary_indices (const FunctionMap &boundary_indicators,
+                                     vector<int> &mapping) const;
+
                                     /**
                                      *  @name Cell iterator functions
                                      */
@@ -1001,6 +1108,23 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      */
     unsigned int n_dofs () const;
 
+                                    /**
+                                     * Return the number of degrees of freedom
+                                     * located on the boundary.
+                                     */
+    unsigned int n_boundary_dofs () const;
+
+                                    /**
+                                     * Return the number of degrees of freedom
+                                     * located on those parts of the boundary
+                                     * which have a boundary indicator listed
+                                     * in the given set. The reason that a
+                                     * #map# rather than a #set# is used is the
+                                     * same as descibed in the section on the
+                                     * #make_boundary_sparsity_pattern# function.
+                                     */
+    unsigned int n_boundary_dofs (const FunctionMap &boundary_indicators) const;
+
                                     /**
                                      * Return a constant reference to the
                                      * selected finite element object.
@@ -1047,6 +1171,10 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                     /**
                                      * Exception
                                      */
+    DeclException0 (ExcInvalidBoundaryIndicator);
+                                    /**
+                                     * Exception
+                                     */
     DeclException0 (ExcInternalError);
                                     /**
                                      * Exception
index 82e76641d08a7a38e19f70bcc0867bc2d571ba52..42c5fea0a5efcc396882dda3d941594082b2e0d7 100644 (file)
@@ -102,7 +102,7 @@ class ProblemBase {
                                      * See the general documentation of this
                                      * class for more detail.
                                      */
-    typedef map<unsigned char,Function<dim>*> FunctionMap;
+    typedef map<unsigned char,const Function<dim>*> FunctionMap;
                                     /**
                                      * Typdedef an iterator which assembles
                                      * matrices and vectors.
index 88d929eff9ffdb542af89ab42ce3dec2b597a9d2..4ab4f49b455534a83456442014101bc6a75f5e02 100644 (file)
@@ -153,7 +153,7 @@ class KellyErrorEstimator {
                                      * for each boundary indicator, which is
                                      * guaranteed by the #map# data type.
                                      */
-    typedef map<unsigned char,Function<dim>*> FunctionMap;
+    typedef map<unsigned char,const Function<dim>*> FunctionMap;
     
     void estimate_error (const DoFHandler<dim>    &dof,
                         const Quadrature<dim-1>  &quadrature,
index 501f29cd63211bdeebeef815bb854a2ab6c1c53b..c61dccc1508bc5c0baa708ba7e011bc52f76b24a 100644 (file)
@@ -8,7 +8,7 @@
 
 #include <base/exceptions.h>
 #include <map>
-
+#include <set>
 
 template <int dim> class Triangulation;
 template <int dim> class DoFHandler;
@@ -80,6 +80,31 @@ class dSMatrix;
  * side.
  *
  *
+ * \subsection{Matrices on the boundary}
+ *
+ * The #create_boundary_mass_matrix# creates the matrix with entries
+ * $m_{ij} = \int_{\Gamma} \phi_i \phi_j dx$, where $\Gamma$ is the union
+ * of boundary parts with indicators contained in a set passed to the
+ * function (i.e. if you want to set up the mass matrix for the parts of
+ * the boundary with indicators zero and 2, you pass the function a set
+ * of #unsigned char#s as parameter #boundary_parts# containing the elements
+ * zero and 2). The $\phi_i$ are the basis functions which have at least
+ * part of their support om $\Gamma$. The mapping between row and column
+ * indices in the mass matrix and the right hand side and the global degree
+ * of freedom numbers of the respective basis functions on the whole domain
+ * is returned as a vector of numbers which has the same size as the dimension
+ * of matrix and right hand side.
+ *
+ * Since in most cases we are not interested in the pure mass matrix on the
+ * boundary, but rather need it to compute the projection of a function to
+ * the boundary, no function is provided to only create the matrix.
+ *
+ * This function needs to get passed a matrix object to hold the resulting sparse
+ * matrix. This object is supposed to be initialized with a suitable sparsity
+ * pattern, which can be created using the
+ * #DoFHandler<>::make_boundary_sparsity_pattern# function.
+ *
+ *
  * \subsection{Right hand sides}
  *
  * In many cases, you will not only want to build the matrix, but also
@@ -102,6 +127,20 @@ class dSMatrix;
 template <int dim>
 class MatrixCreator {
   public:
+                                    /**
+                                     * Declare a data type which denotes a
+                                     * mapping between a boundary indicator
+                                     * and the function denoting the boundary
+                                     * values on this part of the boundary.
+                                     * Only one boundary function may be given
+                                     * for each boundary indicator, which is
+                                     * guaranteed by the #map# data type.
+                                     * 
+                                     * See the general documentation of this
+                                     * class for more detail.
+                                     */
+    typedef map<unsigned char,const Function<dim>*> FunctionMap;
+
                                     /**
                                      * Assemble the mass matrix. If no 
                                      * coefficient is given, it is assumed
@@ -135,6 +174,31 @@ class MatrixCreator {
                                    dVector                  &rhs_vector,
                                    const Function<dim>      *a = 0);
 
+                                    /**
+                                     * Assemble the mass matrix and a right
+                                     * hand side vector along the boundary.
+                                     * If no 
+                                     * coefficient is given, it is assumed
+                                     * to be constant one.
+                                     *
+                                     * The matrix is assumed to already be
+                                     * initialized with a suiting sparsity
+                                     * pattern (the #DoFHandler# provides an
+                                     * appropriate function).
+                                     *
+                                     * See the general doc of this class
+                                     * for more information.
+                                     */
+    static void create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
+                                            const FiniteElement<dim> &fe,
+                                            const Quadrature<dim-1>  &q,
+                                            const Boundary<dim>      &boundary,
+                                            dSMatrix                 &matrix,
+                                            const FunctionMap        &rhs,
+                                            dVector                  &rhs_vector,
+                                            vector<int>              &vec_to_dof_mapping,
+                                            const Function<dim>      *a = 0);
+
                                     /**
                                      * Assemble the Laplace matrix. If no 
                                      * coefficient is given, it is assumed
@@ -184,8 +248,8 @@ class MatrixCreator {
  * The #apply_boundar_values# function inserts boundary conditions of
  * into a system of equations.  To actually do this you have to specify
  * a list of degree of freedom indices along with the value this degree of
- * freedom shall assume. To see how to get such a list, see below in the
- * discussion of the #interpolate_boundary_values# function.
+ * freedom shall assume. To see how to get such a list, see the discussion
+ * of the #VectorTools::interpolate_boundary_values# function.
  *
  * The inclusion into the assemblage process is as follows: when the matrix and
  * vectors are set up, a list of nodes subject to dirichlet bc is made and
@@ -220,7 +284,7 @@ class MatrixCreator {
  * It it not clear whether the deletion of coupling between the boundary degree
  * of freedom and other dofs really forces the corresponding entry in the
  * solution vector to have the right value when using iterative solvers,
- * since their search directions may contains components in the direction
+ * since their search directions may contain components in the direction
  * of the boundary node. For this reason, we perform a very simple line
  * balancing by not setting the main diagonal entry to unity, but rather
  * to the value it had before deleting this line, or to the first nonzero
@@ -238,55 +302,12 @@ class MatrixCreator {
  * the correct value after solving again. This question is an open one as of
  * now and may be answered by future experience.
  *
- *
- * \subsection{Getting a list of boundary values}
- *
- * As discussed above, the #apply_boundary_values# function takes a list
- * of boundary nodes and their values. You can get such a list by interpolation
- * of a boundary function using the #interpolate_boundary_values# function.
- * To use it, you have to
- * specify a list of pairs of boundary indicators (of type #unsigned char#;
- * see the section in the documentation of the \Ref{Triangulation} class for more
- * details) and the according functions denoting the dirichlet boundary values
- * of the nodes on boundary faces with this boundary indicator.
- *
- * Usually, all other boundary conditions, such as inhomogeneous Neumann values
- * or mixed boundary conditions are handled in the weak formulation. No attempt
- * is made to include these into the process of assemblage therefore.
- *
- * Within this function, boundary values are interpolated, i.e. a node is given
- * the point value of the boundary function. In some cases, it may be necessary
- * to use the L2-projection of the boundary function or any other method. For
- * this purpose other functions exist in the #MatrixTools# library (or will
- * exist at least).
- *
- * You should be aware that the boundary function may be evaluated at nodes
- * on the interior of faces. These, however, need not be on the true
- * boundary, but rather are on the approximation of the boundary represented
- * by teh mapping of the unit cell to the real cell. Since this mapping will
- * in most cases not be the exact one at the face, the boundary function is
- * evaluated at points which are not on the boundary and you should make
- * sure that the returned values are reasonable in some sense anyway.
  * 
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim>
 class MatrixTools : public MatrixCreator<dim> {
   public:
-                                    /**
-                                     * Declare a data type which denotes a
-                                     * mapping between a boundary indicator
-                                     * and the function denoting the boundary
-                                     * values on this part of the boundary.
-                                     * Only one boundary function may be given
-                                     * for each boundary indicator, which is
-                                     * guaranteed by the #map# data type.
-                                     * 
-                                     * See the general documentation of this
-                                     * class for more detail.
-                                     */
-    typedef map<unsigned char,Function<dim>*> FunctionMap;
-
                                     /**
                                      * Apply dirichlet boundary conditions
                                      * to the system matrix and vectors
@@ -298,27 +319,6 @@ class MatrixTools : public MatrixCreator<dim> {
                                       dVector               &solution,
                                       dVector               &right_hand_side);
 
-                                    /**
-                                     * Make up the list of node subject
-                                     * to Dirichlet boundary conditions
-                                     * and the values they are to be
-                                     * assigned, by interpolation around
-                                     * the boundary.
-                                     *
-                                     * See the general doc for more
-                                     * information.
-                                     */
-    static void interpolate_boundary_values (const DoFHandler<dim> &dof,
-                                            const FunctionMap     &dirichlet_bc,
-                                            const FiniteElement<dim> &fe,
-                                            const Boundary<dim> &boundary,
-                                            map<int,double>     &boundary_values);
-    
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInvalidBoundaryIndicator);
                                     /**
                                      * Exception
                                      */
index 12641af42bbf3c333b73e12d9b0ac0ee9f4b956b..c5a7ff36dbd42a0ec3a67ca5408eb181474c8460 100644 (file)
@@ -7,6 +7,7 @@
 
 
 #include <base/exceptions.h>
+#include <map>
 
 template <int dim> class DoFHandler;
 template <int dim> class Function;
@@ -81,11 +82,70 @@ enum NormType {
  *   $f_i = \int_\Omega f(x) \phi_i(x) dx$. The solution vector $v$ then is
  *   the projection.
  *
+ *   In order to get proper results, it is necessary to treat boundary
+ *   conditions right. This is done by $L_2$-projection of the trace of the
+ *   given function onto the finite element space restricted to the boundary
+ *   of the domain, then taking this information and using it to eliminate
+ *   the boundary nodes from the mass matrix of the whole domain, using the
+ *   #MatrixTools::apply_boundary_values# function. The projection of the
+ *   trace of the function to the boundary is done with the
+ *   #VectorTools::project_boundary_values# (see below) function. You may
+ *   specify a flag telling the projection that the function has zero boundary
+ *   values, in which case the $L_2$-projection onto the boundary is not
+ *   needed. If it is needed, the #VectorTools::project_boundary_values# is
+ *   called with a map of boundary functions of which all boundary indicators
+ *   from zero to 254 (255 is used for other purposes, see the #Triangulation#
+ *   class documentation) point to the function to be projected. The projection
+ *   to the boundary takes place using a second quadrature formula given to
+ *   the #project# function.
+ *
  *   The solution of the linear system is presently done using a simple CG
  *   method without preconditioning and without multigrid. This is clearly not
  *   too efficient, but sufficient in many cases and simple to implement. This
  *   detail may change in the future.
  *
+ * \item Interpolation of boundary values:
+ *   The #MatrixTools::apply_boundary_values# function takes a list
+ *   of boundary nodes and their values. You can get such a list by interpolation
+ *   of a boundary function using the #interpolate_boundary_values# function.
+ *   To use it, you have to
+ *   specify a list of pairs of boundary indicators (of type #unsigned char#;
+ *   see the section in the documentation of the \Ref{Triangulation} class for more
+ *   details) and the according functions denoting the dirichlet boundary values
+ *   of the nodes on boundary faces with this boundary indicator.
+ *
+ *   Usually, all other boundary conditions, such as inhomogeneous Neumann values
+ *   or mixed boundary conditions are handled in the weak formulation. No attempt
+ *   is made to include these into the process of assemblage therefore.
+ *
+ *   Within this function, boundary values are interpolated, i.e. a node is given
+ *   the point value of the boundary function. In some cases, it may be necessary
+ *   to use the L2-projection of the boundary function or any other method. For
+ *   this purpose other functions refer to the #VectorTools::project_boundary_values#
+ *   function below.
+ *
+ *   You should be aware that the boundary function may be evaluated at nodes
+ *   on the interior of faces. These, however, need not be on the true
+ *   boundary, but rather are on the approximation of the boundary represented
+ *   by the mapping of the unit cell to the real cell. Since this mapping will
+ *   in most cases not be the exact one at the face, the boundary function is
+ *   evaluated at points which are not on the boundary and you should make
+ *   sure that the returned values are reasonable in some sense anyway.
+ *
+ * \item Projection of boundary values:
+ *   The #project_boundary_values# function acts similar to the
+ *   #interpolate_boundary_values# function, apart from the fact that it does
+ *   not get the nodal values of boundary nodes by interpolation but rather
+ *   through the $L_2$-projection of the trace of the function to the boundary.
+ *
+ *   The projection takes place on all boundary parts with boundary indicators
+ *   listed in the map of boundary functions. These boundary parts may or may
+ *   not be contiguous. For these boundary parts, the mass matrix is assembled
+ *   using the #MatrixTools::create_boundary_mass_matrix# function, as well as
+ *   the appropriate right hand side. Then the resulting system of equations is
+ *   solved using a simple CG method (without preconditioning), which is in most
+ *   cases sufficient for the present purpose.
+ *
  * \item Computing errors:
  *   The function #integrate_difference# performs the calculation of the error
  *   between the finite element solution and a given (continuous) reference
@@ -133,7 +193,7 @@ enum NormType {
  *   To get the {\it global} L_1 error, you have to sum up the entries in
  *   #difference#, e.g. using #dVector::l1_norm# function.
  *   For the global L_2 difference, you have to sum up the squares of the
- *   entries and take the root of the sum, e.g. using #dVector::l2_norm.
+ *   entries and take the root of the sum, e.g. using #dVector::l2_norm#.
  *   These two operations represent the
  *   l_1 and l_2 norms of the vectors, but you need not take the absolute
  *   value of each entry, since the cellwise norms are already positive.
@@ -149,6 +209,20 @@ enum NormType {
 template <int dim>
 class VectorTools {
   public:
+                                    /**
+                                     * Declare a data type which denotes a
+                                     * mapping between a boundary indicator
+                                     * and the function denoting the boundary
+                                     * values on this part of the boundary.
+                                     * Only one boundary function may be given
+                                     * for each boundary indicator, which is
+                                     * guaranteed by the #map# data type.
+                                     * 
+                                     * See the general documentation of this
+                                     * class for more detail.
+                                     */
+    typedef map<unsigned char,const Function<dim>*> FunctionMap;
+
                                     /**
                                      * Compute the interpolation of
                                      * #function# at the ansatz points to
@@ -174,10 +248,43 @@ class VectorTools {
                         const ConstraintMatrix   &constraints,
                         const FiniteElement<dim> &fe,
                         const Quadrature<dim>    &q,
+                        const Quadrature<dim-1>  &q_boundary,
                         const Boundary<dim>      &boundary,
                         const Function<dim>      &function,
+                        const bool                has_zero_boundary,
                         dVector                  &vec);
 
+                                    /**
+                                     * Make up the list of node subject
+                                     * to Dirichlet boundary conditions
+                                     * and the values they are to be
+                                     * assigned, by interpolation around
+                                     * the boundary.
+                                     *
+                                     * See the general doc for more
+                                     * information.
+                                     */
+    static void interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                            const FunctionMap     &dirichlet_bc,
+                                            const FiniteElement<dim> &fe,
+                                            const Boundary<dim> &boundary,
+                                            map<int,double>     &boundary_values);
+    
+                                    /**
+                                     * Project #function# to the boundary
+                                     * of the domain, using the given quadrature
+                                     * formula for the faces.
+                                     *
+                                     * See the general documentation of this
+                                     * class for further information.
+                                     */
+    static void project_boundary_values (const DoFHandler<dim>    &dof,
+                                        const FunctionMap        &boundary_functions,
+                                        const FiniteElement<dim> &fe,
+                                        const Quadrature<dim-1>  &q,
+                                        const Boundary<dim>      &boundary,
+                                        map<int,double>          &boundary_values);
+    
                                     /**
                                      * Integrate the difference between
                                      * a finite element function and
@@ -205,6 +312,10 @@ class VectorTools {
                                      * Exception
                                      */
     DeclException0 (ExcInvalidFE);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidBoundaryIndicator);
 };
 
 
index 68326925b7f2f773cebc98183b44cf06bd391162..fbb3c795dcf94555a77a3ee85daf22ba889885b7 100644 (file)
@@ -10,6 +10,7 @@
 #include <fe/fe.h>
 #include <lac/dsmatrix.h>
 #include <map>
+#include <set>
 #include <algorithm>
 
 
@@ -625,6 +626,68 @@ unsigned int DoFHandler<dim>::n_dofs () const {
 
 
 
+unsigned int DoFHandler<1>::n_boundary_dofs () const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <int dim>
+unsigned int DoFHandler<dim>::n_boundary_dofs () const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  
+  set<int> boundary_dofs;
+
+  const unsigned int dofs_per_face = selected_fe->dofs_per_face;
+  vector<int> dofs_on_face(dofs_per_face);
+  active_face_iterator face = begin_active_face (),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (face->at_boundary())
+      {
+       face->get_dof_indices (dofs_on_face);
+       for (unsigned int i=0; i<dofs_per_face; ++i)
+         boundary_dofs.insert(dofs_on_face[i]);
+      };
+  return boundary_dofs.size();
+};    
+
+
+
+unsigned int DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <int dim>
+unsigned int DoFHandler<dim>::n_boundary_dofs (const FunctionMap &boundary_indicators) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (boundary_indicators.find(255) == boundary_indicators.end(),
+         ExcInvalidBoundaryIndicator());
+  
+  set<int> boundary_dofs;
+
+  const unsigned int dofs_per_face = selected_fe->dofs_per_face;
+  vector<int> dofs_on_face(dofs_per_face);
+  active_face_iterator face = begin_active_face (),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (boundary_indicators.find(face->boundary_indicator()) !=
+       boundary_indicators.end())
+      {
+       face->get_dof_indices (dofs_on_face);
+       for (unsigned int i=0; i<dofs_per_face; ++i)
+         boundary_dofs.insert(dofs_on_face[i]);
+      };
+  return boundary_dofs.size();
+};    
+
+
+
 template <int dim>
 const Triangulation<dim> & DoFHandler<dim>::get_tria () const {
   return *tria;
@@ -1059,77 +1122,129 @@ void DoFHandler<2>::make_constraint_matrix (ConstraintMatrix &constraints) const
 
 
 
-void DoFHandler<1>::make_sparsity_pattern (dSMatrixStruct &sparsity) const {
+template <int dim>
+void DoFHandler<dim>::make_sparsity_pattern (dSMatrixStruct &sparsity) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (sparsity.n_rows() == n_dofs(),
+         ExcDifferentDimensions (sparsity.n_rows(), n_dofs()));
+  Assert (sparsity.n_cols() == n_dofs(),
+         ExcDifferentDimensions (sparsity.n_cols(), n_dofs()));
+
+  const unsigned int total_dofs = selected_fe->total_dofs;
+  vector<int> dofs_on_this_cell(total_dofs);
   active_cell_iterator cell = begin_active(),
                       endc = end();
-
-                                  // set up an array which dofs are used
-                                  // on a specific cell
-  unsigned int *dofs_on_this_cell = new unsigned int[selected_fe->total_dofs];
-  
   for (; cell!=endc; ++cell) 
     {
-      unsigned int dof_number=0;
-
-                                      // fill dof indices on vertices
-      for (unsigned int vertex=0; vertex<2; ++vertex)
-       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
-         dofs_on_this_cell[dof_number++] = cell->vertex_dof_index (vertex,d);
-
-                                      // fill dof indices on line
-      for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
-       dofs_on_this_cell[dof_number++] = cell->dof_index (d);
-
+      cell->get_dof_indices (dofs_on_this_cell);
                                       // make sparsity pattern for this cell
-      for (unsigned int i=0; i<selected_fe->total_dofs; ++i)
-       for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
+      for (unsigned int i=0; i<total_dofs; ++i)
+       for (unsigned int j=0; j<total_dofs; ++j)
          sparsity.add (dofs_on_this_cell[i],
                        dofs_on_this_cell[j]);
     };
+};
+
+
 
-  delete[] dofs_on_this_cell;
+void DoFHandler<1>::make_boundary_sparsity_pattern (const vector<int> &,
+                                                   dSMatrixStruct &) const {
+    Assert (selected_fe != 0, ExcNoFESelected());
+    Assert (false, ExcInternalError());
 };
 
 
 
-void DoFHandler<2>::make_sparsity_pattern (dSMatrixStruct &sparsity) const {
-  active_cell_iterator cell = begin_active(),
-                      endc = end();
+template <int dim>
+void DoFHandler<dim>::make_boundary_sparsity_pattern (const vector<int> &dof_to_boundary_mapping,
+                                                     dSMatrixStruct &sparsity) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (dof_to_boundary_mapping.size() == n_dofs(), ExcInternalError());
+  Assert (sparsity.n_rows() == n_boundary_dofs(),
+         ExcDifferentDimensions (sparsity.n_rows(), n_boundary_dofs()));
+  Assert (sparsity.n_cols() == n_boundary_dofs(),
+         ExcDifferentDimensions (sparsity.n_cols(), n_boundary_dofs()));
+  Assert (*max_element(dof_to_boundary_mapping.begin(),
+                      dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1,
+         ExcInternalError());
+
+  const unsigned int total_dofs = selected_fe->dofs_per_face;
+  vector<int> dofs_on_this_face(total_dofs);
+  active_face_iterator face = begin_active_face(),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (face->at_boundary())
+      {
+       face->get_dof_indices (dofs_on_this_face);
 
-                                  // set up an array which dofs are used
-                                  // on a specific cell
-  unsigned int      *dofs_on_this_cell = new unsigned int[selected_fe->total_dofs];
+                                        // make sure all dof indices have a
+                                        // boundary index
+       Assert (*min_element(dofs_on_this_face.begin(),
+                            dofs_on_this_face.end()) >=0,
+               ExcInternalError());
+       
+                                        // make sparsity pattern for this cell
+       for (unsigned int i=0; i<total_dofs; ++i)
+         for (unsigned int j=0; j<total_dofs; ++j)
+           sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
+                         dof_to_boundary_mapping[dofs_on_this_face[j]]);
+      };
+};
 
-  
-  for (; cell!=endc; ++cell) 
-    {
-      int dof_number=0;
 
-                                      // fill dof indices on vertices
-      for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
-       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
-         dofs_on_this_cell[dof_number++] = cell->vertex_dof_index (vertex,d);
 
-      for (unsigned int line=0; line<GeometryInfo<2>::faces_per_cell; ++line)
-       for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
-         dofs_on_this_cell[dof_number++] = cell->line(line)->dof_index (d);
-      
-                                      // fill dof indices on quad
-      for (unsigned int d=0; d<selected_fe->dofs_per_quad; ++d)
-       dofs_on_this_cell[dof_number++] = cell->dof_index (d);
+void DoFHandler<1>::make_boundary_sparsity_pattern (const FunctionMap &,
+                                                   const vector<int> &,
+                                                   dSMatrixStruct &) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcInternalError());
+};
 
-                                      // make sparsity pattern for this cell
-      for (unsigned int i=0; i<selected_fe->total_dofs; ++i)
-       for (unsigned int j=0; j<selected_fe->total_dofs; ++j)
-         sparsity.add (dofs_on_this_cell[i],
-                       dofs_on_this_cell[j]);
-    };
-  
-  delete[] dofs_on_this_cell;
+
+
+template <int dim>
+void DoFHandler<dim>::make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators,
+                                                     const vector<int> &dof_to_boundary_mapping,
+                                                     dSMatrixStruct &sparsity) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (dof_to_boundary_mapping.size() == n_dofs(), ExcInternalError());
+  Assert (boundary_indicators.find(255) == boundary_indicators.end(),
+         ExcInvalidBoundaryIndicator());
+  Assert (sparsity.n_rows() == n_boundary_dofs(boundary_indicators),
+         ExcDifferentDimensions (sparsity.n_rows(), n_boundary_dofs(boundary_indicators)));
+  Assert (sparsity.n_cols() == n_boundary_dofs(boundary_indicators),
+         ExcDifferentDimensions (sparsity.n_cols(), n_boundary_dofs(boundary_indicators)));
+  Assert (*max_element(dof_to_boundary_mapping.begin(),
+                      dof_to_boundary_mapping.end()) == (signed int)sparsity.n_rows()-1,
+         ExcInternalError());
+
+  const unsigned int total_dofs = selected_fe->dofs_per_face;
+  vector<int> dofs_on_this_face(total_dofs);
+  active_face_iterator face = begin_active_face(),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (boundary_indicators.find(face->boundary_indicator()) !=
+       boundary_indicators.end())
+      {
+       face->get_dof_indices (dofs_on_this_face);
+
+                                        // make sure all dof indices have a
+                                        // boundary index
+       Assert (*min_element(dofs_on_this_face.begin(),
+                            dofs_on_this_face.end()) >=0,
+               ExcInternalError());
+                                        // make sparsity pattern for this cell
+       for (unsigned int i=0; i<total_dofs; ++i)
+         for (unsigned int j=0; j<total_dofs; ++j)
+           sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
+                         dof_to_boundary_mapping[dofs_on_this_face[j]]);
+      };
 };
 
 
 
+
+
 template <int dim>
 void DoFHandler<dim>::make_transfer_matrix (const DoFHandler<dim> &transfer_from,
                                            dSMatrixStruct        &transfer_pattern) const {
@@ -1406,6 +1521,22 @@ unsigned int DoFHandler<2>::max_couplings_between_dofs () const {
 
 
 
+unsigned int DoFHandler<1>::max_couplings_between_boundary_dofs () const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcInternalError());
+  return 0;
+};
+
+
+
+unsigned int DoFHandler<2>::max_couplings_between_boundary_dofs () const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  return 3*selected_fe->dofs_per_vertex + 2*selected_fe->dofs_per_line;
+};
+
+
+
+
 unsigned int DoFHandler<1>::max_transfer_entries (const unsigned int max_level_diff) const {
   Assert (max_level_diff<2, ExcOnlyOnelevelTransferImplemented());
   switch (max_level_diff)
@@ -1481,6 +1612,81 @@ void DoFHandler<dim>::distribute_cell_to_dof_vector (const dVector &cell_data,
 
 
 
+void DoFHandler<1>::map_dof_to_boundary_indices (vector<int> &) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void DoFHandler<dim>::map_dof_to_boundary_indices (vector<int> &mapping) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+
+  mapping.clear ();
+  mapping.insert (mapping.end(), n_dofs(), -1);
+
+  const unsigned int dofs_per_face = selected_fe->dofs_per_face;
+  vector<int> dofs_on_face(dofs_per_face);
+  int next_boundary_index = 0;
+  
+  active_face_iterator face = begin_active_face(),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (face->at_boundary()) 
+      {
+       face->get_dof_indices (dofs_on_face);
+       for (unsigned int i=0; i<dofs_per_face; ++i)
+         if (mapping[dofs_on_face[i]] == -1)
+           mapping[dofs_on_face[i]] = next_boundary_index++;
+      };
+
+  Assert (static_cast<unsigned int>(next_boundary_index) == n_boundary_dofs(),
+         ExcInternalError());
+};
+
+
+
+void DoFHandler<1>::map_dof_to_boundary_indices (const FunctionMap &,
+                                                vector<int> &) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void DoFHandler<dim>::map_dof_to_boundary_indices (const FunctionMap &boundary_indicators,
+                                                  vector<int> &mapping) const {
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (boundary_indicators.find(255) == boundary_indicators.end(),
+         ExcInvalidBoundaryIndicator());
+
+  mapping.clear ();
+  mapping.insert (mapping.end(), n_dofs(), -1);
+
+  const unsigned int dofs_per_face = selected_fe->dofs_per_face;
+  vector<int> dofs_on_face(dofs_per_face);
+  int next_boundary_index = 0;
+  
+  active_face_iterator face = begin_active_face(),
+                      endf = end_face();
+  for (; face!=endf; ++face)
+    if (boundary_indicators.find(face->boundary_indicator()) !=
+       boundary_indicators.end())
+      {
+       face->get_dof_indices (dofs_on_face);
+       for (unsigned int i=0; i<dofs_per_face; ++i)
+         if (mapping[dofs_on_face[i]] == -1)
+           mapping[dofs_on_face[i]] = next_boundary_index++;
+      };
+
+  Assert (static_cast<unsigned int>(next_boundary_index) == n_boundary_dofs(),
+         ExcInternalError());
+};
+
+
+
 
 
 void DoFHandler<1>::reserve_space () {
index ebb649a0d0039f8412677c369aaf7e874971f997..42900b3918e8fa89df7bc046a22559d57ca5342e 100644 (file)
@@ -3,6 +3,7 @@
 #include <numerics/assembler.h>
 #include <numerics/base.h>
 #include <numerics/matrices.h>
+#include <numerics/vectors.h>
 #include <grid/dof_constraints.h>
 #include <grid/tria_iterator.h>
 #include <basic/data_io.h>
@@ -118,7 +119,7 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                                   // apply Dirichlet bc as described
                                   // in the docs
   map<int, double> boundary_value_list;
-  MatrixTools<dim>::interpolate_boundary_values (*dof_handler,
+  VectorTools<dim>::interpolate_boundary_values (*dof_handler,
                                                 dirichlet_bc, fe, boundary,
                                                 boundary_value_list);
   MatrixTools<dim>::apply_boundary_values (boundary_value_list,
index c2bd7354ece5036bca0a791b0a64b34a47c452be..36fc572414677cccad8aca3a7ceda390ecec9853 100644 (file)
@@ -74,6 +74,20 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 
 
 
+template <int dim>
+void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
+                                                     const FiniteElement<dim> &fe,
+                                                     const Quadrature<dim-1>  &q,
+                                                     const Boundary<dim>      &boundary,
+                                                     dSMatrix                 &matrix,
+                                                     const FunctionMap        &rhs,
+                                                     dVector                  &rhs_vector,
+                                                     vector<int>              &vec_to_dof_mapping,
+                                                     const Function<dim>      *a) {};
+
+
+
+
 
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
@@ -234,59 +248,6 @@ void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_va
 
 
 
-void
-MatrixTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
-                                            const FunctionMap &,
-                                            const FiniteElement<1> &,
-                                            const Boundary<1> &,
-                                            map<int,double>   &) {
-  Assert (false, ExcNotImplemented());
-};
-
-
-
-
-template <int dim>
-void
-MatrixTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
-                                              const FunctionMap     &dirichlet_bc,
-                                              const FiniteElement<dim> &fe,
-                                              const Boundary<dim>      &boundary,
-                                              map<int,double>   &boundary_values) {
-  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
-         ExcInvalidBoundaryIndicator());
-                                  // use two face iterators, since we need
-                                  // a DoF-iterator for the dof indices, but
-                                  // a Tria-iterator for the fe object
-  DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-                                       endf = dof.end_face();
-  
-  FunctionMap::const_iterator function_ptr;
-
-                                  // field to store the indices of dofs
-  vector<int>         face_dofs (fe.dofs_per_face);
-  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
-  vector<double>      dof_values (fe.dofs_per_face);
-       
-  for (; face!=endf; ++face)
-    if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
-       dirichlet_bc.end()) 
-                                      // face is subject to one of the
-                                      // bc listed in #dirichlet_bc#
-      {
-                                        // get indices, physical location and
-                                        // boundary values of dofs on this
-                                        // face
-       face->get_dof_indices (face_dofs);
-       fe.get_face_ansatz_points (face, boundary, dof_locations);
-       function_ptr->second->value_list (dof_locations, dof_values);
-
-                                        // enter into list
-       for (unsigned int i=0; i<face_dofs.size(); ++i)
-         boundary_values[face_dofs[i]] = dof_values[i];
-      };
-};
-
 
 
 
index bfffb5a645628d8074da8f0ac5af6baa6f3d570e..b395244d4313a9577ef17fe60b481060b381d035 100644 (file)
@@ -71,16 +71,71 @@ void VectorTools<dim>::interpolate (const DoFHandler<dim>    &dof,
 
 
 
+void VectorTools<1>::project (const DoFHandler<1>    &,
+                             const ConstraintMatrix &,
+                             const FiniteElement<1> &,
+                             const Quadrature<1>    &,
+                             const Quadrature<0>    &,
+                             const Boundary<1>      &,
+                             const Function<1>      &,
+                             const bool              ,
+                             dVector                &) {
+                                  // this function should easily be implemented
+                                  // using the template below. However some
+                                  // changes have to be made since faces don't
+                                  // exist in 1D. Maybe integrate the creation of
+                                  // zero boundary values into the
+                                  // project_boundary_values function?
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
 template <int dim>
 void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
-                                 const ConstraintMatrix   &constraints,
-                                 const FiniteElement<dim> &fe,
-                                 const Quadrature<dim>    &q,
-                                 const Boundary<dim>      &boundary,
-                                 const Function<dim>      &function,
-                                 dVector                  &vec) {
-  vec.reinit (dof.n_dofs());
+                               const ConstraintMatrix   &constraints,
+                               const FiniteElement<dim> &fe,
+                               const Quadrature<dim>    &q,
+                               const Quadrature<dim-1>  &q_boundary,
+                               const Boundary<dim>      &boundary,
+                               const Function<dim>      &function,
+                               const bool                has_zero_boundary,
+                               dVector                  &vec) {
+                                  // make up boundary values
+  map<int,double> boundary_values;
+
+  if (has_zero_boundary == false) 
+    {
+                                      // set up a list of boundary functions for
+                                      // the different boundary parts. We want the
+                                      // #function# to hold on all parts of the
+                                      // boundary
+      FunctionMap boundary_functions;
+      for (unsigned char c=0; c<255; ++c)
+       boundary_functions[c] = &function;
+      project_boundary_values (dof, boundary_functions, fe, q_boundary,
+                              boundary, boundary_values);
+    }
+  else
+                                    // no need to project boundary values
+    {
+      DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+                                           endf = dof.end_face();
+      vector<int> face_dof_indices (fe.dofs_per_face);
+      for (; face!=endf; ++face)
+       {
+         face->get_dof_indices (face_dof_indices);
+         for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+                                            // enter zero boundary values
+                                            // for all boundary nodes
+           boundary_values[face_dof_indices[i]] = 0.;
+       };
+    };
   
+      
+                                  // set up mass matrix and right hand side
+  vec.reinit (dof.n_dofs());
   dSMatrixStruct sparsity(dof.n_dofs(),
                          dof.n_dofs(),
                          dof.max_couplings_between_dofs());
@@ -93,9 +148,10 @@ void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
                                          mass_matrix, function, tmp);
 
   constraints.condense (mass_matrix);
-  constraints.condense (tmp);
+  MatrixTools<dim>::apply_boundary_values (boundary_values,
+                                          mass_matrix, vec, tmp);
 
-  int    max_iter  = 4000;
+  int    max_iter  = 1000;
   double tolerance = 1.e-16;
   Control                          control1(max_iter,tolerance);
   PrimitiveVectorMemory<dVector>   memory(tmp.size());
@@ -109,6 +165,103 @@ void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
 
 
 
+void
+VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
+                                            const FunctionMap &,
+                                            const FiniteElement<1> &,
+                                            const Boundary<1> &,
+                                            map<int,double>   &) {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
+template <int dim>
+void
+VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
+                                              const FunctionMap     &dirichlet_bc,
+                                              const FiniteElement<dim> &fe,
+                                              const Boundary<dim>      &boundary,
+                                              map<int,double>   &boundary_values) {
+  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
+         ExcInvalidBoundaryIndicator());
+                                  // use two face iterators, since we need
+                                  // a DoF-iterator for the dof indices, but
+                                  // a Tria-iterator for the fe object
+  DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+                                       endf = dof.end_face();
+  
+  FunctionMap::const_iterator function_ptr;
+
+                                  // field to store the indices of dofs
+  vector<int>         face_dofs (fe.dofs_per_face);
+  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+  vector<double>      dof_values (fe.dofs_per_face);
+       
+  for (; face!=endf; ++face)
+    if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
+       dirichlet_bc.end()) 
+                                      // face is subject to one of the
+                                      // bc listed in #dirichlet_bc#
+      {
+                                        // get indices, physical location and
+                                        // boundary values of dofs on this
+                                        // face
+       face->get_dof_indices (face_dofs);
+       fe.get_face_ansatz_points (face, boundary, dof_locations);
+       function_ptr->second->value_list (dof_locations, dof_values);
+
+                                        // enter into list
+       for (unsigned int i=0; i<face_dofs.size(); ++i)
+         boundary_values[face_dofs[i]] = dof_values[i];
+      };
+};
+
+
+
+template <int dim>
+void
+VectorTools<dim>::project_boundary_values (const DoFHandler<dim>    &dof,
+                                          const FunctionMap        &boundary_functions,
+                                          const FiniteElement<dim> &fe,
+                                          const Quadrature<dim-1>  &q,
+                                          const Boundary<dim>      &boundary,
+                                          map<int,double>   &boundary_values) {
+  vector<int>    dof_to_boundary_mapping;
+  dof.map_dof_to_boundary_indices (dof_to_boundary_mapping);
+  
+                                  // set up sparsity structure
+  dSMatrixStruct sparsity(dof.n_boundary_dofs(boundary_functions),
+                         dof.max_couplings_between_boundary_dofs());
+  dof.make_boundary_sparsity_pattern (boundary_functions, dof_to_boundary_mapping,
+                                     sparsity);
+
+                                  // note: for three or more dimensions, there
+                                  // may be constrained nodes on the boundary
+                                  // in this case the boundary mass matrix has
+                                  // to be condensed and the solution is to
+                                  // be distributed afterwards, which is not
+                                  // yet implemented
+  if (dim<3)
+    sparsity.compress();
+  else
+    Assert (false, ExcNotImplemented());
+  
+
+                                  // make mass matrix and right hand side
+  dSMatrix       mass_matrix(sparsity);
+  dVector        rhs(sparsity.n_rows());
+  
+
+  MatrixTools<dim>::create_boundary_mass_matrix (dof, fe, q, boundary,
+                                                mass_matrix, boundary_functions,
+                                                rhs, dof_to_boundary_mapping);
+  
+};
+
+
+
 
 template <int dim>
 void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.