From 5f171788ad191a11a186cc757a936b8a076b0659 Mon Sep 17 00:00:00 2001 From: richter Date: Tue, 1 Feb 2000 10:35:37 +0000 Subject: [PATCH] Implement multithreaded version of the error estimator. git-svn-id: https://svn.dealii.org/trunk@2302 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/error_estimator.h | 171 ++++- .../source/numerics/error_estimator.cc | 692 +++++++++++------- 2 files changed, 574 insertions(+), 289 deletions(-) diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 69018e41e5..8eeb90c952 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -12,7 +12,21 @@ #include #include - + // 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 /** @@ -172,7 +186,8 @@ template class KellyErrorEstimator { public: - /** + + /** * Declare a data type which denotes a * mapping between a boundary indicator * and the function denoting the boundary @@ -180,9 +195,10 @@ class KellyErrorEstimator { * Only one boundary function may be given * for each boundary indicator, which is * guaranteed by the #map# data type. - */ + */ typedef map*> FunctionMap; + /** * Implementation of the error * estimator described above. You @@ -219,14 +235,28 @@ class KellyErrorEstimator { * 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 &dof, const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const vector &component_mask = vector(), - const Function *coefficients = 0); + const vector &component_mask_ = vector(), + const Function *coefficients = 0, + unsigned int n_threads=N_THREADS); /** * Exception @@ -245,8 +275,12 @@ class KellyErrorEstimator { */ DeclException0 (ExcInvalidBoundaryFunction); + + private: - /** + + + /** * Declare a data type to represent the * mapping between faces and integrated * jumps of gradients. See the general @@ -254,12 +288,107 @@ class KellyErrorEstimator { */ typedef map::face_iterator,double> FaceIntegrals; + /** * Redeclare an active cell iterator. * This is simply for convenience. */ typedef DoFHandler::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 &dof; + const Quadrature &quadrature; + const FunctionMap &neumann_bc; + const Vector &solution; + vector component_mask; + const Function *coefficients; + unsigned int n_threads; + + DoFHandler::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 > > 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 > > > psi; + /** + * The same vector for a neighbor cell + */ + vector< vector > > > neighbor_psi; + /** + * The normal vectors of the finite + * element function on one face + */ + vector< vector > > normal_vectors; + + vector< vector > coefficient_values1; + vector< vector > > coefficient_values; + + vector< vector > JxW_values; + + /** + * A constructor of the + * class Data. All variables are + * passed as references. + */ + Data(const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const Vector &solution, + vector component_mask_, + const Function *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 @@ -278,17 +407,16 @@ class KellyErrorEstimator { * 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 &fe_face_values_cell, - FEFaceValues &fe_face_values_neighbor, - FaceIntegrals &face_integrals, - const Vector&solution, - const vector &component_mask, - const Function *coefficient); - + FEFaceValues &fe_face_values_neighbor); + + /** * The same applies as for the function * above, except that integration is @@ -297,15 +425,12 @@ class KellyErrorEstimator { * 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 &fe_face_values, - FESubfaceValues &fe_subface_values, - FaceIntegrals &face_integrals, - const Vector &solution, - const vector &component_mask, - const Function *coefficient); + FESubfaceValues &fe_subface_values); }; @@ -316,3 +441,5 @@ class KellyErrorEstimator { /* end of #ifndef __error_estimator_H */ #endif /*---------------------------- error_estimator.h ---------------------------*/ + + diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 0a9be5f7dd..f4ca76d363 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -17,6 +17,108 @@ #include #include #include +#include + + // if multithreaded include + // ThreadManager +#ifdef DEAL_II_USE_MT +#include +#endif + +#if deal_II_dimension == 1 + +template <> +KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, + const Quadrature<0> &, + const FunctionMap &, + const Vector &, + vector , + const Function<1> *, + unsigned int ): + dof(* static_cast *> (0)), + quadrature(* static_cast *> (0)), + neumann_bc(* static_cast (0)), + solution(* static_cast *> (0)) +{ + Assert (false, ExcInternalError()); +} + +#else + +template +KellyErrorEstimator::Data::Data(const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const Vector &solution, + vector component_mask_, + const Function *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(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 +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> &, @@ -34,7 +142,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, const Vector &solution, Vector &error, const vector &component_mask_, - const Function<1> *coefficient) + const Function<1> *coefficient, + unsigned int) { const unsigned int n_components = dof.get_fe().n_components(); @@ -174,56 +283,14 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, }; }; -#endif + // #if deal_II_dimension !=1 +#else template -void KellyErrorEstimator::estimate (const DoFHandler &dof, - const Quadrature &quadrature, - const FunctionMap &neumann_bc, - const Vector &solution, - Vector &error, - const vector &component_mask_, - const Function *coefficient) +void * KellyErrorEstimator::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 component_mask ((component_mask_.size() == 0) ? - vector(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 @@ -235,124 +302,235 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, // need not compute all values on the // neighbor cells, so using two objects // gives us a performance gain). - FEFaceValues fe_face_values_cell (dof.get_fe(), quadrature, + FEFaceValues 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 fe_face_values_neighbor (dof.get_fe(), quadrature, update_gradients); - FESubfaceValues fe_subface_values (dof.get_fe(), quadrature, update_gradients); + FEFaceValues fe_face_values_neighbor (data.dof.get_fe(), + data.quadrature, + update_gradients); + FESubfaceValues 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::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::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::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 +void KellyErrorEstimator::estimate (const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const Vector &solution, + Vector &error, + const vector &component_mask, + const Function *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::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::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 + FunData; + + // One struct of this type for every thread + vector + fun_data (data.n_threads, + FunData (data,0,&KellyErrorEstimator::estimate_some)); + + + // get start cells for each thread + for (unsigned int l=0;l::faces_per_cell; ++face_no) + for (unsigned int face_no=0; face_no::faces_per_cell; + ++face_no) { - Assert (face_integrals.find(cell->face(face_no)) != - face_integrals.end(), - ExcInternalError()); - error(present_cell) += (face_integrals[cell->face(face_no)] * + 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 &, - const vector &, - const Function<1> *) { + FEFaceValues<1> &) +{ Assert (false, ExcInternalError()); }; @@ -360,15 +538,13 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iter 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 &, - const vector &, - const Function<1> *) { + FESubfaceValues<1> &) +{ Assert (false, ExcInternalError()); }; @@ -378,37 +554,22 @@ integrate_over_irregular_face (const active_cell_iterator &, template void KellyErrorEstimator:: -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 &fe_face_values_cell, - FEFaceValues &fe_face_values_neighbor, - FaceIntegrals &face_integrals, - const Vector &solution, - const vector &component_mask, - const Function *coefficient) + FEFaceValues &fe_face_values_neighbor) { - const unsigned int n_components = component_mask.size(); const DoFHandler::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 > > psi(n_q_points, vector >(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 @@ -435,18 +596,15 @@ integrate_over_regular_face (const active_cell_iterator &cell, // 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 > > neighbor_psi (n_q_points, - vector >(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 > phi(n_q_points, vector(n_components)); - const vector > &normal_vectors(fe_face_values_cell. - get_normal_vectors()); - - for (unsigned int component=0; componentn_components == 1) + if (data.coefficients->n_components == 1) { - vector coefficient_values (n_q_points); - coefficient->value_list (fe_face_values_cell.get_quadrature_points(), - coefficient_values); - for (unsigned int component=0; componentvalue_list (fe_face_values_cell.get_quadrature_points(), + data.coefficient_values1[this_thread]); + for (unsigned int component=0; component > coefficient_values (n_q_points, - Vector(n_components)); - coefficient->vector_value_list (fe_face_values_cell.get_quadrature_points(), - coefficient_values); - for (unsigned int component=0; componentvector_value_list (fe_face_values_cell.get_quadrature_points(), + data.coefficient_values[this_thread]); + for (unsigned int component=0; componentat_boundary() == true) // neumann boundary face. compute @@ -511,20 +667,21 @@ integrate_over_regular_face (const active_cell_iterator &cell, // 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 > g(n_q_points, Vector(n_components)); - neumann_bc.find(boundary_indicator)->second + + vector > g(data.n_q_points, Vector(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 &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 void KellyErrorEstimator:: -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 &fe_face_values, - FESubfaceValues &fe_subface_values, - FaceIntegrals &face_integrals, - const Vector &solution, - const vector &component_mask, - const Function *coefficient) + FESubfaceValues &fe_subface_values) { - const unsigned int n_components = component_mask.size(); - const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); @@ -580,8 +733,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // finite element, and the first // index the number of the // quadrature point - vector > > psi(n_q_points, vector >(n_components)); - + // store which number #cell# has in the // list of neighbors of #neighbor# unsigned int neighbor_neighbor; @@ -606,26 +758,26 @@ integrate_over_irregular_face (const active_cell_iterator &cell, Assert (!neighbor->child(GeometryInfo:: 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 > > neighbor_psi (n_q_points, - vector >(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 > phi(n_q_points, vector(n_components)); - const vector > &normal_vectors(fe_face_values. - get_normal_vectors()); - - for (unsigned int component=0; componentn_components == 1) + if (data.coefficients->n_components == 1) { - vector coefficient_values (n_q_points); - coefficient->value_list (fe_face_values.get_quadrature_points(), - coefficient_values); - for (unsigned int component=0; componentvalue_list (fe_face_values.get_quadrature_points(), + data.coefficient_values1[this_thread]); + for (unsigned int component=0; component > coefficient_values (n_q_points, - Vector(n_components)); - coefficient->vector_value_list (fe_face_values.get_quadrature_points(), - coefficient_values); - for (unsigned int component=0; componentvector_value_list (fe_face_values.get_quadrature_points(), + data.coefficient_values[this_thread]); + for (unsigned int component=0; component &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; componentface(neighbor_neighbor)] = face_integral; + data.face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral; }; @@ -704,16 +857,21 @@ integrate_over_irregular_face (const active_cell_iterator &cell, for (unsigned int subface_no=0; subface_no::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; -- 2.39.5