// $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
* 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 <typename InputVector>
static void estimate (const Mapping<dim> &mapping,
const std::vector<bool> &component_mask = std::vector<bool>(),
const Function<dim> *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
const std::vector<bool> &component_mask = std::vector<bool>(),
const Function<dim> *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
const std::vector<bool> &component_mask = std::vector<bool>(),
const Function<dim> *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
const std::vector<bool> &component_mask = std::vector<bool>(),
const Function<dim> *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);
/**
std::vector<double> 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.
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);
};
* 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
const std::vector<bool> &component_mask = std::vector<bool>(),
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
const std::vector<bool> &component_mask = std::vector<bool>(),
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
const std::vector<bool> &component_mask = std::vector<bool>(),
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
const std::vector<bool> &component_mask = std::vector<bool>(),
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);
/**
// $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
const std::vector<bool> &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<const InputVector *> solutions (1, &solution);
std::vector<Vector<float>*> 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);
}
const std::vector<bool> &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);
}
const std::vector<bool> &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);
}
const std::vector<bool> &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();
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<n_solution_vectors; ++n)
(*errors[n])(cell_index) = 0;
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
const std::vector<bool> &component_mask,
const Function<dim> *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<const InputVector *> solutions (1, &solution);
std::vector<Vector<float>*> 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);
}
const std::vector<bool> &component_mask,
const Function<dim> *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<dim> 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);
}
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
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
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<GeometryInfo<dim>::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;
}
// 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;
}
const std::vector<bool> &component_mask_,
const Function<dim> *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();
= 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
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<GeometryInfo<dim>::faces_per_cell;
const std::vector<bool> &component_mask,
const Function<dim> *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<dim> 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);
}
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
const unsigned int , \
+ const unsigned int , \
const unsigned int); \
\
template \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
const unsigned int , \
+ const unsigned int , \
const unsigned int); \
\
template \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
const unsigned int , \
+ const unsigned int , \
const unsigned int); \
\
template \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
const unsigned int , \
+ const unsigned int , \
const unsigned int)
INSTANTIATE(Vector<double>);
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ Improved: The <code>KellyErrorEstimator</code> 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.
+ <br>
+ (WB 2006/1/4)
+ </p>
+
<li> <p>
Improved: A new <code>TriaAccessor::ExcCellHasNoChildren</code>
exception will be raised if the <code
--- /dev/null
+//---------------------------- error_estimator_02.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- error_estimator_02.cc ---------------------------
+
+// same as error_estimator_02, but check for material_id instead of
+// subdomain_id
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q.h>
+#include <numerics/vectors.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+ public:
+ MySquareFunction () : Function<dim>(2) {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const
+ { return (component+1)*p.square(); };
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ { values(0) = value(p,0);
+ values(1) = value(p,1); };
+};
+
+
+
+template <int dim>
+Quadrature<dim-1> &
+get_q_face (Function<dim>&)
+{
+ static QGauss4<dim-1> q;
+ return q;
+}
+
+
+Quadrature<0> &
+get_q_face (Function<1>&)
+{
+ Quadrature<0> *q = 0;
+ return *q;
+}
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &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<steps[dim]; ++i)
+ {
+ typename Triangulation<dim>::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<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end ();
+ for (; cell!=endc; ++cell)
+ {
+ unsigned int material = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ if (cell->center()(d) > 0)
+ material |= (1<<d);
+ Assert (material < (1<<dim), ExcInternalError());
+
+ cell->set_material_id (material);
+ }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Functions::CosineFunction<dim> function;
+
+ Triangulation<dim> tria;
+ make_mesh (tria);
+
+ FE_Q<dim> element(3);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(element);
+
+ MappingQ<dim> mapping(3);
+ Quadrature<dim-1> &q_face = get_q_face(function);
+
+ std::map<unsigned char,const Function<dim>*> neumann_bc;
+ neumann_bc[0] = &function;
+
+ Vector<double> v (dof.n_dofs());
+ VectorTools::interpolate (mapping, dof, function, v);
+
+ Vector<float> error1 (tria.n_active_cells());
+ Vector<float> 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<dim>::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<error1.size(); ++i)
+ deallog << error1(i) << std::endl;
+
+ // then do the same with different
+ // material ids and add up the result
+ for (unsigned int material=0; material<(1<<dim); ++material)
+ {
+ deallog << "Material id=" << material << std::endl;
+
+ Vector<float> this_error (tria.n_active_cells());
+ KellyErrorEstimator<dim>::estimate (mapping, dof, q_face, neumann_bc,
+ v, this_error,
+ std::vector<bool>(), 0,
+ multithread_info.n_default_threads,
+ deal_II_numbers::invalid_unsigned_int,
+ material);
+ this_error *= scaling_factor;
+
+ // copy the result into error2. since
+ // every invokation of the kelly
+ // estimator should only operate on
+ // the cells of one material,
+ // whenever there is something in
+ // this_error, the corresponding
+ // entry in error2 should still be
+ // empty
+ for (unsigned int i=0; i<this_error.size(); ++i)
+ {
+ deallog << i << ' ' << this_error(i) << std::endl;
+
+ Assert ((this_error(i)==0) || (error2(i)==0),
+ ExcInternalError());
+ if (this_error(i) != 0)
+ error2(i) = this_error(i);
+ }
+ }
+
+ // now compare the results of the two
+ // computations
+ Assert (error1 == error2, ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile ("error_estimator_02/output");
+ logfile.precision (2);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}