From 45af105020aba6364c171d5dd175af7dc29c97cd Mon Sep 17 00:00:00 2001 From: guido Date: Thu, 25 May 2000 21:30:55 +0000 Subject: [PATCH] VectorTools::integrate_difference: 1. Fixed computation of non-squared norms 2. weight functions with more components ComponentSelectFunction for selecting components git-svn-id: https://svn.dealii.org/trunk@2947 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/multigrid/multigrid.h | 4 +- deal.II/deal.II/include/numerics/vectors.h | 35 +-- deal.II/deal.II/source/numerics/vectors.cc | 222 +++++++++++++----- 3 files changed, 181 insertions(+), 80 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index 1098c4fb96..73f18a8f87 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -16,11 +16,11 @@ #include #include -#include -#include #include #include #include +#include +#include #include diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index b5d64ff3f5..347e24d60c 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -461,20 +461,27 @@ class VectorTools * * The additional argument * #weight# allows to evaluate - * weighted norms. This is useful - * for weighting the error of - * different parts differently. A - * special use is to have - * #weight=0# in some parts of - * the domain, e.g. at the - * location of a shock and - * #weight=1# elsewhere. This - * allows convergence tests in - * smooth parts of in general - * discontinuous solutions. By - * default, no weighting function - * is given, i.e. weight=1 in the - * whole domain. + * weighted norms. The weight + * function may be + * one-dimensional, establishing + * a weight in the domain. It + * also may have as many + * components as the finite + * element function: Then, + * different components get + * different weights. This can + * be applied for instant with + * the characteristic function of + * a subset of the domain or a + * weight function selecting only + * some components of the + * solution. The weight function + * is expected to be positive, + * but negative values are not + * filtered. By default, no + * weighting function is given, + * i.e. weight=1 in the whole + * domain. * * It is assumed that the number * of components of the function diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 8c3fb2a204..f5728bf737 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -802,6 +802,12 @@ VectorTools::integrate_difference (const DoFHandler &dof, Assert( !((n_components == 1) && (norm == mean)), ExcNotUseful()); + 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()); UpdateFlags update_flags = UpdateFlags (update_q_points | @@ -819,14 +825,13 @@ VectorTools::integrate_difference (const DoFHandler &dof, vector > > function_grads (n_q_points, vector >(n_components)); vector weight_values (n_q_points); + vector > weight_vectors (n_q_points, n_components); vector > psi_values (n_q_points, Vector(n_components)); vector > > psi_grads (n_q_points, vector >(n_components)); vector psi_scalar (n_q_points); - vector psi_square (n_q_points); - // tmp vector when we use the // Function functions for // scalar functions @@ -842,6 +847,16 @@ VectorTools::integrate_difference (const DoFHandler &dof, // initialize for this cell fe_values.reinit (cell); + 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); + } + switch (norm) { case mean: @@ -858,15 +873,15 @@ VectorTools::integrate_difference (const DoFHandler &dof, // the function really has only // one component if (fe_is_system) - exact_solution.vector_value_list (fe_values.get_quadrature_points(), - psi_values); + exact_solution.vector_value_list (fe_values.get_quadrature_points(), + psi_values); else { exact_solution.value_list (fe_values.get_quadrature_points(), tmp_values); for (unsigned int i=0; i &dof, for (unsigned int q=0; qn_components > 1) + for (unsigned int q=0; qn_components == 1 + for (unsigned int q=0; qn_components > 1) + for (unsigned int q=0; qn_components == 1 + for (unsigned int q=0; qvalue_list(fe_values.get_quadrature_points(), - weight_values); - for (unsigned int q=0; qn_components > 1) + for (unsigned int q=0; qn_components == 1 + for (unsigned int q=0; qn_components > 1) + for (unsigned int q=0; qn_components == 1 + for (unsigned int q=0; q &dof, psi_grads[q][k] -= function_grads[q][k]; // take square of integrand - fill_n (psi_square.begin(), n_q_points, 0.0); - for (unsigned int k=0; kvalue_list(fe_values.get_quadrature_points(), - weight_values); + for (unsigned int k=0; kn_components > 1) + for (unsigned int q=0; qn_components == 1 + for (unsigned int q=0; q