From c3f987f3d2ab1cbd0f4a91391777fb06a99d18ac Mon Sep 17 00:00:00 2001 From: heltai Date: Fri, 10 Apr 2009 15:55:14 +0000 Subject: [PATCH] Fixed integrate difference to work with codimone grids git-svn-id: https://svn.dealii.org/trunk@18578 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 36 ++++---- .../include/numerics/vectors.templates.h | 82 +++++++++---------- .../deal.II/source/numerics/vectors.inst.in | 50 +++++++++++ deal.II/doc/news/changes.h | 7 ++ 4 files changed, 116 insertions(+), 59 deletions(-) diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index e63d0b2e5f..59e7d255ef 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -1336,7 +1336,7 @@ class VectorTools template static void create_point_source_vector(const hp::MappingCollection &mapping, const hp::DoFHandler &dof, - const Point &p, + const Point &p, Vector &rhs_vector); /** @@ -1351,7 +1351,7 @@ class VectorTools */ template static void create_point_source_vector(const hp::DoFHandler &dof, - const Point &p, + const Point &p, Vector &rhs_vector); /** @@ -1368,7 +1368,7 @@ class VectorTools static void create_boundary_right_hand_side (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, const std::set &boundary_indicators = std::set()); @@ -1381,7 +1381,7 @@ class VectorTools template static void create_boundary_right_hand_side (const DoFHandler &dof, const Quadrature &q, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, const std::set &boundary_indicators = std::set()); @@ -1393,7 +1393,7 @@ class VectorTools static void create_boundary_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, const std::set &boundary_indicators = std::set()); @@ -1410,7 +1410,7 @@ class VectorTools template static void create_boundary_right_hand_side (const hp::DoFHandler &dof, const hp::QCollection &q, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, const std::set &boundary_indicators = std::set()); @@ -1491,11 +1491,11 @@ class VectorTools static void integrate_difference (const Mapping &mapping, const DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const Quadrature &q, const NormType &norm, - const Function *weight=0, + const Function *weight=0, const double exponent = 2.); /** @@ -1506,22 +1506,22 @@ class VectorTools template static void integrate_difference (const DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const Quadrature &q, const NormType &norm, - const Function *weight=0, + const Function *weight=0, const double exponent = 2.); template static void integrate_difference (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const hp::QCollection &q, const NormType &norm, - const Function *weight=0, + const Function *weight=0, const double exponent = 2.); /** @@ -1532,11 +1532,11 @@ class VectorTools template static void integrate_difference (const hp::DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const hp::QCollection &q, const NormType &norm, - const Function *weight=0, + const Function *weight=0, const double exponent = 2.); /** @@ -1558,9 +1558,9 @@ class VectorTools template static void point_difference (const DoFHandler& dof, const InVector& fe_function, - const Function& exact_solution, + const Function& exact_solution, Vector& difference, - const Point& point); + const Point& point); /** * Point error evaluation. Find @@ -1583,9 +1583,9 @@ class VectorTools static void point_difference (const Mapping &mapping, const DoFHandler& dof, const InVector& fe_function, - const Function& exact_solution, + const Function& exact_solution, Vector& difference, - const Point& point); + const Point& point); /** * Evaluate a possibly diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 9ac535b381..1a883eb954 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -847,7 +847,7 @@ void VectorTools::create_point_source_vector (const DoFHandler template void VectorTools::create_point_source_vector (const hp::MappingCollection &mapping, const hp::DoFHandler &dof_handler, - const Point &p, + const Point &p, Vector &rhs_vector) { Assert (rhs_vector.size() == dof_handler.n_dofs(), @@ -857,7 +857,7 @@ void VectorTools::create_point_source_vector (const hp::MappingCollection::active_cell_iterator, Point > + std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof_handler, p); @@ -880,7 +880,7 @@ void VectorTools::create_point_source_vector (const hp::MappingCollection void VectorTools::create_point_source_vector (const hp::DoFHandler &dof_handler, - const Point &p, + const Point &p, Vector &rhs_vector) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -897,7 +897,7 @@ void VectorTools::create_boundary_right_hand_side (const Mapping &mapping, const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, + const Function &rhs_function, Vector &rhs_vector, const std::set &boundary_indicators) { @@ -1019,7 +1019,7 @@ void VectorTools::create_boundary_right_hand_side (const Mapping &, const DoFHandler &, const Quadrature &, - const Function &, + const Function &, Vector &, const std::set &) { @@ -1032,7 +1032,7 @@ template void VectorTools::create_boundary_right_hand_side (const DoFHandler &dof_handler, const Quadrature &quadrature, - const Function &rhs_function, + const Function &rhs_function, Vector &rhs_vector, const std::set &boundary_indicators) { @@ -1053,7 +1053,7 @@ void VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &mapping, const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, + const Function &rhs_function, Vector &rhs_vector, const std::set &boundary_indicators) { @@ -1186,7 +1186,7 @@ void VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &, const hp::DoFHandler &, const hp::QCollection &, - const Function &, + const Function &, Vector &, const std::set &) { @@ -1199,7 +1199,7 @@ template void VectorTools::create_boundary_right_hand_side (const hp::DoFHandler &dof_handler, const hp::QCollection &quadrature, - const Function &rhs_function, + const Function &rhs_function, Vector &rhs_vector, const std::set &boundary_indicators) { @@ -1519,7 +1519,7 @@ interpolate_boundary_values (const Mapping // face face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); - const std::vector > &dof_locations + const std::vector > &dof_locations = fe_values.get_quadrature_points (); if (fe_is_system) @@ -2032,7 +2032,7 @@ VectorTools::interpolate_boundary_values // of dofs on this face face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); - const std::vector > &dof_locations + const std::vector > &dof_locations = fe_values.get_quadrature_points (); if (fe_is_system) @@ -2930,7 +2930,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH const FiniteElement &fe = dof_handler.get_fe(); std::vector face_dofs (fe.dofs_per_face); - std::vector > dof_locations (fe.dofs_per_face); + std::vector > dof_locations (fe.dofs_per_face); // have a map that stores normal vectors // for each vector-dof tuple we want to @@ -3442,11 +3442,11 @@ namespace internal do_integrate_difference (const dealii::hp::MappingCollection &mapping, const DH &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const dealii::hp::QCollection &q, const dealii::VectorTools::NormType &norm, - const Function *weight, + const Function *weight, const double exponent_1) { // we mark the "exponent" parameter @@ -3502,15 +3502,15 @@ namespace internal break; } - dealii::hp::FECollection fe_collection (dof.get_fe()); - dealii::hp::FEValues x_fe_values(mapping, fe_collection, q, update_flags); + dealii::hp::FECollection fe_collection (dof.get_fe()); + dealii::hp::FEValues x_fe_values(mapping, fe_collection, q, update_flags); const unsigned int max_n_q_points = q.max_n_quadrature_points (); std::vector< Vector > function_values (max_n_q_points, Vector(n_components)); - std::vector > > - function_grads (max_n_q_points, std::vector >(n_components)); + std::vector > > + function_grads (max_n_q_points, std::vector >(n_components)); std::vector weight_values (max_n_q_points); @@ -3519,16 +3519,16 @@ namespace internal std::vector > psi_values (max_n_q_points, Vector(n_components)); - std::vector > > - psi_grads (max_n_q_points, std::vector >(n_components)); + std::vector > > + psi_grads (max_n_q_points, std::vector >(n_components)); std::vector psi_scalar (max_n_q_points); // tmp vector when we use the - // Function functions for + // Function functions for // scalar functions std::vector tmp_values (max_n_q_points); - std::vector > tmp_gradients (max_n_q_points); + std::vector > tmp_gradients (max_n_q_points); // loop over all cells typename DH::active_cell_iterator cell = dof.begin_active(), @@ -3539,7 +3539,7 @@ namespace internal // initialize for this cell x_fe_values.reinit (cell); - const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); + const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); const unsigned int n_q_points = fe_values.n_quadrature_points; // resize all out scratch @@ -3549,7 +3549,7 @@ namespace internal function_values.resize (n_q_points, Vector(n_components)); function_grads.resize (n_q_points, - std::vector >(n_components)); + std::vector >(n_components)); weight_values.resize (n_q_points); weight_vectors.resize (n_q_points, @@ -3558,7 +3558,7 @@ namespace internal psi_values.resize (n_q_points, Vector(n_components)); psi_grads.resize (n_q_points, - std::vector >(n_components)); + std::vector >(n_components)); psi_scalar.resize (n_q_points); tmp_values.resize (n_q_points); @@ -3791,11 +3791,11 @@ void VectorTools::integrate_difference (const Mapping &mapping, const DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const Quadrature &q, const NormType &norm, - const Function *weight, + const Function *weight, const double exponent) { internal::VectorTools @@ -3810,16 +3810,16 @@ template void VectorTools::integrate_difference (const DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const Quadrature &q, const NormType &norm, - const Function *weight, + const Function *weight, const double exponent) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); internal::VectorTools - ::do_integrate_difference(hp::StaticMappingQ1::mapping_collection, + ::do_integrate_difference(hp::StaticMappingQ1::mapping_collection, dof, fe_function, exact_solution, difference, hp::QCollection(q), norm, weight, exponent); @@ -3832,11 +3832,11 @@ void VectorTools::integrate_difference (const dealii::hp::MappingCollection &mapping, const dealii::hp::DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const dealii::hp::QCollection &q, const NormType &norm, - const Function *weight, + const Function *weight, const double exponent) { internal::VectorTools @@ -3851,11 +3851,11 @@ template void VectorTools::integrate_difference (const dealii::hp::DoFHandler &dof, const InVector &fe_function, - const Function &exact_solution, + const Function &exact_solution, OutVector &difference, const dealii::hp::QCollection &q, const NormType &norm, - const Function *weight, + const Function *weight, const double exponent) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -3872,9 +3872,9 @@ template void VectorTools::point_difference (const DoFHandler &dof, const InVector &fe_function, - const Function &exact_function, + const Function &exact_function, Vector &difference, - const Point &point) + const Point &point) { point_difference(StaticMappingQ1::mapping, dof, @@ -3890,9 +3890,9 @@ void VectorTools::point_difference (const Mapping &mapping, const DoFHandler &dof, const InVector &fe_function, - const Function &exact_function, + const Function &exact_function, Vector &difference, - const Point &point) + const Point &point) { const FiniteElement& fe = dof.get_fe(); @@ -3902,7 +3902,7 @@ VectorTools::point_difference (const Mapping &mapping, // first find the cell in which this point // is, initialize a quadrature rule with // it, and then a FEValues object - const std::pair::active_cell_iterator, Point > + const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, @@ -3973,7 +3973,7 @@ VectorTools::point_value (const Mapping &mapping, // first find the cell in which this point // is, initialize a quadrature rule with // it, and then a FEValues object - const std::pair::active_cell_iterator, Point > + const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); @@ -4011,7 +4011,7 @@ VectorTools::point_value (const Mapping &mapping, // first find the cell in which this point // is, initialize a quadrature rule with // it, and then a FEValues object - const std::pair::active_cell_iterator, Point > + const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, diff --git a/deal.II/deal.II/source/numerics/vectors.inst.in b/deal.II/deal.II/source/numerics/vectors.inst.in index 13600d5797..3fbe9e4adf 100644 --- a/deal.II/deal.II/source/numerics/vectors.inst.in +++ b/deal.II/deal.II/source/numerics/vectors.inst.in @@ -67,6 +67,56 @@ for (VEC : SERIAL_VECTORS) const NormType&, const Function*, const double); + +#if deal_II_dimension != 3 + template + void VectorTools::integrate_difference, deal_II_dimension+1> + (const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); + + template + void VectorTools::integrate_difference, deal_II_dimension+1 > + (const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); + + template + void VectorTools::integrate_difference, deal_II_dimension+1> + (const Mapping&, + const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); + + template + void VectorTools::integrate_difference, deal_II_dimension+1 > + (const Mapping&, + const DoFHandler&, + const VEC&, + const Function&, + Vector&, + const Quadrature&, + const NormType&, + const Function*, + const double); +#endif + + template void VectorTools::integrate_difference (const DoFHandler&, diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 14d0bcd89a..851560f9a6 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -790,6 +790,13 @@ inconvenience this causes.

deal.II

    +
  1. +

    + Fixed: Integrate Difference now works also on codimension one grids. +
    + (Luca Heltai 2009/04/10) +

    +

  2. Changed: The ConstraintMatrix class can now also handle inhomogeneous -- 2.39.5