From 20efffd1f02c9d68b555a5eb23dea1a7750b5d7b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 27 Apr 1998 18:04:29 +0000 Subject: [PATCH] Add first version of simple error estimator. git-svn-id: https://svn.dealii.org/trunk@206 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/numerics/base.cc | 10 +- .../source/numerics/error_estimator.cc | 222 +++++++++++++++++- 2 files changed, 222 insertions(+), 10 deletions(-) diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 8b054662c7..cf4fad227e 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -78,7 +78,7 @@ void ProblemBase::assemble (const Equation &equation, const Quadrature &quadrature, const FiniteElement &fe, const UpdateFlags update_flags, - const DirichletBC &dirichlet_bc, + const FunctionMap &dirichlet_bc, const Boundary &boundary) { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); @@ -369,7 +369,7 @@ template void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, dVector &solution, dVector &right_hand_side, - const DirichletBC &dirichlet_bc, + const FunctionMap &dirichlet_bc, const FiniteElement &fe, const Boundary &boundary) { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); @@ -471,7 +471,7 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, void -ProblemBase<1>::make_boundary_value_list (const DirichletBC &, +ProblemBase<1>::make_boundary_value_list (const FunctionMap &, const FiniteElement<1> &, const Boundary<1> &, map &) const { @@ -484,7 +484,7 @@ ProblemBase<1>::make_boundary_value_list (const DirichletBC &, template void -ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, +ProblemBase::make_boundary_value_list (const FunctionMap &dirichlet_bc, const FiniteElement &fe, const Boundary &boundary, map &boundary_values) const { @@ -496,7 +496,7 @@ ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_ DoFHandler::active_face_iterator face = dof_handler->begin_active_face(), endf = dof_handler->end_face(); - DirichletBC::const_iterator function_ptr; + FunctionMap::const_iterator function_ptr; // field to store the indices of dofs // initialize once to get the size right diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 3e2cc7ae48..26250d2fe0 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -1,13 +1,32 @@ /* $Id$ */ -#include +#include +#include +#include +#include +#include #include +#include +#include #include +#include +#include + + + +inline double sqr (const double x) { + return x*x; +}; + void KellyErrorEstimator<1>::estimate_error (const DoFHandler<1> &, + const Quadrature<0> &, + const FiniteElement<1> &, + const Boundary<1> &, + const FunctionMap &, const dVector &, dVector &) const { Assert(false, ExcNotImplemented()); @@ -16,9 +35,202 @@ void KellyErrorEstimator<1>::estimate_error (const DoFHandler<1> &, template -void KellyErrorEstimator::estimate_error (const DoFHandler &dof, - const dVector &solution, - dVector &error) const { - ; +void KellyErrorEstimator::estimate_error (const DoFHandler &dof, + const Quadrature &quadrature, + const FiniteElement &fe, + const Boundary &boundary, + const FunctionMap &neumann_bc, + const dVector &solution, + dVector &error) const { + Assert (neumann_bc.find(255) == neumann_bc.end(), + ExcInvalidBoundaryIndicator()); + + // reserve one slot for each cell and set + // it to zero + error.reinit (dof.get_tria().n_active_cells()); + + // number of integration points per face + const unsigned int n_q_points = quadrature.n_quadrature_points; + // number of dofs per cell + const unsigned int n_dofs = fe.total_dofs; + + // make up a fe face values object for the + // restriction of the finite element function + // to a face, for the present cell and its + // neighbor. + FEFaceValues fe_face_values_cell (fe, quadrature, + UpdateFlags(update_gradients | update_JxW_values | + update_jacobians | update_q_points | + update_normal_vectors)); + FEFaceValues fe_face_values_neighbor (fe, quadrature, + UpdateFlags(update_gradients)); + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // loop over all cells + for (unsigned int present_cell=0; cell!=endc; ++cell, ++present_cell) + // loop over all faces of this cell + for (unsigned int face_no=0; face_no<2*dim; ++face_no) + { + const unsigned char boundary_indicator = cell->face(face_no)->boundary_indicator(); + if ((boundary_indicator != 255) && + neumann_bc.find(boundary_indicator)==neumann_bc.end()) + // this face is part of the boundary + // but not of the neumann boundary + // -> nothing to do + continue; + + + // initialize data of the restriction + // of this cell to the present face + fe_face_values_cell.reinit (cell, face_no, fe, boundary); + + // set up a vector of the gradients + // of the finite element function + // on this cell at the quadrature + // points + // + // let psi be a short name for + // [grad u_h] + vector > psi(n_q_points); + + // get a list of the values of + // the solution for the ansatz + // functions on this cell + vector dof_values; + cell->get_dof_values (solution, dof_values); + + // get a list of the gradients of + // the ansatz functions on this + // cell at the quadrature points + const vector > > &shape_grads(fe_face_values_cell.get_shape_grads()); + + // compute the gradients of the solution + // at the quadrature points by summing + // over the ansatz functions. + for (unsigned int j=0; jneighbor(face_no).state() == valid, + ExcInternalError()); + unsigned int neighbor_neighbor; + DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); + + // find which number the current + // face has relative to the neighboring + // cell + for (neighbor_neighbor=0; neighbor_neighbor<2*dim; ++neighbor_neighbor) + if (neighbor->neighbor(neighbor_neighbor) == cell) + break; + + Assert (neighbor_neighborget_dof_values (solution, dof_values); + // get a list of the gradients of the + // + const vector > > &neighbor_grads (fe_face_values_cell. + get_shape_grads()); + // subtract the gradients of the + // solution on the neigbor cell + // at the quadrature points from + // those of the present cell + for (unsigned int j=0; j phi(n_q_points,0); + const vector > &normal_vectors(fe_face_values_cell.get_normal_vectors()); + for (unsigned int point=0; point g(n_q_points); + neumann_bc.find(boundary_indicator)->second + ->value_list (fe_face_values_cell.get_quadrature_points(), + g); + + for (unsigned int point=0; point; +template class KellyErrorEstimator<2>; -- 2.39.5