template <int dim> class Function;
-/**
- * Denote which norm/integral is to be computed. The following possibilities
- * are implemented:
- * \begin{itemize}
- * \item #mean#: the function or difference of functions is integrated
- * on each cell.
- * \item #L1_norm#: the absolute value of the function is integrated.
- * \item #L2_norm#: the square of the function is integrated on each
- * cell; afterwards the root is taken of this value.
- * \end{itemize}
- */
-enum NormType {
- mean,
- L1_norm,
- L2_norm,
- Linfty_norm,
- H1_seminorm,
- H1_norm
-};
-
-
/**
* 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.
- *
- * 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
- * matrix and vectors are changed accordingly. This is done by deleting all
- * entries in the matrix in the line of this degree of freedom, setting the
- * main diagonal entry to one and the right hand side element to the
- * boundary value at this node. This forces this node's value to be as specified.
- * To decouple the remaining linear system of equations and to make the system
- * symmetric again (at least if it was before), one Gauss elimination
- * step is performed with this line, by adding this (now almost empty) line to
- * all other lines which couple with the given degree of freedom and thus
- * eliminating all coupling between this degree of freedom and others. Now
- * also the column consists only of zeroes, apart from the main diagonal entry.
+ * To actually apply the boundary conditions, use is made of the
+ * #MatrixTools::apply_boundary_values# function and by interpolation of
+ * the boundary_values using the #MatrixTool::interpolate_boundary_values#
+ * function. See there for more information.
*
- * It seems as if we had to make clear not to overwrite the lines of other
- * boundary nodes when doing the Gauss elimination step. However, since we
- * reset the right hand side when passing such a node, it is not a problem
- * to change the right hand side values of other boundary nodes not yet
- * processed. It would be a problem to change those entries of nodes already
- * processed, but since the matrix entry of the present column on the row
- * of an already processed node is zero, the Gauss step does not change
- * the right hand side. We need therefore not take special care of other
- * boundary nodes.
- *
- * To make solving faster, we preset the solution vector with the right boundary
- * values. Since boundary nodes can never be hanging nodes, and since all other
- * entries of the solution vector are zero, we need not condense the solution
- * vector if the condensation process is done in-place. If done by copying
- * matrix and vectors to smaller ones, it would also be necessary to condense
- * the solution vector to preserve the preset boundary values.
- *
- * 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
- * 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
- * main diagonal entry if it is zero from a previous Gauss elimination
- * step. Of course we have to change
- * the right hand side appropriately. This is not a very good
- * strategy, but it at least should give the main diagonal entry a value
- * in the right order of dimension, which makes the solving process a bit
- * more stable. A refined algorithm would set the entry to the mean of the
- * other diagonal entries, but this seems to be too expensive.
- *
- * Because of the mentioned question, whether or not a preset solution value
- * which does not couple with other degrees of freedom remains its value or
- * not during solving iteratively, it may or may not be necessary to set
- * the correct value after solving again. This question is an open one as of
- * now and may be answered by future experience.
- *
- * At present, 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.
- * This can be done by overloading the virtual function
- * #make_boundary_value_list# which must return a list of boundary dofs
- * and their corresponding values.
- *
- * 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.
- *
- *
- * \subsection{Computing errors}
- *
- * The function #integrate_difference# performs the calculation of the error
- * between the finite element solution and a given (continuous) reference
- * function in different norms. The integration is performed using a given
- * quadrature formulae and assumes that the given finite element objects equals
- * that used for the computation of the solution.
- *
- * The result ist stored in a vector (named #difference#), where each entry
- * equals the given norm of the difference on one cell. The order of entries
- * is the same as a #cell_iterator# takes when started with #begin_active# and
- * promoted with the #++# operator.
- *
- * You can use the #distribute_cell_to_dof_vector# function of the #DoFHandler#
- * class to convert cell based data to a data vector with values on the degrees
- * of freedom, which can then be attached to a #DataOut# object to be printed.
- *
- * Presently, there is the possibility to compute the following values from the
- * difference, on each cell: #mean#, #L1_norm#, #L2_norm#, #Linfty_norm#,
- * #H1_seminorm#.
- * For the mean difference value, the reference function minus the numerical
- * solution is computed, not the other way round.
- *
- * The infinity norm of the difference on a given cell returns the maximum
- * absolute value of the difference at the quadrature points given by the
- * quadrature formula parameter. This will in some cases not be too good
- * an approximation, since for example the Gauss quadrature formulae do
- * not evaluate the difference at the end or corner points of the cells.
- * You may want to chose a quadrature formula with more quadrature points
- * or one with another distribution of the quadrature points in this case.
- * You should also take into account the superconvergence properties of finite
- * elements in some points: for example in 1D, the standard finite element
- * method is a collocation method and should return the exact value at nodal
- * points. Therefore, the trapezoidal rule should always return a vanishing
- * L-infinity error. Conversely, in 2D the maximum L-infinity error should
- * be located at the vertices or at the center of the cell, which would make
- * it plausible to use the Simpson quadrature rule. On the other hand, there
- * may be superconvergence at Gauss integration points. These examples are not
- * intended as a rule of thumb, rather they are though to illustrate that the
- * use of the wrong quadrature formula may show a significantly wrong result
- * and care should be taken to chose the right formula.
- *
- * The $H_1$ seminorm is the $L_2$ norm of the gradient of the difference. The
- * full $H_1$ norm is the sum of the seminorm and the $L_2$ norm.
- *
- * 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.
- * 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.
- *
- * To get the global mean difference, simply sum up the elements as above.
- * To get the L_\infty norm, take the maximum of the vector elements, e.g.
- * using the #dVector::linfty_norm# function.
- *
- * For the global $H_1$ norm and seminorm, the same rule applies as for the
- * $L_2$ norm: compute the $l_2$ norm of the cell error vector.
*/
template <int dim>
class ProblemBase {
const Boundary<dim> &boundary = StraightBoundary<dim>());
/**
- * Solve the system of equations.
+ * Solve the system of equations. This uses
+ * a simple CG method.
*/
virtual void solve ();
- /**
- * Integrate the difference between
- * the solution computed before and
- * the reference solution, which
- * is given as a continuous function
- * object.
- *
- * See the general documentation of this
- * class for more information.
- */
- void integrate_difference (const Function<dim> &exact_solution,
- dVector &difference,
- const Quadrature<dim> &q,
- const FiniteElement<dim> &fe,
- const NormType &norm,
- const Boundary<dim> &boundary=StraightBoundary<dim>()) const;
-
/**
* Initialize the #DataOut# object with
* the grid and DoF handler used for this
*/
virtual pair<char*,char*> get_solution_name () const;
- /**
- * Make up the list of node subject
- * to Dirichlet boundary conditions
- * and the values they are to be
- * assigned.
- *
- * See the general doc for more
- * information.
- */
- virtual void make_boundary_value_list (const FunctionMap &dirichlet_bc,
- const FiniteElement<dim> &fe,
- const Boundary<dim> &boundary,
- map<int,double> &boundary_values) const;
-
/**
* Exception
*/
/**
* Exception
*/
- DeclException0 (ExcNoMemory);
- /**
- * Exception
- */
- DeclException0 (ExcInvalidFE);
- /**
- * Exception
- */
- DeclException0 (ExcNotImplemented);
- /**
- * Exception
- */
- DeclException0 (ExcInvalidBoundaryIndicator);
- /**
- * Exception
- */
DeclException0 (ExcNoTriaSelected);
protected:
*/
ConstraintMatrix constraints;
- /**
- * Apply dirichlet boundary conditions
- * to the system matrix and vectors
- * as described in the general
- * documentation.
- */
- void apply_dirichlet_bc (dSMatrix &matrix,
- dVector &solution,
- dVector &right_hand_side,
- const FunctionMap &dirichlet_bc,
- const FiniteElement<dim> &fe,
- const Boundary<dim> &boundary);
-
friend class Assembler<dim>;
};
#include <base/exceptions.h>
+#include <map>
+
template <int dim> class Triangulation;
template <int dim> class DoFHandler;
+/**
+ * Provide a collection of functions operating on matrices. These include
+ * the application of boundary conditions to a linear system of equations
+ * and others.
+ *
+ *
+ * \subsection{Boundary conditions}
+ *
+ * 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.
+ *
+ * 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
+ * matrix and vectors are changed accordingly. This is done by deleting all
+ * entries in the matrix in the line of this degree of freedom, setting the
+ * main diagonal entry to one and the right hand side element to the
+ * boundary value at this node. This forces this node's value to be as specified.
+ * To decouple the remaining linear system of equations and to make the system
+ * symmetric again (at least if it was before), one Gauss elimination
+ * step is performed with this line, by adding this (now almost empty) line to
+ * all other lines which couple with the given degree of freedom and thus
+ * eliminating all coupling between this degree of freedom and others. Now
+ * also the column consists only of zeroes, apart from the main diagonal entry.
+ *
+ * It seems as if we had to make clear not to overwrite the lines of other
+ * boundary nodes when doing the Gauss elimination step. However, since we
+ * reset the right hand side when passing such a node, it is not a problem
+ * to change the right hand side values of other boundary nodes not yet
+ * processed. It would be a problem to change those entries of nodes already
+ * processed, but since the matrix entry of the present column on the row
+ * of an already processed node is zero, the Gauss step does not change
+ * the right hand side. We need therefore not take special care of other
+ * boundary nodes.
+ *
+ * To make solving faster, we preset the solution vector with the right boundary
+ * values. Since boundary nodes can never be hanging nodes, and since all other
+ * entries of the solution vector are zero, we need not condense the solution
+ * vector if the condensation process is done in-place. If done by copying
+ * matrix and vectors to smaller ones, it would also be necessary to condense
+ * the solution vector to preserve the preset boundary values.
+ *
+ * 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
+ * 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
+ * main diagonal entry if it is zero from a previous Gauss elimination
+ * step. Of course we have to change
+ * the right hand side appropriately. This is not a very good
+ * strategy, but it at least should give the main diagonal entry a value
+ * in the right order of dimension, which makes the solving process a bit
+ * more stable. A refined algorithm would set the entry to the mean of the
+ * other diagonal entries, but this seems to be too expensive.
+ *
+ * Because of the mentioned question, whether or not a preset solution value
+ * which does not couple with other degrees of freedom remains its value or
+ * not during solving iteratively, it may or may not be necessary to set
+ * 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
+ * as described in the general
+ * documentation.
+ */
+ static void apply_boundary_values (const map<int,double> &boundary_values,
+ dSMatrix &matrix,
+ 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
+ */
+ DeclException0 (ExcNotImplemented);
+};
+
+
+
+
/**
* Equation class to be passed to the #Assembler# if you want to make up the
* mass matrix for your problem. The mass matrix is the matrix with
+
+
+
/*---------------------------- matrices.h ---------------------------*/
/* end of #ifndef __matrices_H */
#endif
/*---------------------------- vectors.h ---------------------------*/
+
+#include <base/exceptions.h>
+
template <int dim> class DoFHandler;
template <int dim> class Function;
template <int dim> class Quadrature;
-template <int dim> class Boundary;
template <int dim> class FiniteElement;
+template <int dim> class Boundary;
+template <int dim> class StraightBoundary;
class ConstraintMatrix;
class dVector;
+/**
+ * Denote which norm/integral is to be computed. The following possibilities
+ * are implemented:
+ * \begin{itemize}
+ * \item #mean#: the function or difference of functions is integrated
+ * on each cell.
+ * \item #L1_norm#: the absolute value of the function is integrated.
+ * \item #L2_norm#: the square of the function is integrated on each
+ * cell; afterwards the root is taken of this value.
+ * \end{itemize}
+ */
+enum NormType {
+ mean,
+ L1_norm,
+ L2_norm,
+ Linfty_norm,
+ H1_seminorm,
+ H1_norm
+};
+
+
+
/**
- * Provide a class which assembles some standard vectors. Among these are
+ * Provide a class which offers some operations on vectors. Amoung these are
+ * assemblage of standard vectors, integration of the difference of a
+ * finite element solution and a continuous function,
* interpolations and projections of continuous functions to the finite
* element space and other operations.
*
* 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 Computing errors:
+ * The function #integrate_difference# performs the calculation of the error
+ * between the finite element solution and a given (continuous) reference
+ * function in different norms. The integration is performed using a given
+ * quadrature formulae and assumes that the given finite element objects equals
+ * that used for the computation of the solution.
+ *
+ * The result ist stored in a vector (named #difference#), where each entry
+ * equals the given norm of the difference on one cell. The order of entries
+ * is the same as a #cell_iterator# takes when started with #begin_active# and
+ * promoted with the #++# operator.
+ *
+ * You can use the #distribute_cell_to_dof_vector# function of the #DoFHandler#
+ * class to convert cell based data to a data vector with values on the degrees
+ * of freedom, which can then be attached to a #DataOut# object to be printed.
+ *
+ * Presently, there is the possibility to compute the following values from the
+ * difference, on each cell: #mean#, #L1_norm#, #L2_norm#, #Linfty_norm#,
+ * #H1_seminorm#.
+ * For the mean difference value, the reference function minus the numerical
+ * solution is computed, not the other way round.
+ *
+ * The infinity norm of the difference on a given cell returns the maximum
+ * absolute value of the difference at the quadrature points given by the
+ * quadrature formula parameter. This will in some cases not be too good
+ * an approximation, since for example the Gauss quadrature formulae do
+ * not evaluate the difference at the end or corner points of the cells.
+ * You may want to chose a quadrature formula with more quadrature points
+ * or one with another distribution of the quadrature points in this case.
+ * You should also take into account the superconvergence properties of finite
+ * elements in some points: for example in 1D, the standard finite element
+ * method is a collocation method and should return the exact value at nodal
+ * points. Therefore, the trapezoidal rule should always return a vanishing
+ * L-infinity error. Conversely, in 2D the maximum L-infinity error should
+ * be located at the vertices or at the center of the cell, which would make
+ * it plausible to use the Simpson quadrature rule. On the other hand, there
+ * may be superconvergence at Gauss integration points. These examples are not
+ * intended as a rule of thumb, rather they are though to illustrate that the
+ * use of the wrong quadrature formula may show a significantly wrong result
+ * a nd care should be taken to chose the right formula.
+ *
+ * The $H_1$ seminorm is the $L_2$ norm of the gradient of the difference. The
+ * full $H_1$ norm is the sum of the seminorm and the $L_2$ norm.
+ *
+ * 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.
+ * 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.
+ *
+ * To get the global mean difference, simply sum up the elements as above.
+ * To get the L_\infty norm, take the maximum of the vector elements, e.g.
+ * using the #dVector::linfty_norm# function.
+ *
+ * For the global $H_1$ norm and seminorm, the same rule applies as for the
+ * $L_2$ norm: compute the $l_2$ norm of the cell error vector.
* \end{itemize}
*/
template <int dim>
-class VectorCreator {
+class VectorTools {
public:
/**
* Compute the interpolation of
const Boundary<dim> &boundary,
const Function<dim> &function,
dVector &vec);
+
+ /**
+ * Integrate the difference between
+ * a finite element function and
+ * the reference function, which
+ * is given as a continuous function
+ * object.
+ *
+ * See the general documentation of this
+ * class for more information.
+ */
+ static void integrate_difference (const DoFHandler<dim> &dof,
+ const dVector &fe_function,
+ const Function<dim> &exact_solution,
+ dVector &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary=StraightBoundary<dim>());
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotImplemented);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidFE);
};
#include <numerics/assembler.h>
#include <numerics/base.h>
+#include <numerics/matrices.h>
#include <grid/dof_constraints.h>
#include <grid/tria_iterator.h>
#include <basic/data_io.h>
-inline double sqr (const double x) {
- return x*x;
-};
-
-
-template <int dim>
-inline double sqr_point (const Point<dim> &p) {
- return p.square();
-};
-
-
-
template <int dim>
ProblemBase<dim>::ProblemBase () :
// apply Dirichlet bc as described
// in the docs
- apply_dirichlet_bc (system_matrix, solution,
- right_hand_side,
- dirichlet_bc, fe, boundary);
+ map<int, double> boundary_value_list;
+ MatrixTools<dim>::interpolate_boundary_values (*dof_handler,
+ dirichlet_bc, fe, boundary,
+ boundary_value_list);
+ MatrixTools<dim>::apply_boundary_values (boundary_value_list,
+ system_matrix, solution,
+ right_hand_side);
};
-template <int dim>
-void ProblemBase<dim>::integrate_difference (const Function<dim> &exact_solution,
- dVector &difference,
- const Quadrature<dim> &q,
- const FiniteElement<dim> &fe,
- const NormType &norm,
- const Boundary<dim> &boundary) const {
- Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
- Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
-
- difference.reinit (tria->n_active_cells());
-
- UpdateFlags update_flags = UpdateFlags (update_q_points |
- update_jacobians |
- update_JxW_values);
- if ((norm==H1_seminorm) || (norm==H1_norm))
- update_flags = UpdateFlags (update_flags | update_gradients);
- FEValues<dim> fe_values(fe, q, update_flags);
-
- // loop over all cells
- DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active(),
- endc = dof_handler->end();
- for (unsigned int index=0; cell != endc; ++cell, ++index)
- {
- double diff=0;
- // initialize for this cell
- fe_values.reinit (cell, fe, boundary);
-
- switch (norm)
- {
- case mean:
- case L1_norm:
- case L2_norm:
- case Linfty_norm:
- case H1_norm:
- {
- // we need the finite element
- // function \psi at the different
- // integration points. Compute
- // it like this:
- // \psi(x_j)=\sum_i v_i \phi_i(x_j)
- // with v_i the nodal values of the
- // solution and \phi_i(x_j) the
- // matrix of the ansatz function
- // values at the integration point
- // x_j. Then the vector
- // of the \psi(x_j) is v*Phi with
- // v being the vector of nodal
- // values on this cell and Phi
- // the matrix.
- //
- // we then need the difference:
- // reference_function(x_j)-\psi_j
- // and assign that to the vector
- // \psi.
- const unsigned int n_q_points = q.n_quadrature_points;
- vector<double> psi (n_q_points);
-
- // in praxi: first compute
- // exact solution vector
- exact_solution.value_list (fe_values.get_quadrature_points(),
- psi);
- // then subtract finite element
- // solution
- if (true)
- {
- vector<double> function_values (n_q_points, 0);
- fe_values.get_function_values (solution, function_values);
-
- transform (psi.begin(), psi.end(),
- function_values.begin(),
- psi.begin(),
- minus<double>());
- };
-
-
- // for L1_norm and Linfty_norm:
- // take absolute
- // value, for the L2_norm take
- // square of psi
- switch (norm)
- {
- case mean:
- break;
- case L1_norm:
- case Linfty_norm:
- transform (psi.begin(), psi.end(),
- psi.begin(), ptr_fun(fabs));
- break;
- case L2_norm:
- case H1_norm:
- transform (psi.begin(), psi.end(),
- psi.begin(), ptr_fun(sqr));
- break;
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // ok, now we have the integrand,
- // let's compute the integral,
- // which is
- // sum_j psi_j JxW_j
- // (or |psi_j| or |psi_j|^2
- switch (norm)
- {
- case mean:
- case L1_norm:
- diff = inner_product (psi.begin(), psi.end(),
- fe_values.get_JxW_values().begin(),
- 0.0);
- break;
- case L2_norm:
- case H1_norm:
- diff = sqrt(inner_product (psi.begin(), psi.end(),
- fe_values.get_JxW_values().begin(),
- 0.0));
- break;
- case Linfty_norm:
- diff = *max_element (psi.begin(), psi.end());
- break;
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // note: the H1_norm uses the result
- // of the L2_norm and control goes
- // over to the next case statement!
- if (norm != H1_norm)
- break;
- };
-
- case H1_seminorm:
- {
- // note: the computation of the
- // H1_norm starts at the previous
- // case statement, but continues
- // here!
-
- // for H1_norm: re-square L2_norm.
- diff = sqr(diff);
-
- // same procedure as above, but now
- // psi is a vector of gradients
- const unsigned int n_q_points = q.n_quadrature_points;
- vector<Point<dim> > psi (n_q_points);
-
- // in praxi: first compute
- // exact solution vector
- exact_solution.gradient_list (fe_values.get_quadrature_points(),
- psi);
-
- // then subtract finite element
- // solution
- if (true)
- {
- vector<Point<dim> > function_grads (n_q_points, Point<dim>());
- fe_values.get_function_grads (solution, function_grads);
-
- transform (psi.begin(), psi.end(),
- function_grads.begin(),
- psi.begin(),
- minus<Point<dim> >());
- };
- // take square of integrand
- vector<double> psi_square (psi.size(), 0.0);
- for (unsigned int i=0; i<n_q_points; ++i)
- psi_square[i] = sqr_point(psi[i]);
-
- // add seminorm to L_2 norm or
- // to zero
- diff += inner_product (psi_square.begin(), psi_square.end(),
- fe_values.get_JxW_values().begin(),
- 0.0);
- diff = sqrt(diff);
-
- break;
- };
-
- default:
- Assert (false, ExcNotImplemented());
- };
-
-
- // append result of this cell
- // to the end of the vector
- difference(index) = diff;
- };
-};
-
-
template <int dim>
void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
-template <int dim>
-void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
- dVector &solution,
- dVector &right_hand_side,
- const FunctionMap &dirichlet_bc,
- const FiniteElement<dim> &fe,
- const Boundary<dim> &boundary) {
- Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
- Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
- ExcInvalidBoundaryIndicator());
-
- // first make up a list of dofs subject
- // to any boundary condition and which
- // value they take; if a node occurs
- // with two bc (e.g. a corner node, with
- // the lines in 2D being subject to
- // different bc's), the last value is taken
- map<int,double> boundary_values;
- make_boundary_value_list (dirichlet_bc, fe, boundary, boundary_values);
-
- map<int,double>::const_iterator dof, endd;
- const unsigned int n_dofs = matrix.m();
- const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern();
- const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
- const int *sparsity_colnums = sparsity.get_column_numbers();
-
- for (dof=boundary_values.begin(), endd=boundary_values.end(); dof != endd; ++dof)
- {
- // for each boundary dof:
-
- // set entries of this line
- // to zero
- for (unsigned int j=sparsity_rowstart[dof->first];
- j<sparsity_rowstart[dof->first+1]; ++j)
- if (sparsity_colnums[j] != dof->first)
- // if not main diagonal entry
- matrix.global_entry(j) = 0.;
-
- // set right hand side to
- // wanted value: if main diagonal
- // entry nonzero, don't touch it
- // and scale rhs accordingly. If
- // zero, take the first main
- // diagonal entry we can find, or
- // one if no nonzero main diagonal
- // element exists. Normally, however,
- // the main diagonal entry should
- // not be zero.
- //
- // store the new rhs entry to make
- // the gauss step more efficient
- double new_rhs;
- if (matrix.diag_element(dof->first) != 0.0)
- new_rhs = right_hand_side(dof->first)
- = dof->second * matrix.diag_element(dof->first);
- else
- {
- double first_diagonal_entry = 1;
- for (unsigned int i=0; i<n_dofs; ++i)
- if (matrix.diag_element(i) != 0)
- {
- first_diagonal_entry = matrix.diag_element(i);
- break;
- };
-
- matrix.set(dof->first, dof->first,
- first_diagonal_entry);
- new_rhs = right_hand_side(dof->first)
- = dof->second * first_diagonal_entry;
- };
-
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const double diagonal_entry = matrix.diag_element(dof->first);
-
- // do the Gauss step
- for (unsigned int row=0; row<n_dofs; ++row)
- for (unsigned int j=sparsity_rowstart[row];
- j<sparsity_rowstart[row+1]; ++j)
- if ((sparsity_colnums[j] == (signed int)dof->first) &&
- ((signed int)row != dof->first))
- // this line has an entry
- // in the regarding column
- // but this is not the main
- // diagonal entry
- {
- // correct right hand side
- right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
- new_rhs;
-
- // set matrix entry to zero
- matrix.global_entry(j) = 0.;
- };
-
-
- // preset solution vector
- solution(dof->first) = dof->second;
- };
-};
-
-
-
-
-
-void
-ProblemBase<1>::make_boundary_value_list (const FunctionMap &,
- const FiniteElement<1> &,
- const Boundary<1> &,
- map<int,double> &) const {
- Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
- Assert (false, ExcNotImplemented());
-};
-
-
-
-
-template <int dim>
-void
-ProblemBase<dim>::make_boundary_value_list (const FunctionMap &dirichlet_bc,
- const FiniteElement<dim> &fe,
- const Boundary<dim> &boundary,
- map<int,double> &boundary_values) const {
- Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
-
- // 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_handler->begin_active_face(),
- endf = dof_handler->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];
- };
-};
-
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/quadrature.h>
+#include <fe/fe.h>
#include <fe/fe_values.h>
#include <numerics/matrices.h>
#include <numerics/assembler.h>
+template <int dim>
+void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_values,
+ dSMatrix &matrix,
+ dVector &solution,
+ dVector &right_hand_side) {
+
+ map<int,double>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ const unsigned int n_dofs = matrix.m();
+ const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern();
+ const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
+ const int *sparsity_colnums = sparsity.get_column_numbers();
+
+ for (; dof != endd; ++dof)
+ {
+ // for each boundary dof:
+
+ // set entries of this line
+ // to zero
+ for (unsigned int j=sparsity_rowstart[dof->first];
+ j<sparsity_rowstart[dof->first+1]; ++j)
+ if (sparsity_colnums[j] != dof->first)
+ // if not main diagonal entry
+ matrix.global_entry(j) = 0.;
+
+ // set right hand side to
+ // wanted value: if main diagonal
+ // entry nonzero, don't touch it
+ // and scale rhs accordingly. If
+ // zero, take the first main
+ // diagonal entry we can find, or
+ // one if no nonzero main diagonal
+ // element exists. Normally, however,
+ // the main diagonal entry should
+ // not be zero.
+ //
+ // store the new rhs entry to make
+ // the gauss step more efficient
+ double new_rhs;
+ if (matrix.diag_element(dof->first) != 0.0)
+ new_rhs = right_hand_side(dof->first)
+ = dof->second * matrix.diag_element(dof->first);
+ else
+ {
+ double first_diagonal_entry = 1;
+ for (unsigned int i=0; i<n_dofs; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ first_diagonal_entry = matrix.diag_element(i);
+ break;
+ };
+
+ matrix.set(dof->first, dof->first,
+ first_diagonal_entry);
+ new_rhs = right_hand_side(dof->first)
+ = dof->second * first_diagonal_entry;
+ };
+
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const double diagonal_entry = matrix.diag_element(dof->first);
+
+ // do the Gauss step
+ for (unsigned int row=0; row<n_dofs; ++row)
+ for (unsigned int j=sparsity_rowstart[row];
+ j<sparsity_rowstart[row+1]; ++j)
+ if ((sparsity_colnums[j] == (signed int)dof->first) &&
+ ((signed int)row != dof->first))
+ // this line has an entry
+ // in the regarding column
+ // but this is not the main
+ // diagonal entry
+ {
+ // correct right hand side
+ right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
+ new_rhs;
+
+ // set matrix entry to zero
+ matrix.global_entry(j) = 0.;
+ };
+
+
+ // preset solution vector
+ solution(dof->first) = dof->second;
+ };
+};
+
+
+
+
+
+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];
+ };
+};
+
+
+
+
+
+
template <int dim>
MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
const Function<dim> * const a) :
template class MatrixCreator<1>;
template class MatrixCreator<2>;
+template class MatrixTools<1>;
+template class MatrixTools<2>;
template class MassMatrix<1>;
template class MassMatrix<2>;
template class LaplaceMatrix<1>;
#include <grid/tria_iterator.h>
#include <grid/dof_constraints.h>
#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/quadrature.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <lac/dvector.h>
#include "../../../mia/vectormemory.h"
#include "../../../mia/cg.h"
+#include <numeric>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+inline double sqr (const double x) {
+ return x*x;
+};
+
+
+template <int dim>
+inline double sqr_point (const Point<dim> &p) {
+ return p.square();
+};
+
+
+
+
+
template <int dim>
-void VectorCreator<dim>::interpolate (const DoFHandler<dim> &dof,
+void VectorTools<dim>::interpolate (const DoFHandler<dim> &dof,
const FiniteElement<dim> &fe,
const Boundary<dim> &boundary,
const Function<dim> &function,
template <int dim>
-void VectorCreator<dim>::project (const DoFHandler<dim> &dof,
+void VectorTools<dim>::project (const DoFHandler<dim> &dof,
const ConstraintMatrix &constraints,
const FiniteElement<dim> &fe,
const Quadrature<dim> &q,
-template VectorCreator<1>;
-template VectorCreator<2>;
+template <int dim>
+void VectorTools<dim>::integrate_difference (const DoFHandler<dim> &dof,
+ const dVector &fe_function,
+ const Function<dim> &exact_solution,
+ dVector &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary) {
+ Assert (fe == dof.get_selected_fe(), ExcInvalidFE());
+
+ difference.reinit (dof.get_tria().n_active_cells());
+
+ UpdateFlags update_flags = UpdateFlags (update_q_points |
+ update_jacobians |
+ update_JxW_values);
+ if ((norm==H1_seminorm) || (norm==H1_norm))
+ update_flags = UpdateFlags (update_flags | update_gradients);
+ FEValues<dim> fe_values(fe, q, update_flags);
+
+ // loop over all cells
+ DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+ endc = dof.end();
+ for (unsigned int index=0; cell != endc; ++cell, ++index)
+ {
+ double diff=0;
+ // initialize for this cell
+ fe_values.reinit (cell, fe, boundary);
+
+ switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ case L2_norm:
+ case Linfty_norm:
+ case H1_norm:
+ {
+ // we need the finite element
+ // function \psi at the different
+ // integration points. Compute
+ // it like this:
+ // \psi(x_j)=\sum_i v_i \phi_i(x_j)
+ // with v_i the nodal values of the
+ // fe_function and \phi_i(x_j) the
+ // matrix of the ansatz function
+ // values at the integration point
+ // x_j. Then the vector
+ // of the \psi(x_j) is v*Phi with
+ // v being the vector of nodal
+ // values on this cell and Phi
+ // the matrix.
+ //
+ // we then need the difference:
+ // reference_function(x_j)-\psi_j
+ // and assign that to the vector
+ // \psi.
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<double> psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+ exact_solution.value_list (fe_values.get_quadrature_points(),
+ psi);
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector<double> function_values (n_q_points, 0);
+ fe_values.get_function_values (fe_function, function_values);
+
+ transform (psi.begin(), psi.end(),
+ function_values.begin(),
+ psi.begin(),
+ minus<double>());
+ };
+
+ // for L1_norm and Linfty_norm:
+ // take absolute
+ // value, for the L2_norm take
+ // square of psi
+ switch (norm)
+ {
+ case mean:
+ break;
+ case L1_norm:
+ case Linfty_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(fabs));
+ break;
+ case L2_norm:
+ case H1_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(sqr));
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // ok, now we have the integrand,
+ // let's compute the integral,
+ // which is
+ // sum_j psi_j JxW_j
+ // (or |psi_j| or |psi_j|^2
+ switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ diff = inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ break;
+ case L2_norm:
+ case H1_norm:
+ diff = sqrt(inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0));
+ break;
+ case Linfty_norm:
+ diff = *max_element (psi.begin(), psi.end());
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // note: the H1_norm uses the result
+ // of the L2_norm and control goes
+ // over to the next case statement!
+ if (norm != H1_norm)
+ break;
+ };
+
+ case H1_seminorm:
+ {
+ // note: the computation of the
+ // H1_norm starts at the previous
+ // case statement, but continues
+ // here!
+
+ // for H1_norm: re-square L2_norm.
+ diff = sqr(diff);
+
+ // same procedure as above, but now
+ // psi is a vector of gradients
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<Point<dim> > psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+ exact_solution.gradient_list (fe_values.get_quadrature_points(),
+ psi);
+
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+ fe_values.get_function_grads (fe_function, function_grads);
+
+ transform (psi.begin(), psi.end(),
+ function_grads.begin(),
+ psi.begin(),
+ minus<Point<dim> >());
+ };
+ // take square of integrand
+ vector<double> psi_square (psi.size(), 0.0);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ psi_square[i] = sqr_point(psi[i]);
+
+ // add seminorm to L_2 norm or
+ // to zero
+ diff += inner_product (psi_square.begin(), psi_square.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ diff = sqrt(diff);
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+
+ // append result of this cell
+ // to the end of the vector
+ difference(index) = diff;
+ };
+};
+
+
+template VectorTools<1>;
+template VectorTools<2>;