From 222528b56d0fbc3d646e01f55abf0bd315428740 Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 6 May 1999 18:32:35 +0000 Subject: [PATCH] Make the Kelly error estimator aware of problems with finite elements composed of more than only one component. git-svn-id: https://svn.dealii.org/trunk@1290 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/error_estimator.h | 30 ++++++++- .../source/numerics/error_estimator.cc | 63 +++++++++++++++++-- 2 files changed, 84 insertions(+), 9 deletions(-) diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 0d09f52c55..144aed5e80 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -165,13 +165,25 @@ class KellyErrorEstimator { * coefficient, but there is a default * value which denotes the constant * coefficient with value one. + * + * You must give the component if the + * finite element in use by the #dof# + * object has more than one component. + * This number shall be between zero + * and the number of components within + * the finite element. If the finite + * element has only one component, + * then the parameter selecting the + * component shall be zero, which is + * also the default value. */ static void estimate (const DoFHandler &dof, const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const Function *coefficient = 0); + const Function *coefficient = 0, + const unsigned int selected_component = 0); /** * Exception @@ -185,7 +197,15 @@ class KellyErrorEstimator { * Exception */ DeclException0 (ExcInvalidBoundaryIndicator); - + /** + * Exception + */ + DeclException2 (ExcInvalidComponent, + int, int, + << "The component you gave (" << arg1 << ") " + << "is invalid with respect to the number " + << "of components in the finite element " + << "(" << arg2 << ")"); private: /** * Declare a data type to represent the @@ -227,6 +247,8 @@ class KellyErrorEstimator { FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, const Vector&solution, + const unsigned int n_components, + const unsigned int selected_component, const Function *coefficient); /** @@ -244,7 +266,9 @@ class KellyErrorEstimator { FESubfaceValues &fe_subface_values, FaceIntegrals &face_integrals, const Vector &solution, - const Function *coefficient); + const unsigned int n_components, + const unsigned int selected_component, + const Function *coefficient); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index c72052848c..ff5989b1d6 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -33,7 +33,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &, const FunctionMap &, const Vector &, Vector &, - const Function<1> *) { + const Function<1> *, + const unsigned int) { Assert(false, ExcNotImplemented()); }; @@ -46,10 +47,13 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const Function *coefficient) + const Function *coefficient, + const unsigned int selected_component) { Assert (neumann_bc.find(255) == neumann_bc.end(), ExcInvalidBoundaryIndicator()); + Assert (selected_component < dof.get_fe().n_components, + ExcInvalidComponent (selected_component, dof.get_fe().n_components)); // create a map of integrals indexed by // the corresponding face. In this map @@ -143,6 +147,8 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, fe_face_values_neighbor, face_integrals, solution, + dof.get_fe().n_components, + selected_component, coefficient); else // otherwise we need to do some @@ -154,6 +160,8 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, fe_face_values_cell, fe_subface_values, face_integrals, solution, + dof.get_fe().n_components, + selected_component, coefficient); }; @@ -195,6 +203,8 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iter FEFaceValues<1> &, FaceIntegrals &, const Vector &, + const unsigned int , + const unsigned int , const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -210,6 +220,8 @@ integrate_over_irregular_face (const active_cell_iterator &, FESubfaceValues<1> &, FaceIntegrals &, const Vector &, + const unsigned int , + const unsigned int , const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -228,6 +240,8 @@ integrate_over_regular_face (const active_cell_iterator &cell, FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, const Vector &solution, + const unsigned int n_components, + const unsigned int selected_component, const Function *coefficient) { const DoFHandler::face_iterator face = cell->face(face_no); @@ -243,7 +257,17 @@ integrate_over_regular_face (const active_cell_iterator &cell, // let psi be a short name for // [a grad u_h] vector > psi(n_q_points); - fe_face_values_cell.get_function_grads (solution, psi); + if (n_components == 1) + fe_face_values_cell.get_function_grads (solution, psi); + else + { + vector > > tmp (n_components, + psi); + fe_face_values_cell.get_function_grads (solution, tmp); + + psi.swap (tmp[selected_component]); + }; + // now compute over the other side of @@ -277,7 +301,16 @@ integrate_over_regular_face (const active_cell_iterator &cell, // the finite element solution // restricted to the neighbor cell vector > neighbor_psi (n_q_points); - fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); + if (n_components == 1) + fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); + else + { + vector > > tmp (n_components, + neighbor_psi); + fe_face_values_neighbor.get_function_grads (solution, tmp); + neighbor_psi.swap (tmp[selected_component]); + }; + // compute the jump in the gradients transform (psi.begin(), psi.end(), @@ -380,6 +413,8 @@ integrate_over_irregular_face (const active_cell_iterator &cell, FESubfaceValues &fe_subface_values, FaceIntegrals &face_integrals, const Vector &solution, + const unsigned int n_components, + const unsigned int selected_component, const Function *coefficient) { const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); Assert (neighbor.state() == valid, ExcInternalError()); @@ -423,14 +458,30 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // 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); + if (n_components == 1) + fe_subface_values.get_function_grads (solution, psi); + else + { + vector > > tmp (n_components, + psi); + fe_subface_values.get_function_grads (solution, tmp); + psi.swap (tmp[selected_component]); + }; // restrict the finite element on the // neighbor cell to the common #subface#. // store the gradient in #neighbor_psi# vector > neighbor_psi (n_q_points); fe_face_values.reinit (neighbor_child, neighbor_neighbor); - fe_face_values.get_function_grads (solution, neighbor_psi); + if (n_components == 1) + fe_face_values.get_function_grads (solution, neighbor_psi); + else + { + vector > > tmp (n_components, + neighbor_psi); + fe_face_values.get_function_grads (solution, tmp); + neighbor_psi.swap (tmp[selected_component]); + }; // compute the jump in the gradients transform (psi.begin(), psi.end(), -- 2.39.5