* 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>
/*---------------------------- dof.h ---------------------------*/
#include <vector>
+#include <map>
#include <base/exceptions.h>
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;
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
*/
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
*/
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
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
*/
*/
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.
/**
* Exception
*/
+ DeclException0 (ExcInvalidBoundaryIndicator);
+ /**
+ * Exception
+ */
DeclException0 (ExcInternalError);
/**
* Exception
* 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.
* 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,
#include <base/exceptions.h>
#include <map>
-
+#include <set>
template <int dim> class Triangulation;
template <int dim> class DoFHandler;
* 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
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
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
* 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
* 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
* 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
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
*/
#include <base/exceptions.h>
+#include <map>
template <int dim> class DoFHandler;
template <int dim> class Function;
* $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
* 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.
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
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
* Exception
*/
DeclException0 (ExcInvalidFE);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidBoundaryIndicator);
};
#include <fe/fe.h>
#include <lac/dsmatrix.h>
#include <map>
+#include <set>
#include <algorithm>
+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;
-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 {
+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)
+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 () {
#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>
// 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,
+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,
-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];
- };
-};
-
+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());
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());
+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,