From 78d4007ece60befced560e01c479525249c8488c Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 7 Nov 2010 12:06:13 +0000
Subject: [PATCH] implement an orthogonal basis of the tangent spaces for all
dimensions (not just 2)
[[Imported from SVN: r6501]]
---
dune/gfe/unitvector.hh | 36 +++++++++++++++++++++++++++++++-----
1 file changed, 31 insertions(+), 5 deletions(-)
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index ee4da9df..d71aaf09 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -293,12 +293,38 @@ public:
Dune::FieldMatrix<double,N-1,N> orthonormalFrame() const {
Dune::FieldMatrix<double,N-1,N> result;
+
+ // Coordinates of the stereographic projection
+ Dune::FieldVector<double,N-1> X;
- if (N==2) {
- result[0][0] = -data_[1];
- result[0][1] = data_[0];
- } else
- DUNE_THROW(Dune::NotImplemented, "orthonormalFrame for N!=2!");
+ if (data_[N-1] <= 0) {
+
+ // Stereographic projection from the north pole onto R^{N-1}
+ for (size_t i=0; i<N-1; i++)
+ X[i] = data_[i] / (1-data_[N-1]);
+
+ } else {
+
+ // Stereographic projection from the south pole onto R^{N-1}
+ for (size_t i=0; i<N-1; i++)
+ X[i] = data_[i] / (1+data_[N-1]);
+
+ }
+
+ double RSquared = X.two_norm2();
+
+ for (size_t i=0; i<N-1; i++)
+ for (size_t j=0; j<N-1; j++)
+ // Note: the matrix is the transpose of the one in the paper
+ result[j][i] = 2*(i==j)*(1+RSquared) - 4*X[i]*X[j];
+
+ for (size_t j=0; j<N-1; j++)
+ result[j][N-1] = 4*X[j];
+
+ // Upper hemisphere: adapt formulas so it is the stereographic projection from the south pole
+ if (data_[N-1] > 0)
+ for (size_t j=0; j<N-1; j++)
+ result[j][N-1] *= -1;
return result;
}
--
GitLab