From 0c2e60a6634ad8423711666ca0f5c80c9df13d19 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 19 May 1998 13:51:14 +0000 Subject: [PATCH] Move some functions from ProblemBase to more useful locations. git-svn-id: https://svn.dealii.org/trunk@315 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/base.h | 217 +----------- deal.II/deal.II/include/numerics/matrices.h | 160 +++++++++ deal.II/deal.II/include/numerics/vectors.h | 121 ++++++- deal.II/deal.II/source/numerics/base.cc | 370 +------------------- deal.II/deal.II/source/numerics/matrices.cc | 153 ++++++++ deal.II/deal.II/source/numerics/vectors.cc | 221 +++++++++++- 6 files changed, 662 insertions(+), 580 deletions(-) diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index f914af9865..82e76641d0 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -25,27 +25,6 @@ template class StraightBoundary; template 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 -}; - - /** @@ -102,136 +81,11 @@ enum NormType { * 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 class ProblemBase { @@ -311,27 +165,11 @@ class ProblemBase { const Boundary &boundary = StraightBoundary()); /** - * 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 &exact_solution, - dVector &difference, - const Quadrature &q, - const FiniteElement &fe, - const NormType &norm, - const Boundary &boundary=StraightBoundary()) const; - /** * Initialize the #DataOut# object with * the grid and DoF handler used for this @@ -362,20 +200,6 @@ class ProblemBase { */ virtual pair 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 &fe, - const Boundary &boundary, - map &boundary_values) const; - /** * Exception */ @@ -383,22 +207,6 @@ class ProblemBase { /** * Exception */ - DeclException0 (ExcNoMemory); - /** - * Exception - */ - DeclException0 (ExcInvalidFE); - /** - * Exception - */ - DeclException0 (ExcNotImplemented); - /** - * Exception - */ - DeclException0 (ExcInvalidBoundaryIndicator); - /** - * Exception - */ DeclException0 (ExcNoTriaSelected); protected: @@ -439,19 +247,6 @@ class ProblemBase { */ 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 &fe, - const Boundary &boundary); - friend class Assembler; }; diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index de9b94d9db..501f29cd63 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -7,6 +7,8 @@ #include +#include + template class Triangulation; template class DoFHandler; @@ -171,6 +173,161 @@ class MatrixCreator { +/** + * 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 +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 + * as described in the general + * documentation. + */ + static void apply_boundary_values (const map &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 &dof, + const FunctionMap &dirichlet_bc, + const FiniteElement &fe, + const Boundary &boundary, + map &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 @@ -377,6 +534,9 @@ class LaplaceMatrix : public Equation { + + + /*---------------------------- matrices.h ---------------------------*/ /* end of #ifndef __matrices_H */ #endif diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index c821ddd9d2..12641af42b 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -5,18 +5,46 @@ /*---------------------------- vectors.h ---------------------------*/ + +#include + template class DoFHandler; template class Function; template class Quadrature; -template class Boundary; template class FiniteElement; +template class Boundary; +template 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. * @@ -57,10 +85,69 @@ class dVector; * 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 -class VectorCreator { +class VectorTools { public: /** * Compute the interpolation of @@ -90,6 +177,34 @@ class VectorCreator { const Boundary &boundary, const Function &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 &dof, + const dVector &fe_function, + const Function &exact_solution, + dVector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &boundary=StraightBoundary()); + + /** + * Exception + */ + DeclException0 (ExcNotImplemented); + /** + * Exception + */ + DeclException0 (ExcInvalidFE); }; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 1a0618c8e8..ebb649a0d0 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -20,18 +21,6 @@ -inline double sqr (const double x) { - return x*x; -}; - - -template -inline double sqr_point (const Point &p) { - return p.square(); -}; - - - template ProblemBase::ProblemBase () : @@ -128,9 +117,13 @@ void ProblemBase::assemble (const Equation &equation, // apply Dirichlet bc as described // in the docs - apply_dirichlet_bc (system_matrix, solution, - right_hand_side, - dirichlet_bc, fe, boundary); + map boundary_value_list; + MatrixTools::interpolate_boundary_values (*dof_handler, + dirichlet_bc, fe, boundary, + boundary_value_list); + MatrixTools::apply_boundary_values (boundary_value_list, + system_matrix, solution, + right_hand_side); }; @@ -156,196 +149,6 @@ void ProblemBase::solve () { -template -void ProblemBase::integrate_difference (const Function &exact_solution, - dVector &difference, - const Quadrature &q, - const FiniteElement &fe, - const NormType &norm, - const Boundary &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 fe_values(fe, q, update_flags); - - // loop over all cells - DoFHandler::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 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 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()); - }; - - - // 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 > 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 > function_grads (n_q_points, Point()); - fe_values.get_function_grads (solution, function_grads); - - transform (psi.begin(), psi.end(), - function_grads.begin(), - psi.begin(), - minus >()); - }; - // take square of integrand - vector psi_square (psi.size(), 0.0); - for (unsigned int i=0; i void ProblemBase::fill_data (DataOut &out) const { @@ -369,163 +172,6 @@ pair ProblemBase::get_solution_name () const { -template -void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, - dVector &solution, - dVector &right_hand_side, - const FunctionMap &dirichlet_bc, - const FiniteElement &fe, - const Boundary &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 boundary_values; - make_boundary_value_list (dirichlet_bc, fe, boundary, boundary_values); - - map::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]; - jfirst+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; ifirst, 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; rowfirst) && - ((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 &) const { - Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); - Assert (false, ExcNotImplemented()); -}; - - - - -template -void -ProblemBase::make_boundary_value_list (const FunctionMap &dirichlet_bc, - const FiniteElement &fe, - const Boundary &boundary, - map &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::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 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 #include #include +#include #include #include #include @@ -145,6 +146,156 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, +template +void MatrixTools::apply_boundary_values (const map &boundary_values, + dSMatrix &matrix, + dVector &solution, + dVector &right_hand_side) { + + map::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]; + jfirst+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; ifirst, 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; rowfirst) && + ((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 &) { + 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 MassMatrix::MassMatrix (const Function * const rhs, const Function * const a) : @@ -367,6 +518,8 @@ void LaplaceMatrix::assemble (dVector &rhs, 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>; diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 061aa87c06..a9136e875a 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include #include #include @@ -16,9 +18,30 @@ #include "../../../mia/vectormemory.h" #include "../../../mia/cg.h" +#include +#include +#include + + + + +inline double sqr (const double x) { + return x*x; +}; + + +template +inline double sqr_point (const Point &p) { + return p.square(); +}; + + + + + template -void VectorCreator::interpolate (const DoFHandler &dof, +void VectorTools::interpolate (const DoFHandler &dof, const FiniteElement &fe, const Boundary &boundary, const Function &function, @@ -49,7 +72,7 @@ void VectorCreator::interpolate (const DoFHandler &dof, template -void VectorCreator::project (const DoFHandler &dof, +void VectorTools::project (const DoFHandler &dof, const ConstraintMatrix &constraints, const FiniteElement &fe, const Quadrature &q, @@ -87,5 +110,195 @@ void VectorCreator::project (const DoFHandler &dof, -template VectorCreator<1>; -template VectorCreator<2>; +template +void VectorTools::integrate_difference (const DoFHandler &dof, + const dVector &fe_function, + const Function &exact_solution, + dVector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &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 fe_values(fe, q, update_flags); + + // loop over all cells + DoFHandler::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 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 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()); + }; + + // 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 > 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 > function_grads (n_q_points, Point()); + fe_values.get_function_grads (fe_function, function_grads); + + transform (psi.begin(), psi.end(), + function_grads.begin(), + psi.begin(), + minus >()); + }; + // take square of integrand + vector psi_square (psi.size(), 0.0); + for (unsigned int i=0; i; +template VectorTools<2>; -- 2.39.5