From: Wolfgang Bangerth Date: Tue, 5 May 1998 14:15:59 +0000 (+0000) Subject: Cleanup and code restructuring. X-Git-Tag: v8.0.0~23000 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3feb99fafe7f9f4ac93bbba1ffab565f29a48eb4;p=dealii.git Cleanup and code restructuring. git-svn-id: https://svn.dealii.org/trunk@254 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index a4a04e6d7a..f1af6816fb 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -156,6 +156,52 @@ class KellyErrorEstimator { * Exception */ DeclException0 (ExcInvalidBoundaryIndicator); + + private: + /** + * Declare a data type to represent the + * mapping between faces and integrated + * jumps of gradients. See the general + * doc of this class for more information. + */ + typedef map::face_iterator,double> FaceIntegrals; + + /** + * Redeclare an active cell iterator. + * This is simply for convenience. + */ + typedef DoFHandler::active_cell_iterator active_cell_iterator; + + /** + * Actually do the computation on a face + * which has no hanging nodes (it is + * regular), i.e. + * either on the other side there is + * nirvana (face is at boundary), or + * the other side's refinement level + * is the same as that of this side, + * then handle the integration of + * these both cases together. + * + * The meaning of the parameters becomes + * clear when looking at the source + * code. This function is only + * externalized from #estimate_error# + * to avoid ending up with a function + * of 500 lines of code. + */ + static void integrate_over_regular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const FiniteElement &fe, + const Boundary &boundary, + const FunctionMap &neumann_bc, + const unsigned int n_q_points, + FEFaceValues &fe_face_values_cell, + FEFaceValues &fe_face_values_neighbor, + FaceIntegrals &face_integrals, + const dVector &solution); + + static void integrate_over_irregular_face (); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index d041c84ded..f50b209fa6 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -57,7 +57,7 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, // loop over the cells and collect the // conrtibutions of the different faces // of the cell. - map::face_iterator, double> face_integrals; + FaceIntegrals face_integrals; // number of integration points per face const unsigned int n_q_points = quadrature.n_quadrature_points; @@ -65,22 +65,25 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, // 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. + // neighbor. In principle we would only need + // one at a time, but this way we can + // have more fine grained access to what + // values really need to be computed (we + // need not compute all values on the + // neighbor cells, so using two objects + // gives us a performance gain). 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 | update_jacobians)); - // loop variables - DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); // loop over all cells - for (; cell!=endc; ++cell) + const active_cell_iterator endc = dof.end(); + for (active_cell_iterator cell = dof.begin_active(); cell!=endc; ++cell) // loop over all faces of this cell for (unsigned int face_no=0; face_no<2*dim; ++face_no) { @@ -115,129 +118,28 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, 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); - fe_face_values_cell.get_function_grads (solution, psi); - - - // now compute over the other side of - // the face - if (boundary_indicator == 255) - // internal face; integrate jump - // of gradient across this face - { - Assert (cell->neighbor(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_neighbor > neighbor_psi (n_q_points); - fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); - - // compute the jump in the gradients - transform (psi.begin(), psi.end(), - neighbor_psi.begin(), - psi.begin(), - minus >()); - }; - - - - - // now psi contains the following: - // - for an internal face, psi=[grad u] - // - for a neumann boundary face, - // psi=grad u - // each component being the - // mentioned value at one of the - // quadrature points - - // next we have to multiply this with - // the normal vector. Since we have - // taken the difference of gradients - // for internal faces, we may chose - // the normal vector of one cell, - // taking that of the neighbor - // would only change the sign. We take - // the outward normal. - // - // let phi be the name of the integrand - vector 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; pointface(face_no)] - = (inner_product (phi.begin(), phi.end(), - fe_face_values_cell.get_JxW_values().begin(), - 0.0) * - cell->diameter() / 24); + if (cell->face(face_no)->has_children() == false) + // if the face is a regular one, i.e. + // either on the other side there is + // nirvana (face is at boundary), or + // the other side's refinement level + // is the same as that of this side, + // then handle the integration of + // these both cases together + integrate_over_regular_face (cell, face_no, fe, + boundary, neumann_bc, + n_q_points, + fe_face_values_cell, + fe_face_values_neighbor, + face_integrals, + solution); + else + // otherwise we need to do some + // special computations which do + // not fit into the framework of + // the above function + integrate_over_irregular_face (); }; @@ -249,7 +151,7 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, error.reinit (dof.get_tria().n_active_cells()); unsigned int present_cell=0; - for (cell=dof.begin_active(); cell!=endc; ++cell, ++present_cell) + for (active_cell_iterator cell=dof.begin_active(); cell!=endc; ++cell, ++present_cell) { // loop over all faces of this cell for (unsigned int face_no=0; face_no<2*dim; ++face_no) @@ -268,6 +170,174 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, +void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iterator &, + const unsigned int , + const FiniteElement<1> &, + const Boundary<1> &, + const FunctionMap &, + const unsigned int , + FEFaceValues<1> &, + FEFaceValues<1> &, + FaceIntegrals &, + const dVector &) { + Assert (false, ExcInternalError()); +}; + + + + +template +void KellyErrorEstimator:: +integrate_over_regular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const FiniteElement &fe, + const Boundary &boundary, + const FunctionMap &neumann_bc, + const unsigned int n_q_points, + FEFaceValues &fe_face_values_cell, + FEFaceValues &fe_face_values_neighbor, + FaceIntegrals &face_integrals, + const dVector &solution) { + DoFHandler::face_iterator face = cell->face(face_no); + + // 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); + fe_face_values_cell.get_function_grads (solution, psi); + + + // now compute over the other side of + // the face + if (face->at_boundary() == false) + // internal face; integrate jump + // of gradient across this face + { + Assert (cell->neighbor(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_neighbor > neighbor_psi (n_q_points); + fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); + + // compute the jump in the gradients + transform (psi.begin(), psi.end(), + neighbor_psi.begin(), + psi.begin(), + minus >()); + }; + + + + + // now psi contains the following: + // - for an internal face, psi=[grad u] + // - for a neumann boundary face, + // psi=grad u + // each component being the + // mentioned value at one of the + // quadrature points + + // next we have to multiply this with + // the normal vector. Since we have + // taken the difference of gradients + // for internal faces, we may chose + // the normal vector of one cell, + // taking that of the neighbor + // would only change the sign. We take + // the outward normal. + // + // let phi be the name of the integrand + vector phi(n_q_points,0); + const vector > &normal_vectors(fe_face_values_cell. + get_normal_vectors()); + + for (unsigned int point=0; pointat_boundary() == true) + // neumann boundary face. compute + // difference between normal + // derivative and boundary function + { + const unsigned char boundary_indicator = face->boundary_indicator(); + + Assert (neumann_bc.find(boundary_indicator) != neumann_bc.end(), + ExcInternalError ()); + // get the values of the boundary + // function at the quadrature + // points + vector 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; pointdiameter() / 24); +}; + + + + +template +void KellyErrorEstimator:: +integrate_over_irregular_face () { + Assert (false, ExcInternalError()); +}; + + + // explicit instantiations template class KellyErrorEstimator<1>;