* Multithreading is only
* implemented in two and three
* dimensions.
+ *
+ * The @p subdomain_id parameter
+ * indicates whether we shall compute
+ * indicators for all cells (in case its
+ * value is the default,
+ * <tt>deal_II_numbers::invalid_unsigned_int</tt>),
+ * or only for the cells belonging to a
+ * certain subdomain with the given
+ * indicator. The latter case is used for
+ * parallel computations where all
+ * processor nodes have the global grid
+ * stored, and could well compute all the
+ * indicators for all cells themselves,
+ * but where it is more efficient to have
+ * each process compute only indicators
+ * for the cells it owns, and have them
+ * exchange the resulting information
+ * afterwards. This is in particular true
+ * for the case where meshes are very
+ * large and computing indicators for @em
+ * every cells is too expensive, while
+ * computing indicators for only local
+ * cells is acceptable. Note that if you
+ * only ask for the indicators of a
+ * certain subdomain to be computed, you
+ * must nevertheless make sure that this
+ * function has access to the correct
+ * node values of @em all degrees of
+ * freedom. This is since the function
+ * needs to access DoF values on
+ * neighboring cells as well, even if
+ * they belong to a different subdomain.
*/
template <typename InputVector>
static void estimate (const Mapping<dim> &mapping,
Vector<float> &error,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Calls the @p{estimate}
Vector<float> &error,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Same function as above, but
std::vector<Vector<float>*> &errors,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Calls the @p{estimate}
std::vector<Vector<float>*> &errors,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
*/
std::vector<double> JxW_values;
+ /**
+ * The subdomain ids we are to care
+ * for.
+ */
+ const unsigned int subdomain_id;
+
/**
* Constructor.
*/
PerThreadData (const unsigned int n_solution_vectors,
const unsigned int n_components,
- const unsigned int n_q_points);
+ const unsigned int n_q_points,
+ const unsigned int subdomain_id);
};
* entries, or an empty
* bit-vector.
*
- * The estimator supports
- * multithreading and splits the
- * cells to
+ * The estimator supports multithreading
+ * and splits the cells to
* @p{multithread_info.n_default_threads}
- * (default) threads. The number
- * of threads to be used in
- * multithreaded mode can be set
- * with the last parameter of the
- * error estimator.
- * Multithreading is only
- * implemented in two and three
- * dimensions.
+ * (default) threads. The number of
+ * threads to be used in multithreaded
+ * mode can be set with the last
+ * parameter of the error estimator.
+ * Multithreading is not presently
+ * implemented for 1d, but we retain the
+ * respective parameter for compatibility
+ * with the function signature in the
+ * general case.
*/
template <typename InputVector>
static void estimate (const Mapping<1> &mapping,
Vector<float> &error,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Calls the @p{estimate}
Vector<float> &error,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Same function as above, but
std::vector<Vector<float>*> &errors,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
* Calls the @p{estimate}
std::vector<Vector<float>*> &errors,
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 n_threads = multithread_info.n_default_threads,
+ const unsigned int subdomain_id = deal_II_numbers::invalid_unsigned_int);
/**
}
+template <int dim>
+static
+void advance_by_n (typename DoFHandler<dim>::active_cell_iterator &cell,
+ const unsigned int n)
+{
+ for (unsigned int t=0;
+ ((t<n) && (cell!=cell->get_dof_handler().end()));
+ ++t, ++cell);
+}
+
+
+
#if deal_II_dimension == 1
template <typename InputVector>
-void KellyErrorEstimator<1>::estimate (const Mapping<1> &mapping,
- const DoFHandler<1> &dof_handler,
- const Quadrature<0> &quadrature,
- const FunctionMap<1>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<1> *coefficients,
- const unsigned int n_threads)
+void
+KellyErrorEstimator<1>::
+estimate (const Mapping<1> &mapping,
+ const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ component_mask, coefficients, n_threads, subdomain_id);
}
template <typename InputVector>
-void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler,
- const Quadrature<0> &quadrature,
- const FunctionMap<1>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<1> *coefficients,
- const unsigned int n_threads)
+void
+KellyErrorEstimator<1>::
+estimate (const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ error, component_mask, coefficients, n_threads, subdomain_id);
}
template <typename InputVector>
-void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler,
- const Quadrature<0> &quadrature,
- const FunctionMap<1>::type &neumann_bc,
- const std::vector<const InputVector*> &solutions,
- std::vector<Vector<float>*> &errors,
- const std::vector<bool> &component_mask,
- const Function<1> *coefficients,
- const unsigned int n_threads)
+void
+KellyErrorEstimator<1>::
+estimate (const DoFHandler<1> &dof_handler,
+ const Quadrature<0> &quadrature,
+ const FunctionMap<1>::type &neumann_bc,
+ const std::vector<const InputVector*> &solutions,
+ std::vector<Vector<float>*> &errors,
+ const std::vector<bool> &component_mask,
+ const Function<1> *coefficients,
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ errors, component_mask, coefficients, n_threads, subdomain_id);
}
std::vector<Vector<float>*> &errors,
const std::vector<bool> &component_mask_,
const Function<1> *coefficient,
- const unsigned int)
+ const unsigned int,
+ const unsigned int subdomain_id)
{
const unsigned int n_components = dof_handler.get_fe().n_components();
const unsigned int n_solution_vectors = solutions.size();
for (unsigned int c=0; c<n_components; ++c)
coefficient_values(c) = 1;
- // loop over all
- // cells. note that the error
- // indicator is only a sum over the
- // two contributions from the two
- // vertices of each cell.
QTrapez<1> quadrature;
- FEValues<1> fe_values (mapping, dof_handler.get_fe(), quadrature, update_gradients);
+ FEValues<1> fe_values (mapping, dof_handler.get_fe(), quadrature,
+ update_gradients);
+
+ // loop over all cells and do something on
+ // the cells which we're told to work
+ // on. note that the error indicator is
+ // only a sum over the two contributions
+ // from the two vertices of each cell.
DoFHandler<1>::active_cell_iterator cell = dof_handler.begin_active();
- for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index)
- {
- for (unsigned int n=0; n<n_solution_vectors; ++n)
- (*errors[n])(cell_index) = 0;
+ 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))
+ {
+ for (unsigned int n=0; n<n_solution_vectors; ++n)
+ (*errors[n])(cell_index) = 0;
- // loop over the two points bounding
- // this line. n==0 is left point,
- // n==1 is right point
- for (unsigned int n=0; n<2; ++n)
- {
- // find right active neighbor
- DoFHandler<1>::cell_iterator neighbor = cell->neighbor(n);
- if (neighbor.state() == IteratorState::valid)
- while (neighbor->has_children())
- neighbor = neighbor->child(n==0 ? 1 : 0);
+ // loop over the two points bounding
+ // this line. n==0 is left point,
+ // n==1 is right point
+ for (unsigned int n=0; n<2; ++n)
+ {
+ // find left or right active
+ // neighbor
+ DoFHandler<1>::cell_iterator neighbor = cell->neighbor(n);
+ if (neighbor.state() == IteratorState::valid)
+ while (neighbor->has_children())
+ neighbor = neighbor->child(n==0 ? 1 : 0);
- // now get the gradients on the
- // both sides of the point
- fe_values.reinit (cell);
-
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- fe_values.get_function_grads (*solutions[s], gradients_here[s]);
-
- if (neighbor.state() == IteratorState::valid)
- {
- fe_values.reinit (neighbor);
-
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- fe_values.get_function_grads (*solutions[s],
- gradients_neighbor[s]);
-
- // extract the
- // gradients of all the
- // components. [0]
- // means: x-derivative,
- // which is the only
- // one here
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- for (unsigned int c=0; c<n_components; ++c)
- grad_neighbor[s](c) = gradients_neighbor[s][n==0 ? 1 : 0][c][0];
- }
- else
- if (neumann_bc.find(n) != neumann_bc.end())
- // if Neumann b.c., then fill
- // the gradients field which
- // will be used later on.
- {
- if (n_components==1)
- {
- const double
- v = neumann_bc.find(n)->second->value(cell->vertex(0));
+ // now get the gradients on the
+ // both sides of the point
+ fe_values.reinit (cell);
+
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ fe_values.get_function_grads (*solutions[s], gradients_here[s]);
+
+ if (neighbor.state() == IteratorState::valid)
+ {
+ fe_values.reinit (neighbor);
+
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ fe_values.get_function_grads (*solutions[s],
+ gradients_neighbor[s]);
+
+ // extract the
+ // gradients of all the
+ // components. [0]
+ // means: x-derivative,
+ // which is the only
+ // one here
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ for (unsigned int c=0; c<n_components; ++c)
+ grad_neighbor[s](c)
+ = gradients_neighbor[s][n==0 ? 1 : 0][c][0];
+ }
+ else
+ if (neumann_bc.find(n) != neumann_bc.end())
+ // if Neumann b.c., then fill
+ // the gradients field which
+ // will be used later on.
+ {
+ if (n_components==1)
+ {
+ const double
+ v = neumann_bc.find(n)->second->value(cell->vertex(0));
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- grad_neighbor[s](0) = v;
- }
- else
- {
- Vector<double> v(n_components);
- neumann_bc.find(n)->second->vector_value(cell->vertex(0), v);
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ grad_neighbor[s](0) = v;
+ }
+ else
+ {
+ Vector<double> v(n_components);
+ neumann_bc.find(n)->second->vector_value(cell->vertex(0), v);
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- grad_neighbor[s] = v;
- };
- }
- else
- // fill with zeroes.
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- grad_neighbor[s].clear ();
-
- // if there is a
- // coefficient, then
- // evaluate it at the
- // present position. if
- // there is none, reuse the
- // preset values.
- if (coefficient != 0)
- {
- if (coefficient->n_components == 1)
- {
- const double c_value = coefficient->value (cell->vertex(n));
- for (unsigned int c=0; c<n_components; ++c)
- coefficient_values(c) = c_value;
- }
- else
- coefficient->vector_value(cell->vertex(n),
- coefficient_values);
- };
-
-
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- for (unsigned int component=0; component<n_components; ++component)
- if (component_mask[component] == true)
- {
- // get gradient
- // here. [0] means
- // x-derivative
- // (there is no
- // other component
- // in 1d)
- const double grad_here = gradients_here[s][n][component][0];
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ grad_neighbor[s] = v;
+ };
+ }
+ else
+ // fill with zeroes.
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ grad_neighbor[s].clear ();
+
+ // if there is a
+ // coefficient, then
+ // evaluate it at the
+ // present position. if
+ // there is none, reuse the
+ // preset values.
+ if (coefficient != 0)
+ {
+ if (coefficient->n_components == 1)
+ {
+ const double c_value = coefficient->value (cell->vertex(n));
+ for (unsigned int c=0; c<n_components; ++c)
+ coefficient_values(c) = c_value;
+ }
+ else
+ coefficient->vector_value(cell->vertex(n),
+ coefficient_values);
+ }
+
+
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ for (unsigned int component=0; component<n_components; ++component)
+ if (component_mask[component] == true)
+ {
+ // get gradient
+ // here. [0] means
+ // x-derivative
+ // (there is no
+ // other component
+ // in 1d)
+ const double grad_here = gradients_here[s][n][component][0];
- const double jump = ((grad_here - grad_neighbor[s](component)) *
- coefficient_values(component));
- (*errors[s])(cell_index) += jump*jump * cell->diameter();
- };
- };
+ const double jump = ((grad_here - grad_neighbor[s](component)) *
+ coefficient_values(component));
+ (*errors[s])(cell_index) += jump*jump * cell->diameter();
+ }
+ }
- for (unsigned int s=0; s<n_solution_vectors; ++s)
- (*errors[s])(cell_index) = std::sqrt((*errors[s])(cell_index));
- };
+ for (unsigned int s=0; s<n_solution_vectors; ++s)
+ (*errors[s])(cell_index) = std::sqrt((*errors[s])(cell_index));
+ }
}
KellyErrorEstimator<dim>::PerThreadData::
PerThreadData (const unsigned int n_solution_vectors,
const unsigned int n_components,
- const unsigned int n_q_points)
+ const unsigned int n_q_points,
+ const unsigned int subdomain_id)
+ :
+ subdomain_id (subdomain_id)
{
// Init the size of a lot of vectors
// needed in the calculations once
phi[i][qp].resize(n_components);
psi[i][qp].resize(n_components);
neighbor_psi[i][qp].resize(n_components);
- };
- };
+ }
+ }
for (unsigned int qp=0;qp<n_q_points;++qp)
coefficient_values[qp].reinit(n_components);
// calls dimension dependent functions
template <int dim>
template <typename InputVector>
-void KellyErrorEstimator<dim>::estimate (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim-1> &quadrature,
- const typename FunctionMap<dim>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<dim> *coefficients,
- const unsigned int n_threads)
+void
+KellyErrorEstimator<dim>::
+estimate (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim-1> &quadrature,
+ const typename FunctionMap<dim>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<dim> *coefficients,
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ component_mask, coefficients, n_threads, subdomain_id);
}
template <int dim>
template <typename InputVector>
-void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim> &dof_handler,
- const Quadrature<dim-1> &quadrature,
- const typename FunctionMap<dim>::type &neumann_bc,
- const InputVector &solution,
- Vector<float> &error,
- const std::vector<bool> &component_mask,
- const Function<dim> *coefficients,
- const unsigned int n_threads)
+void
+KellyErrorEstimator<dim>::
+estimate (const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim-1> &quadrature,
+ const typename FunctionMap<dim>::type &neumann_bc,
+ const InputVector &solution,
+ Vector<float> &error,
+ const std::vector<bool> &component_mask,
+ const Function<dim> *coefficients,
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ error, component_mask, coefficients, n_threads, subdomain_id);
}
PerThreadData &per_thread_data)
{
const unsigned int n_solution_vectors = solutions.size();
+
+ const unsigned int subdomain_id = per_thread_data.subdomain_id;
// make up a fe face values object
// for the restriction of the
// loop over all cells for this thread
- // the iteration of cell is done at the end
- for (; cell!=dof_handler.end(); )
+ for (; cell!=dof_handler.end();
+ advance_by_n<dim>(cell,this_thread.second))
{
// loop over all faces of this cell
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
// make sure we do work
// only once: this face
// if the neighboring cell is less
- // refined than the present one, then
- // do nothing since we integrate
- // over the subfaces when we visit
- // the coarse cells.
+ // refined than the present one,
+ // then do nothing since we
+ // integrate over the subfaces when
+ // we visit the coarse cells.
if (cell->at_boundary(face_no) == false)
if (cell->neighbor(face_no)->level() < cell->level())
continue;
- // if this face is part of the boundary
- // but not of the neumann boundary
- // -> nothing to do. However, to make
- // things easier when summing up the
- // contributions of the faces of cells,
- // we enter this face into the list
- // of faces with contribution zero.
+ // if this face is part of the
+ // boundary but not of the neumann
+ // boundary -> nothing to
+ // do. However, to make things
+ // easier when summing up the
+ // contributions of the faces of
+ // cells, we enter this face into
+ // the list of faces with
+ // contribution zero.
const unsigned char boundary_indicator
= cell->face(face_no)->boundary_indicator();
if (cell->face(face_no)->at_boundary()
for (unsigned int n=0; n<n_solution_vectors; ++n)
face_integrals[cell->face(face_no)][n] = 0;
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)))
+ {
+ // ok, cell is unwanted, but
+ // maybe its neighbor behind
+ // the face we presently work
+ // on? oh is there a face at
+ // all?
+ if (cell->at_boundary(face_no))
+ continue;
+
+ bool care_for_cell = false;
+ if (cell->neighbor(face_no)->has_children() == false)
+ care_for_cell |= (cell->neighbor(face_no)->subdomain_id()
+ == subdomain_id);
+ 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)
+ {
+ care_for_cell = true;
+ break;
+ }
+ }
+
+ // so if none of the neighbors
+ // cares for this subdomain
+ // either, then try next face
+ if (care_for_cell == false)
+ continue;
+ }
+
+
+ // so now we know that we care for
+ // this face, let's do something
+ // about it:
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
+ // 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 (dof_handler, quadrature,
neumann_bc, solutions, component_mask,
cell, face_no,
fe_face_values_cell,
fe_subface_values);
- };
-
- // go to next cell for this
- // thread. note that the cells
- // for each of the threads are
- // interleaved.
- for (unsigned int t=0;
- ((t<this_thread.second) && (cell!=dof_handler.end()));
- ++t, ++cell);
- };
+ }
+ }
}
template <int dim>
template <typename InputVector>
void
-KellyErrorEstimator<dim>::estimate (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim-1> &quadrature,
- const typename FunctionMap<dim>::type &neumann_bc,
- const std::vector<const InputVector *> &solutions,
- std::vector<Vector<float>*> &errors,
- const std::vector<bool> &component_mask_,
- const Function<dim> *coefficients,
- const unsigned int n_threads_)
+KellyErrorEstimator<dim>::
+estimate (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim-1> &quadrature,
+ const typename FunctionMap<dim>::type &neumann_bc,
+ const std::vector<const InputVector *> &solutions,
+ std::vector<Vector<float>*> &errors,
+ const std::vector<bool> &component_mask_,
+ const Function<dim> *coefficients,
+ const unsigned int n_threads_,
+ const unsigned int subdomain_id)
{
const unsigned int n_components = dof_handler.get_fe().n_components();
for (typename FunctionMap<dim>::type::const_iterator i=neumann_bc.begin();
i!=neumann_bc.end(); ++i)
- Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
+ Assert (i->second->n_components == n_components,
+ ExcInvalidBoundaryFunction());
Assert ((component_mask_.size() == 0) ||
(component_mask_.size() == n_components), ExcInvalidComponentMask());
std::vector<bool>(n_components, true) :
component_mask_);
Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
- Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
+ Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
ExcInvalidComponentMask());
// if NOT multithreaded, set
// in multithreaded mode. Negative
// value indicates that the face
// has not yet been processed.
- std::vector<double> default_face_integrals (n_solution_vectors, -10e20);
+ const double invalid_double = -10e20;
+ std::vector<double> default_face_integrals (n_solution_vectors,
+ invalid_double);
FaceIntegrals face_integrals;
- for (active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ for (active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
face_integrals[cell->face(face_no)] = default_face_integrals;
// components
std::vector<PerThreadData*> data_structures (n_threads);
for (unsigned int i=0; i<n_threads; ++i)
- data_structures[i] = new PerThreadData (solutions.size(),
- dof_handler.get_fe().n_components(),
- quadrature.n_quadrature_points);
+ data_structures[i]
+ = new PerThreadData (solutions.size(),
+ dof_handler.get_fe().n_components(),
+ quadrature.n_quadrature_points,
+ subdomain_id);
// split all cells into threads if
// multithreading is used and run
{
delete data_structures[i];
data_structures[i] = 0;
- };
+ }
// finally add up the contributions of the
};
unsigned int present_cell=0;
-
+
+ // now walk over all cells and collect
+ // information from the faces. only do
+ // something if this is a cell we care for
+ // based on the subdomain id
for (active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end();
++cell, ++present_cell)
- {
- // loop over all faces of this cell
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- {
- Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
- ExcInternalError());
- const double factor = cell->diameter() / 24;
+ if ((subdomain_id == deal_II_numbers::invalid_unsigned_int)
+ ||
+ (cell->subdomain_id() == subdomain_id))
+ {
+ // loop over all faces of this cell
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ Assert(face_integrals.find(cell->face(face_no))
+ != face_integrals.end(),
+ ExcInternalError());
+ const double factor = cell->diameter() / 24;
- for (unsigned int n=0; n<n_solution_vectors; ++n)
- (*errors[n])(present_cell) += (face_integrals[cell->face(face_no)][n] *
- factor);
- };
-
- for (unsigned int n=0; n<n_solution_vectors; ++n)
- (*errors[n])(present_cell) = std::sqrt((*errors[n])(present_cell));
- };
+ for (unsigned int n=0; n<n_solution_vectors; ++n)
+ {
+ // make sure that we have
+ // written a meaningful value
+ // into this slot
+ Assert (face_integrals[cell->face(face_no)][n] >= 0,
+ ExcInternalError());
+
+ (*errors[n])(present_cell)
+ += (face_integrals[cell->face(face_no)][n] * factor);
+ }
+ }
+
+ for (unsigned int n=0; n<n_solution_vectors; ++n)
+ (*errors[n])(present_cell) = std::sqrt((*errors[n])(present_cell));
+ }
}
std::vector<Vector<float>*> &errors,
const std::vector<bool> &component_mask,
const Function<dim> *coefficients,
- const unsigned int n_threads)
+ const unsigned int n_threads,
+ const unsigned int subdomain_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);
+ errors, component_mask, coefficients, n_threads, subdomain_id);
}
Vector<float> &, \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
- const unsigned int ); \
+ const unsigned int , \
+ const unsigned int); \
\
template \
void \
Vector<float> &, \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
- const unsigned int ); \
+ const unsigned int , \
+ const unsigned int); \
\
template \
void \
std::vector<Vector<float>*> &, \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
- const unsigned int ); \
+ const unsigned int , \
+ const unsigned int); \
\
template \
void \
std::vector<Vector<float>*> &, \
const std::vector<bool> &, \
const Function<deal_II_dimension> *, \
- const unsigned int )
+ const unsigned int , \
+ const unsigned int)
INSTANTIATE(Vector<double>);
INSTANTIATE(Vector<float>);
projects / dealii-svn.git / commitdiff