From 92432bc28353c9887cb8e77927c0fb4a4b4f6323 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 27 Mar 1998 15:32:42 +0000 Subject: [PATCH] Support for faces of cells and starting support for Dirichlet boundary conditions (untested). git-svn-id: https://svn.dealii.org/trunk@99 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/function.h | 64 ++++++- .../deal.II/Attic/examples/poisson/poisson.cc | 3 + deal.II/deal.II/include/dofs/dof_accessor.h | 56 ++++--- deal.II/deal.II/include/dofs/dof_handler.h | 124 +++++++++++++- deal.II/deal.II/include/grid/tria.h | 135 +++++++++++++-- deal.II/deal.II/include/numerics/base.h | 47 +++++- deal.II/deal.II/source/dofs/dof_accessor.cc | 79 ++++----- .../deal.II/source/dofs/dof_constraints.cc | 142 +++++++++------- deal.II/deal.II/source/dofs/dof_handler.cc | 158 ++++++++++++++++++ deal.II/deal.II/source/grid/tria.cc | 155 +++++++++++++++++ deal.II/deal.II/source/grid/tria_accessor.cc | 4 +- deal.II/deal.II/source/numerics/base.cc | 130 +++++++++++++- tests/big-tests/poisson/poisson.cc | 3 + 13 files changed, 941 insertions(+), 159 deletions(-) diff --git a/deal.II/base/include/base/function.h b/deal.II/base/include/base/function.h index aac5e22336..5ecd440b72 100644 --- a/deal.II/base/include/base/function.h +++ b/deal.II/base/include/base/function.h @@ -5,19 +5,75 @@ /*---------------------------- function.h ---------------------------*/ -//forward declaration -template class Point; -template class vector; +#include +#include +#include + +/** + This class is a model for a continuous function. It returns the value + of the function at a given point through the #operator ()# member function, + which is virtual. It also has a function to return a whole list of function + values at different points to reduce the overhead of the virtual function + calls; this function is preset to successively call the function returning + one value at a time. + + There are other functions return the gradient of the function at one or + several points. You only have to overload those functions you need; the + functions returning several values at a time will call those returning + only one value, while those ones will throw an exception when called but + not overloaded. + + Unless only called a very small number of times, you should overload + both those functions returning only one value as well as those returning + a whole array, since the cost of evaluation of a point value is often + less than the virtual function call itself. +*/ template class Function { public: - virtual double operator () (const Point &p) const = 0; + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; + + /** + * Set #values# to the point values + * of the function at the #points#. + * It is assumed that #values# be + * empty. + */ virtual void value_list (const vector > &points, vector &values) const; + + /** + * Return the gradient of the function + * at the given point. + */ + virtual Point gradient (const Point &p) const; + + /** + * Set #gradients# to the gradients of + * the function at the #points#. + * It is assumed that #values# be + * empty. + */ + virtual void gradient_list (const vector > &points, + vector > &gradients) const; + + + /** + * Exception + */ + DeclException0 (ExcPureFunctionCalled); + /** + * Exception + */ + DeclException0 (ExcVectorNotEmpty); }; diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index d7761bbc59..a1f6444408 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -9,12 +9,15 @@ #include #include #include +#include #include #include #include #include #include + +#include #include #include extern "C" { diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index 63df424b39..925090e51b 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -55,6 +55,10 @@ class DoFAccessor { << "Invalid index " << arg1 << ", index must be between " << arg2 << " and " << arg3 << "."); + /** + * Exception + */ + DeclException0 (ExcVectorNotEmpty); protected: /** @@ -164,6 +168,18 @@ class DoFLineAccessor : public BaseClass, public DoFAccessor { const unsigned int i, const int index) const; + /** + * Return the indices of the dofs of this + * line in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, + * dofs on line. + * + * The dof indices are appended to the end + * of the vector. It is assumed that + * the vector be empty beforehand. + */ + void get_dof_indices (vector &dof_indices) const; + /** * Return the #i#th child as a DoF line * iterator. This function is needed since @@ -237,6 +253,19 @@ class DoFQuadAccessor : public BaseClass, public DoFAccessor { const unsigned int i, const int index) const; + /** + * Return the indices of the dofs of this + * quad 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 dof indices are appended to the end + * of the vector. It is assumed that + * the vector be empty beforehand. + */ + void get_dof_indices (vector &dof_indices) const; + /** * Return a pointer to the #i#th line * bounding this #Quad#. @@ -308,11 +337,11 @@ class DoFSubstructAccessor<1> : public DoFLineAccessor<1,CellAccessor<1> > { 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 + // what a face_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; + typedef void * face_iterator; }; @@ -333,11 +362,11 @@ class DoFSubstructAccessor<2> : public DoFQuadAccessor<2,CellAccessor<2> > { 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 + // what a face_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; + typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > face_iterator; }; @@ -393,22 +422,9 @@ class DoFCellAccessor : public DoFSubstructAccessor { * This function is not implemented in 1D, * and maps to DoFQuadAccessor::line in 2D. */ - typename DoFSubstructAccessor::substruct_iterator + typename DoFSubstructAccessor::face_iterator face (const unsigned int i) const; - /** - * Return the indices of 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 dof indices are appended to the end - * of the vector. It is assumed that - * the vector be empty beforehand. - */ - void get_dof_indices (vector &dof_indices) const; - /** * Return the value of the given vector * restricted to the dofs of this @@ -424,10 +440,6 @@ class DoFCellAccessor : public DoFSubstructAccessor { void get_dof_values (const dVector &values, vector &dof_values) const; - /** - * Exception - */ - DeclException0 (ExcVectorNotEmpty); /** * Exception */ diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 360b6f434c..20aebc1c9f 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -153,9 +153,9 @@ class DoFDimensionInfo<1> { 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; + typedef void * raw_face_iterator; + typedef void * face_iterator; + typedef void * active_face_iterator; }; @@ -182,9 +182,9 @@ class DoFDimensionInfo<2> { 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; + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; }; @@ -436,9 +436,9 @@ 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; + typedef typename DoFDimensionInfo::raw_face_iterator raw_face_iterator; + typedef typename DoFDimensionInfo::face_iterator face_iterator; + typedef typename DoFDimensionInfo::active_face_iterator active_face_iterator; /** @@ -703,6 +703,108 @@ class DoFHandler : public DoFDimensionInfo { /*---------------------------------------*/ + /** + * @name Face iterator functions + */ + /*@{*/ + /** + * Return iterator to the first face, used + * or not, on level #level#. If a level + * has no faces, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 2D and #begin_raw_quad# in 3D. + */ + raw_face_iterator begin_raw_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first used face + * on level #level#. + * + * This function calls #begin_line# + * in 2D and #begin_quad# in 3D. + */ + face_iterator begin_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * face on level #level#. + * + * This function calls #begin_active_line# + * in 2D and #begin_active_quad# in 3D. + */ + active_face_iterator begin_active_face(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 2D and #end_quad# in 3D. + */ + raw_face_iterator end_face () const; + + /** + * Return an iterator pointing to the + * last face, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face () const; + + /** + * Return an iterator pointing to the last + * face of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used face. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face () const; + + /** + * Return an iterator pointing to the last + * used face on level #level#. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active face. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face () const; + + /** + * Return an iterator pointing to the last + * active face on level #level#. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face (const unsigned int level) const; + //@} + + + /*---------------------------------------*/ + /** * @name Line iterator functions */ @@ -928,6 +1030,10 @@ class DoFHandler : public DoFDimensionInfo { DeclException1 (ExcMatrixHasWrongSize, int, << "The matrix has the wrong dimension " << arg1); + /** + * Exception + */ + DeclException0 (ExcFunctionNotUseful); protected: /** diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 2685ec001e..97e3e13103 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -463,7 +463,7 @@ class TriaDimensionInfo; #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 + The same applies for the #face_iterator# types, since lines have no substructures apart from vertices, which are handled in a different way, however. */ @@ -481,9 +481,9 @@ class TriaDimensionInfo<1> { 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; + typedef void * raw_face_iterator; + typedef void * face_iterator; + typedef void * active_face_iterator; }; @@ -508,9 +508,9 @@ class TriaDimensionInfo<1> { 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; + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; \end{verbatim} */ class TriaDimensionInfo<2> { @@ -527,9 +527,9 @@ class TriaDimensionInfo<2> { 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; + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; }; @@ -601,9 +601,9 @@ class TriaDimensionInfo<2> { 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; + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; \end{verbatim} By using the cell iterators, you can write code nearly independent of @@ -1023,9 +1023,9 @@ 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; + typedef typename TriaDimensionInfo::raw_face_iterator raw_face_iterator; + typedef typename TriaDimensionInfo::face_iterator face_iterator; + typedef typename TriaDimensionInfo::active_face_iterator active_face_iterator; /** * Create a triangulation and create @@ -1338,6 +1338,109 @@ class Triangulation : public TriaDimensionInfo { //@} /*---------------------------------------*/ + /*---------------------------------------*/ + + /** + * @name Face iterator functions + */ + /*@{*/ + /** + * Return iterator to the first face, used + * or not, on level #level#. If a level + * has no faces, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 2D and #begin_raw_quad# in 3D. + */ + raw_face_iterator begin_raw_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first used face + * on level #level#. + * + * This function calls #begin_line# + * in 2D and #begin_quad# in 3D. + */ + face_iterator begin_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * face on level #level#. + * + * This function calls #begin_active_line# + * in 2D and #begin_active_quad# in 3D. + */ + active_face_iterator begin_active_face(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 2D and #end_quad# in 3D. + */ + raw_face_iterator end_face () const; + + /** + * Return an iterator pointing to the + * last face, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face () const; + + /** + * Return an iterator pointing to the last + * face of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used face. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face () const; + + /** + * Return an iterator pointing to the last + * used face on level #level#. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active face. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face () const; + + /** + * Return an iterator pointing to the last + * active face on level #level#. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face (const unsigned int level) const; + //@} + + + /*---------------------------------------*/ /** * @name Line iterator functions diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 17a5873dc2..572239610d 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -18,6 +18,8 @@ template class FiniteElement; template class Quadrature; template class DataOut; template class Function; +template class Equation; +template class Assembler; template , class Alloc = alloc> class map; @@ -110,6 +112,16 @@ enum NormType { 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 @@ -123,7 +135,9 @@ enum NormType { 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 + 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 @@ -135,6 +149,13 @@ enum NormType { 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. {\bf Computing errors} @@ -194,7 +215,7 @@ class ProblemBase { * See the general documentation of this * class for more detail. */ - typedef map > DirichletBC; + typedef map* > DirichletBC; /** * Constructor. Use this triangulation and @@ -282,6 +303,18 @@ class ProblemBase { * want anything else. */ 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 DirichletBC &dirichlet_bc, + map &boundary_values) const; /** * Exception @@ -299,6 +332,10 @@ class ProblemBase { * Exception */ DeclException0 (ExcNotImplemented); + /** + * Exception + */ + DeclException0 (ExcInvalidBoundaryIndicator); protected: /** @@ -344,9 +381,9 @@ class ProblemBase { * as described in the general * documentation. */ - void apply_dirichlet_bc (dSMatrix &matrix, - dVector &solution, - dVector &right_hand_side, + void apply_dirichlet_bc (dSMatrix &matrix, + dVector &solution, + dVector &right_hand_side, const DirichletBC &dirichlet_bc); friend class Assembler; diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index 17e7817ad8..8de61c1d6e 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -83,6 +83,23 @@ void DoFLineAccessor::set_vertex_dof_index (const unsigned int ve +template +void +DoFLineAccessor::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_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)); + for (unsigned int d=0; dget_selected_fe().dofs_per_line; ++d) + dof_indices.push_back (dof_index(d)); +}; + + + template TriaIterator > DoFLineAccessor::child (const unsigned int i) const { @@ -177,6 +194,27 @@ void DoFQuadAccessor::set_vertex_dof_index (const unsigned int ve +template +void +DoFQuadAccessor::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_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)); +}; + + + + template TriaIterator > > DoFQuadAccessor::line (const unsigned int i) const { @@ -257,7 +295,7 @@ DoFCellAccessor::child (const unsigned int i) const { -DoFSubstructAccessor<1>::substruct_iterator +DoFSubstructAccessor<1>::face_iterator DoFCellAccessor<1>::face (const unsigned int) const { Assert (false, ExcNotUsefulForThisDimension()); return 0; @@ -265,50 +303,13 @@ DoFCellAccessor<1>::face (const unsigned int) const { -DoFSubstructAccessor<2>::substruct_iterator +DoFSubstructAccessor<2>::face_iterator DoFCellAccessor<2>::face (const unsigned int i) const { return line(i); }; -void -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_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)); - for (unsigned int d=0; dget_selected_fe().dofs_per_line; ++d) - dof_indices.push_back (dof_index(d)); -}; - - - - -void -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_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 { diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 1ec514f779..ba2d1eaf33 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -87,7 +87,7 @@ void ConstraintMatrix::clear () { void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, dSMatrixStruct &condensed) const { Assert (sorted == true, ExcMatrixNotClosed()); - Assert (uncondensed.compressed == true, ExcMatrixNotClosed()); + Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed()); Assert (uncondensed.n_rows() == uncondensed.n_cols(), ExcMatrixNotSquare()); @@ -116,21 +116,23 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, next_constraint = lines.begin(); - for (int row=0; row::const_iterator c = lines.begin(); - while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c; + while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j]) + ++c; for (unsigned int q=0; q!=(*c).entries.size(); ++q) condensed.add (new_line[row], new_line[(*c).entries[q].first]); @@ -138,13 +140,14 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, else // line must be distributed { - for (int j=uncondensed.rowstart[row]; j::const_iterator c = lines.begin(); - while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c; + while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j]) ++c; for (unsigned int p=0; p!=(*c).entries.size(); ++p) for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q) @@ -171,7 +174,7 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const { Assert (sorted == true, ExcMatrixNotClosed()); - Assert (sparsity.compressed == false, ExcMatrixIsClosed()); + Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed()); Assert (sparsity.n_rows() == sparsity.n_cols(), ExcMatrixNotSquare()); @@ -199,46 +202,52 @@ void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const { for (int row=0; row::const_iterator c = lines.begin(); - while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c; + while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c; for (unsigned int q=0; q!=(*c).entries.size(); ++q) // distribute to rows with // appropriate weight condensed.add (new_line[row], new_line[(*c).entries[q].first], - uncondensed.val[j] * (*c).entries[q].second); + uncondensed.global_entry(j) * (*c).entries[q].second); } else // line must be distributed { - for (int j=uncondensed_struct.rowstart[row]; j::const_iterator c = lines.begin(); - while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c; + while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c; for (unsigned int p=0; p!=(*c).entries.size(); ++p) for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q) condensed.add (new_line[(*next_constraint).entries[q].first], new_line[(*c).entries[p].first], - uncondensed.val[j] * + uncondensed.global_entry(j) * (*next_constraint).entries[q].second * (*c).entries[p].second); }; @@ -347,7 +359,7 @@ void ConstraintMatrix::condense (dSMatrix &uncondensed) const { const dSMatrixStruct &sparsity = uncondensed.get_sparsity_pattern (); Assert (sorted == true, ExcMatrixNotClosed()); - Assert (sparsity.compressed == true, ExcMatrixNotClosed()); + Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed()); Assert (sparsity.n_rows() == sparsity.n_cols(), ExcMatrixNotSquare()); @@ -378,67 +390,75 @@ void ConstraintMatrix::condense (dSMatrix &uncondensed) const { if (distribute[row] == -1) // regular line. loop over cols - for (int j=sparsity.rowstart[row]; j::last_active (const unsigned int level) const { //------------------------------------------------------------------ + + +DoFDimensionInfo<1>::raw_face_iterator +DoFHandler<1>::begin_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::face_iterator +DoFHandler<1>::begin_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::active_face_iterator +DoFHandler<1>::begin_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::raw_face_iterator +DoFHandler<1>::end_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::raw_face_iterator +DoFHandler<1>::last_raw_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::raw_face_iterator +DoFHandler<1>::last_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::face_iterator +DoFHandler<1>::last_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::face_iterator +DoFHandler<1>::last_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::active_face_iterator +DoFHandler<1>::last_active_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +DoFDimensionInfo<1>::active_face_iterator +DoFHandler<1>::last_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + + + +DoFDimensionInfo<2>::raw_face_iterator +DoFHandler<2>::begin_raw_face (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +DoFDimensionInfo<2>::face_iterator +DoFHandler<2>::begin_face (const unsigned int level) const { + return begin_line (level); +}; + + + +DoFDimensionInfo<2>::active_face_iterator +DoFHandler<2>::begin_active_face (const unsigned int level) const { + return begin_active_line (level); +}; + + + +DoFDimensionInfo<2>::raw_face_iterator +DoFHandler<2>::end_face () const { + return end_line (); +}; + + + +DoFDimensionInfo<2>::raw_face_iterator +DoFHandler<2>::last_raw_face () const { + return last_raw_line (); +}; + + + +DoFDimensionInfo<2>::raw_face_iterator +DoFHandler<2>::last_raw_face (const unsigned int level) const { + return last_raw_line (level); +}; + + + +DoFDimensionInfo<2>::face_iterator +DoFHandler<2>::last_face () const { + return last_line (); +}; + + + +DoFDimensionInfo<2>::face_iterator +DoFHandler<2>::last_face (const unsigned int level) const { + return last_line (level); +}; + + + +DoFDimensionInfo<2>::active_face_iterator +DoFHandler<2>::last_active_face () const { + return last_active_line (); +}; + + + +DoFDimensionInfo<2>::active_face_iterator +DoFHandler<2>::last_active_face (const unsigned int level) const { + return last_active_line (level); +}; + + + +// --------------------------------------------------------------- + + + template typename DoFHandler::raw_line_iterator DoFHandler::begin_raw_line (const unsigned int level) const { diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index fe955a0547..6130774957 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -956,6 +956,161 @@ Triangulation<2>::last_active (const unsigned int level) const { + + +TriaDimensionInfo<1>::raw_face_iterator +Triangulation<1>::begin_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::face_iterator +Triangulation<1>::begin_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::active_face_iterator +Triangulation<1>::begin_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::raw_face_iterator +Triangulation<1>::end_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::raw_face_iterator +Triangulation<1>::last_raw_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::raw_face_iterator +Triangulation<1>::last_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::face_iterator +Triangulation<1>::last_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::face_iterator +Triangulation<1>::last_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::active_face_iterator +Triangulation<1>::last_active_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +TriaDimensionInfo<1>::active_face_iterator +Triangulation<1>::last_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + + +TriaDimensionInfo<2>::raw_face_iterator +Triangulation<2>::begin_raw_face (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +TriaDimensionInfo<2>::face_iterator +Triangulation<2>::begin_face (const unsigned int level) const { + return begin_line (level); +}; + + + +TriaDimensionInfo<2>::active_face_iterator +Triangulation<2>::begin_active_face (const unsigned int level) const { + return begin_active_line (level); +}; + + + +TriaDimensionInfo<2>::raw_face_iterator +Triangulation<2>::end_face () const { + return end_line (); +}; + + + +TriaDimensionInfo<2>::raw_face_iterator +Triangulation<2>::last_raw_face () const { + return last_raw_line (); +}; + + + +TriaDimensionInfo<2>::raw_face_iterator +Triangulation<2>::last_raw_face (const unsigned int level) const { + return last_raw_line (level); +}; + + + +TriaDimensionInfo<2>::face_iterator +Triangulation<2>::last_face () const { + return last_line (); +}; + + + +TriaDimensionInfo<2>::face_iterator +Triangulation<2>::last_face (const unsigned int level) const { + return last_line (level); +}; + + + +TriaDimensionInfo<2>::active_face_iterator +Triangulation<2>::last_active_face () const { + return last_active_line (); +}; + + + +TriaDimensionInfo<2>::active_face_iterator +Triangulation<2>::last_active_face (const unsigned int level) const { + return last_active_line (level); +}; + + + + + template typename TriaDimensionInfo::raw_line_iterator Triangulation::begin_raw_line (unsigned int level) const { diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index a6228c505c..ff74f023cd 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -529,7 +529,7 @@ void CellAccessor<1>::set_material_id (const unsigned char mat_id) const { -Triangulation<1>::substruct_iterator +Triangulation<1>::face_iterator CellAccessor<1>::face (const unsigned int) const { Assert (false, ExcNotUsefulForThisDimension()); return 0; @@ -563,7 +563,7 @@ void CellAccessor<2>::set_material_id (const unsigned char mat_id) const { -Triangulation<2>::substruct_iterator +Triangulation<2>::face_iterator CellAccessor<2>::face (const unsigned int i) const { return line(i); }; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 08e75e5239..84c946b708 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -11,6 +11,7 @@ #include "../../../mia/vectormemory.h" #include "../../../mia/cg.h" +#include #include @@ -273,7 +274,134 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, dVector &solution, dVector &right_hand_side, const DirichletBC &dirichlet_bc) { -// l; + 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, boundary_values); + + map::const_iterator dof, endd; + const unsigned int n_dofs = (unsigned int)matrix.m(); + const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern(); + const int *sparsity_rowstart = sparsity.get_rowstart_indices(), + *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 (int j=sparsity_rowstart[(*dof).first]; + j::make_boundary_value_list (const DirichletBC &, + map &) const { + Assert (false, ExcNotImplemented()); +}; + + + + +template +void +ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, + map &boundary_values) const { + DoFHandler::active_face_iterator face = dof_handler->begin_active_face(), + endf = dof_handler->end_face(); + DirichletBC::const_iterator function_ptr; + 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 + vector face_dofs; + vector > dof_locations; + vector dof_values; + + face->get_dof_indices (face_dofs); +//physical location + (*function_ptr).second->value_list (dof_locations, dof_values); + + // enter into list + for (unsigned int i=0; i #include #include +#include #include #include #include #include #include + +#include #include #include extern "C" { -- 2.39.5