Skip to content
Snippets Groups Projects
Commit bbc8b54f authored by Oliver Sander's avatar Oliver Sander Committed by sander@FU-BERLIN.DE
Browse files

make the number type a parameter, but default to 'double' if no type is given 

[[Imported from SVN: r8145]]
parent 0705eb2d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/realtuple.hh
+ 21
− 21
View file @ bbc8b54f
......@@ -11,29 +11,29 @@
Currently this class only exists for testing purposes.
*/
template <int N>
template <int N, class T=double>
class RealTuple
{
public:
typedef double ctype;
typedef T ctype;
/** \brief The type used for global coordinates */
typedef Dune::FieldVector<double,N> CoordinateType;
typedef Dune::FieldVector<T,N> CoordinateType;
/** \brief Dimension of the manifold formed by unit vectors */
static const int dim = N;
typedef Dune::FieldVector<double,N> EmbeddedTangentVector;
typedef Dune::FieldVector<T,N> EmbeddedTangentVector;
typedef Dune::FieldVector<double,N> TangentVector;
typedef Dune::FieldVector<T,N> TangentVector;
/** \brief Default constructor */
RealTuple()
{}
/** \brief Construction from a scalar */
RealTuple(double v)
RealTuple(T v)
{
data_ = v;
}
......@@ -44,11 +44,11 @@ public:
{}
/** \brief Constructor from FieldVector*/
RealTuple(const Dune::FieldVector<double,N>& other)
RealTuple(const Dune::FieldVector<T,N>& other)
: data_(other)
{}
RealTuple& operator=(const Dune::FieldVector<double,N>& other) {
RealTuple& operator=(const Dune::FieldVector<T,N>& other) {
data_ = other;
return *this;
}
......@@ -61,7 +61,7 @@ public:
/** \brief Geodesic distance between two points
Simply the Euclidean distance */
static double distance(const RealTuple& a, const RealTuple& b) {
static T distance(const RealTuple& a, const RealTuple& b) {
return (a.data_ - b.data_).two_norm();
}
......@@ -80,9 +80,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
Dune::FieldMatrix<double,N,N> result;
Dune::FieldMatrix<T,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = 2*(i==j);
......@@ -94,9 +94,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const RealTuple& a, const RealTuple& b) {
static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const RealTuple& a, const RealTuple& b) {
Dune::FieldMatrix<double,N,N> result;
Dune::FieldMatrix<T,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = -2*(i==j);
......@@ -108,16 +108,16 @@ public:
The result is the constant zero-tensor.
*/
static Tensor3<double,N,N,N> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
return Tensor3<double,N,N,N>(0);
static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
return Tensor3<T,N,N,N>(0);
}
/** \brief Compute the mixed third derivative \partial d^3 / \partial da db^2
The result is the constant zero-tensor.
*/
static Tensor3<double,N,N,N> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const RealTuple& a, const RealTuple& b) {
return Tensor3<double,N,N,N>(0);
static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const RealTuple& a, const RealTuple& b) {
return Tensor3<T,N,N,N>(0);
}
/** \brief Project tangent vector of R^n onto the tangent space */
......@@ -126,7 +126,7 @@ public:
}
/** \brief The global coordinates, if you really want them */
const Dune::FieldVector<double,N>& globalCoordinates() const {
const Dune::FieldVector<T,N>& globalCoordinates() const {
return data_;
}
......@@ -134,9 +134,9 @@ public:
In general this frame field, may of course not be continuous, but for RealTuples it is.
*/
Dune::FieldMatrix<double,N,N> orthonormalFrame() const {
Dune::FieldMatrix<T,N,N> orthonormalFrame() const {
Dune::FieldMatrix<double,N,N> result;
Dune::FieldMatrix<T,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
......@@ -152,7 +152,7 @@ public:
private:
Dune::FieldVector<double,N> data_;
Dune::FieldVector<T,N> data_;
};
......
dune/gfe/unitvector.hh
+ 55
− 54
View file @ bbc8b54f
......@@ -9,18 +9,19 @@
/** \brief A unit vector in R^N
\tparam N Dimension of the embedding space
\tparam T The type used for individual coordinates
*/
template <int N>
template <int N, class T=double>
class UnitVector
{
/** \brief Computes sin(x) / x without getting unstable for small x */
static double sinc(const double& x) {
static T sinc(const T& x) {
return (x < 1e-4) ? 1 - (x*x/6) : std::sin(x)/x;
}
/** \brief Compute the derivative of arccos^2 without getting unstable for x close to 1 */
static double derivativeOfArcCosSquared(const double& x) {
const double eps = 1e-4;
static T derivativeOfArcCosSquared(const T& x) {
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
} else if (x < -1+eps) { // The function is not differentiable
......@@ -30,8 +31,8 @@ class UnitVector
}
/** \brief Compute the second derivative of arccos^2 without getting unstable for x close to 1 */
static double secondDerivativeOfArcCosSquared(const double& x) {
const double eps = 1e-4;
static T secondDerivativeOfArcCosSquared(const T& x) {
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return 2.0/3 - 8*(x-1)/15;
} else if (x < -1+eps) { // The function is not differentiable
......@@ -41,14 +42,14 @@ class UnitVector
}
/** \brief Compute the third derivative of arccos^2 without getting unstable for x close to 1 */
static double thirdDerivativeOfArcCosSquared(const double& x) {
const double eps = 1e-4;
static T thirdDerivativeOfArcCosSquared(const T& x) {
const T eps = 1e-4;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
return -8.0/15 + 24*(x-1)/35;
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else {
double d = 1-x*x;
T d = 1-x*x;
return 6*x/(d*d) - 6*x*x*std::acos(x)/(d*d*std::sqrt(d)) - 2*std::acos(x)/(d*std::sqrt(d));
}
}
......@@ -56,10 +57,10 @@ class UnitVector
public:
/** \brief The type used for coordinates */
typedef double ctype;
typedef T ctype;
/** \brief The type used for global coordinates */
typedef Dune::FieldVector<double,N> CoordinateType;
typedef Dune::FieldVector<T,N> CoordinateType;
/** \brief Dimension of the manifold formed by unit vectors */
static const int dim = N-1;
......@@ -67,30 +68,30 @@ public:
/** \brief Dimension of the Euclidean space the manifold is embedded in */
static const int embeddedDim = N;
typedef Dune::FieldVector<double,N-1> TangentVector;
typedef Dune::FieldVector<T,N-1> TangentVector;
typedef Dune::FieldVector<double,N> EmbeddedTangentVector;
typedef Dune::FieldVector<T,N> EmbeddedTangentVector;
/** \brief Default constructor */
UnitVector()
{}
/** \brief Constructor from a vector. The vector gets normalized */
UnitVector(const Dune::FieldVector<double,N>& vector)
UnitVector(const Dune::FieldVector<T,N>& vector)
: data_(vector)
{
data_ /= data_.two_norm();
}
/** \brief Constructor from an array. The array gets normalized */
UnitVector(const Dune::array<double,N>& vector)
UnitVector(const Dune::array<T,N>& vector)
{
for (int i=0; i<N; i++)
data_[i] = vector[i];
data_ /= data_.two_norm();
}
UnitVector<N>& operator=(const Dune::FieldVector<double,N>& vector)
UnitVector<N>& operator=(const Dune::FieldVector<T,N>& vector)
{
data_ = vector;
data_ /= data_.two_norm();
......@@ -100,7 +101,7 @@ public:
/** \brief The exponential map */
static UnitVector exp(const UnitVector& p, const TangentVector& v) {
Dune::FieldMatrix<double,N-1,N> frame = p.orthonormalFrame();
Dune::FieldMatrix<T,N-1,N> frame = p.orthonormalFrame();
EmbeddedTangentVector ev;
frame.mtv(v,ev);
......@@ -113,7 +114,7 @@ public:
assert( std::abs(p.data_*v) < 1e-5 );
const double norm = v.two_norm();
const T norm = v.two_norm();
UnitVector result = p;
result.data_ *= std::cos(norm);
result.data_.axpy(sinc(norm), v);
......@@ -121,12 +122,12 @@ public:
}
/** \brief Length of the great arc connecting the two points */
static double distance(const UnitVector& a, const UnitVector& b) {
static T distance(const UnitVector& a, const UnitVector& b) {
// Not nice: we are in a class for unit vectors, but the class is actually
// supposed to handle perturbations of unit vectors as well. Therefore
// we normalize here.
double x = a.data_ * b.data_/a.data_.two_norm()/b.data_.two_norm();
T x = a.data_ * b.data_/a.data_.two_norm()/b.data_.two_norm();
// paranoia: if the argument is just eps larger than 1 acos returns NaN
x = std::min(x,1.0);
......@@ -139,7 +140,7 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& a, const UnitVector& b) {
double x = a.data_ * b.data_;
T x = a.data_ * b.data_;
EmbeddedTangentVector result = a.data_;
......@@ -158,13 +159,13 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
double sp = p.data_ * q.data_;
T sp = p.data_ * q.data_;
Dune::FieldVector<double,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
Dune::FieldMatrix<double,N,N> A;
Dune::FieldMatrix<T,N,N> A;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
A[i][j] = pProjected[i]*pProjected[j];
......@@ -172,7 +173,7 @@ public:
A *= secondDerivativeOfArcCosSquared(sp);
// Compute matrix B (see notes)
Dune::FieldMatrix<double,N,N> Pq;
Dune::FieldMatrix<T,N,N> Pq;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
Pq[i][j] = (i==j) - q.data_[i]*q.data_[j];
......@@ -187,33 +188,33 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const UnitVector& a, const UnitVector& b) {
static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const UnitVector& a, const UnitVector& b) {
double sp = a.data_ * b.data_;
T sp = a.data_ * b.data_;
// Compute vector A (see notes)
Dune::FieldMatrix<double,1,N> row;
Dune::FieldMatrix<T,1,N> row;
row[0] = b.projectOntoTangentSpace(a.globalCoordinates());
Dune::FieldVector<double,N> tmp = a.projectOntoTangentSpace(b.globalCoordinates());
Dune::FieldMatrix<double,N,1> column;
Dune::FieldVector<T,N> tmp = a.projectOntoTangentSpace(b.globalCoordinates());
Dune::FieldMatrix<T,N,1> column;
for (int i=0; i<N; i++) // turn row vector into column vector
column[i] = tmp[i];
Dune::FieldMatrix<double,N,N> A;
Dune::FieldMatrix<T,N,N> A;
// A = row * column
Dune::FMatrixHelp::multMatrix(column,row,A);
A *= secondDerivativeOfArcCosSquared(sp);
// Compute matrix B (see notes)
Dune::FieldMatrix<double,N,N> Pp, Pq;
Dune::FieldMatrix<T,N,N> Pp, Pq;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++) {
Pp[i][j] = (i==j) - a.data_[i]*a.data_[j];
Pq[i][j] = (i==j) - b.data_[i]*b.data_[j];
}
Dune::FieldMatrix<double,N,N> B;
Dune::FieldMatrix<T,N,N> B;
Dune::FMatrixHelp::multMatrix(Pp,Pq,B);
// Bring it all together
......@@ -227,19 +228,19 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Tensor3<double,N,N,N> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
Tensor3<double,N,N,N> result;
Tensor3<T,N,N,N> result;
double sp = p.data_ * q.data_;
T sp = p.data_ * q.data_;
// The projection matrix onto the tangent space at p and q
Dune::FieldMatrix<double,N,N> Pq;
Dune::FieldMatrix<T,N,N> Pq;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
Pq[i][j] = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
Dune::FieldVector<double,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
......@@ -262,24 +263,24 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
static Tensor3<double,N,N,N> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const UnitVector& p, const UnitVector& q) {
static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const UnitVector& p, const UnitVector& q) {
Tensor3<double,N,N,N> result;
Tensor3<T,N,N,N> result;
double sp = p.data_ * q.data_;
T sp = p.data_ * q.data_;
// The projection matrix onto the tangent space at p and q
Dune::FieldMatrix<double,N,N> Pp, Pq;
Dune::FieldMatrix<T,N,N> Pp, Pq;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++) {
Pp[i][j] = (i==j) - p.globalCoordinates()[i]*p.globalCoordinates()[j];
Pq[i][j] = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
}
Dune::FieldVector<double,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
Dune::FieldVector<double,N> qProjected = p.projectOntoTangentSpace(q.globalCoordinates());
Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
Dune::FieldVector<T,N> qProjected = p.projectOntoTangentSpace(q.globalCoordinates());
Tensor3<double,N,N,N> derivativeOfPqOTimesPq;
Tensor3<T,N,N,N> derivativeOfPqOTimesPq;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
for (int k=0; k<N; k++) {
......@@ -288,10 +289,10 @@ public:
derivativeOfPqOTimesPq[i][j][k] += Pp[i][l] * (Pq[j][l]*pProjected[k] + pProjected[j]*Pq[k][l]);
}
result = thirdDerivativeOfArcCosSquared(sp) * Tensor3<double,N,N,N>::product(qProjected,pProjected,pProjected)
result = thirdDerivativeOfArcCosSquared(sp) * Tensor3<T,N,N,N>::product(qProjected,pProjected,pProjected)
+ secondDerivativeOfArcCosSquared(sp) * derivativeOfPqOTimesPq
- secondDerivativeOfArcCosSquared(sp) * sp * Tensor3<double,N,N,N>::product(qProjected,Pq)
- derivativeOfArcCosSquared(sp) * Tensor3<double,N,N,N>::product(qProjected,Pq);
- secondDerivativeOfArcCosSquared(sp) * sp * Tensor3<T,N,N,N>::product(qProjected,Pq)
- derivativeOfArcCosSquared(sp) * Tensor3<T,N,N,N>::product(qProjected,Pq);
return result;
}
......@@ -313,12 +314,12 @@ public:
This basis is of course not globally continuous.
*/
Dune::FieldMatrix<double,N-1,N> orthonormalFrame() const {
Dune::FieldMatrix<T,N-1,N> orthonormalFrame() const {
Dune::FieldMatrix<double,N-1,N> result;
Dune::FieldMatrix<T,N-1,N> result;
// Coordinates of the stereographic projection
Dune::FieldVector<double,N-1> X;
Dune::FieldVector<T,N-1> X;
if (data_[N-1] <= 0) {
......@@ -334,7 +335,7 @@ public:
}
double RSquared = X.two_norm2();
T RSquared = X.two_norm2();
for (size_t i=0; i<N-1; i++)
for (size_t j=0; j<N-1; j++)
......@@ -365,7 +366,7 @@ public:
private:
Dune::FieldVector<double,N> data_;
Dune::FieldVector<T,N> data_;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment