From: Wolfgang Bangerth Date: Wed, 20 May 1998 17:09:45 +0000 (+0000) Subject: First part of huge changes to allow projection of the trace of a function to the... X-Git-Tag: v8.0.0~22930 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a49bad903ff1171e517dedc98f9916493497128c;p=dealii.git First part of huge changes to allow projection of the trace of a function to the boundary or parts thereof, to allow the creation of mass matrices on the boundary, use boundary projection in the projection of functions to the ansatz space, etc. git-svn-id: https://svn.dealii.org/trunk@324 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/function.h b/deal.II/base/include/base/function.h index 8c530205eb..4cdfe88e4d 100644 --- a/deal.II/base/include/base/function.h +++ b/deal.II/base/include/base/function.h @@ -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 diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 7d3cd74743..bac5a62d56 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -5,6 +5,7 @@ /*---------------------------- dof.h ---------------------------*/ #include +#include #include @@ -25,8 +26,8 @@ template class TriaIterator; template class TriaActiveIterator; template class Triangulation; - template class FiniteElementBase; +template class Function; class dVector; class dSMatrix; @@ -447,6 +448,17 @@ class DoFHandler : public DoFDimensionInfo { typedef typename DoFDimensionInfo::face_iterator face_iterator; typedef typename DoFDimensionInfo::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*> FunctionMap; + /** * Constructor. Take #tria# as the @@ -531,6 +543,57 @@ class DoFHandler : public DoFDimensionInfo { */ 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 &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 &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 { */ 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 { 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 &mapping) const; + + /** + * DOC FIXME: algorithm + */ + void map_dof_to_boundary_indices (const FunctionMap &boundary_indicators, + vector &mapping) const; + /** * @name Cell iterator functions */ @@ -1001,6 +1108,23 @@ class DoFHandler : public DoFDimensionInfo { */ 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 { /** * Exception */ + DeclException0 (ExcInvalidBoundaryIndicator); + /** + * Exception + */ DeclException0 (ExcInternalError); /** * Exception diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 82e76641d0..42c5fea0a5 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -102,7 +102,7 @@ class ProblemBase { * See the general documentation of this * class for more detail. */ - typedef map*> FunctionMap; + typedef map*> FunctionMap; /** * Typdedef an iterator which assembles * matrices and vectors. diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 88d929eff9..4ab4f49b45 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -153,7 +153,7 @@ class KellyErrorEstimator { * for each boundary indicator, which is * guaranteed by the #map# data type. */ - typedef map*> FunctionMap; + typedef map*> FunctionMap; void estimate_error (const DoFHandler &dof, const Quadrature &quadrature, diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 501f29cd63..c61dccc150 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -8,7 +8,7 @@ #include #include - +#include template class Triangulation; template 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 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*> FunctionMap; + /** * Assemble the mass matrix. If no * coefficient is given, it is assumed @@ -135,6 +174,31 @@ class MatrixCreator { dVector &rhs_vector, const Function *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 &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const FunctionMap &rhs, + dVector &rhs_vector, + vector &vec_to_dof_mapping, + const Function *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 class MatrixTools : public 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*> FunctionMap; - /** * Apply dirichlet boundary conditions * to the system matrix and vectors @@ -298,27 +319,6 @@ class MatrixTools : public MatrixCreator { 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 &dof, - const FunctionMap &dirichlet_bc, - const FiniteElement &fe, - const Boundary &boundary, - map &boundary_values); - - - /** - * Exception - */ - DeclException0 (ExcInvalidBoundaryIndicator); /** * Exception */ diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 12641af42b..c5a7ff36db 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -7,6 +7,7 @@ #include +#include template class DoFHandler; template 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 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*> FunctionMap; + /** * Compute the interpolation of * #function# at the ansatz points to @@ -174,10 +248,43 @@ class VectorTools { const ConstraintMatrix &constraints, const FiniteElement &fe, const Quadrature &q, + const Quadrature &q_boundary, const Boundary &boundary, const Function &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 &dof, + const FunctionMap &dirichlet_bc, + const FiniteElement &fe, + const Boundary &boundary, + map &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 &dof, + const FunctionMap &boundary_functions, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + map &boundary_values); + /** * Integrate the difference between * a finite element function and @@ -205,6 +312,10 @@ class VectorTools { * Exception */ DeclException0 (ExcInvalidFE); + /** + * Exception + */ + DeclException0 (ExcInvalidBoundaryIndicator); }; diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index 68326925b7..fbb3c795dc 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -10,6 +10,7 @@ #include #include #include +#include #include @@ -625,6 +626,68 @@ unsigned int DoFHandler::n_dofs () const { +unsigned int DoFHandler<1>::n_boundary_dofs () const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template +unsigned int DoFHandler::n_boundary_dofs () const { + Assert (selected_fe != 0, ExcNoFESelected()); + + set boundary_dofs; + + const unsigned int dofs_per_face = selected_fe->dofs_per_face; + vector 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::n_boundary_dofs (const FunctionMap &) const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (false, ExcNotImplemented()); + return 0; +}; + + +template +unsigned int DoFHandler::n_boundary_dofs (const FunctionMap &boundary_indicators) const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (boundary_indicators.find(255) == boundary_indicators.end(), + ExcInvalidBoundaryIndicator()); + + set boundary_dofs; + + const unsigned int dofs_per_face = selected_fe->dofs_per_face; + vector 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 const Triangulation & DoFHandler::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 +void DoFHandler::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 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; ddofs_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; ddofs_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; itotal_dofs; ++i) - for (unsigned int j=0; jtotal_dofs; ++j) + for (unsigned int i=0; i::make_boundary_sparsity_pattern (const vector &, + 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 +void DoFHandler::make_boundary_sparsity_pattern (const vector &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 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::vertices_per_cell; ++vertex) - for (unsigned int d=0; ddofs_per_vertex; ++d) - dofs_on_this_cell[dof_number++] = cell->vertex_dof_index (vertex,d); - for (unsigned int line=0; line::faces_per_cell; ++line) - for (unsigned int d=0; ddofs_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; ddofs_per_quad; ++d) - dofs_on_this_cell[dof_number++] = cell->dof_index (d); +void DoFHandler<1>::make_boundary_sparsity_pattern (const FunctionMap &, + const vector &, + dSMatrixStruct &) const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (false, ExcInternalError()); +}; - // make sparsity pattern for this cell - for (unsigned int i=0; itotal_dofs; ++i) - for (unsigned int j=0; jtotal_dofs; ++j) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - }; - - delete[] dofs_on_this_cell; + + +template +void DoFHandler::make_boundary_sparsity_pattern (const FunctionMap &boundary_indicators, + const vector &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 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 void DoFHandler::make_transfer_matrix (const DoFHandler &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::distribute_cell_to_dof_vector (const dVector &cell_data, +void DoFHandler<1>::map_dof_to_boundary_indices (vector &) const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (false, ExcNotImplemented()); +}; + + + +template +void DoFHandler::map_dof_to_boundary_indices (vector &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 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(next_boundary_index) == n_boundary_dofs(), + ExcInternalError()); +}; + + + +void DoFHandler<1>::map_dof_to_boundary_indices (const FunctionMap &, + vector &) const { + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (false, ExcNotImplemented()); +}; + + + +template +void DoFHandler::map_dof_to_boundary_indices (const FunctionMap &boundary_indicators, + vector &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 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(next_boundary_index) == n_boundary_dofs(), + ExcInternalError()); +}; + + + void DoFHandler<1>::reserve_space () { diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index ebb649a0d0..42900b3918 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -118,7 +119,7 @@ void ProblemBase::assemble (const Equation &equation, // apply Dirichlet bc as described // in the docs map boundary_value_list; - MatrixTools::interpolate_boundary_values (*dof_handler, + VectorTools::interpolate_boundary_values (*dof_handler, dirichlet_bc, fe, boundary, boundary_value_list); MatrixTools::apply_boundary_values (boundary_value_list, diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index c2bd7354ec..36fc572414 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -74,6 +74,20 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, +template +void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + dSMatrix &matrix, + const FunctionMap &rhs, + dVector &rhs_vector, + vector &vec_to_dof_mapping, + const Function *a) {}; + + + + template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, @@ -234,59 +248,6 @@ void MatrixTools::apply_boundary_values (const map &boundary_va -void -MatrixTools<1>::interpolate_boundary_values (const DoFHandler<1> &, - const FunctionMap &, - const FiniteElement<1> &, - const Boundary<1> &, - map &) { - Assert (false, ExcNotImplemented()); -}; - - - - -template -void -MatrixTools::interpolate_boundary_values (const DoFHandler &dof, - const FunctionMap &dirichlet_bc, - const FiniteElement &fe, - const Boundary &boundary, - map &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::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 face_dofs (fe.dofs_per_face); - vector > dof_locations (face_dofs.size(), Point()); - vector 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::interpolate (const DoFHandler &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 void VectorTools::project (const DoFHandler &dof, - const ConstraintMatrix &constraints, - const FiniteElement &fe, - const Quadrature &q, - const Boundary &boundary, - const Function &function, - dVector &vec) { - vec.reinit (dof.n_dofs()); + const ConstraintMatrix &constraints, + const FiniteElement &fe, + const Quadrature &q, + const Quadrature &q_boundary, + const Boundary &boundary, + const Function &function, + const bool has_zero_boundary, + dVector &vec) { + // make up boundary values + map 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::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); + vector face_dof_indices (fe.dofs_per_face); + for (; face!=endf; ++face) + { + face->get_dof_indices (face_dof_indices); + for (unsigned int i=0; i::project (const DoFHandler &dof, mass_matrix, function, tmp); constraints.condense (mass_matrix); - constraints.condense (tmp); + MatrixTools::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 memory(tmp.size()); @@ -109,6 +165,103 @@ void VectorTools::project (const DoFHandler &dof, +void +VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &, + const FunctionMap &, + const FiniteElement<1> &, + const Boundary<1> &, + map &) { + Assert (false, ExcNotImplemented()); +}; + + + + +template +void +VectorTools::interpolate_boundary_values (const DoFHandler &dof, + const FunctionMap &dirichlet_bc, + const FiniteElement &fe, + const Boundary &boundary, + map &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::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 face_dofs (fe.dofs_per_face); + vector > dof_locations (face_dofs.size(), Point()); + vector 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 +void +VectorTools::project_boundary_values (const DoFHandler &dof, + const FunctionMap &boundary_functions, + const FiniteElement &fe, + const Quadrature &q, + const Boundary &boundary, + map &boundary_values) { + vector 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::create_boundary_mass_matrix (dof, fe, q, boundary, + mass_matrix, boundary_functions, + rhs, dof_to_boundary_mapping); + +}; + + + template void VectorTools::integrate_difference (const DoFHandler &dof,