From 47629428e77465edde9eba8565f7e793476531a8 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 16 May 2015 00:55:55 +0200 Subject: [PATCH] Added template Number to DerivativeForms. --- doc/news/changes.h | 6 + include/deal.II/base/derivative_form.h | 158 +++++++++--------- tests/base/derivative_forms_01.cc | 77 +++++++++ ...erivative_forms_01.with_trilinos=on.output | 55 ++++++ 4 files changed, 216 insertions(+), 80 deletions(-) create mode 100644 tests/base/derivative_forms_01.cc create mode 100644 tests/base/derivative_forms_01.with_trilinos=on.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 32fa7b998c..1d9097ce29 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -377,6 +377,12 @@ inconvenience this causes.
    +
  1. New: DerivativeForm() now takes an additional optional + template argument specifying the type, similarly to Tensor() classes. +
    + (Luca Heltai, 2015/05/16) +
    2. +
  2. New: Utilities::MPI::min() functions.
    (Timo Heister, 2015/05/12) diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index e6a4b57af1..69d3820fab 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -29,11 +29,11 @@ DEAL_II_NAMESPACE_OPEN * R}^{\text{spacedim}}$, the second derivative a bilinear map from ${\mathbb * R}^{\text{dim}} \times {\mathbb R}^{\text{dim}}$ to ${\mathbb * R}^{\text{spacedim}}$ and so on. In deal.II we represent these derivaties - * using objects of type DerivativeForm<1,dim,spacedim>, - * DerivativeForm<2,dim,spacedim> and so on. - * @author Sebastian Pauletti, 2011 + * using objects of type DerivativeForm<1,dim,spacedim,Number>, + * DerivativeForm<2,dim,spacedim,Number> and so on. + * @author Sebastian Pauletti, 2011, Luca Heltai, 2015 */ -template +template class DerivativeForm { public: @@ -45,50 +45,50 @@ public: /** * Constructor from a second order tensor. */ - DerivativeForm (const Tensor<2,dim> &); + DerivativeForm (const Tensor<2,dim,Number> &); /** * Read-Write access operator. */ - Tensor &operator [] (const unsigned int i); + Tensor &operator [] (const unsigned int i); /** * Read-only access operator. */ - const Tensor &operator [] (const unsigned int i) const; + const Tensor &operator [] (const unsigned int i) const; /** * Assignment operator. */ - DerivativeForm &operator = (const DerivativeForm &); + DerivativeForm &operator = (const DerivativeForm &); /** * Assignment operator. */ - DerivativeForm &operator = (const Tensor<2,dim> &); + DerivativeForm &operator = (const Tensor<2,dim, Number> &); /** * Assignment operator. */ - DerivativeForm &operator = (const Tensor<1,dim> &); + DerivativeForm &operator = (const Tensor<1,dim, Number> &); /** - * Converts a DerivativeForm <1,dim, dim> to Tensor<2,dim>. If the + * Converts a DerivativeForm <1,dim, dim> to Tensor<2,dim,Number>. If the * derivative is the Jacobian of F, then Tensor[i] = grad(F^i). */ - operator Tensor<2,dim>() const; + operator Tensor<2,dim,Number>() const; /** - * Converts a DerivativeForm <1, dim, 1> to Tensor<1,dim>. + * Converts a DerivativeForm <1, dim, 1> to Tensor<1,dim,Number>. */ - operator Tensor<1,dim>() const; + operator Tensor<1,dim,Number>() const; /** * Return the transpose of a rectangular DerivativeForm, that is to say * viewed as a two dimensional matrix. */ - DerivativeForm<1, spacedim, dim> transpose () const; + DerivativeForm<1, spacedim, dim, Number> transpose () const; /** @@ -104,7 +104,7 @@ public: * covariant matrix, namely $DF*G^{-1}$, where $G = DF^{t}*DF$. If $DF$ is * square, covariant from gives $DF^{-t}$. */ - DerivativeForm<1, dim, spacedim> covariant_form() const; + DerivativeForm<1, dim, spacedim, Number> covariant_form() const; @@ -127,14 +127,14 @@ private: /** * Auxiliary function that computes (*this) * T^{t} */ - DerivativeForm<1, dim, spacedim> times_T_t (Tensor<2,dim> T) const; + DerivativeForm<1, dim, spacedim, Number> times_T_t (Tensor<2,dim,Number> T) const; private: /** * Array of tensors holding the subelements. */ - Tensor tensor[spacedim]; + Tensor tensor[spacedim]; }; @@ -147,9 +147,9 @@ private: -template +template inline -DerivativeForm::DerivativeForm () +DerivativeForm::DerivativeForm () { // default constructor. not specifying an initializer list calls // the default constructor of the subobjects, which initialize them @@ -158,9 +158,9 @@ DerivativeForm::DerivativeForm () -template +template inline -DerivativeForm::DerivativeForm(const Tensor<2,dim> &T) +DerivativeForm::DerivativeForm(const Tensor<2,dim,Number> &T) { Assert( (dim == spacedim) && (order==1), ExcMessage("Only allowed for square tensors.")); @@ -171,11 +171,11 @@ DerivativeForm::DerivativeForm(const Tensor<2,dim> &T) -template +template inline -DerivativeForm & -DerivativeForm:: -operator = (const DerivativeForm &ta) +DerivativeForm & +DerivativeForm:: +operator = (const DerivativeForm &ta) { for (unsigned int j=0; j &ta) -template +template inline -DerivativeForm &DerivativeForm:: -operator = (const Tensor<2,dim> &ta) +DerivativeForm &DerivativeForm:: +operator = (const Tensor<2,dim,Number> &ta) { Assert( (dim == spacedim) && (order==1), ExcMessage("Only allowed for square tensors.")); @@ -201,10 +201,10 @@ operator = (const Tensor<2,dim> &ta) -template +template inline -DerivativeForm &DerivativeForm:: -operator = (const Tensor<1,dim> &T) +DerivativeForm &DerivativeForm:: +operator = (const Tensor<1,dim,Number> &T) { Assert( (1 == spacedim) && (order==1), ExcMessage("Only allowed for spacedim==1 and order==1.")); @@ -217,9 +217,9 @@ operator = (const Tensor<1,dim> &T) -template +template inline -Tensor &DerivativeForm:: +Tensor &DerivativeForm:: operator[] (const unsigned int i) { Assert (i +template inline -const Tensor &DerivativeForm:: +const Tensor &DerivativeForm:: operator[] (const unsigned int i) const { Assert (i +template inline -DerivativeForm::operator Tensor<1,dim>() const +DerivativeForm::operator Tensor<1,dim,Number>() const { Assert( (1 == spacedim) && (order==1), ExcMessage("Only allowed for spacedim==1.")); @@ -254,16 +254,14 @@ DerivativeForm::operator Tensor<1,dim>() const -template +template inline -DerivativeForm::operator Tensor<2,dim>() const +DerivativeForm::operator Tensor<2,dim,Number>() const { Assert( (dim == spacedim) && (order==1), ExcMessage("Only allowed for square tensors.")); - - - Tensor<2,dim> t; + Tensor<2,dim,Number> t; if ((dim == spacedim) && (order==1)) for (unsigned int j=0; j::operator Tensor<2,dim>() const -template +template inline -DerivativeForm<1,spacedim,dim> -DerivativeForm:: +DerivativeForm<1,spacedim,dim,Number> +DerivativeForm:: transpose () const { Assert(order==1, ExcMessage("Only for rectangular DerivativeForm.")); - DerivativeForm<1,spacedim,dim> tt; + DerivativeForm<1,spacedim,dim,Number> tt; for (unsigned int i=0; i +template inline -DerivativeForm<1, dim, spacedim> -DerivativeForm::times_T_t (Tensor<2,dim> T) const +DerivativeForm<1, dim, spacedim,Number> +DerivativeForm::times_T_t (Tensor<2,dim,Number> T) const { Assert( order==1, ExcMessage("Only for order == 1.")); - DerivativeForm<1,dim, spacedim> dest; + DerivativeForm<1,dim, spacedim,Number> dest; for (unsigned int i=0; i::times_T_t (Tensor<2,dim> T) const } -template +template inline double -DerivativeForm::determinant () const +DerivativeForm::determinant () const { Assert( order==1, ExcMessage("Only for order == 1.")); if (dim == spacedim) { - Tensor<2,dim> T = (Tensor<2,dim>)( (*this) ); + Tensor<2,dim,Number> T = (Tensor<2,dim,Number>)( (*this) ); return dealii::determinant(T); } else { Assert( spacedim>dim, ExcMessage("Only for spacedim>dim.")); DerivativeForm<1,spacedim,dim> DF_t = this->transpose(); - Tensor<2, dim> G; //First fundamental form + Tensor<2,dim,Number> G; //First fundamental form for (unsigned int i=0; i::determinant () const -template +template inline -DerivativeForm<1,dim,spacedim> -DerivativeForm::covariant_form() const +DerivativeForm<1,dim,spacedim,Number> +DerivativeForm::covariant_form() const { if (dim == spacedim) { - Tensor<2,dim> DF_t (dealii::transpose(invert( (Tensor<2,dim>)(*this) ))); + Tensor<2,dim,Number> DF_t (dealii::transpose(invert( (Tensor<2,dim,Number>)(*this) ))); DerivativeForm<1,dim, spacedim> result = DF_t; return (result); } @@ -353,7 +351,7 @@ DerivativeForm::covariant_form() const { DerivativeForm<1,spacedim,dim> DF_t = this->transpose(); - Tensor<2, dim> G; //First fundamental form + Tensor<2,dim,Number> G; //First fundamental form for (unsigned int i=0; i::covariant_form() const } -template +template inline std::size_t -DerivativeForm::memory_consumption () +DerivativeForm::memory_consumption () { - return sizeof(DerivativeForm); + return sizeof(DerivativeForm); } #endif // DOXYGEN @@ -387,13 +385,13 @@ DerivativeForm::memory_consumption () * @relates DerivativeForm * @author Sebastian Pauletti, 2011 */ -template +template inline -Tensor<1, spacedim> -apply_transformation (const DerivativeForm<1,dim,spacedim> &DF, - const Tensor<1,dim> &T) +Tensor<1,spacedim,Number> +apply_transformation (const DerivativeForm<1,dim,spacedim,Number> &DF, + const Tensor<1,dim,Number> &T) { - Tensor<1, spacedim> dest; + Tensor<1,spacedim,Number> dest; for (unsigned int i=0; i &DF, * @author Sebastian Pauletti, 2011 */ //rank=2 -template +template inline DerivativeForm<1, spacedim, dim> -apply_transformation (const DerivativeForm<1,dim,spacedim> &DF, - const Tensor<2,dim> &T) +apply_transformation (const DerivativeForm<1,dim,spacedim,Number> &DF, + const Tensor<2,dim,Number> &T) { DerivativeForm<1, spacedim, dim> dest; @@ -428,13 +426,13 @@ apply_transformation (const DerivativeForm<1,dim,spacedim> &DF, * @relates DerivativeForm * @author Sebastian Pauletti, 2011 */ -template +template inline -Tensor<2, spacedim> -apply_transformation (const DerivativeForm<1,dim,spacedim> &DF1, - const DerivativeForm<1,dim,spacedim> &DF2) +Tensor<2,spacedim,Number> +apply_transformation (const DerivativeForm<1,dim,spacedim,Number> &DF1, + const DerivativeForm<1,dim,spacedim,Number> &DF2) { - Tensor<2, spacedim> dest; + Tensor<2,spacedim,Number> dest; for (unsigned int i=0; i &DF1, * @relates DerivativeForm * @author Sebastian Pauletti, 2011 */ -template +template inline -DerivativeForm<1,spacedim,dim> -transpose (const DerivativeForm<1,dim,spacedim> &DF) +DerivativeForm<1,spacedim,dim,Number> +transpose (const DerivativeForm<1,dim,spacedim,Number> &DF) { - DerivativeForm<1,spacedim,dim> tt; + DerivativeForm<1,spacedim,dim,Number> tt; tt = DF.transpose(); return tt; } diff --git a/tests/base/derivative_forms_01.cc b/tests/base/derivative_forms_01.cc new file mode 100644 index 0000000000..1f9994d95d --- /dev/null +++ b/tests/base/derivative_forms_01.cc @@ -0,0 +1,77 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#include "../tests.h" +#include +#include "Sacado.hpp" + +// Compute the derivative of F: R^dim -> R^spacedim +// We construct the function F_i, i=[0,spacedim) +// and we compute the first and second derivatives using +// Sacado + +typedef typename Sacado::Fad::DFad Sdouble; +typedef typename Sacado::Fad::DFad SSdouble; + +template +void test() { + Tensor<1,dim,SSdouble> x; + for(unsigned int j=0; j F; + for(unsigned int i=0; i dF; + for(unsigned int i=0; i ddF; + for(unsigned int i=0; i(); + test<1,2>(); + test<1,3>(); + test<2,2>(); + test<2,3>(); + test<3,3>(); +} diff --git a/tests/base/derivative_forms_01.with_trilinos=on.output b/tests/base/derivative_forms_01.with_trilinos=on.output new file mode 100644 index 0000000000..0fa161312f --- /dev/null +++ b/tests/base/derivative_forms_01.with_trilinos=on.output @@ -0,0 +1,55 @@ + +DEAL::dim = 1, spacedim = 1 +DEAL::x : 1.57080 [ 1.00000 ] [ 1.00000 [ ] ] +DEAL::F[0] : 1.00000 [ 6.12323e-17 ] [ 6.12323e-17 [ -1.00000 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 ] +DEAL::ddF[0] : -1.00000 +DEAL::dim = 1, spacedim = 2 +DEAL::x : 1.57080 [ 1.00000 ] [ 1.00000 [ ] ] +DEAL::F[0] : 1.00000 [ 6.12323e-17 ] [ 6.12323e-17 [ -1.00000 ] ] +DEAL::F[1] : 1.22465e-16 [ -1.00000 ] [ -1.00000 [ -1.22465e-16 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 ] +DEAL::dF[1] : -1.00000 [ -1.22465e-16 ] +DEAL::ddF[0] : -1.00000 +DEAL::ddF[1] : -1.22465e-16 +DEAL::dim = 1, spacedim = 3 +DEAL::x : 1.57080 [ 1.00000 ] [ 1.00000 [ ] ] +DEAL::F[0] : 1.00000 [ 6.12323e-17 ] [ 6.12323e-17 [ -1.00000 ] ] +DEAL::F[1] : 1.22465e-16 [ -1.00000 ] [ -1.00000 [ -1.22465e-16 ] ] +DEAL::F[2] : -1.00000 [ -1.83697e-16 ] [ -1.83697e-16 [ 1.00000 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 ] +DEAL::dF[1] : -1.00000 [ -1.22465e-16 ] +DEAL::dF[2] : -1.83697e-16 [ 1.00000 ] +DEAL::ddF[0] : -1.00000 +DEAL::ddF[1] : -1.22465e-16 +DEAL::ddF[2] : 1.00000 +DEAL::dim = 2, spacedim = 2 +DEAL::x : 1.57080 [ 1.00000 0.00000 ] [ 1.00000 [ ] 0.00000 [ ] ] 3.14159 [ 0.00000 1.00000 ] [ 0.00000 [ ] 1.00000 [ ] ] +DEAL::F[0] : 1.00000 [ 6.12323e-17 -1.00000 ] [ 6.12323e-17 [ -1.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 ] ] +DEAL::F[1] : -1.00000 [ -1.00000 -1.83697e-16 ] [ -1.00000 [ -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 ] +DEAL::dF[1] : -1.00000 [ -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 ] +DEAL::ddF[0] : -1.00000 0.00000 0.00000 -1.22465e-16 +DEAL::ddF[1] : -1.22465e-16 0.00000 0.00000 1.00000 +DEAL::dim = 2, spacedim = 3 +DEAL::x : 1.57080 [ 1.00000 0.00000 ] [ 1.00000 [ ] 0.00000 [ ] ] 3.14159 [ 0.00000 1.00000 ] [ 0.00000 [ ] 1.00000 [ ] ] +DEAL::F[0] : 1.00000 [ 6.12323e-17 -1.00000 ] [ 6.12323e-17 [ -1.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 ] ] +DEAL::F[1] : -1.00000 [ -1.00000 -1.83697e-16 ] [ -1.00000 [ -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 ] ] +DEAL::F[2] : -1.00000 [ -1.83697e-16 1.00000 ] [ -1.83697e-16 [ 1.00000 0.00000 ] 1.00000 [ 0.00000 2.44929e-16 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 ] +DEAL::dF[1] : -1.00000 [ -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 ] +DEAL::dF[2] : -1.83697e-16 [ 1.00000 0.00000 ] 1.00000 [ 0.00000 2.44929e-16 ] +DEAL::ddF[0] : -1.00000 0.00000 0.00000 -1.22465e-16 +DEAL::ddF[1] : -1.22465e-16 0.00000 0.00000 1.00000 +DEAL::ddF[2] : 1.00000 0.00000 0.00000 2.44929e-16 +DEAL::dim = 3, spacedim = 3 +DEAL::x : 1.57080 [ 1.00000 0.00000 0.00000 ] [ 1.00000 [ ] 0.00000 [ ] 0.00000 [ ] ] 3.14159 [ 0.00000 1.00000 0.00000 ] [ 0.00000 [ ] 1.00000 [ ] 0.00000 [ ] ] 4.71239 [ 0.00000 0.00000 1.00000 ] [ 0.00000 [ ] 0.00000 [ ] 1.00000 [ ] ] +DEAL::F[0] : 2.22045e-16 [ 6.12323e-17 -1.00000 -1.83697e-16 ] [ 6.12323e-17 [ -1.00000 0.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 0.00000 1.00000 ] ] +DEAL::F[1] : -1.00000 [ -1.00000 -1.83697e-16 1.00000 ] [ -1.00000 [ -1.22465e-16 0.00000 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 0.00000 ] 1.00000 [ 0.00000 0.00000 2.44929e-16 ] ] +DEAL::F[2] : -2.22045e-16 [ -1.83697e-16 1.00000 3.06162e-16 ] [ -1.83697e-16 [ 1.00000 0.00000 0.00000 ] 1.00000 [ 0.00000 2.44929e-16 0.00000 ] 3.06162e-16 [ 0.00000 0.00000 -1.00000 ] ] +DEAL::dF[0] : 6.12323e-17 [ -1.00000 0.00000 0.00000 ] -1.00000 [ 0.00000 -1.22465e-16 0.00000 ] -1.83697e-16 [ 0.00000 0.00000 1.00000 ] +DEAL::dF[1] : -1.00000 [ -1.22465e-16 0.00000 0.00000 ] -1.83697e-16 [ 0.00000 1.00000 0.00000 ] 1.00000 [ 0.00000 0.00000 2.44929e-16 ] +DEAL::dF[2] : -1.83697e-16 [ 1.00000 0.00000 0.00000 ] 1.00000 [ 0.00000 2.44929e-16 0.00000 ] 3.06162e-16 [ 0.00000 0.00000 -1.00000 ] +DEAL::ddF[0] : -1.00000 0.00000 0.00000 0.00000 -1.22465e-16 0.00000 0.00000 0.00000 1.00000 +DEAL::ddF[1] : -1.22465e-16 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 2.44929e-16 +DEAL::ddF[2] : 1.00000 0.00000 0.00000 0.00000 2.44929e-16 0.00000 0.00000 0.00000 -1.00000 -- 2.39.5