From aa4d82856186d9d8f611f864b86d2a1e3d3d2818 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 16 Dec 2010 21:25:41 +0000 Subject: [PATCH] Specialize functions for 3d. git-svn-id: https://svn.dealii.org/trunk@22984 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-38/step-38.cc | 77 ++++++++++++++++------------- 1 file changed, 42 insertions(+), 35 deletions(-) diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 3d64e00288..e5d78044a6 100755 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -59,28 +59,33 @@ class Solution : public Function }; -template -double Solution::value (const Point &p, - const unsigned int) const +template <> +double +Solution<3>::value (const Point<3> &p, + const unsigned int) const { - return sin(numbers::PI * p(0))*cos(numbers::PI * p(1))*exp(p(2)); + return (std::sin(numbers::PI * p(0)) * + std::cos(numbers::PI * p(1))*exp(p(2))); } -template -Tensor<1,dim> Solution::gradient (const Point &p, - const unsigned int) const + +template <> +Tensor<1,3> +Solution<3>::gradient (const Point<3> &p, + const unsigned int) const { - double dPi = numbers::PI; + using numbers::PI; - Tensor<1,dim> return_value; + Tensor<1,3> return_value; - return_value[0] = dPi *cos(dPi * p(0))*cos(dPi * p(1))*exp(p(2)); - return_value[1] = -dPi *sin(dPi * p(0))*sin(dPi * p(1))*exp(p(2)); - return_value[2] = sin(dPi * p(0))*cos(dPi * p(1))*exp(p(2)); + return_value[0] = PI *cos(PI * p(0))*cos(PI * p(1))*exp(p(2)); + return_value[1] = -PI *sin(PI * p(0))*sin(PI * p(1))*exp(p(2)); + return_value[2] = sin(PI * p(0))*cos(PI * p(1))*exp(p(2)); return return_value; } + template class RightHandSide : public Function { @@ -91,48 +96,50 @@ class RightHandSide : public Function const unsigned int component = 0) const; }; -template -double RightHandSide::value (const Point &p, - const unsigned int comp) const + +template <> +double +RightHandSide<3>::value (const Point<3> &p, + const unsigned int comp) const { - Assert(dim == 3, ExcInternalError()); - - double dPi = numbers::PI; + using numbers::PI; // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu) - Tensor<2,dim> hessian; + Tensor<2,3> hessian; - hessian[0][0] = -dPi*dPi*sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); - hessian[1][1] = -dPi*dPi*sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); - hessian[2][2] = sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); + hessian[0][0] = -PI*PI*sin(PI*p(0))*cos(PI*p(1))*exp(p(2)); + hessian[1][1] = -PI*PI*sin(PI*p(0))*cos(PI*p(1))*exp(p(2)); + hessian[2][2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2)); - hessian[0][1] = -dPi*dPi*cos(dPi*p(0))*sin(dPi*p(1))*exp(p(2)); - hessian[1][0] = -dPi*dPi*cos(dPi*p(0))*sin(dPi*p(1))*exp(p(2)); + hessian[0][1] = -PI*PI*cos(PI*p(0))*sin(PI*p(1))*exp(p(2)); + hessian[1][0] = -PI*PI*cos(PI*p(0))*sin(PI*p(1))*exp(p(2)); - hessian[0][2] = dPi*cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); - hessian[2][0] = dPi*cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); + hessian[0][2] = PI*cos(PI*p(0))*cos(PI*p(1))*exp(p(2)); + hessian[2][0] = PI*cos(PI*p(0))*cos(PI*p(1))*exp(p(2)); - hessian[1][2] = -dPi*sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2)); - hessian[2][1] = -dPi*sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2)); + hessian[1][2] = -PI*sin(PI*p(0))*sin(PI*p(1))*exp(p(2)); + hessian[2][1] = -PI*sin(PI*p(0))*sin(PI*p(1))*exp(p(2)); - Tensor<1,dim> gradient; - gradient[0] = dPi * cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); - gradient[1] = - dPi * sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2)); - gradient[2] = sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2)); + Tensor<1,3> gradient; + gradient[0] = PI * cos(PI*p(0))*cos(PI*p(1))*exp(p(2)); + gradient[1] = - PI * sin(PI*p(0))*sin(PI*p(1))*exp(p(2)); + gradient[2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2)); double curvature; - Point normal; + Point<3> normal; double dLength; - curvature = dim-1; + curvature = 3-1; // dim-1 dLength = sqrt(p(0)*p(0)+p(1)*p(1)+p(2)*p(2)); normal[0] = p(0)/dLength; normal[1] = p(1)/dLength; normal[2] = p(2)/dLength; - return -trace(hessian) + (hessian * normal) * normal + (gradient * normal)*curvature; + return (-trace(hessian) + + (hessian * normal) * normal + + (gradient * normal) * curvature); } -- 2.39.5