<h3>Specific improvements</h3>
<ol>
+<li> Fixed: The function FEValuesViews::SymmetricTensor::divergence had a bug.
+This is now fixed.
+<br>
+(Wolfgang Bangerth, Feifei Cheng, Venkat Vallala 2011/04/21)
+
<li> 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.
}
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
// 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;
--- /dev/null
+//----------------------------------------------------------------------
+// $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 <base/logstream.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+ const FiniteElement<dim>& fe)
+{
+ deallog << "FE=" << fe.get_name()
+ << std::endl;
+
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ Vector<double> fe_function(dof.n_dofs());
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ fe_function(i) = i+1;
+
+ const QGauss<dim> quadrature(2);
+ FEValues<dim> 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<Tensor<1,dim> > 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<unsigned int> local_dof_indices (fe.dofs_per_cell);
+ dof.begin_active()->get_dof_indices (local_dof_indices);
+
+ for (unsigned int q=0; q<quadrature.size(); ++q)
+ {
+ Tensor<1,dim> div_alt;
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ div_alt += fe_values[extractor].divergence (i,q) *
+ fe_function(local_dof_indices[i]);
+
+ deallog << "q_point=" << q << std::endl
+ << " method 1: " << divergences[q] << std::endl
+ << " method 2: " << div_alt << std::endl
+ << std::endl;
+ Assert ((divergences[q] - div_alt).norm() <= divergences[q].norm(),
+ ExcInternalError());
+ }
+}
+
+
+
+template<int dim>
+void test_hyper_sphere()
+{
+ Triangulation<dim> tr;
+ GridGenerator::hyper_ball(tr);
+
+ static const HyperBallBoundary<dim> boundary;
+ tr.set_boundary (0, boundary);
+
+ FESystem<dim> fe (FE_Q<dim>(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>();
+}
--- /dev/null
+
+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::
--- /dev/null
+
+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
projects / dealii.git / commitdiff