From: Wolfgang Bangerth Date: Fri, 22 Apr 2011 04:07:07 +0000 (+0000) Subject: Fix a bug in FEValuesViews::SymmetricTensor::divergence. At least I believe so. X-Git-Tag: v8.0.0~4137 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=86c7de8ce3f5893ac40b02fac8c2dcd700ae114b;p=dealii.git Fix a bug in FEValuesViews::SymmetricTensor::divergence. At least I believe so. git-svn-id: https://svn.dealii.org/trunk@23620 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ee575ccf9c..f8844005cb 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,11 @@ should be fixed now.

Specific improvements

    +
  1. Fixed: The function FEValuesViews::SymmetricTensor::divergence had a bug. +This is now fixed. +
    +(Wolfgang Bangerth, Feifei Cheng, Venkat Vallala 2011/04/21) +
  2. Fixed: Under some conditions, FEFaceValues applied to an FESystem element that contained elements of type FE_Nothing would receive an erroneous exception. This is now fixed. diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h index 1fcf713bef..448699d111 100644 --- a/deal.II/include/deal.II/fe/fe_values.h +++ b/deal.II/include/deal.II/fe/fe_values.h @@ -4285,9 +4285,9 @@ namespace FEValuesViews } else if (snc != -1) { - // have a single non-zero component + // we have a single non-zero component // when the symmetric tensor is - // repsresented in unrolled form. + // represented in unrolled form. // this implies we potentially have // two non-zero components when // represented in component form! we @@ -4299,43 +4299,48 @@ namespace FEValuesViews // is a first order tensor. // // assume the second-order tensor is - // A with componets A_{ij}. then + // A with components A_{ij}. then // A_{ij} = A_{ji} and there is only // one (if diagonal) or two non-zero // entries in the tensorial // representation. define the // divergence as: - // b_i := \dfrac{\partial A_{ij}}{\partial x_j}. + // b_i := \dfrac{\partial phi_{ij}}{\partial x_j}. + // (which is incidentally also + // b_j := \dfrac{\partial phi_{ij}}{\partial x_i}). + // In both cases, a sum is implied. // - // Now, knowing the row ii and - // collumn jj of the non-zero entry - // we compute the divergence as - // b_ii = \dfrac{\partial A_{ij}}{\partial x_jj} (no sum) - // and if ii =! jj (not on a diagonal) - // b_jj = \dfrac{\partial A_{ij}}{\partial x_ii} (no sum) - - divergence_type return_value; - - // non-zero index in unrolled format + // Now, we know the nonzero component + // in unrolled form: it is indicated + // by 'snc'. we can figure out which + // tensor components belong to this: const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; - const unsigned int ii = value_type::unrolled_to_component_indices(comp)[0]; const unsigned int jj = value_type::unrolled_to_component_indices(comp)[1]; - // value of the non-zero tensor - // component - const double A_ij = fe_values.shape_values(snc,q_point); - - // the gradient of the non-zero shape - // function + // given the form of the divergence + // above, if ii=jj there is only a + // single nonzero component of the + // full tensor and the gradient + // equals + // b_ii := \dfrac{\partial phi_{ii,ii}}{\partial x_ii}. + // all other entries of 'b' are zero + // + // on the other hand, if ii!=jj, then + // there are two nonzero entries in + // the full tensor and + // b_ii := \dfrac{\partial phi_{ii,jj}}{\partial x_ii}. + // b_jj := \dfrac{\partial phi_{ii,jj}}{\partial x_jj}. + // again, all other entries of 'b' are + // zero const Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point]; - return_value[ii] = A_ij * phi_grad[jj]; + divergence_type return_value; + return_value[ii] = phi_grad[jj]; - // if we are not on a diagonal if (ii != jj) - return_value[jj] = A_ij * phi_grad[ii]; + return_value[jj] = phi_grad[ii]; return return_value; diff --git a/tests/deal.II/fe_values_view_23.cc b/tests/deal.II/fe_values_view_23.cc new file mode 100644 index 0000000000..37e4fd1c09 --- /dev/null +++ b/tests/deal.II/fe_values_view_23.cc @@ -0,0 +1,110 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2007, 2008, 2010, 2011 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. +// +//---------------------------------------------------------------------- + + +// there was a bug in getting the divergence of shape functions for the +// SymmetricTensor extractors. test that it is fixed by comparing with +// get_function_divergences + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + deallog << "FE=" << fe.get_name() + << std::endl; + + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + Vector fe_function(dof.n_dofs()); + for (unsigned int i=0; i quadrature(2); + FEValues fe_values (fe, quadrature, + update_values | update_gradients); + fe_values.reinit (dof.begin_active()); + + // let the FEValues object compute the + // divergences at quadrature points + std::vector > divergences (quadrature.size()); + FEValuesExtractors::SymmetricTensor<2> extractor(0); + fe_values[extractor] + .get_function_divergences (fe_function, divergences); + + // now do the same "by hand" + std::vector local_dof_indices (fe.dofs_per_cell); + dof.begin_active()->get_dof_indices (local_dof_indices); + + for (unsigned int q=0; q div_alt; + for (unsigned int i=0; i +void test_hyper_sphere() +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const HyperBallBoundary boundary; + tr.set_boundary (0, boundary); + + FESystem fe (FE_Q(1), + SymmetricTensor<2,dim>::n_independent_components); + test(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("fe_values_view_23/output"); + deallog << std::setprecision (3); + + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-7); + + test_hyper_sphere<2>(); + test_hyper_sphere<3>(); +} diff --git a/tests/deal.II/fe_values_view_23/cmp/generic b/tests/deal.II/fe_values_view_23/cmp/generic new file mode 100644 index 0000000000..62bab9f620 --- /dev/null +++ b/tests/deal.II/fe_values_view_23/cmp/generic @@ -0,0 +1,51 @@ + +DEAL::FE=FESystem<2>[FE_Q<2>(1)^3] +DEAL::q_point=0 +DEAL:: method 1: 15.5 15.5 +DEAL:: method 2: 15.5 15.5 +DEAL:: +DEAL::q_point=1 +DEAL:: method 1: 18.3 18.3 +DEAL:: method 2: 18.3 18.3 +DEAL:: +DEAL::q_point=2 +DEAL:: method 1: 16.2 16.2 +DEAL:: method 2: 16.2 16.2 +DEAL:: +DEAL::q_point=3 +DEAL:: method 1: 20.7 20.7 +DEAL:: method 2: 20.7 20.7 +DEAL:: +DEAL::FE=FESystem<3>[FE_Q<3>(1)^6] +DEAL::q_point=0 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=1 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=2 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=3 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=4 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=5 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=6 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: +DEAL::q_point=7 +DEAL:: method 1: 99.4 99.4 99.4 +DEAL:: method 2: 99.4 99.4 99.4 +DEAL:: diff --git a/tests/deal.II/fe_values_view_23/fe_values_view_09/cmp/generic b/tests/deal.II/fe_values_view_23/fe_values_view_09/cmp/generic new file mode 100644 index 0000000000..977609c8f4 --- /dev/null +++ b/tests/deal.II/fe_values_view_23/fe_values_view_09/cmp/generic @@ -0,0 +1,91 @@ + +DEAL::FE=FESystem<2>[FE_Q<2>(1)-FE_RaviartThomas<2>(1)-FE_Nedelec<2>(1)] +DEAL::component=0 +DEAL::0.807 4.36 +DEAL::0.807 5.29 +DEAL::1.31 4.07 +DEAL::1.31 5.59 +DEAL::component=1 +DEAL::18.2 352. +DEAL::-18.2 337. +DEAL::333. 180. +DEAL::-295. 296. +DEAL::component=2 +DEAL::43.2 -31.7 +DEAL::43.2 118. +DEAL::125. -61.1 +DEAL::125. -40.8 +DEAL::component=3 +DEAL::2.58e-16 5.84 +DEAL::2.58e-16 5.84 +DEAL::3.58e-16 9.52 +DEAL::3.58e-16 9.52 +DEAL::component=4 +DEAL::5.84 -6.75 +DEAL::5.84 6.75 +DEAL::9.52 -11.0 +DEAL::9.52 11.0 +DEAL::FE=FESystem<3>[FE_Q<3>(1)-FE_RaviartThomas<3>(1)-FE_Nedelec<3>(1)] +DEAL::component=0 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::2.37 4.73 9.46 +DEAL::component=1 +DEAL::8.18e-13 -275. -138. +DEAL::-5.71e-13 -275. -138. +DEAL::-367. -275. 4.71e+03 +DEAL::367. -275. 4.92e+03 +DEAL::-184. 4.57e+03 -138. +DEAL::184. 4.79e+03 -138. +DEAL::4.80e+03 4.57e+03 4.71e+03 +DEAL::-4.38e+03 4.79e+03 4.92e+03 +DEAL::component=2 +DEAL::-275. 1.30e-12 -138. +DEAL::-275. -459. 4.80e+03 +DEAL::-275. -1.40e-12 -138. +DEAL::-275. 459. 5.02e+03 +DEAL::4.67e+03 -91.8 -138. +DEAL::4.67e+03 3.52e+03 4.80e+03 +DEAL::4.88e+03 91.8 -138. +DEAL::4.88e+03 -3.09e+03 5.02e+03 +DEAL::component=3 +DEAL::-275. -138. 1.36e-12 +DEAL::-275. 4.90e+03 -367. +DEAL::4.76e+03 -138. -184. +DEAL::4.76e+03 4.90e+03 2.23e+03 +DEAL::-275. -138. -3.64e-12 +DEAL::-275. 5.11e+03 367. +DEAL::4.97e+03 -138. 184. +DEAL::4.97e+03 5.11e+03 -1.81e+03 +DEAL::component=4 +DEAL::-7.70e-16 5.60 22.4 +DEAL::-7.70e-16 5.60 22.4 +DEAL::-6.69e-16 5.60 22.4 +DEAL::-6.69e-16 5.60 22.4 +DEAL::1.58e-16 5.60 22.4 +DEAL::1.58e-16 5.60 22.4 +DEAL::9.52e-16 5.60 22.4 +DEAL::9.52e-16 5.60 22.4 +DEAL::component=5 +DEAL::5.60 5.78e-16 22.4 +DEAL::5.60 -1.39e-15 22.4 +DEAL::5.60 5.78e-16 22.4 +DEAL::5.60 -1.39e-15 22.4 +DEAL::5.60 -6.60e-16 22.4 +DEAL::5.60 3.53e-16 22.4 +DEAL::5.60 -6.60e-16 22.4 +DEAL::5.60 3.53e-16 22.4 +DEAL::component=6 +DEAL::5.60 11.2 -1.38e-15 +DEAL::5.60 11.2 1.96e-16 +DEAL::5.60 11.2 1.23e-17 +DEAL::5.60 11.2 1.50e-15 +DEAL::5.60 11.2 -1.38e-15 +DEAL::5.60 11.2 1.96e-16 +DEAL::5.60 11.2 1.23e-17 +DEAL::5.60 11.2 1.50e-15