Repository navigation
Expand file tree
/
Copy pathanalytic_soln.hpp
More file actions
116 lines (100 loc) · 3.46 KB
/
Copy pathanalytic_soln.hpp
File metadata and controls
116 lines (100 loc) · 3.46 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#ifndef _analytic_soln_hpp_
#define _analytic_soln_hpp_
//@HEADER
// ************************************************************************
//
// MiniFE: Simple Finite Element Assembly and Solve
// Copyright (2006-2013) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can 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.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
//
// ************************************************************************
//@HEADER
#include <cmath>
#ifndef MINIFE_SCALAR
#define MINIFE_SCALAR double;
#endif
namespace miniFE {
typedef MINIFE_SCALAR Scalar;
// The 'soln' function below computes the analytic solution for
// steady state temperature in a brick-shaped domain (formally called
// a rectangular parallelepiped). The inputs to the function are
// the x,y,z coordinates of the point at which temperature is to be
// computed, and the number of terms p,q in the series expansion.
//
// The equations used for the temperature solution are equations 9 and 10
// in section 6.2 of Carslaw & Jaeger, "Conduction of Heat in Solids".
//
// The paralellepiped being used is defined by this domain:
// 0 <= x <= 1.0
// 0 <= y <= 1.0
// 0 <= z <= 1.0
//
// With boundary conditions prescribing the temperature to be 1.0 on
// the x==1.0 face, and 0.0 on all other faces.
//
// Thus, in the equations from Carslaw & Jaeger, the following constants
// are used:
//
// a == b == c == 1.0 (the extents of the domain)
// v1 == 0.0 (temperature at x == 0.0)
// v2 == 1.0 (temperature at x == 1.0)
//
const Scalar PI = 3.141592653589793238462;
const Scalar PI_SQR = PI*PI;
const Scalar term0 = 16.0/(PI_SQR);
inline Scalar fcn_l(int p, int q)
{
return std::sqrt((2*p+1)*(2*p+1)*PI_SQR + (2*q+1)*(2*q+1)*PI_SQR);
}
inline Scalar fcn(int n, Scalar u)
{
return (2*n+1)*PI*u;
}
inline Scalar soln(Scalar x, Scalar y, Scalar z, int max_p, int max_q)
{
Scalar sum = 0;
for(int p=0; p<=max_p; ++p) {
const Scalar p21y = fcn(p, y);
const Scalar sin_py = std::sin(p21y)/(2*p+1);
for(int q=0; q<=max_q; ++q) {
const Scalar q21z = fcn(q, z);
const Scalar sin_qz = std::sin(q21z)/(2*q+1);
const Scalar l = fcn_l(p, q);
const Scalar sinh1 = std::sinh(l*x);
const Scalar sinh2 = std::sinh(l);
const Scalar tmp = (sinh1*sin_py)*(sin_qz/sinh2);
//if the scalar l gets too big, sinh(l) becomes inf.
//if that happens, tmp is a NaN.
//crude check for NaN:
//if tmp != tmp, tmp is NaN
if (tmp == tmp) {
sum += tmp;
}
else {
//if we got a NaN, break out of this inner loop and go to
//the next iteration of the outer loop.
break;
}
}
}
return term0*sum;
}
}//namespace miniFE
#endif /* _analytic_soln_hpp_ */