From: bangerth Date: Fri, 28 Jul 2006 21:00:42 +0000 (+0000) Subject: First steps to adjusting Kelly to hp DoFHandler X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=53082efebf57a9c7cdb062528524d19695ad9b69;p=dealii-svn.git First steps to adjusting Kelly to hp DoFHandler git-svn-id: https://svn.dealii.org/trunk@13479 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 5cf95ca20f..fda3d3dab0 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -18,14 +18,17 @@ #include #include #include -#include #include -#include #include -template class FEFaceValues; -template class FESubfaceValues; +namespace hp +{ + template class QCollection; + template class MappingCollection; + template class FEFaceValues; + template class FESubfaceValues; +} @@ -332,7 +335,7 @@ class KellyErrorEstimator /** * Same function as above, but * accepts more than one solution - * vectors and returns one error + * vector and returns one error * vector for each solution * vector. For the reason of * existence of this function, @@ -431,14 +434,18 @@ class KellyErrorEstimator * general documentation of this * class for more information. */ - typedef typename std::map::face_iterator,std::vector > FaceIntegrals; + typedef + typename std::map::face_iterator,std::vector > + FaceIntegrals; /** * Redeclare an active cell iterator. * This is simply for convenience. */ - typedef typename DoFHandler::active_cell_iterator active_cell_iterator; + typedef + typename DoFHandler::active_cell_iterator + active_cell_iterator; /** @@ -592,9 +599,9 @@ class KellyErrorEstimator * seperatly. */ template - static void estimate_some (const Mapping &mapping, + static void estimate_some (const hp::MappingCollection &mapping, const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const typename FunctionMap::type &neumann_bc, const std::vector &solutions, const std::vector &component_mask, @@ -627,7 +634,7 @@ class KellyErrorEstimator static void integrate_over_regular_face (const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const typename FunctionMap::type &neumann_bc, const std::vector &solutions, const std::vector &component_mask, @@ -636,8 +643,8 @@ class KellyErrorEstimator PerThreadData &per_thread_data, const active_cell_iterator &cell, const unsigned int face_no, - FEFaceValues &fe_face_values_cell, - FEFaceValues &fe_face_values_neighbor); + hp::FEFaceValues &fe_face_values_cell, + hp::FEFaceValues &fe_face_values_neighbor); /** @@ -654,7 +661,7 @@ class KellyErrorEstimator static void integrate_over_irregular_face (const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const std::vector &solutions, const std::vector &component_mask, const Function *coefficients, @@ -662,8 +669,8 @@ class KellyErrorEstimator PerThreadData &per_thread_data, const active_cell_iterator &cell, const unsigned int face_no, - FEFaceValues &fe_face_values, - FESubfaceValues &fe_subface_values); + hp::FEFaceValues &fe_face_values, + hp::FESubfaceValues &fe_subface_values); /** * By the resolution of Defect diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 2822785992..d1a62edb52 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -25,8 +25,11 @@ #include #include #include +#include #include #include +#include +#include #include #include @@ -455,9 +458,9 @@ estimate (const DoFHandler &dof_handler, template template void KellyErrorEstimator:: -estimate_some (const Mapping &mapping, +estimate_some (const hp::MappingCollection &mapping, const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const typename FunctionMap::type &neumann_bc, const std::vector &solutions, const std::vector &component_mask, @@ -470,6 +473,8 @@ estimate_some (const Mapping &mapping, const unsigned int subdomain_id = per_thread_data.subdomain_id; const unsigned int material_id = per_thread_data.material_id; + + const hp::FECollection fe (dof_handler.get_fe()); // make up a fe face values object // for the restriction of the @@ -488,24 +493,24 @@ estimate_some (const Mapping &mapping, // in debug mode, make sure that // some data matches, so compute // quadrature points always - FEFaceValues fe_face_values_cell (mapping, - dof_handler.get_fe(), - quadrature, - UpdateFlags( - update_gradients | - update_JxW_values | - ((!neumann_bc.empty() || - (coefficients != 0)) ? - update_q_points : 0) | - update_normal_vectors)); - FEFaceValues fe_face_values_neighbor (mapping, - dof_handler.get_fe(), + hp::FEFaceValues fe_face_values_cell (mapping, + fe, + quadrature, + UpdateFlags( + update_gradients | + update_JxW_values | + ((!neumann_bc.empty() || + (coefficients != 0)) ? + update_q_points : 0) | + update_normal_vectors)); + hp::FEFaceValues fe_face_values_neighbor (mapping, + fe, quadrature, update_gradients); - FESubfaceValues fe_subface_values (mapping, - dof_handler.get_fe(), - quadrature, - update_gradients); + hp::FESubfaceValues fe_subface_values (mapping, + fe, + quadrature, + update_gradients); active_cell_iterator cell = dof_handler.begin_active(); @@ -817,9 +822,9 @@ estimate (const Mapping &mapping, // work around another nasty bug in // icc7 Threads::ThreadGroup<> threads; - void (*estimate_some_ptr) (const Mapping &, + void (*estimate_some_ptr) (const hp::MappingCollection &, const DoFHandler &, - const Quadrature &, + const hp::QCollection &, const typename FunctionMap::type &, const std::vector &, const std::vector &, @@ -828,11 +833,17 @@ estimate (const Mapping &mapping, FaceIntegrals &, PerThreadData &) = &KellyErrorEstimator::template estimate_some; + + hp::MappingCollection mapping_collection; + mapping_collection.push_back (mapping); + + hp::QCollection face_quadratures; + face_quadratures.push_back (quadrature); for (unsigned int i=0; i template void KellyErrorEstimator:: integrate_over_regular_face (const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const typename FunctionMap::type &neumann_bc, const std::vector &solutions, const std::vector &component_mask, @@ -941,11 +952,11 @@ integrate_over_regular_face (const DoFHandler &dof_handler, PerThreadData &per_thread_data, const active_cell_iterator &cell, const unsigned int face_no, - FEFaceValues &fe_face_values_cell, - FEFaceValues &fe_face_values_neighbor) + hp::FEFaceValues &fe_face_values_cell, + hp::FEFaceValues &fe_face_values_neighbor) { const typename DoFHandler::face_iterator face = cell->face(face_no); - const unsigned int n_q_points = quadrature.n_quadrature_points, + const unsigned int n_q_points = quadrature[cell->active_fe_index()].n_quadrature_points, n_components = dof_handler.get_fe().n_components(), n_solution_vectors = solutions.size(); @@ -957,7 +968,8 @@ integrate_over_regular_face (const DoFHandler &dof_handler, // get gradients of the finite element // function on this cell for (unsigned int n=0; n &dof_handler, // get gradients on neighbor cell for (unsigned int n=0; n &dof_handler, // would only change the sign. We take // the outward normal. - per_thread_data.normal_vectors=fe_face_values_cell.get_normal_vectors(); + per_thread_data.normal_vectors = + fe_face_values_cell.get_present_fe_values().get_normal_vectors(); for (unsigned int n=0; n &dof_handler, if (coefficients->n_components == 1) { - coefficients->value_list (fe_face_values_cell.get_quadrature_points(), + coefficients->value_list (fe_face_values_cell.get_present_fe_values() + .get_quadrature_points(), per_thread_data.coefficient_values1); for (unsigned int n=0; n &dof_handler, else // vector-valued coefficient { - coefficients->vector_value_list (fe_face_values_cell.get_quadrature_points(), + coefficients->vector_value_list (fe_face_values_cell.get_present_fe_values() + .get_quadrature_points(), per_thread_data.coefficient_values); for (unsigned int n=0; n &dof_handler, { std::vector g(n_q_points); neumann_bc.find(boundary_indicator)->second - ->value_list (fe_face_values_cell.get_quadrature_points(), g); + ->value_list (fe_face_values_cell.get_present_fe_values() + .get_quadrature_points(), g); for (unsigned int n=0; n &dof_handler, { std::vector > g(n_q_points, Vector(n_components)); neumann_bc.find(boundary_indicator)->second - ->vector_value_list (fe_face_values_cell.get_quadrature_points(), + ->vector_value_list (fe_face_values_cell.get_present_fe_values() + .get_quadrature_points(), g); for (unsigned int n=0; n &dof_handler, // mentioned value at one of the // quadrature points - per_thread_data.JxW_values = fe_face_values_cell.get_JxW_values(); + per_thread_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 @@ -1125,7 +1145,7 @@ template template void KellyErrorEstimator:: integrate_over_irregular_face (const DoFHandler &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &quadrature, const std::vector &solutions, const std::vector &component_mask, const Function *coefficients, @@ -1133,11 +1153,11 @@ integrate_over_irregular_face (const DoFHandler &dof_handler, PerThreadData &per_thread_data, const active_cell_iterator &cell, const unsigned int face_no, - FEFaceValues &fe_face_values, - FESubfaceValues &fe_subface_values) + hp::FEFaceValues &fe_face_values, + hp::FESubfaceValues &fe_subface_values) { const typename DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); - const unsigned int n_q_points = quadrature.n_quadrature_points, + const unsigned int n_q_points = quadrature[cell->active_fe_index()].n_quadrature_points, n_components = dof_handler.get_fe().n_components(), n_solution_vectors = solutions.size(); const typename DoFHandler::face_iterator @@ -1188,13 +1208,15 @@ integrate_over_irregular_face (const DoFHandler &dof_handler, // store the gradient of the // solution in psi for (unsigned int n=0; n &dof_handler, // // let phi be the name of the integrand - per_thread_data.normal_vectors=fe_face_values.get_normal_vectors(); + per_thread_data.normal_vectors + = fe_face_values.get_present_fe_values().get_normal_vectors(); for (unsigned int n=0; n &dof_handler, // scalar coefficient if (coefficients->n_components == 1) { - coefficients->value_list (fe_face_values.get_quadrature_points(), + coefficients->value_list (fe_face_values.get_present_fe_values() + .get_quadrature_points(), per_thread_data.coefficient_values1); for (unsigned int n=0; n &dof_handler, else // vector-valued coefficient { - coefficients->vector_value_list (fe_face_values.get_quadrature_points(), + coefficients->vector_value_list (fe_face_values.get_present_fe_values() + .get_quadrature_points(), per_thread_data.coefficient_values); for (unsigned int n=0; n &dof_handler, // neighbor cell, while the // latter is on the refined // face of the big cell here - per_thread_data.JxW_values = fe_face_values.get_JxW_values(); + per_thread_data.JxW_values + = fe_face_values.get_present_fe_values().get_JxW_values(); // take the square of the phi[i] // for integration, and sum up