From: Wolfgang Bangerth Date: Fri, 3 Apr 1998 17:44:30 +0000 (+0000) Subject: Full support for Dirichlet BC and bug fixes. X-Git-Tag: v8.0.0~23120 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4b269082a483e98c88ab68c94ca54265049d9320;p=dealii.git Full support for Dirichlet BC and bug fixes. git-svn-id: https://svn.dealii.org/trunk@134 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 20aebc1c9f..1d39fb5de8 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -1095,7 +1095,14 @@ class DoFHandler : public DoFDimensionInfo { /** * Store a copy of the finite element given - * latest for the distribution of dofs. + * latest for the distribution of dofs. In + * fact, since the FE class itself has not + * much functionality, this object only + * stores numbers such as the number of + * dofs per line etc. Calling any of the + * more specific functions will result in + * an error about calling a pure virtual + * function. */ FiniteElement *selected_fe; diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 592dbc5ba3..f1188b0b49 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -35,7 +35,7 @@ template class Quadrature; the real cell at the quadrature points and so on. The Jacobian matrix is defined to be - $$ J_{ij} = {d\xi_i \over d\x_j} $$ + $$ J_{ij} = {d\xi_i \over dx_j} $$ which is the form needed to compute the gradient on the real cell from the gradient on the unit cell. If we want to transform the area element $dx dy$ from the real to the unit cell, we have to take the determinant of @@ -392,7 +392,7 @@ class FiniteElementBase { * should really be pure, but then we could * not make copies of a finite element * object even if we did not intend to use - * this function. Therefore, we ommit the + * this function. Therefore, we omit the * #=0# signature and implement this function * by throwing an exception. * #p# is a point on the reference element. @@ -406,7 +406,7 @@ class FiniteElementBase { * should really be pure, but then we could * not make copies of a finite element * object even if we did not intend to use - * this function. Therefore, we ommit the + * this function. Therefore, we omit the * #=0# signature and implement this function * by throwing an exception. * #p# is a point on the reference element, @@ -472,6 +472,10 @@ class FiniteElementBase { * higher dimensions, it depends on the * present fe and needs reimplementation * by the user. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, @@ -481,6 +485,21 @@ class FiniteElementBase { const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const; + + /** + * Return the ansatz points this FE has + * on a face if a cell would have the + * given face as a side. This function is + * needed for higher order elements, if + * we want to use curved boundary + * approximations. + * + * The function assumes that the fields + * already have the right number of + * elements. + */ + virtual void face_ansatz_points (const Triangulation::face_iterator &face, + vector > &ansatz_points) const; /** * Comparison operator. We also check for @@ -512,7 +531,13 @@ class FiniteElementBase { * Exception */ DeclException0 (ExcNotImplemented); - + /** + * Exception + */ + DeclException2 (ExcWrongFieldDimension, + int, int, + << "The field has not the assumed dimension " << arg2 + << ", but has " << arg1 << " elements."); protected: /** * Have #N=2^dim# matrices keeping the @@ -598,6 +623,12 @@ class FiniteElement; elements which are no more valid. Consequence: make sure the copy constructor is correct. + + This class should really be a pure one, with functions having the #=0# + signature. However, some instances of this class need to be floating around + anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values + of #dof_per*# available), so we do not make it pure but rather implement + those functions which should in fact be pure to throw an error. */ class FiniteElement<1> : public FiniteElementBase<1> { public: @@ -670,6 +701,10 @@ class FiniteElement<1> : public FiniteElementBase<1> { * distance to the origin. The standard * implementation distributes the dofs on * the line equidistantly. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell, const vector > &unit_points, @@ -679,6 +714,20 @@ class FiniteElement<1> : public FiniteElementBase<1> { const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const; + + /** + * Return the ansatz points this FE has + * on a face if a cell would have the + * given face as a side. This function is + * needed for higher order elements, if + * we want to use curved boundary + * approximations. + * + * Question: is this function useful in 1D? + * At present it is not implemented. + */ + virtual void face_ansatz_points (const Triangulation<1>::face_iterator &face, + vector > &ansatz_points) const; }; @@ -719,6 +768,12 @@ class FiniteElement<1> : public FiniteElementBase<1> { If you want to extend this class (not by derivation, but by adding new elements), see \Ref{FiniteElement<1>} + + This class should really be a pure one, with functions having the #=0# + signature. However, some instances of this class need to be floating around + anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values + of #dof_per*# available), so we do not make it pure but rather implement + those functions which should in fact be pure to throw an error. */ class FiniteElement<2> : public FiniteElementBase<2> { public: @@ -796,6 +851,10 @@ class FiniteElement<2> : public FiniteElementBase<2> { * function is therefore not implemented * by the FE<2> base class, but is made * pure virtual. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation<2>::cell_iterator &cell, const vector > &unit_points, @@ -805,6 +864,21 @@ class FiniteElement<2> : public FiniteElementBase<2> { const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const; + + /** + * Return the ansatz points this FE has + * on a face if a cell would have the + * given face as a side. This function is + * needed for higher order elements, if + * we want to use curved boundary + * approximations. + * + * The function assumes that the fields + * already have the right number of + * elements. + */ + virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face, + vector > &ansatz_points) const; }; diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index c6171f5268..e4da5ae11c 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -68,6 +68,10 @@ class FELinear : public FiniteElement { * function is therefore not implemented * by the FE<2> base class, but is made * pure virtual. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, @@ -77,6 +81,21 @@ class FELinear : public FiniteElement { const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const; + + /** + * Return the ansatz points this FE has + * on a face if a cell would have the + * given face as a side. Since we have no + * degrees of freedom on the faces for + * the linear ansatz, the ansatz points are + * simply the vertices of the face. + * + * The function assumes that the fields + * already have the right number of + * elements. + */ + virtual void face_ansatz_points (const Triangulation::face_iterator &face, + vector > &ansatz_points) const; }; @@ -133,6 +152,10 @@ class FEQuadratic : public FiniteElement { * function is therefore not implemented * by the FE<2> base class, but is made * pure virtual. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, @@ -198,6 +221,10 @@ class FECubic : public FiniteElement { * function is therefore not implemented * by the FE<2> base class, but is made * pure virtual. + * + * The function assumes that the fields + * already have the right number of + * elements. */ virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 7f0422227a..4a716bb56d 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -331,8 +331,9 @@ class ProblemBase { * See the general doc for more * information. */ - virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc, - map &boundary_values) const; + virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc, + const FiniteElement &fe, + map &boundary_values) const; /** * Exception @@ -406,7 +407,8 @@ class ProblemBase { void apply_dirichlet_bc (dSMatrix &matrix, dVector &solution, dVector &right_hand_side, - const DirichletBC &dirichlet_bc); + const DirichletBC &dirichlet_bc, + const FiniteElement &fe); friend class Assembler; }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 8543071a75..c4907af43f 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -37,6 +37,7 @@ FEValues::FEValues (const FiniteElement &fe, JxW_values(quadrature.n_quadrature_points, 0), quadrature_points(quadrature.n_quadrature_points, Point()), unit_quadrature_points(quadrature.n_quadrature_points, Point()), + ansatz_points (fe.total_dofs, Point()), jacobi_matrices (quadrature.n_quadrature_points, dFMatrix(dim,dim)), update_flags (update_flags) @@ -253,6 +254,14 @@ void FiniteElementBase::fill_fe_values (const typename Triangulation:: +template +void FiniteElementBase::face_ansatz_points (const typename Triangulation::face_iterator &, + vector > &) const { + Assert (false, ExcPureFunctionCalled()); +}; + + + /*------------------------------- FiniteElement ----------------------*/ @@ -272,6 +281,14 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + + // local mesh width const double h=(cell->vertex(1)(0) - cell->vertex(0)(0)); @@ -306,6 +323,13 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce +void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &, + vector > &) const { + // is this function useful in 1D? + Assert (false, ExcPureFunctionCalled()); +}; + + bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const { return ((dofs_per_vertex == f.dofs_per_vertex) && (dofs_per_line == f.dofs_per_line) && @@ -328,6 +352,14 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, +void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &, + vector > &) const { + // is this function useful in 1D? + Assert (false, ExcPureFunctionCalled()); +}; + + + /*------------------------------- Explicit Instantiations -------------*/ diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 93cb12965e..0604e9e7a2 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -5,6 +5,8 @@ #include + + FELinear<1>::FELinear () : FiniteElement<1> (1, 0) { @@ -93,6 +95,12 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, +void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &, + vector > &) const { + Assert (false, ExcPureFunctionCalled()); +}; + + FELinear<2>::FELinear () : @@ -239,8 +247,15 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, const bool compute_ansatz_points, vector > &q_points, const bool compute_q_points) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + const unsigned int dim=2; - const unsigned int n_vertices=4; + const unsigned int n_vertices=(1< vertices[n_vertices]; @@ -318,6 +333,20 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, +template +void FELinear::face_ansatz_points (const typename Triangulation::face_iterator &face, + vector > &ansatz_points) const { + Assert (ansatz_points.size() == (1<<(dim-1)), + typename FiniteElementBase::ExcWrongFieldDimension (ansatz_points.size(), + 1<<(dim-1))); + + for (unsigned int vertex=0; vertex<(1<<(dim-1)); ++vertex) + ansatz_points[vertex] = face->vertex(vertex); +}; + + + + @@ -383,13 +412,20 @@ FEQuadratic::shape_grad (const unsigned int i, void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, - const vector > &, - vector &, + const vector > &unit_points, + vector &jacobians, const bool, - vector > &, + vector > &ansatz_points, const bool, - vector > &, + vector > &q_points, const bool) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + Assert (false, typename FiniteElementBase<2>::ExcNotImplemented()); }; @@ -452,13 +488,20 @@ FECubic::shape_grad (const unsigned int i, void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, - const vector > &, - vector &, + const vector > &unit_points, + vector &jacobians, const bool, - vector > &, + vector > &ansatz_points, const bool, - vector > &, + vector > &q_points, const bool) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + Assert (false, typename FiniteElementBase<2>::ExcNotImplemented()); }; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index c3c3edbb94..7b555ba916 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -113,7 +113,8 @@ 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); + right_hand_side, + dirichlet_bc, fe); // condense system matrix in-place constraints.condense (system_matrix); @@ -307,11 +308,12 @@ template void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, dVector &solution, dVector &right_hand_side, - const DirichletBC &dirichlet_bc) { + const DirichletBC &dirichlet_bc, + const FiniteElement &fe) { 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 @@ -319,7 +321,7 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, // 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); + make_boundary_value_list (dirichlet_bc, fe, boundary_values); map::const_iterator dof, endd; const unsigned int n_dofs = (unsigned int)matrix.m(); @@ -346,10 +348,16 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, // zero, take the first main // diagonal entry we can find, or // one if no nonzero main diagonal - // element exists. + // 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) - right_hand_side((*dof).first) = (*dof).second * - matrix.diag_element((*dof).first); + new_rhs = right_hand_side((*dof).first) + = (*dof).second * matrix.diag_element((*dof).first); else { double first_diagonal_entry = 1; @@ -362,19 +370,19 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, matrix.set((*dof).first, (*dof).first, first_diagonal_entry); - right_hand_side((*dof).first) = (*dof).second * first_diagonal_entry; + new_rhs = right_hand_side((*dof).first) + = (*dof).second * first_diagonal_entry; }; - + // store the only nonzero entry // of this line for the Gauss // elimination step const double diagonal_entry = matrix.diag_element((*dof).first); - // do the Gauss step - for (unsigned int row=0; row::apply_dirichlet_bc (dSMatrix &matrix, { // correct right hand side right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry * - (*dof).second; + new_rhs; // set matrix entry to zero matrix.global_entry(j) = 0.; @@ -402,6 +410,7 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, void ProblemBase<1>::make_boundary_value_list (const DirichletBC &, + const FiniteElement<1> &, map &) const { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); Assert (false, ExcNotImplemented()); @@ -413,13 +422,28 @@ ProblemBase<1>::make_boundary_value_list (const DirichletBC &, template void ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, + const FiniteElement &fe, 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(); + Triangulation::active_face_iterator tface = tria->begin_active_face(); + DirichletBC::const_iterator function_ptr; - for (; face!=endf; ++face) + + // field to store the indices of dofs + // initialize once to get the size right + // for the following fields. + vector face_dofs; + face->get_dof_indices (face_dofs); + vector > dof_locations (face_dofs.size(), Point()); + vector dof_values; + + for (; face!=endf; ++face, ++tface) if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) != dirichlet_bc.end()) // face is subject to one of the @@ -428,17 +452,10 @@ ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, // get indices, physical location and // boundary values of dofs on this // face - vector face_dofs; - vector > dof_locations; - vector dof_values; - + face_dofs.erase (face_dofs.begin(), face_dofs.end()); + dof_values.erase (dof_values.begin(), dof_values.end()); face->get_dof_indices (face_dofs); - -//physical location -// Assert (false, ExcNotImplemented()); - dof_locations.insert (dof_locations.begin(), - face_dofs.size(), Point()); - + fe.face_ansatz_points (tface, dof_locations); (*function_ptr).second->value_list (dof_locations, dof_values); // enter into list