From: kronbichler Date: Thu, 22 Aug 2013 02:05:32 +0000 (+0000) Subject: Split up error_estimator into two files, improve the generic face estimate part so... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b2afbba7d490438d920f69fc673b4093e259e13;p=dealii-svn.git Split up error_estimator into two files, improve the generic face estimate part so it compiles quicker, too. git-svn-id: https://svn.dealii.org/trunk@30400 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/dofs/dof_accessor.templates.h b/deal.II/include/deal.II/dofs/dof_accessor.templates.h index 8006dba0d0..989591f843 100644 --- a/deal.II/include/deal.II/dofs/dof_accessor.templates.h +++ b/deal.II/include/deal.II/dofs/dof_accessor.templates.h @@ -3492,6 +3492,7 @@ namespace internal template +inline TriaIterator > DoFCellAccessor::face (const unsigned int i) const { @@ -3519,6 +3520,7 @@ DoFCellAccessor::get_dof_indices (std::vector & } + template inline void DoFCellAccessor::get_mg_dof_indices (std::vector &dof_indices) const @@ -3526,6 +3528,8 @@ void DoFCellAccessor::get_mg_dof_indices (std::vector::get_mg_dof_indices (this->level (), dof_indices); } + + template inline void DoFCellAccessor::set_mg_dof_indices (const std::vector &dof_indices) @@ -3534,6 +3538,7 @@ void DoFCellAccessor::set_mg_dof_indices (const std::vector inline void DoFCellAccessor::get_active_or_mg_dof_indices (std::vector &dof_indices) const @@ -3548,6 +3553,7 @@ void DoFCellAccessor::get_active_or_mg_dof_indices (std::vector template +inline void DoFCellAccessor::get_dof_values (const InputVector &values, Vector &local_values) const @@ -3564,6 +3570,7 @@ DoFCellAccessor::get_dof_values (const InputVector &values, template template +inline void DoFCellAccessor::get_dof_values (const InputVector &values, ForwardIterator local_values_begin, @@ -3579,6 +3586,7 @@ DoFCellAccessor::get_dof_values (const InputVector &values, template template +inline void DoFCellAccessor::get_dof_values (const ConstraintMatrix &constraints, const InputVector &values, @@ -3596,6 +3604,7 @@ DoFCellAccessor::get_dof_values (const ConstraintMatrix &constraints, template template +inline void DoFCellAccessor::set_dof_values (const Vector &local_values, OutputVector &values) const diff --git a/deal.II/include/deal.II/numerics/error_estimator.h b/deal.II/include/deal.II/numerics/error_estimator.h index f47c98b30b..6f15d97c2e 100644 --- a/deal.II/include/deal.II/numerics/error_estimator.h +++ b/deal.II/include/deal.II/numerics/error_estimator.h @@ -234,114 +234,59 @@ class KellyErrorEstimator { public: /** - * Implementation of the error - * estimator described above. You - * may give a coefficient, but - * there is a default value which - * denotes the constant - * coefficient with value - * one. The coefficient function - * may either be a scalar one, in - * which case it is used for all - * components of the finite - * element, or a vector-valued - * one with as many components as - * there are in the finite - * element; in the latter case, - * each component is weighted by - * the respective component in - * the coefficient. + * Implementation of the error estimator described above. You may give a + * coefficient, but there is a default value which denotes the constant + * coefficient with value one. The coefficient function may either be a + * scalar one, in which case it is used for all components of the finite + * element, or a vector-valued one with as many components as there are in + * the finite element; in the latter case, each component is weighted by 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 - * (see @ref GlossComponentMask) - * 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 + * (see @ref GlossComponentMask) 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. * - * The @p subdomain_id parameter - * indicates whether we shall compute - * indicators for all cells (in case its - * value is the default, - * numbers::invalid_unsigned_int), - * 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 cell 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. + * The @p subdomain_id parameter indicates whether we shall compute + * indicators for all cells (in case its value is the default, + * numbers::invalid_unsigned_int), 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 cell + * 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. * - * The @p material_id parameter has a similar - * meaning: if not set to its default value - * (which is numbers::invalid_material_id), - * it means that indicators will only be - * computed for cells with this particular - * material id. + * The @p material_id parameter has a similar meaning: if not set to its + * default value (which is numbers::invalid_material_id), it means that + * indicators will only be computed for cells with this particular material + * id. * - * The @p n_threads parameter used to - * indicate the number of threads to be - * used to compute the error - * estimator. This parameter is now - * ignored, with the number of threads - * determined automatically. The - * parameter is retained for - * compatibility with old versions of the - * library. + * The @p n_threads parameter used to indicate the number of threads to be + * used to compute the error estimator. This parameter is now ignored, with + * the number of threads determined automatically. The parameter is retained + * for compatibility with old versions of the library. * - * @note If the DoFHandler object - * given as an argument to this - * function builds on a - * parallel::distributed::Triangulation, - * this function skips - * computations on all cells that - * are not locally owned. In that - * case, the only valid value for - * the subdomain_id argument - * (besides the invalid value) is - * the subdomain id that is - * associated with the currently + * @note If the DoFHandler object given as an argument to this function + * builds on a parallel::distributed::Triangulation, this function skips + * computations on all cells that are not locally owned. In that case, the + * only valid value for the subdomain_id argument (besides the invalid + * value) is the subdomain id that is associated with the currently * processor, as reported by * parallel::distributed::Triangulation::locally_owned_subdomain(). Even - * though nothing is computed on - * cells that we don't locally - * own, the error indicator - * vector must still have a - * length equal to the number of - * active cell in the mesh as - * reported by + * though nothing is computed on cells that we don't locally own, the error + * indicator vector must still have a length equal to the number of active + * cell in the mesh as reported by * parallel::distributed::Triangulation::n_locally_owned_active_cells(). */ template @@ -358,8 +303,7 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Calls the @p estimate - * function, see above, with + * Calls the @p estimate function, see above, with * mapping=MappingQ1@(). */ template @@ -375,31 +319,17 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Same function as above, but - * accepts more than one solution - * vector and returns one error - * vector for each solution - * vector. For the reason of - * existence of this function, - * see the general documentation - * of this class. + * Same function as above, but accepts more than one solution vector and + * returns one error vector for each solution vector. For the reason of + * existence of this function, see the general documentation of this class. * - * Since we do not want to force - * the user of this function to - * copy around their solution - * vectors, the vector of - * solution vectors takes - * pointers to the solutions, - * rather than being a vector of - * vectors. This makes it simpler - * to have the solution vectors - * somewhere in memory, rather - * than to have them collected - * somewhere special. (Note that - * it is not possible to - * construct of vector of - * references, so we had to use a - * vector of pointers.) + * Since we do not want to force the user of this function to copy around + * their solution vectors, the vector of solution vectors takes pointers to + * the solutions, rather than being a vector of vectors. This makes it + * simpler to have the solution vectors somewhere in memory, rather than to + * have them collected somewhere special. (Note that it is not possible to + * construct of vector of references, so we had to use a vector of + * pointers.) */ template static void estimate (const Mapping &mapping, @@ -415,8 +345,7 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Calls the @p estimate - * function, see above, with + * Calls the @p estimate function, see above, with * mapping=MappingQ1@(). */ template @@ -433,10 +362,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const Mapping &mapping, @@ -453,10 +380,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const DH &dof, @@ -472,10 +397,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const Mapping &mapping, @@ -492,10 +415,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const DH &dof, @@ -561,51 +482,28 @@ class KellyErrorEstimator<1,spacedim> { public: /** - * Implementation of the error - * estimator described above. You - * may give a coefficient, but - * there is a default value which - * denotes the constant - * coefficient with value - * one. The coefficient function - * may either be a scalar one, in - * which case it is used for all - * components of the finite - * element, or a vector-valued - * one with as many components as - * there are in the finite - * element; in the latter case, - * each component is weighted by - * the respective component in - * the coefficient. + * Implementation of the error estimator described above. You may give a + * coefficient, but there is a default value which denotes the constant + * coefficient with value one. The coefficient function may either be a + * scalar one, in which case it is used for all components of the finite + * element, or a vector-valued one with as many components as there are in + * the finite element; in the latter case, each component is weighted by 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. All the other parameters - * are as in the general case used for 2d - * and higher. + * 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 - * 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 not presently - * implemented for 1d, but we retain the - * respective parameter for compatibility - * with the function signature in the - * general case. + * The estimator supports multithreading and splits the cells to + * 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 not presently + * implemented for 1d, but we retain the respective parameter for + * compatibility with the function signature in the general case. */ template static void estimate (const Mapping<1,spacedim> &mapping, @@ -621,8 +519,7 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Calls the @p estimate - * function, see above, with + * Calls the @p estimate function, see above, with * mapping=MappingQ1<1>(). */ template @@ -638,31 +535,17 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Same function as above, but - * accepts more than one solution - * vectors and returns one error - * vector for each solution - * vector. For the reason of - * existence of this function, - * see the general documentation - * of this class. + * Same function as above, but accepts more than one solution vectors and + * returns one error vector for each solution vector. For the reason of + * existence of this function, see the general documentation of this class. * - * Since we do not want to force - * the user of this function to - * copy around their solution - * vectors, the vector of - * solution vectors takes - * pointers to the solutions, - * rather than being a vector of - * vectors. This makes it simpler - * to have the solution vectors - * somewhere in memory, rather - * than to have them collected - * somewhere special. (Note that - * it is not possible to - * construct of vector of - * references, so we had to use a - * vector of pointers.) + * Since we do not want to force the user of this function to copy around + * their solution vectors, the vector of solution vectors takes pointers to + * the solutions, rather than being a vector of vectors. This makes it + * simpler to have the solution vectors somewhere in memory, rather than to + * have them collected somewhere special. (Note that it is not possible to + * construct of vector of references, so we had to use a vector of + * pointers.) */ template static void estimate (const Mapping<1,spacedim> &mapping, @@ -678,8 +561,7 @@ public: const types::material_id material_id = numbers::invalid_material_id); /** - * Calls the @p estimate - * function, see above, with + * Calls the @p estimate function, see above, with * mapping=MappingQ1<1>(). */ template @@ -696,10 +578,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const Mapping<1,spacedim> &mapping, @@ -716,10 +596,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const DH &dof, @@ -735,10 +613,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const Mapping<1,spacedim> &mapping, @@ -755,10 +631,8 @@ public: /** - * Equivalent to the set of functions - * above, except that this one takes a - * quadrature collection for hp finite - * element dof handlers. + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp finite element dof handlers. */ template static void estimate (const DH &dof, diff --git a/deal.II/source/dofs/dof_accessor.cc b/deal.II/source/dofs/dof_accessor.cc index 59b8c48334..53b67ac4d1 100644 --- a/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/source/dofs/dof_accessor.cc @@ -117,8 +117,8 @@ get_interpolated_dof_values (const InputVector &values, typename BaseClass::ExcVectorDoesNotMatch()); if (!this->has_children()) - // if this cell has no children: simply - // return the exact values on this cell + // if this cell has no children: simply return the exact values on this + // cell this->get_dof_values (values, interpolated_values); else // otherwise clobber them from the children @@ -128,80 +128,45 @@ get_interpolated_dof_values (const InputVector &values, interpolated_values = 0; - // later on we will have to - // push the values interpolated - // from the child to the mother - // cell into the output - // vector. unfortunately, there - // are two types of elements: - // ones where you add up the - // contributions from the - // different child cells, and - // ones where you overwrite. + // later on we will have to push the values interpolated from the child + // to the mother cell into the output vector. unfortunately, there are + // two types of elements: ones where you add up the contributions from + // the different child cells, and ones where you overwrite. // - // an example for the first is - // piecewise constant (and - // discontinuous) elements, - // where we build the value on - // the coarse cell by averaging - // the values from the cell - // (i.e. by adding up a - // fraction of the values of - // their values) + // an example for the first is piecewise constant (and discontinuous) + // elements, where we build the value on the coarse cell by averaging + // the values from the cell (i.e. by adding up a fraction of the values + // of their values) // - // an example for the latter - // are the usual continuous - // elements. the value on a - // vertex of a coarse cell must - // there be the same, - // irrespective of the adjacent - // cell we are presently on. so - // we always overwrite. in - // fact, we must, since we - // cannot know in advance how - // many neighbors there will - // be, so there is no way to - // compute the average with - // fixed factors + // an example for the latter are the usual continuous elements. the + // value on a vertex of a coarse cell must there be the same, + // irrespective of the adjacent cell we are presently on. so we always + // overwrite. in fact, we must, since we cannot know in advance how many + // neighbors there will be, so there is no way to compute the average + // with fixed factors // - // so we have to find out to - // which type this element - // belongs. the difficulty is: - // the finite element may be a - // composed one, so we can only - // hope to do this for each - // shape function - // individually. in fact, there - // are even weird finite - // elements (for example the - // Raviart-Thomas element) - // which have shape functions - // that are additive (interior - // ones) and others that are - // overwriting (face degrees of - // freedom that need to be - // continuous across the - // face). to avoid checking - // this over and over again, we - // do this once now and cache - // the results + // so we have to find out to which type this element belongs. the + // difficulty is: the finite element may be a composed one, so we can + // only hope to do this for each shape function individually. in fact, + // there are even weird finite elements (for example the Raviart-Thomas + // element) which have shape functions that are additive (interior ones) + // and others that are overwriting (face degrees of freedom that need to + // be continuous across the face). to avoid checking this over and over + // again, we do this once now and cache the results std::vector restriction_is_additive (dofs_per_cell); for (unsigned int i=0; in_children(); ++child) { - // get the values from the present - // child, if necessary by + // get the values from the present child, if necessary by // interpolation itself this->child(child)->get_interpolated_dof_values (values, tmp1); - // interpolate these to the mother - // cell + // interpolate these to the mother cell fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1); - // and add up or set them - // in the output vector + // and add up or set them in the output vector for (unsigned int i=0; i &local_values, typename BaseClass::ExcVectorDoesNotMatch()); if (!this->has_children()) - // if this cell has no children: simply - // set the values on this cell + // if this cell has no children: simply set the values on this cell this->set_dof_values (local_values, values); else // otherwise distribute them to the children @@ -245,12 +209,8 @@ set_dof_values_by_interpolation (const Vector &local_values, Assert (this->child(child)->get_fe().dofs_per_cell == dofs_per_cell, ExcNotImplemented()); - // prolong the given data - // to the present - // cell. FullMatrix only - // wants us to call vmult - // if the matrix size is - // actually non-zero, so + // prolong the given data to the present cell. FullMatrix only wants + // us to call vmult if the matrix size is actually non-zero, so // check that case if (tmp.size() > 0) { diff --git a/deal.II/source/numerics/CMakeLists.txt b/deal.II/source/numerics/CMakeLists.txt index db675fcd7a..7ab970c447 100644 --- a/deal.II/source/numerics/CMakeLists.txt +++ b/deal.II/source/numerics/CMakeLists.txt @@ -29,6 +29,7 @@ SET(_src derivative_approximation.cc dof_output_operator.cc error_estimator.cc + error_estimator_1d.cc fe_field_function.cc matrix_tools.cc histogram.cc @@ -55,6 +56,7 @@ SET(_inst derivative_approximation.inst.in dof_output_operator.inst.in error_estimator.inst.in + error_estimator_1d.inst.in fe_field_function.inst.in matrix_tools.inst.in point_value_history.inst.in diff --git a/deal.II/source/numerics/error_estimator.cc b/deal.II/source/numerics/error_estimator.cc index 46771ecc59..3fb633c4a5 100644 --- a/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/source/numerics/error_estimator.cc @@ -51,24 +51,14 @@ DEAL_II_NAMESPACE_OPEN namespace { - inline - double sqr (const double x) - { - return x*x; - } - - template inline void advance_by_n (CellIterator &cell, const unsigned int n) { - // store a pointer to the end - // iterator, since we can't get at - // it any more once cell is already - // the end iterator (in that case - // dereferencing cell-> triggers an - // assertion) + // store a pointer to the end iterator, since we can't get at it any more + // once cell is already the end iterator (in that case dereferencing + // cell-> triggers an assertion) const CellIterator endc = cell->get_dof_handler().end(); for (unsigned int t=0; ((t struct ParallelData @@ -148,14 +109,12 @@ namespace internal const dealii::hp::FECollection finite_element; /** - * The quadrature formulas to be used for - * the faces. + * The quadrature formulas to be used for the faces. */ const dealii::hp::QCollection face_quadratures; /** - * FEFaceValues objects to integrate over - * the faces of the current and + * FEFaceValues objects to integrate over the faces of the current and * potentially of neighbor cells. */ dealii::hp::FEFaceValues fe_face_values_cell; @@ -163,38 +122,22 @@ namespace internal dealii::hp::FESubfaceValues fe_subface_values; /** - * A vector to store the jump - * of the normal vectors in - * the quadrature points for - * each of the solution - * vectors (i.e. a temporary - * value). This vector is not - * allocated inside the - * functions that use it, but - * rather globally, since - * memory allocation is slow, - * in particular in presence - * of multiple threads where - * synchronisation makes - * things even slower. + * A vector to store the jump of the normal vectors in the quadrature + * points for each of the solution vectors (i.e. a temporary + * value). This vector is not allocated inside the functions that use + * it, but rather globally, since memory allocation is slow, in + * particular in presence of multiple threads where synchronisation + * makes things even slower. */ std::vector > > phi; /** - * A vector for the gradients of - * the finite element function - * on one cell + * 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 third index be the - * component of the finite - * element, and the second - * index the number of the - * quadrature point. The - * first index denotes the - * index of the solution - * vector. + * Let psi be a short name for a grad u_h, where the third + * index be the component of the finite element, and the second index + * the number of the quadrature point. The first index denotes the index + * of the solution vector. */ std::vector > > > psi; @@ -204,43 +147,35 @@ namespace internal std::vector > > > neighbor_psi; /** - * The normal vectors of the finite - * element function on one face + * The normal vectors of the finite element function on one face */ std::vector > normal_vectors; /** - * Two arrays needed for the - * values of coefficients in - * the jumps, if they are - * given. + * Two arrays needed for the values of coefficients in the jumps, if + * they are given. */ std::vector coefficient_values1; std::vector > coefficient_values; /** - * Array for the products of - * Jacobian determinants and - * weights of quadraturs - * points. + * Array for the products of Jacobian determinants and weights of + * quadraturs points. */ std::vector JxW_values; /** - * The subdomain id we are to care - * for. + * The subdomain id we are to care for. */ const types::subdomain_id subdomain_id; /** - * The material id we are to care - * for. + * The material id we are to care for. */ const types::material_id material_id; /** - * Some more references to input data to - * the KellyErrorEstimator::estimate() - * function. + * Some more references to input data to the + * KellyErrorEstimator::estimate() function. */ const typename FunctionMap::type *neumann_bc; const ComponentMask component_mask; @@ -262,10 +197,9 @@ namespace internal const Function *coefficients); /** - * Resize the arrays so that they fit the - * number of quadrature points associated - * with the given finite element index - * into the hp collections. + * Resize the arrays so that they fit the number of quadrature points + * associated with the given finite element index into the hp + * collections. */ void resize (const unsigned int active_fe_index); }; @@ -363,11 +297,9 @@ namespace internal /** - * Copy data from the - * local_face_integrals map of a single - * ParallelData object into a global such - * map. This is the copier stage of a - * WorkStream pipeline. + * Copy data from the local_face_integrals map of a single ParallelData + * object into a global such map. This is the copier stage of a WorkStream + * pipeline. */ template void @@ -375,16 +307,14 @@ namespace internal std::map > &face_integrals) { - // now copy locally computed elements - // into the global map + // now copy locally computed elements into the global map for (typename std::map >::const_iterator p=local_face_integrals.begin(); p!=local_face_integrals.end(); ++p) { - // double check that the - // element does not already - // exists in the global map + // double check that the element does not already exists in the + // global map Assert (face_integrals.find (p->first) == face_integrals.end(), ExcInternalError()); @@ -400,105 +330,40 @@ namespace internal /** - * Actually do the computation on - * a face which has no hanging - * nodes (it is regular), 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. + * Actually do the computation based on the evaluated gradients in + * ParallelData. */ - template - void - integrate_over_regular_face (const std::vector &solutions, - ParallelData ¶llel_data, - std::map > &local_face_integrals, - const typename DH::active_cell_iterator &cell, - const unsigned int face_no, - dealii::hp::FEFaceValues &fe_face_values_cell, - dealii::hp::FEFaceValues &fe_face_values_neighbor) + template + std::vector + integrate_over_face (ParallelData ¶llel_data, + const typename DH::face_iterator &face, + dealii::hp::FEFaceValues &fe_face_values_cell) { const unsigned int dim = DH::dimension; - const typename DH::face_iterator face = cell->face(face_no); - const unsigned int n_q_points = parallel_data.face_quadratures[cell->active_fe_index()].size(), + const unsigned int n_q_points = parallel_data.psi[0].size(), n_components = parallel_data.finite_element.n_components(), - n_solution_vectors = solutions.size(); - - - // initialize data of the restriction - // of this cell to the present face - fe_face_values_cell.reinit (cell, face_no, - cell->active_fe_index()); + n_solution_vectors = parallel_data.psi.size(); - // get gradients of the finite element - // function on this cell - for (unsigned int n=0; nat_boundary() == false) - // internal face; integrate jump - // of gradient across this face { - Assert (cell->neighbor(face_no).state() == IteratorState::valid, - ExcInternalError()); - - const typename DH::active_cell_iterator neighbor = cell->neighbor(face_no); - - // find which number the - // current face has relative to - // the neighboring cell - const unsigned int neighbor_neighbor - = cell->neighbor_of_neighbor (face_no); - Assert (neighbor_neighbor::faces_per_cell, - ExcInternalError()); - - // get restriction of finite element - // function of @p{neighbor} to the - // common face. in the hp case, use the - // quadrature formula that matches the - // one we would use for the present - // cell - fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor, - cell->active_fe_index()); - - // get gradients on neighbor cell + // compute the jump in the gradients for (unsigned int n=0; nn_components == 1) { - parallel_data.coefficients ->value_list (fe_face_values_cell.get_present_fe_values() .get_quadrature_points(), @@ -545,18 +409,15 @@ namespace internal if (face->at_boundary() == true) - // neumann boundary face. compute - // difference between normal - // derivative and boundary function + // neumann boundary face. compute difference between normal derivative + // and boundary function { const types::boundary_id boundary_indicator = face->boundary_indicator(); Assert (parallel_data.neumann_bc->find(boundary_indicator) != parallel_data.neumann_bc->end(), ExcInternalError ()); - // get the values of the boundary - // function at the quadrature - // points + // get the values of the boundary function at the quadrature points if (n_components == 1) { std::vector g(n_q_points); @@ -587,26 +448,96 @@ namespace internal // now phi contains the following: // - for an internal face, phi=[a du/dn] - // - for a neumann boundary face, - // phi=a du/dn-g - // each component being the - // mentioned value at one of the - // quadrature points + // - for a neumann boundary face, phi=a du/dn-g + // each component being the mentioned value at one of the quadrature + // points parallel_data.JxW_values = fe_face_values_cell.get_present_fe_values().get_JxW_values(); - // take the square of the phi[i] - // for integration, and sum up + // take the square of the phi[i] for integration, and sum up std::vector face_integral (n_solution_vectors, 0); for (unsigned int n=0; n(parallel_data.phi[n][p][component]) * parallel_data.JxW_values[p]; - local_face_integrals[face] = face_integral; + return face_integral; + } + + + + /** + * Actually do the computation on a face which has no hanging nodes (it is + * regular), 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. + */ + template + void + integrate_over_regular_face (const std::vector &solutions, + ParallelData ¶llel_data, + std::map > &local_face_integrals, + const typename DH::active_cell_iterator &cell, + const unsigned int face_no, + dealii::hp::FEFaceValues &fe_face_values_cell, + dealii::hp::FEFaceValues &fe_face_values_neighbor) + { + const unsigned int dim = DH::dimension; + + const typename DH::face_iterator face = cell->face(face_no); + const unsigned int n_q_points = parallel_data.face_quadratures[cell->active_fe_index()].size(), + n_components = parallel_data.finite_element.n_components(), + n_solution_vectors = solutions.size(); + + + // initialize data of the restriction + // of this cell to the present face + fe_face_values_cell.reinit (cell, face_no, + cell->active_fe_index()); + + // get gradients of the finite element + // function on this cell + for (unsigned int n=0; nat_boundary() == false) + // internal face; integrate jump of gradient across this face + { + Assert (cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + + const typename DH::active_cell_iterator neighbor = cell->neighbor(face_no); + + // find which number the current face has relative to the + // neighboring cell + const unsigned int neighbor_neighbor + = cell->neighbor_of_neighbor (face_no); + Assert (neighbor_neighbor::faces_per_cell, + ExcInternalError()); + + // get restriction of finite element function of @p{neighbor} to the + // common face. in the hp case, use the quadrature formula that + // matches the one we would use for the present cell + fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor, + cell->active_fe_index()); + + // get gradients on neighbor cell + for (unsigned int n=0; n void - integrate_over_irregular_face (const std::vector &solutions, - ParallelData ¶llel_data, + integrate_over_irregular_face (const std::vector &solutions, + ParallelData ¶llel_data, std::map > &local_face_integrals, const typename DH::active_cell_iterator &cell, const unsigned int face_no, @@ -643,20 +574,15 @@ namespace internal Assert (neighbor.state() == IteratorState::valid, ExcInternalError()); Assert (face->has_children(), ExcInternalError()); - // set up a vector of the gradients - // of the finite element function - // on this cell at the quadrature - // points + + // 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 - - // store which number @p{cell} has - // in the list of neighbors of + // 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 + + // store which number @p{cell} has in the list of neighbors of // @p{neighbor} const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no); @@ -666,138 +592,38 @@ namespace internal // loop over all subfaces for (unsigned int subface_no=0; subface_non_children(); ++subface_no) { - // get an iterator pointing to the - // cell behind the present subface + // get an iterator pointing to the cell behind the present subface const typename DH::active_cell_iterator neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); Assert (!neighbor_child->has_children(), ExcInternalError()); - // restrict the finite element - // on the present cell to the - // subface + // restrict the finite element on the present cell to the subface fe_subface_values.reinit (cell, face_no, subface_no, cell->active_fe_index()); - // restrict the finite element - // on the neighbor cell to the - // common @p{subface}. + // restrict the finite element on the neighbor cell to the common + // @p{subface}. fe_face_values.reinit (neighbor_child, neighbor_neighbor, cell->active_fe_index()); - // store the gradient of the - // solution in psi + // store the gradient of the solution in psi for (unsigned int n=0; nn_components == 1) - { - parallel_data.coefficients - ->value_list (fe_face_values.get_present_fe_values() - .get_quadrature_points(), - parallel_data.coefficient_values1); - for (unsigned int n=0; nvector_value_list (fe_face_values.get_present_fe_values() - .get_quadrature_points(), - parallel_data.coefficient_values); - for (unsigned int n=0; n face_integral (n_solution_vectors, 0); - for (unsigned int n=0; nface(neighbor_neighbor)] - = face_integral; + // call generic evaluate function + local_face_integrals[neighbor_child->face(neighbor_neighbor)] = + integrate_over_face (parallel_data, face, fe_face_values); } - - // finally loop over all subfaces to - // collect the contributions of the - // subfaces and store them with the - // mother face + // finally loop over all subfaces to collect the contributions of the + // subfaces and store them with the mother face std::vector sum (n_solution_vectors, 0); for (unsigned int subface_no=0; subface_non_children(); ++subface_no) { @@ -816,29 +642,26 @@ namespace internal /** - * Computate the error on the faces of a - * single cell. + * Computate the error on the faces of a single cell. * - * This function is only needed - * in two or three dimensions. - * The error estimator in one - * dimension is implemented - * seperatly. + * This function is only needed in two or three dimensions. The error + * estimator in one dimension is implemented separately. */ - template + template void estimate_one_cell (const typename DH::active_cell_iterator &cell, ParallelData ¶llel_data, std::map > &local_face_integrals, const std::vector &solutions) { + const unsigned int dim = DH::dimension; + const unsigned int spacedim = DH::dimension; const unsigned int n_solution_vectors = solutions.size(); const types::subdomain_id subdomain_id = parallel_data.subdomain_id; const unsigned int material_id = parallel_data.material_id; - // empty our own copy of the local face - // integrals + // empty our own copy of the local face integrals local_face_integrals.clear(); // loop over all faces of this cell @@ -848,20 +671,11 @@ namespace internal const typename DH::face_iterator face=cell->face(face_no); - // make sure we do work - // only once: this face - // may either be regular - // or irregular. if it is - // regular and has a - // neighbor, then we - // visit the face twice, - // once from every - // side. let the one with - // the lower index do the - // work. if it is at the - // boundary, or if the - // face is irregular, - // then do the work below + // make sure we do work only once: this face may either be regular + // or irregular. if it is regular and has a neighbor, then we visit + // the face twice, once from every side. let the one with the lower + // index do the work. if it is at the boundary, or if the face is + // irregular, then do the work below if ((face->has_children() == false) && !cell->at_boundary(face_no) && (!cell->neighbor_is_coarser(face_no) && @@ -870,24 +684,17 @@ namespace internal cell->neighbor(face_no)->level() < cell->level())))) 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 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 (face->at_boundary() == false) if (cell->neighbor_is_coarser(face_no)) 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. if (face->at_boundary() && (parallel_data.neumann_bc->find(face->boundary_indicator()) == @@ -898,12 +705,9 @@ namespace internal continue; } - // 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 + // 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 == numbers::invalid_subdomain_id) || @@ -913,11 +717,8 @@ namespace internal || (cell->material_id() == material_id))) ) { - // ok, cell is unwanted, but - // maybe its neighbor behind - // the face we presently work - // on? oh is there a face at - // all? + // ok, cell is unwanted, but maybe its neighbor behind the face + // we presently work on? oh is there a face at all? if (face->at_boundary()) continue; @@ -950,32 +751,24 @@ namespace internal } } - // so if none of the neighbors - // cares for this subdomain or - // material either, then try - // next face + // so if none of the neighbors cares for this subdomain or + // material 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. first re-size the - // arrays we may use to the correct + // so now we know that we care for this face, let's do something + // about it. first re-size the arrays we may use to the correct // size: parallel_data.resize (cell->active_fe_index()); // then do the actual integration if (face->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 + // 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 (solutions, parallel_data, local_face_integrals, @@ -984,10 +777,8 @@ namespace internal parallel_data.fe_face_values_neighbor); else - // otherwise we need to do some - // special computations which do - // not fit into the framework of - // the above function + // otherwise we need to do some special computations which do not + // fit into the framework of the above function integrate_over_irregular_face (solutions, parallel_data, local_face_integrals, @@ -1003,418 +794,6 @@ namespace internal - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const Mapping<1,spacedim> &mapping, - const DH &dof_handler, - const Quadrature<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const InputVector &solution, - Vector &error, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - // just pass on to the other function - const std::vector solutions (1, &solution); - std::vector*> errors (1, &error); - estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const DH &dof_handler, - const Quadrature<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const InputVector &solution, - Vector &error, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const DH &dof_handler, - const Quadrature<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const std::vector &solutions, - std::vector*> &errors, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const Mapping<1,spacedim> &mapping, - const DH &dof_handler, - const hp::QCollection<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const InputVector &solution, - Vector &error, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - // just pass on to the other function - const std::vector solutions (1, &solution); - std::vector*> errors (1, &error); - estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const DH &dof_handler, - const hp::QCollection<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const InputVector &solution, - Vector &error, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - - -template -template -void -KellyErrorEstimator<1,spacedim>:: -estimate (const DH &dof_handler, - const hp::QCollection<0> &quadrature, - const typename FunctionMap::type &neumann_bc, - const std::vector &solutions, - std::vector*> &errors, - const ComponentMask &component_mask, - const Function *coefficients, - const unsigned int n_threads, - const types::subdomain_id subdomain_id, - const types::material_id material_id) -{ - estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id, material_id); -} - - - - -template -template -void KellyErrorEstimator<1,spacedim>:: -estimate (const Mapping<1,spacedim> &/*mapping*/, - const DH &/*dof_handler*/, - const hp::QCollection<0> &, - const typename FunctionMap::type &/*neumann_bc*/, - const std::vector &/*solutions*/, - std::vector*> &/*errors*/, - const ComponentMask &/*component_mask_*/, - const Function */*coefficient*/, - const unsigned int, - const types::subdomain_id /*subdomain_id*/, - const types::material_id /*material_id*/) -{ - Assert (false, ExcInternalError()); -} - - - -template -template -void KellyErrorEstimator<1,spacedim>:: -estimate (const Mapping<1,spacedim> &mapping, - const DH &dof_handler, - const Quadrature<0> &, - const typename FunctionMap::type &neumann_bc, - const std::vector &solutions, - std::vector*> &errors, - const ComponentMask &component_mask, - const Function *coefficient, - const unsigned int, - const types::subdomain_id subdomain_id_, - const types::material_id material_id) -{ -#ifdef DEAL_II_WITH_P4EST - if (dynamic_cast*> - (&dof_handler.get_tria()) - != 0) - Assert ((subdomain_id_ == numbers::invalid_subdomain_id) - || - (subdomain_id_ == - dynamic_cast&> - (dof_handler.get_tria()).locally_owned_subdomain()), - ExcMessage ("For parallel distributed triangulations, the only " - "valid subdomain_id that can be passed here is the " - "one that corresponds to the locally owned subdomain id.")); - - const types::subdomain_id subdomain_id - = ((dynamic_cast*> - (&dof_handler.get_tria()) - != 0) - ? - dynamic_cast&> - (dof_handler.get_tria()).locally_owned_subdomain() - : - subdomain_id_); -#else - const types::subdomain_id subdomain_id - = subdomain_id_; -#endif - - const unsigned int n_components = dof_handler.get_fe().n_components(); - const unsigned int n_solution_vectors = solutions.size(); - - // sanity checks - Assert (neumann_bc.find(numbers::internal_face_boundary_id) == neumann_bc.end(), - ExcInvalidBoundaryIndicator()); - - for (typename FunctionMap::type::const_iterator i=neumann_bc.begin(); - i!=neumann_bc.end(); ++i) - Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); - - Assert (component_mask.represents_n_components(n_components), - ExcInvalidComponentMask()); - Assert (component_mask.n_selected_components(n_components) > 0, - ExcInvalidComponentMask()); - - Assert ((coefficient == 0) || - (coefficient->n_components == n_components) || - (coefficient->n_components == 1), - ExcInvalidCoefficient()); - - Assert (solutions.size() > 0, - ExcNoSolutions()); - Assert (solutions.size() == errors.size(), - ExcIncompatibleNumberOfElements(solutions.size(), errors.size())); - for (unsigned int n=0; nsize() == dof_handler.n_dofs(), - ExcInvalidSolutionVector()); - - Assert ((coefficient == 0) || - (coefficient->n_components == n_components) || - (coefficient->n_components == 1), - ExcInvalidCoefficient()); - - for (typename FunctionMap::type::const_iterator i=neumann_bc.begin(); - i!=neumann_bc.end(); ++i) - Assert (i->second->n_components == n_components, - ExcInvalidBoundaryFunction()); - - // reserve one slot for each cell and set - // it to zero - for (unsigned int n=0; n > > > - gradients_here (n_solution_vectors, - std::vector > >(2, std::vector >(n_components))); - std::vector > > > - gradients_neighbor (gradients_here); - std::vector > - grad_neighbor (n_solution_vectors, Vector(n_components)); - - // reserve some space for - // coefficient values at one point. - // if there is no coefficient, then - // we fill it by unity once and for - // all and don't set it any more - Vector coefficient_values (n_components); - if (coefficient == 0) - for (unsigned int c=0; c quadrature; - const hp::QCollection<1> q_collection(quadrature); - - const hp::FECollection<1,spacedim> fe (dof_handler.get_fe()); - - hp::MappingCollection<1,spacedim> mapping_collection; - mapping_collection.push_back (mapping); - - hp::FEValues<1,spacedim> fe_values (mapping_collection, fe, q_collection, - 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. - typename DH::active_cell_iterator cell = dof_handler.begin_active(); - for (unsigned int cell_index=0; cell != dof_handler.end(); - ++cell, ++cell_index) - if (((subdomain_id == numbers::invalid_subdomain_id) - || - (cell->subdomain_id() == subdomain_id)) - && - ((material_id == numbers::invalid_material_id) - || - (cell->material_id() == material_id))) - { - for (unsigned int n=0; nneighbor(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; ssecond->value(cell->vertex(0)); - - for (unsigned int s=0; s v(n_components); - neumann_bc.find(n)->second->vector_value(cell->vertex(0), v); - - for (unsigned int s=0; sn_components == 1) - { - const double c_value = coefficient->value (cell->vertex(n)); - for (unsigned int c=0; cvector_value(cell->vertex(n), - coefficient_values); - } - - - for (unsigned int s=0; sdiameter(); - } - } - - for (unsigned int s=0; s @@ -1581,20 +960,14 @@ estimate (const Mapping &mapping, const unsigned int n_solution_vectors = solutions.size(); - // Map of integrals indexed by - // the corresponding face. In this map - // we store the integrated jump of the - // gradient for each face. - // At the end of the function, we again - // loop over the cells and collect the - // contributions of the different faces - // of the cell. + // Map of integrals indexed by the corresponding face. In this map we store + // the integrated jump of the gradient for each face. At the end of the + // function, we again loop over the cells and collect the contributions of + // the different faces of the cell. std::map > face_integrals; - // all the data needed in the error - // estimator by each of the threads - // is gathered in the following - // stuctures + // all the data needed in the error estimator by each of the threads is + // gathered in the following structures const hp::MappingCollection mapping_collection(mapping); const internal::ParallelData parallel_data (dof_handler.get_fe(), @@ -1612,18 +985,16 @@ estimate (const Mapping &mapping, // now let's work on all those cells: WorkStream::run (dof_handler.begin_active(), static_cast(dof_handler.end()), - std_cxx1x::bind (&internal::estimate_one_cell, + std_cxx1x::bind (&internal::estimate_one_cell, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, std_cxx1x::ref(solutions)), std_cxx1x::bind (&internal::copy_local_to_global, std_cxx1x::_1, std_cxx1x::ref(face_integrals)), parallel_data, sample_local_face_integrals); - // finally add up the contributions of the - // faces for each cell + // finally add up the contributions of the faces for each cell - // reserve one slot for each cell and set - // it to zero + // reserve one slot for each cell and set it to zero for (unsigned int n=0; n &mapping, (*errors[n])(i)=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 + // 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 unsigned int present_cell=0; for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); @@ -1658,9 +1027,8 @@ estimate (const Mapping &mapping, for (unsigned int n=0; nface(face_no)][n] >= 0, ExcInternalError()); diff --git a/deal.II/source/numerics/error_estimator.inst.in b/deal.II/source/numerics/error_estimator.inst.in index 34b06e672b..8ca601fbcf 100644 --- a/deal.II/source/numerics/error_estimator.inst.in +++ b/deal.II/source/numerics/error_estimator.inst.in @@ -24,7 +24,7 @@ template class KellyErrorEstimator; for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; DH : DOFHANDLER_TEMPLATES ) { -#if deal_II_dimension <= deal_II_space_dimension +#if deal_II_dimension != 1 && deal_II_dimension <= deal_II_space_dimension template void diff --git a/deal.II/source/numerics/error_estimator_1d.cc b/deal.II/source/numerics/error_estimator_1d.cc new file mode 100644 index 0000000000..15b65d0bc8 --- /dev/null +++ b/deal.II/source/numerics/error_estimator_1d.cc @@ -0,0 +1,439 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const Mapping<1,spacedim> &mapping, + const DH &dof_handler, + const Quadrature<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const InputVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + // just pass on to the other function + const std::vector solutions (1, &solution); + std::vector*> errors (1, &error); + estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, + component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const DH &dof_handler, + const Quadrature<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const InputVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, + error, component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const DH &dof_handler, + const Quadrature<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const std::vector &solutions, + std::vector*> &errors, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, + errors, component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const Mapping<1,spacedim> &mapping, + const DH &dof_handler, + const hp::QCollection<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const InputVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + // just pass on to the other function + const std::vector solutions (1, &solution); + std::vector*> errors (1, &error); + estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, + component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const DH &dof_handler, + const hp::QCollection<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const InputVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, + error, component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + + +template +template +void +KellyErrorEstimator<1,spacedim>:: +estimate (const DH &dof_handler, + const hp::QCollection<0> &quadrature, + const typename FunctionMap::type &neumann_bc, + const std::vector &solutions, + std::vector*> &errors, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id) +{ + estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, + errors, component_mask, coefficients, n_threads, subdomain_id, material_id); +} + + + + +template +template +void KellyErrorEstimator<1,spacedim>:: +estimate (const Mapping<1,spacedim> &/*mapping*/, + const DH &/*dof_handler*/, + const hp::QCollection<0> &, + const typename FunctionMap::type &/*neumann_bc*/, + const std::vector &/*solutions*/, + std::vector*> &/*errors*/, + const ComponentMask &/*component_mask_*/, + const Function */*coefficient*/, + const unsigned int, + const types::subdomain_id /*subdomain_id*/, + const types::material_id /*material_id*/) +{ + Assert (false, ExcInternalError()); +} + + + +template +template +void KellyErrorEstimator<1,spacedim>:: +estimate (const Mapping<1,spacedim> &mapping, + const DH &dof_handler, + const Quadrature<0> &, + const typename FunctionMap::type &neumann_bc, + const std::vector &solutions, + std::vector*> &errors, + const ComponentMask &component_mask, + const Function *coefficient, + const unsigned int, + const types::subdomain_id subdomain_id_, + const types::material_id material_id) +{ +#ifdef DEAL_II_WITH_P4EST + if (dynamic_cast*> + (&dof_handler.get_tria()) + != 0) + Assert ((subdomain_id_ == numbers::invalid_subdomain_id) + || + (subdomain_id_ == + dynamic_cast&> + (dof_handler.get_tria()).locally_owned_subdomain()), + ExcMessage ("For parallel distributed triangulations, the only " + "valid subdomain_id that can be passed here is the " + "one that corresponds to the locally owned subdomain id.")); + + const types::subdomain_id subdomain_id + = ((dynamic_cast*> + (&dof_handler.get_tria()) + != 0) + ? + dynamic_cast&> + (dof_handler.get_tria()).locally_owned_subdomain() + : + subdomain_id_); +#else + const types::subdomain_id subdomain_id + = subdomain_id_; +#endif + + const unsigned int n_components = dof_handler.get_fe().n_components(); + const unsigned int n_solution_vectors = solutions.size(); + + // sanity checks + Assert (neumann_bc.find(numbers::internal_face_boundary_id) == neumann_bc.end(), + ExcInvalidBoundaryIndicator()); + + for (typename FunctionMap::type::const_iterator i=neumann_bc.begin(); + i!=neumann_bc.end(); ++i) + Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); + + Assert (component_mask.represents_n_components(n_components), + ExcInvalidComponentMask()); + Assert (component_mask.n_selected_components(n_components) > 0, + ExcInvalidComponentMask()); + + Assert ((coefficient == 0) || + (coefficient->n_components == n_components) || + (coefficient->n_components == 1), + ExcInvalidCoefficient()); + + Assert (solutions.size() > 0, + ExcNoSolutions()); + Assert (solutions.size() == errors.size(), + ExcIncompatibleNumberOfElements(solutions.size(), errors.size())); + for (unsigned int n=0; nsize() == dof_handler.n_dofs(), + ExcInvalidSolutionVector()); + + Assert ((coefficient == 0) || + (coefficient->n_components == n_components) || + (coefficient->n_components == 1), + ExcInvalidCoefficient()); + + for (typename FunctionMap::type::const_iterator i=neumann_bc.begin(); + i!=neumann_bc.end(); ++i) + Assert (i->second->n_components == n_components, + ExcInvalidBoundaryFunction()); + + // reserve one slot for each cell and set it to zero + for (unsigned int n=0; n > > > + gradients_here (n_solution_vectors, + std::vector > >(2, std::vector >(n_components))); + std::vector > > > + gradients_neighbor (gradients_here); + std::vector > + grad_neighbor (n_solution_vectors, Vector(n_components)); + + // reserve some space for coefficient values at one point. if there is no + // coefficient, then we fill it by unity once and for all and don't set it + // any more + Vector coefficient_values (n_components); + if (coefficient == 0) + for (unsigned int c=0; c quadrature; + const hp::QCollection<1> q_collection(quadrature); + + const hp::FECollection<1,spacedim> fe (dof_handler.get_fe()); + + hp::MappingCollection<1,spacedim> mapping_collection; + mapping_collection.push_back (mapping); + + hp::FEValues<1,spacedim> fe_values (mapping_collection, fe, q_collection, + 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. + typename DH::active_cell_iterator cell = dof_handler.begin_active(); + for (unsigned int cell_index=0; cell != dof_handler.end(); + ++cell, ++cell_index) + if (((subdomain_id == numbers::invalid_subdomain_id) + || + (cell->subdomain_id() == subdomain_id)) + && + ((material_id == numbers::invalid_material_id) + || + (cell->material_id() == material_id))) + { + for (unsigned int n=0; nneighbor(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; ssecond->value(cell->vertex(0)); + + for (unsigned int s=0; s v(n_components); + neumann_bc.find(n)->second->vector_value(cell->vertex(0), v); + + for (unsigned int s=0; sn_components == 1) + { + const double c_value = coefficient->value (cell->vertex(n)); + for (unsigned int c=0; cvector_value(cell->vertex(n), + coefficient_values); + } + + + for (unsigned int s=0; sdiameter(); + } + } + + for (unsigned int s=0; s:: +estimate > (const Mapping &, + const DH &, + const Quadrature &, + const FunctionMap::type &, + const VEC &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > ( + const DH &, + const Quadrature &, + const FunctionMap::type &, + const VEC &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > (const Mapping &, + const DH &, + const Quadrature &, + const FunctionMap::type &, + const std::vector &, + std::vector *> &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > ( + const DH &, + const Quadrature &, + const FunctionMap::type &, + const std::vector &, + std::vector *> &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > (const Mapping &, + const DH &, + const hp::QCollection &, + const FunctionMap::type &, + const VEC &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > ( + const DH &, + const hp::QCollection &, + const FunctionMap::type &, + const VEC &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > (const Mapping &, + const DH &, + const hp::QCollection &, + const FunctionMap::type &, + const std::vector &, + std::vector *> &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +template +void +KellyErrorEstimator:: +estimate > ( + const DH &, + const hp::QCollection &, + const FunctionMap::type &, + const std::vector &, + std::vector *> &, + const ComponentMask &, + const Function *, + const unsigned int , + const unsigned int , + const types::material_id); + +#endif +}