From 7e7c934f0250aa2c1a066bfefe5db84ed4af4701 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 15 Dec 2006 03:13:24 +0000 Subject: [PATCH] Use the same quadrature formula for the error estimator for both sides of a face, even for the hp case. git-svn-id: https://svn.dealii.org/trunk@14250 0785d39b-7218-0410-832d-ea1e28bc413d --- .../source/numerics/error_estimator.cc | 136 ++++++++++++------ 1 file changed, 95 insertions(+), 41 deletions(-) diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 86925590b6..705c681a8d 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -461,25 +461,61 @@ template KellyErrorEstimator::PerThreadData:: PerThreadData (const unsigned int n_solution_vectors, const unsigned int n_components, - const unsigned int n_q_points, + const unsigned int max_n_q_points, const unsigned int subdomain_id, const unsigned int material_id) : subdomain_id (subdomain_id), material_id (material_id) { - // Init the size of a lot of vectors - // needed in the calculations once - // per thread. - normal_vectors.resize(n_q_points); - coefficient_values1.resize(n_q_points); - coefficient_values.resize(n_q_points); - JxW_values.resize(n_q_points); + // Init the size of a lot of vectors needed + // in the calculations once per thread. we + // will later resize them as necessary, but + // for the moment just reserve the maximal + // memory necessary to avoid later + // re-allocation + normal_vectors.resize(max_n_q_points); + coefficient_values1.resize(max_n_q_points); + coefficient_values.resize(max_n_q_points); + JxW_values.resize(max_n_q_points); phi.resize(n_solution_vectors); psi.resize(n_solution_vectors); neighbor_psi.resize(n_solution_vectors); + for (unsigned int i=0; i +void +KellyErrorEstimator::PerThreadData:: +resize (const unsigned int n_components, + const unsigned int n_q_points) +{ + const unsigned int n_solution_vectors = phi.size(); + + normal_vectors.resize(n_q_points); + coefficient_values1.resize(n_q_points); + coefficient_values.resize(n_q_points); + JxW_values.resize(n_q_points); + for (unsigned int i=0; i &mapping, PerThreadData &per_thread_data) { const unsigned int n_solution_vectors = solutions.size(); - + const unsigned int n_components = dof_handler.get_fe().n_components(); + const unsigned int subdomain_id = per_thread_data.subdomain_id; const unsigned int material_id = per_thread_data.material_id; @@ -812,7 +849,15 @@ estimate_some (const hp::MappingCollection &mapping, // so now we know that we care for // this face, let's do something - // about it: + // about it. first re-size the + // arrays we may use to the correct + // size: + per_thread_data.resize (n_components, + quadrature[cell->active_fe_index()] + .n_quadrature_points); + + + // then do the actual integration if (face->has_children() == false) // if the face is a regular one, // i.e. either on the other side @@ -850,34 +895,13 @@ estimate_some (const hp::MappingCollection &mapping, -template -template -void -KellyErrorEstimator:: -estimate (const Mapping &/*mapping*/, - const DH &/*dof_handler*/, - const hp::QCollection &/*quadrature*/, - const typename FunctionMap::type &/*neumann_bc*/, - const std::vector &/*solutions*/, - std::vector*> &/*errors*/, - const std::vector &/*component_mask_*/, - const Function */*coefficients*/, - const unsigned int /*n_threads_*/, - const unsigned int /*subdomain_id*/, - const unsigned int /*material_id*/) -{ - Assert (false, ExcNotImplemented()); -} - - - template template void KellyErrorEstimator:: estimate (const Mapping &mapping, const DH &dof_handler, - const Quadrature &quadrature, + const hp::QCollection &face_quadratures, const typename FunctionMap::type &neumann_bc, const std::vector &solutions, std::vector*> &errors, @@ -974,7 +998,7 @@ estimate (const Mapping &mapping, data_structures[i] = new PerThreadData (solutions.size(), dof_handler.get_fe().n_components(), - quadrature.n_quadrature_points, + face_quadratures.max_n_quadrature_points(), subdomain_id, material_id); // split all cells into threads if @@ -999,9 +1023,6 @@ estimate (const Mapping &mapping, hp::MappingCollection mapping_collection; mapping_collection.push_back (mapping); - hp::QCollection face_quadratures; - face_quadratures.push_back (quadrature); - for (unsigned int i=0; i &mapping, } + +template +template +void +KellyErrorEstimator:: +estimate (const Mapping &mapping, + const DH &dof_handler, + const Quadrature &quadrature, + const typename FunctionMap::type &neumann_bc, + const std::vector &solutions, + std::vector*> &errors, + const std::vector &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const unsigned int subdomain_id, + const unsigned int material_id) +{ + // forward to the function with the QCollection + estimate (mapping, dof_handler, + hp::QCollection(quadrature), + neumann_bc, solutions, + errors, component_mask, coefficients, + n_threads, subdomain_id, material_id); +} + + template template void KellyErrorEstimator::estimate (const DH &dof_handler, @@ -1146,7 +1193,8 @@ integrate_over_regular_face (const DH &dof_handl // initialize data of the restriction // of this cell to the present face - fe_face_values_cell.reinit (cell, face_no); + fe_face_values_cell.reinit (cell, face_no, + cell->active_fe_index()); // get gradients of the finite element // function on this cell @@ -1175,8 +1223,12 @@ integrate_over_regular_face (const DH &dof_handl // get restriction of finite element // function of @p{neighbor} to the - // common face. - fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor); + // 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; nactive_fe_index()); // restrict the finite element // on the neighbor cell to the // common @p{subface}. - fe_face_values.reinit (neighbor_child, neighbor_neighbor); + fe_face_values.reinit (neighbor_child, neighbor_neighbor, + cell->active_fe_index()); // store the gradient of the // solution in psi -- 2.39.5