From: Wolfgang Bangerth Date: Wed, 6 May 1998 14:12:10 +0000 (+0000) Subject: Finish support for error estimator also on unstructured grids. X-Git-Tag: v8.0.0~22997 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bfd1b5d1dc8b7039a7a1a6f723f4779c91e9d1ee;p=dealii.git Finish support for error estimator also on unstructured grids. git-svn-id: https://svn.dealii.org/trunk@257 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 f1af6816fb..83dc2d8ec6 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -61,6 +61,20 @@ class dVector; sum up the contributions of the faces (which are the integrated square of the jumps) of each cell and take the square root. + We store the contribution of each face in a #map#, as provided by the + C++ standard library, with the iterator pointing to that face being the + key into the map. In fact, we do not store the indicator per face, but + only the integral listed above. When looping the second time over all + cells, we have to sum up the contributions of the faces, multiply them + with $\frac h{24}$ and take the square root. By doing the multiplication + with $h$ in the second loop, we avoid problems to decide with which $h$ + to multiply, that of the cell on the one or that of the cell on the other + side of the face. + + $h$ is taken to be the greatest length of the diagonals of the cell. For + more or less uniform cells without deformed angles, this coincides with + the diameter of the cell. + {\bf Boundary values} @@ -201,7 +215,23 @@ class KellyErrorEstimator { FaceIntegrals &face_integrals, const dVector &solution); - static void integrate_over_irregular_face (); + /** + * The same applies as for the function + * above, except that integration is + * over face #face_no# of #cell#, where + * the respective neighbor is refined, + * so that the integration is a bit more + * complex. + */ + static void integrate_over_irregular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const FiniteElement &fe, + const Boundary &boundary, + const unsigned int n_q_points, + FEFaceValues &fe_face_values, + FESubfaceValues &fe_subface_values, + FaceIntegrals &face_integrals, + const dVector &solution); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index f50b209fa6..e9d7835c3a 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -80,12 +81,15 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, FEFaceValues fe_face_values_neighbor (fe, quadrature, UpdateFlags(update_gradients | update_jacobians)); - + FESubfaceValues fe_subface_values (fe, quadrature, + UpdateFlags(update_gradients | + update_jacobians)); + // loop over all cells 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) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { // if we already visited this // face: do nothing @@ -139,7 +143,11 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, // special computations which do // not fit into the framework of // the above function - integrate_over_irregular_face (); + integrate_over_irregular_face (cell, face_no, fe, boundary, + n_q_points, + fe_face_values_cell, + fe_subface_values, + face_integrals, solution); }; @@ -154,12 +162,13 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, 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) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { Assert (face_integrals.find(cell->face(face_no)) != face_integrals.end(), ExcInternalError()); - error(present_cell) += face_integrals[cell->face(face_no)]; + error(present_cell) += (face_integrals[cell->face(face_no)] * + cell->diameter() / 24); }; error(present_cell) = sqrt(error(present_cell)); @@ -198,7 +207,7 @@ integrate_over_regular_face (const active_cell_iterator &cell, FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, const dVector &solution) { - DoFHandler::face_iterator face = cell->face(face_no); + const DoFHandler::face_iterator face = cell->face(face_no); // initialize data of the restriction // of this cell to the present face @@ -230,18 +239,19 @@ integrate_over_regular_face (const active_cell_iterator &cell, // find which number the current // face has relative to the neighboring // cell - for (neighbor_neighbor=0; neighbor_neighbor<2*dim; ++neighbor_neighbor) + for (neighbor_neighbor=0; neighbor_neighbor::faces_per_cell; + ++neighbor_neighbor) if (neighbor->neighbor(neighbor_neighbor) == cell) break; - Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); // get restriction of finite element // function of #neighbor# to the // common face. fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor, fe, boundary); - + // get a list of the gradients of // the finite element solution // restricted to the neighbor cell @@ -321,19 +331,144 @@ integrate_over_regular_face (const active_cell_iterator &cell, // perform integration by multiplication // with weights and summation. - face_integrals[face] = (inner_product (phi.begin(), phi.end(), - fe_face_values_cell.get_JxW_values().begin(), - 0.0) * - cell->diameter() / 24); + face_integrals[face] = inner_product (phi.begin(), phi.end(), + fe_face_values_cell.get_JxW_values().begin(), + 0.0); }; +void KellyErrorEstimator<1>:: +integrate_over_irregular_face (const active_cell_iterator &, + const unsigned int , + const FiniteElement<1> &, + const Boundary<1> &, + const unsigned int , + FEFaceValues<1> &, + FESubfaceValues<1> &, + FaceIntegrals &, + const dVector &) { + Assert (false, ExcInternalError()); +}; + + + template void KellyErrorEstimator:: -integrate_over_irregular_face () { - Assert (false, ExcInternalError()); +integrate_over_irregular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const FiniteElement &fe, + const Boundary &boundary, + const unsigned int n_q_points, + FEFaceValues &fe_face_values, + FESubfaceValues &fe_subface_values, + FaceIntegrals &face_integrals, + const dVector &solution) { + const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); + Assert (neighbor.state() == valid, ExcInternalError()); + Assert (neighbor->has_children(), ExcInternalError()); + Assert (neighbor->child(0)->has_children()==false, + ExcInternalError()); + // 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); + + // store which number #cell# has in the + // list of neighbors of #neighbor# + unsigned int neighbor_neighbor; + for (neighbor_neighbor=0; neighbor_neighbor::faces_per_cell; + ++neighbor_neighbor) + if (neighbor->neighbor(neighbor_neighbor) == cell) + break; + Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); + + // loop over all subfaces + for (unsigned int subface_no=0; subface_no::subfaces_per_face; + ++subface_no) + { + // get an iterator pointing to the + // cell behind the present subface + const active_cell_iterator neighbor_child + = neighbor->child(GeometryInfo:: + child_cell_on_face(neighbor_neighbor,subface_no)); + Assert (neighbor_child->face(neighbor_neighbor) == + cell->face(face_no)->child(subface_no), + ExcInternalError()); + + // restrict the finite element on the + // present cell to the subface and + // store the gradient of the solution + // in psi + fe_subface_values.reinit (cell, face_no, subface_no, + fe, boundary); + fe_subface_values.get_function_grads (solution, psi); + + // restrict the finite element on the + // neighbor cell to the common #subface#. + // store the gradient in #neighbor_psi# + vector > neighbor_psi (n_q_points); + fe_face_values.reinit (neighbor_child, neighbor_neighbor, + fe, boundary); + fe_face_values.get_function_grads (solution, neighbor_psi); + + // compute the jump in the gradients + transform (psi.begin(), psi.end(), + neighbor_psi.begin(), + psi.begin(), + minus >()); + + // 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. + get_normal_vectors()); + + for (unsigned int point=0; pointface(neighbor_neighbor)] + = inner_product (phi.begin(), phi.end(), + fe_face_values.get_JxW_values(). + begin(), + 0.0); + }; + + + // finally loop over all subfaces to + // collect the contributions of the + // subfaces and store them with the + // mother face + double sum=0; + DoFHandler::face_iterator face = cell->face(face_no); + for (unsigned int subface_no=0; subface_no::subfaces_per_face; + ++subface_no) + { + Assert (face_integrals.find(face->child(subface_no)) != + face_integrals.end(), + ExcInternalError()); + sum += face_integrals[face->child(subface_no)]; + }; + face_integrals[face] = sum; };