#include <dofs/dof_handler.h>
#include <map>
-
+ // if multithreaded set number
+ // of threads to use.
+ // The default number of threads can
+ // be set during the compiling with
+ // the flag N_THREADS.
+#ifdef DEAL_II_USE_MT
+ #ifndef N_THREADS
+ #define N_THREADS 4
+ #endif
+#else
+ #ifdef N_THREADS
+ #undef N_THREADS
+ #endif
+ #define N_THREADS 1
+#endif
/**
template <int dim>
class KellyErrorEstimator {
public:
- /**
+
+ /**
* Declare a data type which denotes a
* mapping between a boundary indicator
* and the function denoting the boundary
* Only one boundary function may be given
* for each boundary indicator, which is
* guaranteed by the #map# data type.
- */
+ */
typedef map<unsigned char,const Function<dim>*> FunctionMap;
+
/**
* Implementation of the error
* estimator described above. You
* bit-vector with all-set
* entries, or an empty
* bit-vector.
+ *
+ * The estimator supports
+ * multithreading and splits
+ * the cells to 4 (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.
*/
+
static void estimate (const DoFHandler<dim> &dof,
const Quadrature<dim-1> &quadrature,
const FunctionMap &neumann_bc,
const Vector<double> &solution,
Vector<float> &error,
- const vector<bool> &component_mask = vector<bool>(),
- const Function<dim> *coefficients = 0);
+ const vector<bool> &component_mask_ = vector<bool>(),
+ const Function<dim> *coefficients = 0,
+ unsigned int n_threads=N_THREADS);
/**
* Exception
*/
DeclException0 (ExcInvalidBoundaryFunction);
+
+
private:
- /**
+
+
+ /**
* Declare a data type to represent the
* mapping between faces and integrated
* jumps of gradients. See the general
*/
typedef map<typename DoFHandler<dim>::face_iterator,double> FaceIntegrals;
+
/**
* Redeclare an active cell iterator.
* This is simply for convenience.
*/
typedef DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+
+ /**
+ * All data needed by the several functions
+ * of the error estimator is gathered in
+ * this struct. It is passed as an reference
+ * to the seperate functions in the object.
+ * Global data is not possible because no
+ * real object is created.
+ */
+ struct Data
+ {
+ const DoFHandler<dim> &dof;
+ const Quadrature<dim-1> &quadrature;
+ const FunctionMap &neumann_bc;
+ const Vector<double> &solution;
+ vector<bool> component_mask;
+ const Function<dim> *coefficients;
+ unsigned int n_threads;
+
+ DoFHandler<dim>::active_cell_iterator endc;
+ unsigned int n_components;
+ unsigned int n_q_points;
+ FaceIntegrals face_integrals;
+
+ /**
+ * A vector to store the jump of the
+ * normal vectors in the quadrature
+ * points.
+ * There is one vector for every
+ * thread used in the estimator.
+ * The allocation of memory has to
+ * be global to enable fast use
+ * of multithreading
+ */
+ vector< vector<vector<double> > > phi;
+ /**
+ * A vector for the gradients of
+ * the finite element function
+ * on one cell
+ *
+ * Let psi be a short name for
+ * #a grad u_h#, where the second
+ * index be the component of the
+ * finite element, and the first
+ * index the number of the
+ * quadrature point.
+ */
+ vector< vector<vector<Tensor<1,dim> > > > psi;
+ /**
+ * The same vector for a neighbor cell
+ */
+ vector< vector<vector<Tensor<1,dim> > > > neighbor_psi;
+ /**
+ * The normal vectors of the finite
+ * element function on one face
+ */
+ vector< vector<Point<dim> > > normal_vectors;
+
+ vector< vector<double> > coefficient_values1;
+ vector< vector<Vector<double> > > coefficient_values;
+
+ vector< vector<double> > JxW_values;
+
+ /**
+ * A constructor of the
+ * class Data. All variables are
+ * passed as references.
+ */
+ Data(const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &quadrature,
+ const FunctionMap &neumann_bc,
+ const Vector<double> &solution,
+ vector<bool> component_mask_,
+ const Function<dim> *coefficients,
+ unsigned int n_threads);
+ };
+
+ /**
+ * Computates the error on all cells
+ * of the domain with the number n,
+ * satisfying
+ * #n=this_thread (mod n_threads)#
+ * This enumeration is choosen to
+ * generate a random distribution
+ * of all cells.
+ *
+ * This function is only needed in
+ * two or three dimensions.
+ * The Errorestimator in one dimension
+ * is implemented seperatly.
+ */
+ static void * estimate_some(Data &data, unsigned int this_thread);
+
/**
* Actually do the computation on a face
* which has no hanging nodes (it is
* to avoid ending up with a function
* of 500 lines of code.
*/
- static void integrate_over_regular_face (const active_cell_iterator &cell,
+
+
+ static void integrate_over_regular_face (Data &data,
+ int this_thread,
+ const active_cell_iterator &cell,
const unsigned int face_no,
- const FunctionMap &neumann_bc,
- const unsigned int n_q_points,
FEFaceValues<dim> &fe_face_values_cell,
- FEFaceValues<dim> &fe_face_values_neighbor,
- FaceIntegrals &face_integrals,
- const Vector<double>&solution,
- const vector<bool> &component_mask,
- const Function<dim> *coefficient);
-
+ FEFaceValues<dim> &fe_face_values_neighbor);
+
+
/**
* The same applies as for the function
* above, except that integration is
* so that the integration is a bit more
* complex.
*/
- static void integrate_over_irregular_face (const active_cell_iterator &cell,
+ static void integrate_over_irregular_face (Data &data,
+ int this_thread,
+ const active_cell_iterator &cell,
const unsigned int face_no,
- const unsigned int n_q_points,
FEFaceValues<dim> &fe_face_values,
- FESubfaceValues<dim> &fe_subface_values,
- FaceIntegrals &face_integrals,
- const Vector<double> &solution,
- const vector<bool> &component_mask,
- const Function<dim> *coefficient);
+ FESubfaceValues<dim> &fe_subface_values);
};
/* end of #ifndef __error_estimator_H */
#endif
/*---------------------------- error_estimator.h ---------------------------*/
+
+
#include <algorithm>
#include <cmath>
#include <vector>
+#include <base/timer.h>
+
+ // if multithreaded include
+ // ThreadManager
+#ifdef DEAL_II_USE_MT
+#include <base/thread_manager.h>
+#endif
+
+#if deal_II_dimension == 1
+
+template <>
+KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &,
+ const Quadrature<0> &,
+ const FunctionMap &,
+ const Vector<double> &,
+ vector<bool> ,
+ const Function<1> *,
+ unsigned int ):
+ dof(* static_cast <const DoFHandler<1> *> (0)),
+ quadrature(* static_cast <const Quadrature<0> *> (0)),
+ neumann_bc(* static_cast <const FunctionMap *> (0)),
+ solution(* static_cast <const Vector<double> *> (0))
+{
+ Assert (false, ExcInternalError());
+}
+
+#else
+
+template <int dim>
+KellyErrorEstimator<dim>::Data::Data(const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &quadrature,
+ const FunctionMap &neumann_bc,
+ const Vector<double> &solution,
+ vector<bool> component_mask_,
+ const Function<dim> *coefficients,
+ unsigned int n_threads):
+ dof(dof),
+ quadrature(quadrature),
+ neumann_bc(neumann_bc),
+ solution(solution),
+ coefficients(coefficients),
+ n_threads(n_threads)
+{
+ n_components = dof.get_fe().n_components();
+
+ // if no mask given: treat all components
+ component_mask = ((component_mask_.size() == 0) ?
+ vector<bool>(n_components, true) :
+ component_mask_);
+
+ Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
+ Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
+ ExcInvalidComponentMask());
+
+ Assert ((coefficients == 0) ||
+ (coefficients->n_components == n_components) ||
+ (coefficients->n_components == 1),
+ ExcInvalidCoefficient());
+
+ Assert (neumann_bc.find(255) == neumann_bc.end(),
+ ExcInvalidBoundaryIndicator());
+
+ for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i)
+ Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
+
+ // the last cell, often needed
+ endc=dof.end();
+
+ n_q_points=quadrature.n_quadrature_points;
+
+ // Init the size of a lot of vectors
+ // needed in the calculations once
+ // per thread.
+ phi.resize(n_threads);
+ psi.resize(n_threads);
+ neighbor_psi.resize(n_threads);
+ normal_vectors.resize(n_threads);
+ coefficient_values1.resize(n_threads);
+ coefficient_values.resize(n_threads);
+ JxW_values.resize(n_threads);
+
+ for (unsigned int t=0;t<n_threads;++t)
+ {
+ phi[t].resize(n_q_points);
+ psi[t].resize(n_q_points);
+ neighbor_psi[t].resize(n_q_points);
+ normal_vectors[t].resize(n_q_points);
+ coefficient_values1[t].resize(n_q_points);
+ coefficient_values[t].resize(n_q_points);
+ JxW_values[t].resize(n_q_points);
+
+ for (unsigned int qp=0;qp<n_q_points;++qp)
+ {
+ phi[t][qp].resize(n_components);
+ psi[t][qp].resize(n_components);
+ neighbor_psi[t][qp].resize(n_components);
+ coefficient_values[t][qp].reinit(n_components);
+ }
+ }
+}
+
+#endif
inline static double sqr (const double x) {
};
-
#if deal_II_dimension == 1
+template <>
+void * KellyErrorEstimator<1>::estimate_some (Data &,const unsigned int )
+{
+ Assert (false, ExcInternalError() );
+ return 0;
+}
+
template <>
void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof,
const Quadrature<0> &,
const Vector<double> &solution,
Vector<float> &error,
const vector<bool> &component_mask_,
- const Function<1> *coefficient)
+ const Function<1> *coefficient,
+ unsigned int)
{
const unsigned int n_components = dof.get_fe().n_components();
};
};
-#endif
+ // #if deal_II_dimension !=1
+#else
template <int dim>
-void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &quadrature,
- const FunctionMap &neumann_bc,
- const Vector<double> &solution,
- Vector<float> &error,
- const vector<bool> &component_mask_,
- const Function<dim> *coefficient)
+void * KellyErrorEstimator<dim>::estimate_some (Data &data, unsigned int this_thread)
{
- const unsigned int n_components = dof.get_fe().n_components();
-
- // if no mask given: treat all components
- vector<bool> component_mask ((component_mask_.size() == 0) ?
- vector<bool>(n_components, true) :
- component_mask_);
- Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
- Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
- ExcInvalidComponentMask());
-
- Assert ((coefficient == 0) ||
- (coefficient->n_components == n_components) ||
- (coefficient->n_components == 1),
- ExcInvalidCoefficient());
-
- Assert (neumann_bc.find(255) == neumann_bc.end(),
- ExcInvalidBoundaryIndicator());
-
- for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i)
- Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
-
- // create a map of integrals indexed by
- // the corresponding face. In this map
- // we store the integrated jump of the
- // gradient for each face. By doing so,
- // we can check whether we have already
- // integrated along this face by testing
- // whether the respective face is already
- // a key in this map.
- // At the end of the function, we again
- // loop over the cells and collect the
- // conrtibutions of the different faces
- // of the cell.
- FaceIntegrals face_integrals;
-
- // number of integration points per face
- const unsigned int n_q_points = quadrature.n_quadrature_points;
// make up a fe face values object for the
// restriction of the finite element function
// need not compute all values on the
// neighbor cells, so using two objects
// gives us a performance gain).
- FEFaceValues<dim> fe_face_values_cell (dof.get_fe(), quadrature,
+ FEFaceValues<dim> fe_face_values_cell (data.dof.get_fe(),
+ data.quadrature,
UpdateFlags(update_gradients |
update_JxW_values |
- ((!neumann_bc.empty() ||
- (coefficient != 0)) ?
+ ((!data.neumann_bc.empty() ||
+ (data.coefficients != 0)) ?
update_q_points : 0) |
update_normal_vectors));
- FEFaceValues<dim> fe_face_values_neighbor (dof.get_fe(), quadrature, update_gradients);
- FESubfaceValues<dim> fe_subface_values (dof.get_fe(), quadrature, update_gradients);
+ FEFaceValues<dim> fe_face_values_neighbor (data.dof.get_fe(),
+ data.quadrature,
+ update_gradients);
+ FESubfaceValues<dim> fe_subface_values (data.dof.get_fe(),
+ data.quadrature,
+ update_gradients);
- // 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<GeometryInfo<dim>::faces_per_cell; ++face_no)
- {
- // if we already visited this
- // face: do nothing
- if (face_integrals.find(cell->face(face_no)) !=
- face_integrals.end())
- continue;
-
- // 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.
- 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.
- const unsigned char boundary_indicator
- = cell->face(face_no)->boundary_indicator();
- if ((boundary_indicator != 255) &&
- neumann_bc.find(boundary_indicator)==neumann_bc.end())
- {
- face_integrals[cell->face(face_no)] = 0;
+
+ DoFHandler<dim>::active_cell_iterator cell=data.dof.begin_active();
+
+ // calculate the start cell for this
+ // thread. the enumeration is choosen
+ // in this strange way to generate a
+ // "random" distribution of the cells.
+ // if the sequence of the iterator would
+ // be used, the threads would take widely
+ // spread times to calculate their cells.
+ for (unsigned int t=0;t<this_thread;++t,++cell);
+
+ // loop over all cells for this thread
+ // the iteration of cell is done at the end
+ for (;cell!=data.endc;)
+ {
+
+ // loop over all faces of this cell
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ // if we already visited this
+ // face: do nothing
+ if (data.face_integrals[cell->face(face_no)] >=0)
continue;
- };
+
+
+ // 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.
+ 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.
+ const unsigned char boundary_indicator
+ = cell->face(face_no)->boundary_indicator();
+ if ((boundary_indicator != 255) &&
+ data.neumann_bc.find(boundary_indicator)==data.neumann_bc.end())
+ {
+ data.face_integrals[cell->face(face_no)] = 0;
+ continue;
+ };
+
+
+
+ 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 (data,
+ this_thread,
+ cell, face_no,
+ fe_face_values_cell,
+ fe_face_values_neighbor);
+
+ else
+ // otherwise we need to do some
+ // special computations which do
+ // not fit into the framework of
+ // the above function
+ integrate_over_irregular_face (data,
+ this_thread,cell, face_no,
+ fe_face_values_cell,
+ fe_subface_values);
+ };
-
- 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, neumann_bc,
- n_q_points,
- fe_face_values_cell,
- fe_face_values_neighbor,
- face_integrals,
+ // next cell in this thread
+ for (unsigned int t=0;((t<data.n_threads)&&(cell!=data.endc));++t,++cell) {};
+ };
+ return 0;
+}
+
+template <int dim>
+void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &quadrature,
+ const FunctionMap &neumann_bc,
+ const Vector<double> &solution,
+ Vector<float> &error,
+ const vector<bool> &component_mask,
+ const Function<dim> *coefficients,
+ unsigned int n_threads)
+{
+ // if NOT multithreaded, set n_threads to one
+#ifndef DEAL_II_USE_MT
+ n_threads = 1;
+#endif
+
+ // all the data needed in the error-
+ // estimator is gathered in this stuct.
+ KellyErrorEstimator<dim>::Data data (dof,
+ quadrature,
+ neumann_bc,
solution,
component_mask,
- coefficient);
- else
- // otherwise we need to do some
- // special computations which do
- // not fit into the framework of
- // the above function
- integrate_over_irregular_face (cell, face_no,
- n_q_points,
- fe_face_values_cell,
- fe_subface_values,
- face_integrals, solution,
- component_mask,
- coefficient);
- };
-
+ coefficients,
+ n_threads);
+ // map of integrals indexed by
+ // the corresponding face. In this map
+ // we store the integrated jump of the
+ // gradient for each face. By doing so,
+ // we can check whether we have already
+ // integrated along this face by testing
+ // whether the respective face is already
+ // a key in this map.
+ // At the end of the function, we again
+ // loop over the cells and collect the
+ // conrtibutions of the different faces
+ // of the cell.
+ // the values for all faces are set to
+ // -10e20. It would cost a lot of time
+ // to syncronisise the initialisation
+ // of the map in multithreaded mode.
+ // negative value indicates that the
+ // face is not calculated.
+
+ for (active_cell_iterator cell=data.dof.begin_active(); cell!=data.endc; ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ data.face_integrals[cell->face(face_no)]=-10e20;
+
+
+ // split all cells into threads
+ // if multithreading is used
+#ifdef DEAL_II_USE_MT
- // finally add up the contributions of the
+ ThreadManager thread_manager;
+
+ // all data needed to start
+ // one thread is gathered
+ // in this struct.
+ typedef ThreadManager::Fun_Data2
+ <class Data &, unsigned int > FunData;
+
+ // One struct of this type for every thread
+ vector<FunData>
+ fun_data (data.n_threads,
+ FunData (data,0,&KellyErrorEstimator::estimate_some));
+
+
+ // get start cells for each thread
+ for (unsigned int l=0;l<data.n_threads;++l)
+ {
+ fun_data[l].arg2=l;
+ };
+
+ // now spawn the threads
+ for (unsigned int i=0;i<data.n_threads; ++i)
+ {
+ thread_manager.spawn(&fun_data[i],THR_SCOPE_SYSTEM | THR_DETACHED);
+ };
+ // wait for all threads to return
+ thread_manager.wait();
+
+ // ... ifdef DEAL_II_USE_MT
+#else
+ // just one thread, calculate
+ // error on all cells
+ KellyErrorEstimator::estimate_some(data,0);
+
+#endif
+
+
+ // finally add up the contributions of the
// faces for each cell
// reserve one slot for each cell and set
// it to zero
- error.reinit (dof.get_tria().n_active_cells());
+ error.reinit (data.dof.get_tria().n_active_cells());
+ for (unsigned int i=0;i<data.dof.get_tria().n_active_cells();++i)
+ error(i)=0;
unsigned int present_cell=0;
- for (active_cell_iterator cell=dof.begin_active(); cell!=endc; ++cell, ++present_cell)
+
+ for (active_cell_iterator cell=data.dof.begin_active();
+ cell!=data.endc;
+ ++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)
+ 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());
- error(present_cell) += (face_integrals[cell->face(face_no)] *
+ Assert(data.face_integrals.find(cell->face(face_no)) != data.face_integrals.end(),
+ ExcInternalError());
+ error(present_cell) += (data.face_integrals[cell->face(face_no)] *
cell->diameter() / 24);
};
-
error(present_cell) = sqrt(error(present_cell));
};
};
+#endif
#if deal_II_dimension == 1
template <>
-void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iterator &,
- const unsigned int ,
- const FunctionMap &,
+void KellyErrorEstimator<1>::integrate_over_regular_face (Data &,
+ int ,
+ const active_cell_iterator &,
const unsigned int ,
FEFaceValues<1> &,
- FEFaceValues<1> &,
- FaceIntegrals &,
- const Vector<double> &,
- const vector<bool> &,
- const Function<1> *) {
+ FEFaceValues<1> &)
+{
Assert (false, ExcInternalError());
};
template <>
void KellyErrorEstimator<1>::
-integrate_over_irregular_face (const active_cell_iterator &,
- const unsigned int ,
+integrate_over_irregular_face (Data &,
+ int,
+ const active_cell_iterator &,
const unsigned int ,
FEFaceValues<1> &,
- FESubfaceValues<1> &,
- FaceIntegrals &,
- const Vector<double> &,
- const vector<bool> &,
- const Function<1> *) {
+ FESubfaceValues<1> &)
+{
Assert (false, ExcInternalError());
};
template <int dim>
void KellyErrorEstimator<dim>::
-integrate_over_regular_face (const active_cell_iterator &cell,
- const unsigned int face_no,
- const FunctionMap &neumann_bc,
- const unsigned int n_q_points,
+integrate_over_regular_face (Data &data,
+ int this_thread,
+ const active_cell_iterator &cell,
+ const unsigned int face_no,
FEFaceValues<dim> &fe_face_values_cell,
- FEFaceValues<dim> &fe_face_values_neighbor,
- FaceIntegrals &face_integrals,
- const Vector<double> &solution,
- const vector<bool> &component_mask,
- const Function<dim> *coefficient)
+ FEFaceValues<dim> &fe_face_values_neighbor)
{
- const unsigned int n_components = component_mask.size();
const DoFHandler<dim>::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);
- // 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
- // [a grad u_h], where the second
- // index be the component of the
- // finite element, and the first
- // index the number of the
- // quadrature point
- vector<vector<Tensor<1,dim> > > psi(n_q_points, vector<Tensor<1,dim> >(n_components));
- fe_face_values_cell.get_function_grads (solution, psi);
+ // get gradients of the finite element
+ // function on this cell
+ fe_face_values_cell.get_function_grads (data.solution, data.psi[this_thread]);
// now compute over the other side of
// the face
// function of #neighbor# to the
// common face.
fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
-
- // get a list of the gradients of
- // the finite element solution
- // restricted to the neighbor cell
- vector<vector<Tensor<1,dim> > > neighbor_psi (n_q_points,
- vector<Tensor<1,dim> >(n_components));
- fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+
+ // get gradients on neighbor cell
+ fe_face_values_neighbor.get_function_grads (data.solution,
+ data.neighbor_psi[this_thread]);
// compute the jump in the gradients
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int p=0; p<n_q_points; ++p)
- psi[p][component] -= neighbor_psi[p][component];
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int p=0; p<data.n_q_points; ++p)
+ data.psi[this_thread][p][component] -= data.neighbor_psi[this_thread][p][component];
};
// taking that of the neighbor
// would only change the sign. We take
// the outward normal.
- //
- // let phi be the name of the integrand
- vector<vector<double> > phi(n_q_points, vector<double>(n_components));
- const vector<Point<dim> > &normal_vectors(fe_face_values_cell.
- get_normal_vectors());
-
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] = psi[point][component]*normal_vectors[point];
-
+
+ data.normal_vectors[this_thread]=fe_face_values_cell.get_normal_vectors();
+
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] = data.psi[this_thread][point][component]*
+ data.normal_vectors[this_thread][point];
+
// if a coefficient was given: use that
// to scale the jump in the gradient
- if (coefficient != 0)
+ if (data.coefficients != 0)
{
// scalar coefficient
- if (coefficient->n_components == 1)
+ if (data.coefficients->n_components == 1)
{
- vector<double> coefficient_values (n_q_points);
- coefficient->value_list (fe_face_values_cell.get_quadrature_points(),
- coefficient_values);
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] *= coefficient_values[point];
+
+ data.coefficients->value_list (fe_face_values_cell.get_quadrature_points(),
+ data.coefficient_values1[this_thread]);
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] *=
+ data.coefficient_values1[this_thread][point];
}
else
- // vector-valued coefficient
+ // vector-valued coefficient
{
- vector<Vector<double> > coefficient_values (n_q_points,
- Vector<double>(n_components));
- coefficient->vector_value_list (fe_face_values_cell.get_quadrature_points(),
- coefficient_values);
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] *= coefficient_values[point](component);
+ data.coefficients->vector_value_list (fe_face_values_cell.get_quadrature_points(),
+ data.coefficient_values[this_thread]);
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] *=
+ data.coefficient_values[this_thread][point](component);
};
};
-
+
if (face->at_boundary() == true)
// neumann boundary face. compute
// derivative and boundary function
{
const unsigned char boundary_indicator = face->boundary_indicator();
-
- Assert (neumann_bc.find(boundary_indicator) != neumann_bc.end(),
+
+ Assert (data.neumann_bc.find(boundary_indicator) != data.neumann_bc.end(),
ExcInternalError ());
// get the values of the boundary
// function at the quadrature
// points
- vector<Vector<double> > g(n_q_points, Vector<double>(n_components));
- neumann_bc.find(boundary_indicator)->second
+
+ vector<Vector<double> > g(data.n_q_points, Vector<double>(data.n_components));
+ data.neumann_bc.find(boundary_indicator)->second
->vector_value_list (fe_face_values_cell.get_quadrature_points(),
g);
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] -= g[point](component);
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] -= g[point](component);
};
// each component being the
// mentioned value at one of the
// quadrature points
- const vector<double> &JxW_values = fe_face_values_cell.get_JxW_values();
+
+ data.JxW_values[this_thread] = fe_face_values_cell.get_JxW_values();
// take the square of the phi[i]
// for integration, and sum up
double face_integral = 0;
- for (unsigned int component=0; component<n_components; ++component)
- if (component_mask[component] == true)
- for (unsigned int p=0; p<n_q_points; ++p)
- face_integral += sqr(phi[p][component]) *
- JxW_values[p];
+ for (unsigned int component=0; component<data.n_components; ++component)
+ if (data.component_mask[component] == true)
+ for (unsigned int p=0; p<data.n_q_points; ++p)
+ face_integral += sqr(data.phi[this_thread][p][component]) *
+ data.JxW_values[this_thread][p];
- face_integrals[face] = face_integral;
+ data.face_integrals[face] = face_integral;
};
template <int dim>
void KellyErrorEstimator<dim>::
-integrate_over_irregular_face (const active_cell_iterator &cell,
+integrate_over_irregular_face (Data &data,
+ int this_thread,
+ const active_cell_iterator &cell,
const unsigned int face_no,
- const unsigned int n_q_points,
FEFaceValues<dim> &fe_face_values,
- FESubfaceValues<dim> &fe_subface_values,
- FaceIntegrals &face_integrals,
- const Vector<double> &solution,
- const vector<bool> &component_mask,
- const Function<dim> *coefficient)
+ FESubfaceValues<dim> &fe_subface_values)
{
- const unsigned int n_components = component_mask.size();
-
const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
Assert (neighbor.state() == valid, ExcInternalError());
Assert (neighbor->has_children(), ExcInternalError());
// finite element, and the first
// index the number of the
// quadrature point
- vector<vector<Tensor<1,dim> > > psi(n_q_points, vector<Tensor<1,dim> >(n_components));
-
+
// store which number #cell# has in the
// list of neighbors of #neighbor#
unsigned int neighbor_neighbor;
Assert (!neighbor->child(GeometryInfo<dim>::
child_cell_on_face(neighbor_neighbor,subface_no))->has_children(),
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_subface_values.get_function_grads (solution, psi);
+ fe_subface_values.get_function_grads (data.solution, data.psi[this_thread]);
// restrict the finite element on the
// neighbor cell to the common #subface#.
// store the gradient in #neighbor_psi#
- vector<vector<Tensor<1,dim> > > neighbor_psi (n_q_points,
- vector<Tensor<1,dim> >(n_components));
+
fe_face_values.reinit (neighbor_child, neighbor_neighbor);
- fe_face_values.get_function_grads (solution, neighbor_psi);
+ fe_face_values.get_function_grads (data.solution, data.neighbor_psi[this_thread]);
// compute the jump in the gradients
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int p=0; p<n_q_points; ++p)
- psi[p][component] -= neighbor_psi[p][component];
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int p=0; p<data.n_q_points; ++p)
+ data.psi[this_thread][p][component] -=
+ data.neighbor_psi[this_thread][p][component];
// note that unlike for the
// case of regular faces
// the outward normal.
//
// let phi be the name of the integrand
- vector<vector<double> > phi(n_q_points, vector<double>(n_components));
- const vector<Point<dim> > &normal_vectors(fe_face_values.
- get_normal_vectors());
-
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] = psi[point][component]*normal_vectors[point];
+
+ data.normal_vectors[this_thread]=fe_face_values.get_normal_vectors();
+
+
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] =
+ data.psi[this_thread][point][component]*
+ data.normal_vectors[this_thread][point];
// if a coefficient was given: use that
// to scale the jump in the gradient
- if (coefficient != 0)
+ if (data.coefficients != 0)
{
// scalar coefficient
- if (coefficient->n_components == 1)
+ if (data.coefficients->n_components == 1)
{
- vector<double> coefficient_values (n_q_points);
- coefficient->value_list (fe_face_values.get_quadrature_points(),
- coefficient_values);
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] *= coefficient_values[point];
+ data.coefficients->value_list (fe_face_values.get_quadrature_points(),
+ data.coefficient_values1[this_thread]);
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] *=
+ data.coefficient_values1[this_thread][point];
}
else
// vector-valued coefficient
{
- vector<Vector<double> > coefficient_values (n_q_points,
- Vector<double>(n_components));
- coefficient->vector_value_list (fe_face_values.get_quadrature_points(),
- coefficient_values);
- for (unsigned int component=0; component<n_components; ++component)
- for (unsigned int point=0; point<n_q_points; ++point)
- phi[point][component] *= coefficient_values[point](component);
+ data.coefficients->vector_value_list (fe_face_values.get_quadrature_points(),
+ data.coefficient_values[this_thread]);
+ for (unsigned int component=0; component<data.n_components; ++component)
+ for (unsigned int point=0; point<data.n_q_points; ++point)
+ data.phi[this_thread][point][component] *=
+ data.coefficient_values[this_thread][point](component);
};
};
-
- const vector<double> &JxW_values = fe_face_values.get_JxW_values();
+
+ data.JxW_values[this_thread] = fe_face_values.get_JxW_values();
// take the square of the phi[i]
// for integration, and sum up
double face_integral = 0;
- for (unsigned int component=0; component<n_components; ++component)
- if (component_mask[component] == true)
- for (unsigned int p=0; p<n_q_points; ++p)
- face_integral += sqr(phi[p][component]) *
- JxW_values[p];
+ for (unsigned int component=0; component<data.n_components; ++component)
+ if (data.component_mask[component] == true)
+ for (unsigned int p=0; p<data.n_q_points; ++p)
+ face_integral += sqr(data.phi[this_thread][p][component]) *
+ data.JxW_values[this_thread][p];
- face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
+ data.face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
};
for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
++subface_no)
{
- Assert (face_integrals.find(face->child(subface_no)) !=
- face_integrals.end(),
+ Assert (data.face_integrals.find(face->child(subface_no)) !=
+ data.face_integrals.end(),
+ ExcInternalError());
+ Assert (data.face_integrals[face->child(subface_no)]>=0,
ExcInternalError());
- sum += face_integrals[face->child(subface_no)];
+ sum += data.face_integrals[face->child(subface_no)];
};
- face_integrals[face] = sum;
+
+ data.face_integrals[face] = sum;
+
};
// explicit instantiations
+
template class KellyErrorEstimator<deal_II_dimension>;
projects / dealii-svn.git / commitdiff