From: Wolfgang Bangerth Date: Wed, 4 Jan 2006 15:09:17 +0000 (+0000) Subject: Allow to also select which cells to work on by using the material id, not only the... X-Git-Tag: v8.0.0~12712 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=39c70cc080048a509c336310e43cbdbd9ef72647;p=dealii.git Allow to also select which cells to work on by using the material id, not only the subdomain id. git-svn-id: https://svn.dealii.org/trunk@11947 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 5ca9f0be90..2c12d93ab2 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -290,6 +290,13 @@ class KellyErrorEstimator * needs to access DoF values on * neighboring cells as well, even if * they belong to a different subdomain. + * + * The @p material_id parameter has a + * similar meaning: if not set to its + * default value, it means that + * indicators will only be computed for + * cells with this particular material + * id. */ template static void estimate (const Mapping &mapping, @@ -301,8 +308,9 @@ class KellyErrorEstimator const std::vector &component_mask = std::vector(), const Function *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); - + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); + /** * Calls the @p estimate * function, see above, with @@ -317,7 +325,8 @@ class KellyErrorEstimator const std::vector &component_mask = std::vector(), const Function *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** * Same function as above, but @@ -356,7 +365,8 @@ class KellyErrorEstimator const std::vector &component_mask = std::vector(), const Function *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** * Calls the @p estimate @@ -372,7 +382,8 @@ class KellyErrorEstimator const std::vector &component_mask = std::vector(), const Function *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** @@ -543,10 +554,15 @@ class KellyErrorEstimator std::vector JxW_values; /** - * The subdomain ids we are to care + * The subdomain id we are to care * for. */ const unsigned int subdomain_id; + /** + * The material id we are to care + * for. + */ + const unsigned int material_id; /** * Constructor. @@ -554,7 +570,8 @@ class KellyErrorEstimator PerThreadData (const unsigned int n_solution_vectors, const unsigned int n_components, const unsigned int n_q_points, - const unsigned int subdomain_id); + const unsigned int subdomain_id, + const unsigned int material_id); }; @@ -699,23 +716,20 @@ class KellyErrorEstimator<1> * the respective component in * the coefficient. * - * You might give a list of - * components you want to - * evaluate, in case the finite - * element used by the - * DoFHandler object is - * vector-valued. You then have - * to set those entries to true - * in the bit-vector - * @p component_mask for which the - * respective component is to be - * used in the error - * estimator. The default is to - * use all components, which is - * done by either providing a - * bit-vector with all-set - * entries, or an empty - * bit-vector. + * You might give a list of components + * you want to evaluate, in case the + * finite element used by the DoFHandler + * object is vector-valued. You then have + * to set those entries to true in the + * bit-vector @p component_mask for which + * the respective component is to be used + * in the error estimator. The default is + * to use all components, which is done + * by either providing a bit-vector with + * all-set entries, or an empty + * bit-vector. All the other parameters + * are as in the general case used for 2d + * and higher. * * The estimator supports multithreading * and splits the cells to @@ -740,7 +754,8 @@ class KellyErrorEstimator<1> const std::vector &component_mask = std::vector(), const Function<1> *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** * Calls the @p estimate @@ -756,7 +771,8 @@ class KellyErrorEstimator<1> const std::vector &component_mask = std::vector(), const Function<1> *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** * Same function as above, but @@ -795,7 +811,8 @@ class KellyErrorEstimator<1> const std::vector &component_mask = std::vector(), const Function<1> *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** * Calls the @p estimate @@ -811,7 +828,8 @@ class KellyErrorEstimator<1> const std::vector &component_mask = std::vector(), const Function<1> *coefficients = 0, const unsigned int n_threads = multithread_info.n_default_threads, - const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int); + const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int, + const unsigned int material_id = deal_II_numbers::invalid_unsigned_int); /** diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 05dbb1f59c..1b56a8111c 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -76,13 +76,14 @@ estimate (const Mapping<1> &mapping, const std::vector &component_mask, const Function<1> *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { // just pass on to the other function const std::vector solutions (1, &solution); std::vector*> errors (1, &error); estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id); + component_mask, coefficients, n_threads, subdomain_id, material_id); } @@ -97,12 +98,13 @@ estimate (const DoFHandler<1> &dof_handler, const std::vector &component_mask, const Function<1> *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); static const MappingQ1<1> mapping; estimate(mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id); + error, component_mask, coefficients, n_threads, subdomain_id, material_id); } @@ -118,12 +120,13 @@ estimate (const DoFHandler<1> &dof_handler, const std::vector &component_mask, const Function<1> *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); static const MappingQ1<1> mapping; estimate(mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id); + errors, component_mask, coefficients, n_threads, subdomain_id, material_id); } @@ -139,7 +142,8 @@ estimate (const Mapping<1> &mapping, const std::vector &component_mask_, const Function<1> *coefficient, const unsigned int, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { const unsigned int n_components = dof_handler.get_fe().n_components(); const unsigned int n_solution_vectors = solutions.size(); @@ -232,9 +236,13 @@ estimate (const Mapping<1> &mapping, DoFHandler<1>::active_cell_iterator cell = dof_handler.begin_active(); for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index) - if ((subdomain_id == deal_II_numbers::invalid_unsigned_int) - || - (cell->subdomain_id() == subdomain_id)) + if (((subdomain_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->subdomain_id() == subdomain_id)) + && + ((material_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->material_id() == material_id))) { for (unsigned int n=0; n::PerThreadData:: PerThreadData (const unsigned int n_solution_vectors, const unsigned int n_components, const unsigned int n_q_points, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) : - subdomain_id (subdomain_id) + subdomain_id (subdomain_id), + material_id (material_id) { // Init the size of a lot of vectors // needed in the calculations once @@ -409,13 +419,14 @@ estimate (const Mapping &mapping, const std::vector &component_mask, const Function *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { // just pass on to the other function const std::vector solutions (1, &solution); std::vector*> errors (1, &error); estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id); + component_mask, coefficients, n_threads, subdomain_id, material_id); } @@ -431,12 +442,14 @@ estimate (const DoFHandler &dof_handler, const std::vector &component_mask, const Function *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); static const MappingQ1 mapping; estimate(mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id); + error, component_mask, coefficients, n_threads, + subdomain_id, material_id); } @@ -459,6 +472,7 @@ estimate_some (const Mapping &mapping, const unsigned int n_solution_vectors = solutions.size(); const unsigned int subdomain_id = per_thread_data.subdomain_id; + const unsigned int material_id = per_thread_data.material_id; // make up a fe face values object // for the restriction of the @@ -589,15 +603,20 @@ estimate_some (const Mapping &mapping, continue; } - // finally: note that we only - // have to do something if either - // the present cell is on the - // subdomain we care for, or if one - // of the neighbors behind the face - // is on the subdomain we care for - if ( ! ((subdomain_id == deal_II_numbers::invalid_unsigned_int) - || - (cell->subdomain_id() == subdomain_id))) + // finally: note that we only have + // to do something if either the + // present cell is on the subdomain + // we care for (and the same for + // material_id), or if one of the + // neighbors behind the face is on + // the subdomain we care for + if ( ! ( ((subdomain_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->subdomain_id() == subdomain_id)) + && + ((material_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->material_id() == material_id))) ) { // ok, cell is unwanted, but // maybe its neighbor behind @@ -609,14 +628,28 @@ estimate_some (const Mapping &mapping, bool care_for_cell = false; if (cell->neighbor(face_no)->has_children() == false) - care_for_cell |= (cell->neighbor(face_no)->subdomain_id() - == subdomain_id); + care_for_cell |= ((cell->neighbor(face_no)->subdomain_id() + == subdomain_id) || + (subdomain_id == deal_II_numbers::invalid_unsigned_int)) + && + ((cell->neighbor(face_no)->material_id() + == material_id) || + (material_id == deal_II_numbers::invalid_unsigned_int)); else { for (unsigned int sf=0; sf::subfaces_per_face; ++sf) - if (cell->neighbor_child_on_subface(face_no,sf) - ->subdomain_id() == subdomain_id) + if (((cell->neighbor_child_on_subface(face_no,sf) + ->subdomain_id() == subdomain_id) + && + (material_id == + deal_II_numbers::invalid_unsigned_int)) + || + ((cell->neighbor_child_on_subface(face_no,sf) + ->material_id() == material_id) + && + (subdomain_id == + deal_II_numbers::invalid_unsigned_int))) { care_for_cell = true; break; @@ -624,8 +657,9 @@ estimate_some (const Mapping &mapping, } // so if none of the neighbors - // cares for this subdomain - // either, then try next face + // cares for this subdomain or + // material either, then try + // next face if (care_for_cell == false) continue; } @@ -684,7 +718,8 @@ estimate (const Mapping &mapping, const std::vector &component_mask_, const Function *coefficients, const unsigned int n_threads_, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { const unsigned int n_components = dof_handler.get_fe().n_components(); @@ -774,7 +809,7 @@ estimate (const Mapping &mapping, = new PerThreadData (solutions.size(), dof_handler.get_fe().n_components(), quadrature.n_quadrature_points, - subdomain_id); + subdomain_id, material_id); // split all cells into threads if // multithreading is used and run @@ -837,9 +872,13 @@ estimate (const Mapping &mapping, for (active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell, ++present_cell) - if ((subdomain_id == deal_II_numbers::invalid_unsigned_int) - || - (cell->subdomain_id() == subdomain_id)) + if ( ((subdomain_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->subdomain_id() == subdomain_id)) + && + ((material_id == deal_II_numbers::invalid_unsigned_int) + || + (cell->material_id() == material_id))) { // loop over all faces of this cell for (unsigned int face_no=0; face_no::faces_per_cell; @@ -879,12 +918,14 @@ void KellyErrorEstimator::estimate (const DoFHandler &do const std::vector &component_mask, const Function *coefficients, const unsigned int n_threads, - const unsigned int subdomain_id) + const unsigned int subdomain_id, + const unsigned int material_id) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); static const MappingQ1 mapping; estimate(mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id); + errors, component_mask, coefficients, n_threads, + subdomain_id, material_id); } @@ -1295,6 +1336,7 @@ estimate (const Mapping &, \ const std::vector &, \ const Function *, \ const unsigned int , \ + const unsigned int , \ const unsigned int); \ \ template \ @@ -1308,6 +1350,7 @@ estimate (const DoFHandler &, \ const std::vector &, \ const Function *, \ const unsigned int , \ + const unsigned int , \ const unsigned int); \ \ template \ @@ -1322,6 +1365,7 @@ estimate (const Mapping &, \ const std::vector &, \ const Function *, \ const unsigned int , \ + const unsigned int , \ const unsigned int); \ \ template \ @@ -1335,6 +1379,7 @@ estimate (const DoFHandler &, \ const std::vector &, \ const Function *, \ const unsigned int , \ + const unsigned int , \ const unsigned int) INSTANTIATE(Vector); diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index fc91f2879b..779ed6ff1a 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -258,6 +258,14 @@ inconvenience this causes.

deal.II

    +
  1. + Improved: The KellyErrorEstimator class now also allows to + select the cells on which it is supposed to work by giving a material + ID, instead of only a subdomain ID. +
    + (WB 2006/1/4) +

    +
  2. Improved: A new TriaAccessor::ExcCellHasNoChildren exception will be raised if the +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +class MySquareFunction : public Function +{ + public: + MySquareFunction () : Function(2) {}; + + virtual double value (const Point &p, + const unsigned int component) const + { return (component+1)*p.square(); }; + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); + values(1) = value(p,1); }; +}; + + + +template +Quadrature & +get_q_face (Function&) +{ + static QGauss4 q; + return q; +} + + +Quadrature<0> & +get_q_face (Function<1>&) +{ + Quadrature<0> *q = 0; + return *q; +} + + + +template +void make_mesh (Triangulation &tria) +{ + + GridGenerator::hyper_cube(tria, -1, 1); + + // refine the mesh in a random way so as to + // generate as many cells with + // hanging nodes as possible + tria.refine_global (4-dim); + const double steps[4] = { /*d=0*/ 0, 7, 3, 3 }; + for (unsigned int i=0; i::active_cell_iterator + cell = tria.begin_active(); + for (unsigned int index=0; cell != tria.end(); ++cell, ++index) + if (index % (3*dim) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement (); + } + + // we now have a number of cells, + // flag them with some material + // ids based on their position, in + // particular we take the quadrant + // (octant) + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end (); + for (; cell!=endc; ++cell) + { + unsigned int material = 0; + for (unsigned int d=0; dcenter()(d) > 0) + material |= (1<set_material_id (material); + } +} + + + + +template +void +check () +{ + Functions::CosineFunction function; + + Triangulation tria; + make_mesh (tria); + + FE_Q element(3); + DoFHandler dof(tria); + dof.distribute_dofs(element); + + MappingQ mapping(3); + Quadrature &q_face = get_q_face(function); + + std::map*> neumann_bc; + neumann_bc[0] = &function; + + Vector v (dof.n_dofs()); + VectorTools::interpolate (mapping, dof, function, v); + + Vector error1 (tria.n_active_cells()); + Vector error2 (tria.n_active_cells()); + + // compute error by looking at all cells at + // once and output this as the base line + // results. scale results so that they show + // up in the output file as a reasonable + // number + KellyErrorEstimator::estimate (mapping, dof, q_face, neumann_bc, + v, error1); + const double scaling_factor = 10000.*error1.size()/error1.l1_norm(); + error1 *= scaling_factor; + + deallog << "Estimated error indicators:" << std::endl; + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +}