From: Wolfgang Bangerth Date: Wed, 25 Mar 1998 17:48:41 +0000 (+0000) Subject: Support for error integration and initial support for dirichlet boundary conditions. X-Git-Tag: v8.0.0~23156 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9748fed859bcdf17c121bdb593fced8ed87c32e9;p=dealii.git Support for error integration and initial support for dirichlet boundary conditions. git-svn-id: https://svn.dealii.org/trunk@95 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/function.h b/deal.II/base/include/base/function.h index a00e86d2a6..aac5e22336 100644 --- a/deal.II/base/include/base/function.h +++ b/deal.II/base/include/base/function.h @@ -7,7 +7,7 @@ //forward declaration template class Point; -template class vector; +template class vector; diff --git a/deal.II/deal.II/Makefile b/deal.II/deal.II/Makefile index ebf512092e..1de1718ad3 100644 --- a/deal.II/deal.II/Makefile +++ b/deal.II/deal.II/Makefile @@ -3,7 +3,7 @@ all: deal.II doc examples -clean: lib-clean doc-clean examples-clean +clean: lib-clean examples-clean doc-clean deal.II: diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index b0d40ed8f4..35ee247952 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -36,15 +36,21 @@ Let ParameterHandler and DataIn/Out throw exceptions. Make more tests Replace function objects by mem_fun, fun_ptr, ... when member function templates are supported. +Invent a way to tell FEValues::reinit which fields to reinit. If you + only want to integrate a FE solution over a cell, transforming the + gradients may be a rather expensive operation! + Check in FEValues::shape_grad whether the gardients have been + reinit'd last time to avoid people taking gradients which have not + been computed. + +implement ProblemBase::distribute_cell_to_dof_vector -DEAL: -make dSMatrixStruct warn if max # of entries per row is overrun! -let dSMatrixStruct have an empty constructor to allow for member - objects being initialized later + +DEAL: why do loops in dFMatrix run backwards? Let an empty dfmatrix assume it has dimensions zero rather than one! diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index aa74ccb7a4..63df424b39 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -7,6 +7,8 @@ template class Triangulation; template class DoFHandler; +class dVector; + #include #include @@ -239,7 +241,8 @@ class DoFQuadAccessor : public BaseClass, public DoFAccessor { * Return a pointer to the #i#th line * bounding this #Quad#. */ - TriaIterator > > line (const unsigned int i) const; + TriaIterator > > + line (const unsigned int i) const; /** * Return the #i#th child as a DoF quad @@ -299,10 +302,17 @@ class DoFSubstructAccessor<1> : public DoFLineAccessor<1,CellAccessor<1> > { * Constructor */ DoFSubstructAccessor (Triangulation<1> *tria, - const int level, - const int index, - const void *local_data) : + const int level, + const int index, + const void *local_data) : DoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a substruct_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef void * substruct_iterator; }; @@ -317,10 +327,17 @@ class DoFSubstructAccessor<2> : public DoFQuadAccessor<2,CellAccessor<2> > { * Constructor */ DoFSubstructAccessor (Triangulation<2> *tria, - const int level, - const int index, - const void *local_data) : + const int level, + const int index, + const void *local_data) : DoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a substruct_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > substruct_iterator; }; @@ -369,6 +386,16 @@ class DoFCellAccessor : public DoFSubstructAccessor { */ TriaIterator > child (const unsigned int) const; + /** + * Return an iterator to the #i#th face + * of this cell. + * + * This function is not implemented in 1D, + * and maps to DoFQuadAccessor::line in 2D. + */ + typename DoFSubstructAccessor::substruct_iterator + face (const unsigned int i) const; + /** * Return the indices of the dofs of this * cell in the standard ordering: dofs @@ -377,10 +404,38 @@ class DoFCellAccessor : public DoFSubstructAccessor { * dofs on quad 0, etc. * * The dof indices are appended to the end - * of the vector. It is not assumed that + * of the vector. It is assumed that * the vector be empty beforehand. */ - void dof_indices (vector &dof_indices) const; + void get_dof_indices (vector &dof_indices) const; + + /** + * Return the value of the given vector + * restricted to the dofs of this + * cell in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * The values are appended to the end + * of the vector. It is assumed that + * the vector be empty beforehand. + */ + void get_dof_values (const dVector &values, + vector &dof_values) const; + + /** + * Exception + */ + DeclException0 (ExcVectorNotEmpty); + /** + * Exception + */ + DeclException0 (ExcVectorDoesNotMatch); + /** + * Exception + */ + DeclException0 (ExcNotUsefulForThisDimension); }; diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index f660cb4ad3..360b6f434c 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -152,6 +152,10 @@ class DoFDimensionInfo<1> { typedef raw_line_iterator raw_cell_iterator; typedef line_iterator cell_iterator; typedef active_line_iterator active_cell_iterator; + + typedef void * raw_substruct_iterator; + typedef void * substruct_iterator; + typedef void * active_substruct_iterator; }; @@ -177,6 +181,10 @@ class DoFDimensionInfo<2> { typedef raw_quad_iterator raw_cell_iterator; typedef quad_iterator cell_iterator; typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_substruct_iterator; + typedef line_iterator substruct_iterator; + typedef active_line_iterator active_substruct_iterator; }; @@ -428,7 +436,11 @@ class DoFHandler : public DoFDimensionInfo { typedef typename DoFDimensionInfo::cell_iterator cell_iterator; typedef typename DoFDimensionInfo::active_cell_iterator active_cell_iterator; + typedef typename DoFDimensionInfo::raw_substruct_iterator raw_substruct_iterator; + typedef typename DoFDimensionInfo::substruct_iterator substruct_iterator; + typedef typename DoFDimensionInfo::active_substruct_iterator active_substruct_iterator; + /** * Constructor. Take #tria# as the * triangulation to work on. diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 7757b1148f..f74d3e8414 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -72,6 +72,15 @@ class FEValues { double shape_value (const unsigned int i, const unsigned int j) const; + /** + * Return a pointer to the matrix holding + * all values of shape functions at all + * integration points, on the present cell. + * For the format of this matrix, see the + * documentation for the matrix itself. + */ + const dFMatrix & get_shape_values () const; + /** * Return the gradient of the #i#th shape * function at the #j# quadrature point. @@ -88,18 +97,40 @@ class FEValues { const Point & shape_grad (const unsigned int i, const unsigned int j) const; + /** + * Return a pointer to the matrix holding + * all gradients of shape functions at all + * integration points, on the present cell. + * For the format of this matrix, see the + * documentation for the matrix itself. + */ + const vector > > & get_shape_grads () const; + /** * Return the position of the #i#th * quadrature point in real space. */ const Point & quadrature_point (const unsigned int i) const; + /** + * Return a pointer to the vector of + * quadrature points. + */ + const vector > & get_quadrature_points () const; + /** * Return the Jacobi determinant times * the weight of the #i#th quadrature * point. */ double JxW (const unsigned int i) const; + + /** + * Return a pointer to the array holding + * the JxW values at the different + * quadrature points. + */ + const vector & get_JxW_values () const; /** * Reinitialize the gradients, Jacobi @@ -641,6 +672,48 @@ class FiniteElement<2> : public FiniteElementBase<2> { + + +/*------------------------ Inline functions -----------------------------------*/ + + + +template +inline +const dFMatrix & FEValues::get_shape_values () const { + return shape_values; +}; + + + + +template +inline +const vector > > & +FEValues::get_shape_grads () const { + return shape_gradients; +}; + + + +template +inline +const vector > & +FEValues::get_quadrature_points () const { + return quadrature_points; +}; + + + +template +inline +const vector & +FEValues::get_JxW_values () const { + return JxW_values; +}; + + + /*---------------------------- fe.h ---------------------------*/ /* end of #ifndef __fe_H */ diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index a7fd173119..2685ec001e 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -462,6 +462,10 @@ class TriaDimensionInfo; #Triangulation<1>#, the #quad_iterator#s are declared as #void *#. Thus these types exist, but are useless and will certainly make any involuntary use visible. + + The same applies for the #substruct_iterator# types, since lines + have no substructures apart from vertices, which are handled in + a different way, however. */ class TriaDimensionInfo<1> { public: @@ -476,6 +480,10 @@ class TriaDimensionInfo<1> { typedef raw_line_iterator raw_cell_iterator; typedef line_iterator cell_iterator; typedef active_line_iterator active_cell_iterator; + + typedef void * raw_substruct_iterator; + typedef void * substruct_iterator; + typedef void * active_substruct_iterator; }; @@ -499,6 +507,10 @@ class TriaDimensionInfo<1> { typedef raw_quad_iterator raw_cell_iterator; typedef quad_iterator cell_iterator; typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_substruct_iterator; + typedef line_iterator substruct_iterator; + typedef active_line_iterator active_substruct_iterator; \end{verbatim} */ class TriaDimensionInfo<2> { @@ -514,6 +526,10 @@ class TriaDimensionInfo<2> { typedef raw_quad_iterator raw_cell_iterator; typedef quad_iterator cell_iterator; typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_substruct_iterator; + typedef line_iterator substruct_iterator; + typedef active_line_iterator active_substruct_iterator; }; @@ -584,10 +600,17 @@ class TriaDimensionInfo<2> { typedef quad_line_iterator raw_cell_iterator; typedef quad_iterator cell_iterator; typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_substruct_iterator; + typedef line_iterator substruct_iterator; + typedef active_line_iterator active_substruct_iterator; \end{verbatim} By using the cell iterators, you can write code nearly independent of - the spatial dimension. + the spatial dimension. The same applies for substructure iterators, + where a substructure is defined as a face of a cell. The face of a + cell is be a vertex in 1D and a line in 2D; however, vertices are + handled in a different way and therefore lines have no faces. The #Triangulation# class offers functions like #begin_active# which gives you an iterator to the first active cell. There are quite a lot of functions @@ -1000,6 +1023,10 @@ class Triangulation : public TriaDimensionInfo { typedef typename TriaDimensionInfo::cell_iterator cell_iterator; typedef typename TriaDimensionInfo::active_cell_iterator active_cell_iterator; + typedef typename TriaDimensionInfo::raw_substruct_iterator raw_substruct_iterator; + typedef typename TriaDimensionInfo::substruct_iterator substruct_iterator; + typedef typename TriaDimensionInfo::active_substruct_iterator active_substruct_iterator; + /** * Create a triangulation and create * the first level of the hierarchy. diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index bff8f3d57f..14f4040a1c 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -731,10 +731,18 @@ class TriaSubstructAccessor<1> : public LineAccessor<1> { * Constructor */ TriaSubstructAccessor (Triangulation<1> *tria, - const int level, - const int index, - const void *local_data) : + const int level, + const int index, + const void *local_data) : LineAccessor<1> (tria,level,index,local_data) {}; + + // do this here, since this way the + // CellAccessor has the possibility to know + // what a substruct_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef void * substruct_iterator; }; @@ -749,10 +757,18 @@ class TriaSubstructAccessor<2> : public QuadAccessor<2> { * Constructor */ TriaSubstructAccessor (Triangulation<2> *tria, - const int level, - const int index, - const void *local_data) : + const int level, + const int index, + const void *local_data) : QuadAccessor<2> (tria,level,index,local_data) {}; + + // do this here, since this way the + // CellAccessor has the possibility to know + // what a substruct_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef TriaIterator<2,LineAccessor<2> > substruct_iterator; }; @@ -853,6 +869,16 @@ class CellAccessor : public TriaSubstructAccessor { */ TriaIterator > child (const unsigned int i) const; + /** + * Return an iterator to the #i#th face + * of this cell. + * + * This function is not implemented in 1D, + * and maps to QuadAccessor::line in 2D. + */ + typename TriaSubstructAccessor::substruct_iterator + face (const unsigned int i) const; + /** * Return the material id of this cell. */ diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 7b2772ef41..17a5873dc2 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -10,6 +10,7 @@ #include + // forward declaration template class Triangulation; template class DoFHandler; @@ -18,6 +19,30 @@ template class Quadrature; template class DataOut; template class Function; +template , class Alloc = alloc> class map; + + + + +/** + 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 +}; + + /** @@ -37,12 +62,15 @@ template class Function; The #assemble# member function does the assemblage of the system matrix and the given number of right hand sides. It does the following steps: \begin{itemize} + \item Initialize solution vector with zero entries. \item Create sparsity pattern of the system matrix and condense it with the constraints induced by hanging nodes. \item Initialize an assembler object. \item Loop over all cells and assemble matrix and vectors using the given quadrature formula and the equation object which contains the weak formulation of the equation. + \item Apply Dirichlet boundary conditions. See the section on boundary + conditions for more details. \item Condense the system matrix and right hand side with the constraints induced by hanging nodes. \end{itemize} @@ -54,10 +82,120 @@ template class Function; which results after having called the #assemble# function is solved. After this, the solution vector is distributed again, i.e. the constrained nodes are given their correct values. + + + {\bf Boundary conditions} + + During assemblage of matrices and right hand side, use is made of dirichlet + boundary conditions (in short: bc) specified to the #assemble# function. You + can 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 this 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 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. 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. + + + {\bf 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 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#. + 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. + + To get the {\it global} L_1 error, you have to sum up the entries in + #difference#, e.g. using the STL function + #accumulate (d.begin(), d.end(), 0)#, if #d# is the difference vector. + For the global L_2 difference, you have to sum up the squares of the + entries and take the root of the sum. 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 STL function #max_element (d.begin(), d.end())#. */ template class ProblemBase { 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 > DirichletBC; + /** * Constructor. Use this triangulation and * degree of freedom object during the @@ -80,16 +218,20 @@ class ProblemBase { * Initiate the process of assemblage of * vectors and system matrix. Use the * given equation object and the given - * quadrature formula. + * quadrature formula. Also use the list + * of dirichlet boundary value functions + * (by default, no dirichlet bc are assumed + * which means that all bc are included + * into the weak formulation). * * For what exactly happens here, refer to * the general doc of this class. */ virtual void assemble (const Equation &equation, const Quadrature &q, - const FiniteElement &fe); - - + const FiniteElement &fe, + const DirichletBC &dirichlet_bc = DirichletBC()); + /** * Solve the system of equations. */ @@ -100,41 +242,16 @@ class ProblemBase { * the solution computed before and * the reference solution, which * is given as a continuous function - * object. The integration is - * performed using the given quadrature - * rule and assumes that the given - * finite element objects equals - * that used for the computation of - * the solution. - * - * The result ist stored in the - * #difference# object, where each - * entry equals the L_1 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 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. + * object. * - * To get the global L_1 error, - * you have to sum up the entries - * in #difference#, e.g. using the STL - * function - * #accumulate (d.begin(), d.end(), 0)#, - * if #d# is the difference vector. - */ - void integrate_L1_difference (const Function &exact_solution, - vector &difference, - const Quadrature &q, - const FiniteElement &fe) const; + * See the general documentation of this + * class for more information. + */ + void integrate_difference (const Function &exact_solution, + vector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm) const; /** * Initialize the #DataOut# object with @@ -178,6 +295,10 @@ class ProblemBase { * Exception */ DeclException0 (ExcInvalidFE); + /** + * Exception + */ + DeclException0 (ExcNotImplemented); protected: /** @@ -216,8 +337,19 @@ class ProblemBase { * hanging nodes. */ 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 DirichletBC &dirichlet_bc); - friend class Assembler; + friend class Assembler; }; diff --git a/deal.II/deal.II/include/numerics/data_io.h b/deal.II/deal.II/include/numerics/data_io.h index 6d35cb4956..e77800bab8 100644 --- a/deal.II/deal.II/include/numerics/data_io.h +++ b/deal.II/deal.II/include/numerics/data_io.h @@ -21,8 +21,6 @@ template class DoFHandler; class dVector; -template class map; -template struct less; diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index cc905ae674..17e7817ad8 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -256,13 +256,29 @@ DoFCellAccessor::child (const unsigned int i) const { + +DoFSubstructAccessor<1>::substruct_iterator +DoFCellAccessor<1>::face (const unsigned int) const { + Assert (false, ExcNotUsefulForThisDimension()); + return 0; +}; + + + +DoFSubstructAccessor<2>::substruct_iterator +DoFCellAccessor<2>::face (const unsigned int i) const { + return line(i); +}; + + + void -DoFCellAccessor<1>::dof_indices (vector &dof_indices) const { +DoFCellAccessor<1>::get_dof_indices (vector &dof_indices) const { Assert (dof_handler != 0, ExcInvalidObject()); Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == 0, ExcVectorNotEmpty()); - dof_indices.reserve (dof_indices.size() + - dof_handler->get_selected_fe().total_dofs); + dof_indices.reserve (dof_handler->get_selected_fe().total_dofs); for (unsigned int vertex=0; vertex<2; ++vertex) for (unsigned int d=0; dget_selected_fe().dofs_per_vertex; ++d) dof_indices.push_back (vertex_dof_index(vertex,d)); @@ -274,23 +290,72 @@ DoFCellAccessor<1>::dof_indices (vector &dof_indices) const { void -DoFCellAccessor<2>::dof_indices (vector &dof_indices) const { +DoFCellAccessor<2>::get_dof_indices (vector &dof_indices) const { Assert (dof_handler != 0, ExcInvalidObject()); Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == 0, ExcVectorNotEmpty()); - dof_indices.reserve (dof_indices.size() + - dof_handler->get_selected_fe().total_dofs); + dof_indices.reserve (dof_handler->get_selected_fe().total_dofs); for (unsigned int vertex=0; vertex<4; ++vertex) for (unsigned int d=0; dget_selected_fe().dofs_per_vertex; ++d) dof_indices.push_back (vertex_dof_index(vertex,d)); for (unsigned int line=0; line<4; ++line) for (unsigned int d=0; dget_selected_fe().dofs_per_line; ++d) dof_indices.push_back (this->line(line)->dof_index(d)); + for (unsigned int d=0; dget_selected_fe().dofs_per_quad; ++d) + dof_indices.push_back (dof_index(d)); +}; + + + + +void +DoFCellAccessor<1>::get_dof_values (const dVector &values, + vector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == 0, ExcVectorNotEmpty()); + Assert ((unsigned int)values.n() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + dof_values.reserve (dof_handler->get_selected_fe().total_dofs); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; dget_selected_fe().dofs_per_vertex; ++d) + dof_values.push_back (values(vertex_dof_index(vertex,d))); + for (unsigned int d=0; dget_selected_fe().dofs_per_line; ++d) + dof_values.push_back (values(dof_index(d))); }; +void +DoFCellAccessor<2>::get_dof_values (const dVector &values, + vector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == 0, ExcVectorNotEmpty()); + Assert ((unsigned int)values.n() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + dof_values.reserve (dof_handler->get_selected_fe().total_dofs); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dget_selected_fe().dofs_per_vertex; ++d) + dof_values.push_back (values(vertex_dof_index(vertex,d))); + for (unsigned int line=0; line<4; ++line) + for (unsigned int d=0; dget_selected_fe().dofs_per_line; ++d) + dof_values.push_back (values(this->line(line)->dof_index(d))); + for (unsigned int d=0; dget_selected_fe().dofs_per_quad; ++d) + dof_values.push_back (values(dof_index(d))); +}; + + + + + + + + // explicit instantiations template class DoFLineAccessor<1,CellAccessor<1> >; diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index 8dec8ca805..643ceb295c 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -1037,8 +1037,8 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // both cells have no children { vector old_dofs, new_dofs; - old_cell->dof_indices (old_dofs); - new_cell->dof_indices (new_dofs); + old_cell->get_dof_indices (old_dofs); + new_cell->get_dof_indices (new_dofs); Assert (old_dofs.size() == selected_fe->total_dofs, ExcInternalError ()); Assert (new_dofs.size() == selected_fe->total_dofs, @@ -1062,7 +1062,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // numbers of old dofs vector old_dof_indices; - old_cell->dof_indices (old_dof_indices); + old_cell->get_dof_indices (old_dof_indices); Assert (old_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1071,7 +1071,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat { // numbers of child dofs vector child_dof_indices; - child[c]->dof_indices (child_dof_indices); + child[c]->get_dof_indices (child_dof_indices); Assert (child_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1094,7 +1094,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // numbers of new dofs vector new_dof_indices; - new_cell->dof_indices(new_dof_indices); + new_cell->get_dof_indices(new_dof_indices); Assert (new_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1103,7 +1103,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat { // numbers of child dofs vector child_dof_indices; - child[c]->dof_indices (child_dof_indices); + child[c]->get_dof_indices (child_dof_indices); Assert (child_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1133,8 +1133,8 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // both cells have no children { vector old_dofs, new_dofs; - old_cell->dof_indices (old_dofs); - new_cell->dof_indices (new_dofs); + old_cell->get_dof_indices (old_dofs); + new_cell->get_dof_indices (new_dofs); Assert (old_dofs.size() == selected_fe->total_dofs, ExcInternalError ()); Assert (new_dofs.size() == selected_fe->total_dofs, @@ -1158,7 +1158,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // numbers of old dofs vector old_dof_indices; - old_cell->dof_indices (old_dof_indices); + old_cell->get_dof_indices (old_dof_indices); Assert (old_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1167,7 +1167,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat { // numbers of child dofs vector child_dof_indices; - child[c]->dof_indices (child_dof_indices); + child[c]->get_dof_indices (child_dof_indices); Assert (child_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1191,7 +1191,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat // numbers of new dofs vector new_dof_indices; - new_cell->dof_indices(new_dof_indices); + new_cell->get_dof_indices(new_dof_indices); Assert (new_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); @@ -1200,7 +1200,7 @@ void DoFHandler::transfer_cell (const typename DoFHandler::cell_iterat { // numbers of child dofs vector child_dof_indices; - child[c]->dof_indices (child_dof_indices); + child[c]->get_dof_indices (child_dof_indices); Assert (child_dof_indices.size() == selected_fe->total_dofs, ExcInternalError ()); diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index c9df0fa519..a6228c505c 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -529,6 +529,14 @@ void CellAccessor<1>::set_material_id (const unsigned char mat_id) const { +Triangulation<1>::substruct_iterator +CellAccessor<1>::face (const unsigned int) const { + Assert (false, ExcNotUsefulForThisDimension()); + return 0; +}; + + + /*------------------------ Functions: CellAccessor<2> -----------------------*/ @@ -554,6 +562,13 @@ void CellAccessor<2>::set_material_id (const unsigned char mat_id) const { }; + +Triangulation<2>::substruct_iterator +CellAccessor<2>::face (const unsigned int i) const { + return line(i); +}; + + /*------------------------ Functions: CellAccessor -----------------------*/ diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index 44f76cedad..218a66cddb 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -96,11 +96,12 @@ void Assembler::assemble (const Equation &equation) { present_level, present_index), fe); + const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs; // clear cell matrix if (assemble_matrix) - for (unsigned int i=0; iget_selected_fe().total_dofs; ++i) - for (unsigned int j=0; jget_selected_fe().total_dofs; ++j) + for (unsigned int i=0; i::assemble (const Equation &equation) { // get indices of dofs vector dofs; - dof_indices (dofs); + get_dof_indices (dofs); // distribute cell matrix if (assemble_matrix) - for (unsigned int i=0; iget_selected_fe().total_dofs; ++i) - for (unsigned int j=0; jget_selected_fe().total_dofs; ++j) + for (unsigned int i=0; iget_selected_fe().total_dofs; ++j) + for (unsigned int j=0; j + + + + +inline double sqr (double x) { + return x*x; +}; @@ -22,8 +30,8 @@ ProblemBase::ProblemBase (Triangulation *tria, DoFHandler *dof) : tria(tria), dof_handler(dof), - system_sparsity(1,1,1), // dummy initialisation, is later reinit'd - system_matrix() // dummy initialisation, is later reinit'd + system_sparsity(1,1,1), // dummy initialisation, is later reinit'd + system_matrix() // dummy initialisation, is later reinit'd { Assert (tria == &dof->get_tria(), ExcDofAndTriaDontMatch()); }; @@ -36,9 +44,10 @@ ProblemBase::~ProblemBase () {}; template -void ProblemBase::assemble (const Equation &equation, - const Quadrature &quadrature, - const FiniteElement &fe) { +void ProblemBase::assemble (const Equation &equation, + const Quadrature &quadrature, + const FiniteElement &fe, + const DirichletBC &dirichlet_bc) { system_sparsity.reinit (dof_handler->n_dofs(), dof_handler->n_dofs(), dof_handler->max_couplings_between_dofs()); @@ -53,6 +62,8 @@ void ProblemBase::assemble (const Equation &equation, // reinite system matrix system_matrix.reinit (system_sparsity); + // reinit solution vector, preset + // with zeroes. solution.reinit (dof_handler->n_dofs()); // create assembler @@ -73,6 +84,11 @@ void ProblemBase::assemble (const Equation &equation, } while ((++assembler).state() == valid); + // apply Dirichlet bc as described + // in the docs + apply_dirichlet_bc (system_matrix, solution, + right_hand_side, dirichlet_bc); + // condense system matrix in-place constraints.condense (system_matrix); @@ -100,28 +116,135 @@ void ProblemBase::solve () { -/* template -void ProblemBase::integrate_L1_difference (const Function &exact_solution, - vector &difference, - const Quadrature &q, - const FiniteElement &fe) const { +void ProblemBase::integrate_difference (const Function &exact_solution, + vector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm) const { Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE()); difference.erase (difference.begin(), difference.end()); difference.reserve (tria->n_cells()); + FEValues fe_values(fe, q); + // loop over all cells - DoFHandler::active_dof_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - for (; cell != end; ++cell) + // (we need an iterator on the triangulation for + // fe_values.reinit, but we also need an iterator + // on the dofhandler; we could generate the first + // out of the second or vice versa each time + // we need them, but conversion seems to be slow. + // therefore, we keep two iterators, hoping that + // incrementing is faster than conversion...) + DoFHandler::active_cell_iterator cell = dof_handler->begin_active(), + endc = dof_handler->end(); + Triangulation::active_cell_iterator tria_cell = dof_handler->get_tria().begin(); + + for (; cell != endc; ++cell, ++tria_cell) { double diff=0; + // initialize for this cell + fe_values.reinit (tria_cell, fe); + + switch (norm) + { + case mean: + case L1_norm: + case L2_norm: + case Linfty_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_dofs = fe.total_dofs; + const dFMatrix & shape_values = fe_values.get_shape_values(); + vector dof_values; + cell->get_dof_values (solution, dof_values); + + vector psi; + + // in praxi: first compute + // exact solution vector + exact_solution.value_list (fe_values.get_quadrature_points(), + psi); + // then subtract finite element + // solution + for (unsigned int j=0; j::fill_data (DataOut &out) const { }; -/* + template pair ProblemBase::get_solution_name () const { - return make_pair("solution", ""); + return pair("solution", ""); +}; + + + + +template +void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, + dVector &solution, + dVector &right_hand_side, + const DirichletBC &dirichlet_bc) { +// l; }; -*/ -