From 8183fa171ba87c2c157709297d7815f814bfdb98 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 5 May 2006 19:27:43 +0000 Subject: [PATCH] Implement VectorTools::integrate_difference for hp as well. On the way, fix a problem where we run into an assertion for W1infty_seminorm. git-svn-id: https://svn.dealii.org/trunk@13071 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 28 + .../include/numerics/vectors.templates.h | 699 ++++++++++-------- .../source/numerics/vectors.instance.h | 43 ++ tests/hp/Makefile | 3 +- tests/hp/integrate_difference.cc | 111 +++ tests/hp/integrate_difference/.cvsignore | 1 + tests/hp/integrate_difference/cmp/generic | 23 + 7 files changed, 612 insertions(+), 296 deletions(-) create mode 100644 tests/hp/integrate_difference.cc create mode 100644 tests/hp/integrate_difference/.cvsignore create mode 100644 tests/hp/integrate_difference/cmp/generic diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 66a4e4190b..f1eed46568 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -34,6 +34,8 @@ template class DoFHandler; namespace hp { template class DoFHandler; + template class MappingCollection; + template class QCollection; } class ConstraintMatrix; @@ -925,6 +927,32 @@ class VectorTools 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, + OutVector &difference, + const hp::QCollection &q, + const NormType &norm, + const Function *weight=0, + const double exponent = 2.); + + /** + * Calls the @p integrate_difference + * function, see above, with + * mapping=MappingQ1@(). + */ + template + static void integrate_difference (const hp::DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const hp::QCollection &q, + const NormType &norm, + const Function *weight=0, + const double exponent = 2.); + /** * Point error evaluation. Find * the first cell containing the diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index a9a6088eb1..32e9174973 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -29,8 +29,10 @@ #include #include #include -#include +#include #include +#include +#include #include #include @@ -1326,313 +1328,377 @@ VectorTools::project_boundary_values (const DoFHandler &dof, -template -void -VectorTools::integrate_difference (const Mapping &mapping, - const DoFHandler &dof, - const InVector &fe_function, - const Function &exact_solution, - OutVector &difference, - const Quadrature &q, - const NormType &norm, - const Function *weight, - const double exponent_1) +namespace internal { - // we mark the "exponent" parameter - // to this function "const" since - // it is strictly incoming, but we - // need to set it to something - // different later on, if - // necessary, so have a read-write - // version of it: - double exponent = exponent_1; + namespace VectorTools + { + template + static + void + do_integrate_difference (const ::hp::MappingCollection &mapping, + const DH &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const ::hp::QCollection &q, + const ::VectorTools::NormType &norm, + const Function *weight, + const double exponent_1) + { + // we mark the "exponent" parameter + // to this function "const" since + // it is strictly incoming, but we + // need to set it to something + // different later on, if + // necessary, so have a read-write + // version of it: + double exponent = exponent_1; - const unsigned int n_q_points = q.n_quadrature_points; - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); + const unsigned int n_components = dof.get_fe().n_components(); + const bool fe_is_system = (n_components != 1); - if (weight!=0) - { - Assert ((weight->n_components==1) || (weight->n_components==n_components), - ExcDimensionMismatch(weight->n_components, n_components)); - } + if (weight!=0) + { + Assert ((weight->n_components==1) || (weight->n_components==n_components), + ExcDimensionMismatch(weight->n_components, n_components)); + } - difference.reinit (dof.get_tria().n_active_cells()); - - switch (norm) - { - case L2_norm: - case H1_seminorm: - case H1_norm: - exponent = 2.; - break; - case L1_norm: - exponent = 1.; - break; - default: - break; - } + difference.reinit (dof.get_tria().n_active_cells()); - UpdateFlags update_flags = UpdateFlags (update_q_points | - update_JxW_values); - switch (norm) - { - case H1_seminorm: - case W1p_seminorm: - case W1infty_seminorm: - update_flags |= UpdateFlags (update_gradients); - break; - case H1_norm: - case W1p_norm: - case W1infty_norm: - update_flags |= UpdateFlags (update_gradients); - // no break! - default: - update_flags |= UpdateFlags (update_values); - break; - } - - FEValues fe_values(mapping, fe, q, update_flags); - - std::vector< Vector > function_values (n_q_points, - Vector(n_components)); - std::vector > > function_grads (n_q_points, - std::vector >(n_components)); - std::vector weight_values (n_q_points); - std::vector > weight_vectors (n_q_points, - Vector(n_components)); - - std::vector > psi_values (n_q_points, - Vector(n_components)); - std::vector > > psi_grads (n_q_points, - std::vector >(n_components)); - std::vector psi_scalar (n_q_points); - // tmp vector when we use the - // Function functions for - // scalar functions - std::vector tmp_values (fe_values.n_quadrature_points); - std::vector > tmp_gradients (fe_values.n_quadrature_points); + switch (norm) + { + case ::VectorTools::L2_norm: + case ::VectorTools::H1_seminorm: + case ::VectorTools::H1_norm: + exponent = 2.; + break; + case ::VectorTools::L1_norm: + exponent = 1.; + break; + default: + break; + } - // loop over all cells - typename DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (unsigned int index=0; cell != endc; ++cell, ++index) - { - double diff=0; - // initialize for this cell - fe_values.reinit (cell); + UpdateFlags update_flags = UpdateFlags (update_q_points | + update_JxW_values); + switch (norm) + { + case ::VectorTools::H1_seminorm: + case ::VectorTools::W1p_seminorm: + case ::VectorTools::W1infty_seminorm: + update_flags |= UpdateFlags (update_gradients); + break; + case ::VectorTools::H1_norm: + case ::VectorTools::W1p_norm: + case ::VectorTools::W1infty_norm: + update_flags |= UpdateFlags (update_gradients); + // no break! + default: + update_flags |= UpdateFlags (update_values); + break; + } + + ::hp::FECollection fe_collection (dof.get_fe()); + ::hp::FEValues x_fe_values(mapping, fe_collection, q, update_flags); + + const unsigned int max_n_q_points = q.max_n_quadrature_points (); - if (weight!=0) - { - if (weight->n_components>1) - weight->vector_value_list (fe_values.get_quadrature_points(), - weight_vectors); - else - { - weight->value_list (fe_values.get_quadrature_points(), - weight_values); - for (unsigned int k=0;k > + function_values (max_n_q_points, Vector(n_components)); + std::vector > > + function_grads (max_n_q_points, std::vector >(n_components)); + + std::vector + weight_values (max_n_q_points); + std::vector > + weight_vectors (max_n_q_points, Vector(n_components)); + + 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_scalar (max_n_q_points); + + // tmp vector when we use the + // Function functions for + // scalar functions + std::vector tmp_values (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(), + endc = dof.end(); + for (unsigned int index=0; cell != endc; ++cell, ++index) + { + double diff=0; + // initialize for this cell + x_fe_values.reinit (cell); + + const 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 + // arrays to the number of + // quadrature points we use + // for the present cell + function_values.resize (n_q_points, + Vector(n_components)); + function_grads.resize (n_q_points, + std::vector >(n_components)); + + weight_values.resize (n_q_points); + weight_vectors.resize (n_q_points, + Vector(n_components)); + + psi_values.resize (n_q_points, + Vector(n_components)); + psi_grads.resize (n_q_points, + std::vector >(n_components)); + psi_scalar.resize (n_q_points); + + tmp_values.resize (n_q_points); + tmp_gradients.resize (n_q_points); + + if (weight!=0) + { + if (weight->n_components>1) + weight->vector_value_list (fe_values.get_quadrature_points(), + weight_vectors); + else + { + weight->value_list (fe_values.get_quadrature_points(), + weight_values); + for (unsigned int k=0;k +void +VectorTools::integrate_difference (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const Quadrature &q, + const NormType &norm, + const Function *weight, + const double exponent) +{ + internal::VectorTools + ::do_integrate_difference (hp::MappingCollection(mapping), + dof, fe_function, exact_solution, + difference, hp::QCollection(q), + norm, weight, exponent); } @@ -1648,9 +1714,52 @@ VectorTools::integrate_difference (const DoFHandler &dof, const double exponent) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - static const MappingQ1 mapping; - integrate_difference(mapping, dof, fe_function, exact_solution, - difference, q, norm, weight, exponent); + internal::VectorTools + ::do_integrate_difference(hp::StaticMappingQ1::mapping_collection, + dof, fe_function, exact_solution, + difference, hp::QCollection(q), + norm, weight, exponent); +} + + + +template +void +VectorTools::integrate_difference (const ::hp::MappingCollection &mapping, + const ::hp::DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const ::hp::QCollection &q, + const NormType &norm, + const Function *weight, + const double exponent) +{ + internal::VectorTools + ::do_integrate_difference (hp::MappingCollection(mapping), + dof, fe_function, exact_solution, + difference, q, + norm, weight, exponent); +} + + +template +void +VectorTools::integrate_difference (const ::hp::DoFHandler &dof, + const InVector &fe_function, + const Function &exact_solution, + OutVector &difference, + const ::hp::QCollection &q, + const NormType &norm, + const Function *weight, + const double exponent) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + internal::VectorTools + ::do_integrate_difference(hp::StaticMappingQ1::mapping_collection, + dof, fe_function, exact_solution, + difference, q, + norm, weight, exponent); } diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h index 7178be51a1..643dd8057b 100644 --- a/deal.II/deal.II/source/numerics/vectors.instance.h +++ b/deal.II/deal.II/source/numerics/vectors.instance.h @@ -80,6 +80,49 @@ void VectorTools::integrate_difference const Function*, const double); +template +void VectorTools::integrate_difference +(const hp::DoFHandler&, + const VEC&, + const Function&, + Vector&, + const hp::QCollection&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const hp::DoFHandler&, + const VEC&, + const Function&, + Vector&, + const hp::QCollection&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const hp::MappingCollection&, + const hp::DoFHandler&, + const VEC&, + const Function&, + Vector&, + const hp::QCollection&, + const NormType&, + const Function*, + const double); +template +void VectorTools::integrate_difference +(const hp::MappingCollection&, + const hp::DoFHandler&, + const VEC&, + const Function&, + Vector&, + const hp::QCollection&, + const NormType&, + const Function*, + const double); + template void VectorTools::point_difference ( const DoFHandler&, diff --git a/tests/hp/Makefile b/tests/hp/Makefile index 59aee90e20..5aad649e49 100644 --- a/tests/hp/Makefile +++ b/tests/hp/Makefile @@ -30,7 +30,8 @@ tests_x = hp_dof_handler \ n_dofs \ random \ step-* \ - n_active_fe_indices + n_active_fe_indices \ + integrate_difference # from above list of regular expressions, generate the real set of # tests diff --git a/tests/hp/integrate_difference.cc b/tests/hp/integrate_difference.cc new file mode 100644 index 0000000000..ab2bb3c37d --- /dev/null +++ b/tests/hp/integrate_difference.cc @@ -0,0 +1,111 @@ +//---------------------------- integrate_difference.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- integrate_difference.cc --------------------------- + + +// call VectorTools::integrate_difference with fe's distributed in the +// same random way as in hp/random + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test () +{ + deallog << "dim=" << dim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global (2); + tria.begin_active()->set_refine_flag (); + tria.execute_coarsening_and_refinement (); + tria.refine_global (4-dim); + + hp::FECollection fe_collection; + hp::QCollection q_collection; + for (unsigned int i=1; i<=4; ++i) + { + fe_collection.push_back(FE_Q (i)); + q_collection.push_back (QGauss (i+2)); + } + + + hp::DoFHandler dof_handler(tria); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + cell->set_active_fe_index (rand() % fe_collection.size()); + + dof_handler.distribute_dofs(fe_collection); + + Vector vec (dof_handler.n_dofs()); + for (unsigned int i=0; i diff (tria.n_active_cells()); + + VectorTools::NormType norms[] = + { + VectorTools::mean, + VectorTools::L1_norm, + VectorTools::L2_norm, + VectorTools::Linfty_norm, + VectorTools::H1_seminorm, + VectorTools::W1p_seminorm + }; + for (unsigned int i=0; i(), + diff, + q_collection, + norms[i]); + deallog << "i=" << i << ", diff=" << diff.l2_norm() << std::endl; + } +} + + +int main () +{ + std::ofstream logfile("integrate_difference/output"); + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1> (); + test<2> (); + test<3> (); + + deallog << "OK" << std::endl; +} diff --git a/tests/hp/integrate_difference/.cvsignore b/tests/hp/integrate_difference/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/integrate_difference/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/hp/integrate_difference/cmp/generic b/tests/hp/integrate_difference/cmp/generic new file mode 100644 index 0000000000..b71900a524 --- /dev/null +++ b/tests/hp/integrate_difference/cmp/generic @@ -0,0 +1,23 @@ + +DEAL::dim=1 +DEAL::i=0, diff=11. +DEAL::i=1, diff=11. +DEAL::i=2, diff=73. +DEAL::i=3, diff=5.3e+02 +DEAL::i=4, diff=3.2e+02 +DEAL::i=5, diff=3.2e+02 +DEAL::dim=2 +DEAL::i=0, diff=97. +DEAL::i=1, diff=97. +DEAL::i=2, diff=1.7e+03 +DEAL::i=3, diff=3.5e+04 +DEAL::i=4, diff=6.3e+03 +DEAL::i=5, diff=6.3e+03 +DEAL::dim=3 +DEAL::i=0, diff=5.5e+02 +DEAL::i=1, diff=5.5e+02 +DEAL::i=2, diff=1.3e+04 +DEAL::i=3, diff=3.5e+05 +DEAL::i=4, diff=3.7e+04 +DEAL::i=5, diff=3.7e+04 +DEAL::OK -- 2.39.5