diff --git a/Numeric/GmshMatrix.cpp b/Numeric/fullMatrix.cpp
similarity index 87%
rename from Numeric/GmshMatrix.cpp
rename to Numeric/fullMatrix.cpp
index 286df546a452f2a30e0d319a7a4c5ea9e7a17540..69928a6f359742a691b557f0ee9de49e18ec2f62 100644
--- a/Numeric/GmshMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -6,7 +6,7 @@
#include <complex>
#include <string.h>
#include "GmshConfig.h"
-#include "GmshMatrix.h"
+#include "fullMatrix.h"
#include "GmshMessage.h"
//#if defined(_MSC_VER)
@@ -39,7 +39,7 @@ extern "C" {
}
template<>
-void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c)
+void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c)
{
int M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1();
@@ -49,8 +49,8 @@ void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c
}
template<>
-void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<double> > &b,
- gmshMatrix<std::complex<double> > &c)
+void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<double> > &b,
+ fullMatrix<std::complex<double> > &c)
{
int M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1();
@@ -60,7 +60,7 @@ void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<doubl
}
template<>
-void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b,
+void fullMatrix<double>::gemm(fullMatrix<double> &a, fullMatrix<double> &b,
double alpha, double beta)
{
int M = size1(), N = size2(), K = a.size2();
@@ -70,8 +70,8 @@ void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b,
}
template<>
-void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &a,
- gmshMatrix<std::complex<double> > &b,
+void fullMatrix<std::complex<double> >::gemm(fullMatrix<std::complex<double> > &a,
+ fullMatrix<std::complex<double> > &b,
std::complex<double> alpha,
std::complex<double> beta)
{
@@ -82,7 +82,7 @@ void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &
}
template<>
-void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y)
+void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y)
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
double alpha = 1., beta = 0.;
@@ -91,8 +91,8 @@ void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y
}
template<>
-void gmshMatrix<std::complex<double> >::mult(const gmshVector<std::complex<double> > &x,
- gmshVector<std::complex<double> > &y)
+void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<double> > &x,
+ fullVector<std::complex<double> > &y)
{
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
std::complex<double> alpha = 1., beta = 0.;
@@ -123,7 +123,7 @@ extern "C" {
}
template<>
-bool gmshMatrix<double>::invertInPlace()
+bool fullMatrix<double>::invertInPlace()
{
int N = size1(), nrhs = N, lda = N, ldb = N, info;
int *ipiv = new int[N];
@@ -148,8 +148,8 @@ bool gmshMatrix<double>::invertInPlace()
template<>
-bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR,
- gmshVector<double> &DI, gmshMatrix<double> &VR,
+bool fullMatrix<double>::eig(fullMatrix<double> &VL, fullVector<double> &DR,
+ fullVector<double> &DI, fullMatrix<double> &VR,
bool sortRealPart)
{
int N = size1(), info;
@@ -193,7 +193,7 @@ bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR,
}
template<>
-bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<double> &result)
+bool fullMatrix<double>::lu_solve(const fullVector<double> &rhs, fullVector<double> &result)
{
int N = size1(), nrhs = 1, lda = N, ldb = N, info;
int *ipiv = new int[N];
@@ -209,7 +209,7 @@ bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<doub
}
template<>
-bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
+bool fullMatrix<double>::invert(fullMatrix<double> &result)
{
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
@@ -231,9 +231,9 @@ bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
}
template<>
-double gmshMatrix<double>::determinant() const
+double fullMatrix<double>::determinant() const
{
- gmshMatrix<double> tmp(*this);
+ fullMatrix<double> tmp(*this);
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
F77NAME(dgetrf)(&M, &N, tmp._data, &lda, ipiv, &info);
@@ -253,12 +253,12 @@ double gmshMatrix<double>::determinant() const
}
template<>
-bool gmshMatrix<double>::svd(gmshMatrix<double> &V, gmshVector<double> &S)
+bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
{
- gmshMatrix<double> VT(V.size2(), V.size1());
+ fullMatrix<double> VT(V.size2(), V.size1());
int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
- gmshVector<double> WORK(LWORK);
+ fullVector<double> WORK(LWORK);
F77NAME(dgesvd)("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
VT._data, &LDVT, WORK._data, &LWORK, &info);
V = VT.transpose();
diff --git a/Numeric/GmshMatrix.h b/Numeric/fullMatrix.h
similarity index 77%
rename from Numeric/GmshMatrix.h
rename to Numeric/fullMatrix.h
index f8604f1312e413349aaf451ec78466c3b9a31ea3..982a86dc17373729d29adfbe5d5ad82c1ebd33b7 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -11,27 +11,27 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
-template <class scalar> class gmshMatrix;
+template <class scalar> class fullMatrix;
template <class scalar>
-class gmshVector
+class fullVector
{
private:
int _r;
scalar *_data;
- friend class gmshMatrix<scalar>;
+ friend class fullMatrix<scalar>;
public:
- gmshVector(int r) : _r(r)
+ fullVector(int r) : _r(r)
{
_data = new scalar[_r];
scale(0.);
}
- gmshVector(const gmshVector<scalar> &other) : _r(other._r)
+ fullVector(const fullVector<scalar> &other) : _r(other._r)
{
_data = new scalar[_r];
for(int i = 0; i < _r; ++i) _data[i] = other._data[i];
}
- ~gmshVector() { if(_data) delete [] _data; }
+ ~fullVector() { if(_data) delete [] _data; }
inline scalar operator () (int i) const
{
return _data[i];
@@ -55,7 +55,7 @@ class gmshVector
for(int i = 0; i < _r; ++i) _data[i] *= s;
}
- inline scalar operator *(const gmshVector<scalar> & other)
+ inline scalar operator *(const fullVector<scalar> & other)
{
scalar s = 0.0;
for(int i = 0; i < _r; ++i) s += _data[i]*other._data[i];
@@ -73,27 +73,27 @@ class gmshVector
};
template <class scalar>
-class gmshMatrix
+class fullMatrix
{
private:
int _r, _c;
scalar *_data;
public:
- gmshMatrix(int r, int c) : _r(r), _c(c)
+ fullMatrix(int r, int c) : _r(r), _c(c)
{
_data = new scalar[_r * _c];
scale(0.);
}
- gmshMatrix(const gmshMatrix<scalar> &other) : _r(other._r), _c(other._c)
+ fullMatrix(const fullMatrix<scalar> &other) : _r(other._r), _c(other._c)
{
_data = new scalar[_r * _c];
gM_memcpy(other);
}
- gmshMatrix() : _r(0), _c(0), _data(0) {}
- ~gmshMatrix() { if(_data) delete [] _data; }
+ fullMatrix() : _r(0), _c(0), _data(0) {}
+ ~fullMatrix() { if(_data) delete [] _data; }
inline int size1() const { return _r; }
inline int size2() const { return _c; }
- gmshMatrix<scalar> & operator = (const gmshMatrix<scalar> &other)
+ fullMatrix<scalar> & operator = (const fullMatrix<scalar> &other)
{
if(this != &other){
_r = other._r;
@@ -103,7 +103,7 @@ class gmshMatrix
}
return *this;
}
- void gM_memcpy(const gmshMatrix<scalar> &other)
+ void gM_memcpy(const fullMatrix<scalar> &other)
{
for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
}
@@ -115,14 +115,14 @@ class gmshMatrix
{
return _data[i + _r * j];
}
- void copy(const gmshMatrix<scalar> &a, int i0, int ni, int j0, int nj,
+ void copy(const fullMatrix<scalar> &a, int i0, int ni, int j0, int nj,
int desti0, int destj0)
{
for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
(*this)(desti, destj) = a(i, j);
}
- void mult(const gmshMatrix<scalar> &b, gmshMatrix<scalar> &c)
+ void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)
#if !defined(HAVE_BLAS)
{
c.scale(0.);
@@ -133,11 +133,11 @@ class gmshMatrix
}
#endif
;
- void gemm(gmshMatrix<scalar> &a, gmshMatrix<scalar> &b,
+ void gemm(fullMatrix<scalar> &a, fullMatrix<scalar> &b,
scalar alpha=1., scalar beta=1.)
#if !defined(HAVE_BLAS)
{
- gmshMatrix<scalar> temp(a.size1(), b.size2());
+ fullMatrix<scalar> temp(a.size1(), b.size2());
a.mult(b, temp);
temp.scale(alpha);
scale(beta);
@@ -160,13 +160,13 @@ class gmshMatrix
{
for(int i = 0; i < _r * _c; ++i) _data[i] += a;
}
- inline void add(const gmshMatrix<scalar> &m)
+ inline void add(const fullMatrix<scalar> &m)
{
for(int i = 0; i < size1(); i++)
for(int j = 0; j < size2(); j++)
(*this)(i, j) += m(i, j);
}
- void mult(const gmshVector<scalar> &x, gmshVector<scalar> &y)
+ void mult(const fullVector<scalar> &x, fullVector<scalar> &y)
#if !defined(HAVE_BLAS)
{
y.scale(0.);
@@ -176,9 +176,9 @@ class gmshMatrix
}
#endif
;
- inline gmshMatrix<scalar> transpose()
+ inline fullMatrix<scalar> transpose()
{
- gmshMatrix<scalar> T(size2(), size1());
+ fullMatrix<scalar> T(size2(), size1());
for(int i = 0; i < size1(); i++)
for(int j = 0; j < size2(); j++)
T(j, i) = (*this)(i, j);
@@ -197,7 +197,7 @@ class gmshMatrix
_data[j + _r * i] = t;
}
}
- bool lu_solve(const gmshVector<scalar> &rhs, gmshVector<scalar> &result)
+ bool lu_solve(const fullVector<scalar> &rhs, fullVector<scalar> &result)
#if !defined(HAVE_LAPACK)
{
Msg::Error("LU factorization requires LAPACK");
@@ -213,10 +213,10 @@ class gmshMatrix
}
#endif
;
- bool eig(gmshMatrix<scalar> &VL, // left eigenvectors
- gmshVector<double> &DR, // Real part of eigenvalues
- gmshVector<double> &DI, // Im part of eigen
- gmshMatrix<scalar> &VR,
+ bool eig(fullMatrix<scalar> &VL, // left eigenvectors
+ fullVector<double> &DR, // Real part of eigenvalues
+ fullVector<double> &DI, // Im part of eigen
+ fullMatrix<scalar> &VR,
bool sortRealPart=false) // if true: sorted from min 'DR' to max 'DR'
#if !defined(HAVE_LAPACK)
{
@@ -225,7 +225,7 @@ class gmshMatrix
}
#endif
;
- bool invert(gmshMatrix<scalar> &result)
+ bool invert(fullMatrix<scalar> &result)
#if !defined(HAVE_LAPACK)
{
Msg::Error("LU factorization requires LAPACK");
@@ -233,11 +233,11 @@ class gmshMatrix
}
#endif
;
- gmshMatrix<scalar> cofactor(int i, int j) const
+ fullMatrix<scalar> cofactor(int i, int j) const
{
int ni = size1();
int nj = size2();
- gmshMatrix<scalar> cof(ni - 1, nj - 1);
+ fullMatrix<scalar> cof(ni - 1, nj - 1);
for(int I = 0; I < ni; I++){
for(int J = 0; J < nj; J++){
if(J != j && I != i)
@@ -254,7 +254,7 @@ class gmshMatrix
}
#endif
;
- bool svd(gmshMatrix<scalar> &V, gmshVector<scalar> &S)
+ bool svd(fullMatrix<scalar> &V, fullVector<scalar> &S)
#if !defined(HAVE_LAPACK)
{
Msg::Error("Singular value decomposition requires LAPACK");
diff --git a/Numeric/functionSpace.cpp b/Numeric/functionSpace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..60dbf13ce84f542c02218038f5792e7e846d2723
--- /dev/null
+++ b/Numeric/functionSpace.cpp
@@ -0,0 +1,672 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Koen Hillewaert
+//
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include "functionSpace.h"
+
+static fullMatrix<double> generate1DMonomials(int order)
+{
+ fullMatrix<double> monomials(order + 1, 1);
+ for (int i = 0; i < order + 1; i++) monomials(i, 0) = i;
+ return monomials;
+}
+
+static fullMatrix<double> generate1DPoints(int order)
+{
+ fullMatrix<double> line(order + 1, 1);
+ line(0, 0) = -1.;
+ line(1, 0) = 1.;
+ double dd = 2. / order;
+ for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
+ return line;
+}
+
+static fullMatrix<double> generatePascalTriangle(int order)
+{
+ fullMatrix<double> monomials((order + 1) * (order + 2) / 2, 2);
+ int index = 0;
+ for (int i = 0; i <= order; i++) {
+ for (int j = 0; j <= i; j++) {
+ monomials(index, 0) = i - j;
+ monomials(index, 1) = j;
+ index++;
+ }
+ }
+ return monomials;
+}
+
+// generate the exterior hull of the Pascal triangle for Serendipity element
+
+static fullMatrix<double> generatePascalSerendipityTriangle(int order)
+{
+ fullMatrix<double> monomials(3 * order, 2);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+
+ int index = 1;
+ for (int i = 1; i <= order; i++) {
+ if (i == order) {
+ for (int j = 0; j <= i; j++) {
+ monomials(index, 0) = i - j;
+ monomials(index, 1) = j;
+ index++;
+ }
+ }
+ else {
+ monomials(index, 0) = i;
+ monomials(index, 1) = 0;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = i;
+ index++;
+ }
+ }
+ return monomials;
+}
+
+// generate the monomials subspace of all monomials of order exactly == p
+
+static fullMatrix<double> generateMonomialSubspace(int dim, int p)
+{
+ fullMatrix<double> monomials;
+
+ switch (dim) {
+ case 1:
+ monomials = fullMatrix<double>(1, 1);
+ monomials(0, 0) = p;
+ break;
+ case 2:
+ monomials = fullMatrix<double>(p + 1, 2);
+ for (int k = 0; k <= p; k++) {
+ monomials(k, 0) = p - k;
+ monomials(k, 1) = k;
+ }
+ break;
+ case 3:
+ monomials = fullMatrix<double>((p + 1) * (p + 2) / 2, 3);
+ int index = 0;
+ for (int i = 0; i <= p; i++) {
+ for (int k = 0; k <= p - i; k++) {
+ monomials(index, 0) = p - i - k;
+ monomials(index, 1) = k;
+ monomials(index, 2) = i;
+ index++;
+ }
+ }
+ break;
+ }
+ return monomials;
+}
+
+// generate external hull of the Pascal tetrahedron
+
+static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
+{
+ int nbMonomials = 4 + 6 * std::max(0, order - 1) +
+ 4 * std::max(0, (order - 2) * (order - 1) / 2);
+ fullMatrix<double> monomials(nbMonomials, 3);
+
+ monomials.set_all(0);
+ int index = 1;
+ for (int p = 1; p < order; p++) {
+ for (int i = 0; i < 3; i++) {
+ int j = (i + 1) % 3;
+ int k = (i + 2) % 3;
+ for (int ii = 0; ii < p; ii++) {
+ monomials(index, i) = p - ii;
+ monomials(index, j) = ii;
+ monomials(index, k) = 0;
+ index++;
+ }
+ }
+ }
+ fullMatrix<double> monomialsMaxOrder = generateMonomialSubspace(3, order);
+ int nbMaxOrder = monomialsMaxOrder.size1();
+ monomials.copy(monomialsMaxOrder, 0, nbMaxOrder, 0, 3, index, 0);
+ return monomials;
+}
+
+// generate Pascal tetrahedron
+
+static fullMatrix<double> generatePascalTetrahedron(int order)
+{
+ int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
+
+ fullMatrix<double> monomials(nbMonomials, 3);
+
+ int index = 0;
+ for (int p = 0; p <= order; p++) {
+ fullMatrix<double> monOrder = generateMonomialSubspace(3, p);
+ int nb = monOrder.size1();
+ monomials.copy(monOrder, 0, nb, 0, 3, index, 0);
+ index += nb;
+ }
+
+ return monomials;
+}
+
+static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
+//static int nbdoftriangleserendip(int order) { return 3 * order; }
+
+//KH : caveat : node coordinates are not yet coherent with node numbering associated
+// to numbering of principal vertices of face !!!!
+
+// uv surface - orientation v0-v2-v1
+static void nodepositionface0(int order, double *u, double *v, double *w)
+{
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= order; v[2]= 0.; w[2] = 0.;
+
+ // edges
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = 0.;
+ v[3 + k] = k + 1;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = k + 1;
+ v[3 + order - 1 + k] = order - 1 - k ;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = 0.;
+ }
+
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3* (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
+ }
+ }
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+// uw surface - orientation v0-v1-v3
+static void nodepositionface1(int order, double *u, double *v, double *w)
+{
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+ u[0] = 0.; v[0]= 0.; w[0] = 0.;
+ u[1] = order; v[1]= 0.; w[1] = 0.;
+ u[2] = 0.; v[2]= 0.; w[2] = order;
+ // edges
+ for (int k = 0; k < (order - 1); k++){
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = 0.;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = 0.;
+ w[3 + order - 1+ k ] = k + 1;
+
+ u[3 + 2 * (order - 1) + k] = 0. ;
+ v[3 + 2 * (order - 1) + k] = 0.;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
+ }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+// vw surface - orientation v0-v3-v2
+static void nodepositionface2(int order, double *u, double *v, double *w)
+{
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= 0.; w[1] = order;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
+ // edges
+ for (int k = 0; k < (order - 1); k++){
+
+ u[3 + k] = 0.;
+ v[3 + k] = 0.;
+ w[3 + k] = k + 1;
+
+ u[3 + order - 1 + k] = 0.;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = order - 1 - k;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = 0.;
+ }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+// uvw surface - orientation v3-v1-v2
+static void nodepositionface3(int order, double *u, double *v, double *w)
+{
+ int ndofT = nbdoftriangle(order);
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+ u[0]= 0.; v[0]= 0.; w[0] = order;
+ u[1]= order; v[1]= 0.; w[1] = 0.;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
+ // edges
+ for (int k = 0; k < (order - 1); k++){
+
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = order - 1 - k;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = 0.;
+
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = k + 1;
+ }
+ if (order > 2){
+ int nbdoftemp = nbdoftriangle(order - 3);
+ nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ &w[3 + 3 * (order - 1)]);
+ for (int k = 0; k < nbdoftemp; k++){
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ }
+ }
+ for (int k = 0; k < ndofT; k++){
+ u[k] = u[k] / order;
+ v[k] = v[k] / order;
+ w[k] = w[k] / order;
+ }
+}
+
+static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
+{
+ int nbPoints =
+ (serendip ?
+ 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
+ (order + 1) * (order + 2) * (order + 3) / 6);
+
+ fullMatrix<double> point(nbPoints, 3);
+
+ double overOrder = (order == 0 ? 1. : 1. / order);
+
+ point(0, 0) = 0.;
+ point(0, 1) = 0.;
+ point(0, 2) = 0.;
+
+ if (order > 0) {
+ point(1, 0) = order;
+ point(1, 1) = 0;
+ point(1, 2) = 0;
+
+ point(2, 0) = 0.;
+ point(2, 1) = order;
+ point(2, 2) = 0.;
+
+ point(3, 0) = 0.;
+ point(3, 1) = 0.;
+ point(3, 2) = order;
+
+ // edges e5 and e6 switched in original version, opposite direction
+ // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
+
+ if (order > 1) {
+ for (int k = 0; k < (order - 1); k++) {
+ point(4 + k, 0) = k + 1;
+ point(4 + order - 1 + k, 0) = order - 1 - k;
+ point(4 + 2 * (order - 1) + k, 0) = 0.;
+ point(4 + 3 * (order - 1) + k, 0) = 0.;
+ // point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
+ // point(4 + 5 * (order - 1) + k, 0) = 0.;
+ point(4 + 4 * (order - 1) + k, 0) = 0.;
+ point(4 + 5 * (order - 1) + k, 0) = k+1;
+
+ point(4 + k, 1) = 0.;
+ point(4 + order - 1 + k, 1) = k + 1;
+ point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 3 * (order - 1) + k, 1) = 0.;
+ // point(4 + 4 * (order - 1) + k, 1) = 0.;
+ // point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 1) = k+1;
+ point(4 + 5 * (order - 1) + k, 1) = 0.;
+
+ point(4 + k, 2) = 0.;
+ point(4 + order - 1 + k, 2) = 0.;
+ point(4 + 2 * (order - 1) + k, 2) = 0.;
+ point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+ }
+
+ if (order > 2) {
+ int ns = 4 + 6 * (order - 1);
+ int nbdofface = nbdoftriangle(order - 3);
+
+ double *u = new double[nbdofface];
+ double *v = new double[nbdofface];
+ double *w = new double[nbdofface];
+
+ nodepositionface0(order - 3, u, v, w);
+
+ // u-v plane
+
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3);
+ }
+
+ ns = ns + nbdofface;
+
+ // u-w plane
+
+ nodepositionface1(order - 3, u, v, w);
+
+ for (int i=0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) ;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ // v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface2(order - 3, u, v, w);
+
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3);
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ // u-v-w plane
+
+ ns = ns + nbdofface;
+
+ nodepositionface3(order - 3, u, v, w);
+
+ for (int i = 0; i < nbdofface; i++){
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 1) = v[i] * (order - 3) + 1.;
+ point(ns + i, 2) = w[i] * (order - 3) + 1.;
+ }
+
+ ns = ns + nbdofface;
+
+ delete [] u;
+ delete [] v;
+ delete [] w;
+
+ if (!serendip && order > 3) {
+
+ fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
+ for (int k = 0; k < interior.size1(); k++) {
+ point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
+ point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
+ point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
+ }
+ }
+ }
+ }
+ }
+
+ point.scale(overOrder);
+ return point;
+}
+
+static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
+{
+ int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
+ fullMatrix<double> point(nbPoints, 2);
+
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+
+ double dd = 1. / order;
+
+ if (order > 0) {
+ point(1, 0) = 1;
+ point(1, 1) = 0;
+ point(2, 0) = 0;
+ point(2, 1) = 1;
+
+ int index = 3;
+
+ if (order > 1) {
+
+ double ksi = 0;
+ double eta = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi += dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ ksi = 1.;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi -= dd;
+ eta += dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ eta = 1.;
+ ksi = 0.;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ eta -= dd;
+ point(index, 0) = ksi;
+ point(index, 1) = eta;
+ }
+
+ if (order > 2 && !serendip) {
+ fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
+ inner.scale(1. - 3. * dd);
+ inner.add(dd);
+ point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+ return point;
+}
+
+static fullMatrix<double> generateLagrangeMonomialCoefficients
+ (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
+{
+ if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
+ return fullMatrix<double>(1, 1);
+ }
+
+ int ndofs = monomial.size1();
+ int dim = monomial.size2();
+
+ fullMatrix<double> Vandermonde(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+ for (int k = 0; k < dim; k++) dd *= pow(point(j, k), monomial(i, k));
+ Vandermonde(i, j) = dd;
+ }
+ }
+
+ fullMatrix<double> coefficient(ndofs, ndofs);
+ Vandermonde.invert(coefficient);
+ return coefficient;
+}
+
+std::map<int, functionSpace> functionSpaces::fs;
+
+const functionSpace &functionSpaces::find(int tag)
+{
+ std::map<int, functionSpace>::const_iterator it = fs.find(tag);
+ if (it != fs.end()) return it->second;
+
+ functionSpace F;
+
+ switch (tag){
+ case MSH_PNT:
+ F.monomials = generate1DMonomials(0);
+ F.points = generate1DPoints(0);
+ break;
+ case MSH_LIN_2 :
+ F.monomials = generate1DMonomials(1);
+ F.points = generate1DPoints(1);
+ break;
+ case MSH_LIN_3 :
+ F.monomials = generate1DMonomials(2);
+ F.points = generate1DPoints(2);
+ break;
+ case MSH_LIN_4:
+ F.monomials = generate1DMonomials(3);
+ F.points = generate1DPoints(3);
+ break;
+ case MSH_LIN_5:
+ F.monomials = generate1DMonomials(4);
+ F.points = generate1DPoints(4);
+ break;
+ case MSH_LIN_6:
+ F.monomials = generate1DMonomials(5);
+ F.points = generate1DPoints(5);
+ break;
+ case MSH_TRI_3 :
+ F.monomials = generatePascalTriangle(1);
+ F.points = gmshGeneratePointsTriangle(1, false);
+ break;
+ case MSH_TRI_6 :
+ F.monomials = generatePascalTriangle(2);
+ F.points = gmshGeneratePointsTriangle(2, false);
+ break;
+ case MSH_TRI_9 :
+ F.monomials = generatePascalSerendipityTriangle(3);
+ F.points = gmshGeneratePointsTriangle(3, true);
+ break;
+ case MSH_TRI_10 :
+ F.monomials = generatePascalTriangle(3);
+ F.points = gmshGeneratePointsTriangle(3, false);
+ break;
+ case MSH_TRI_12 :
+ F.monomials = generatePascalSerendipityTriangle(4);
+ F.points = gmshGeneratePointsTriangle(4, true);
+ break;
+ case MSH_TRI_15 :
+ F.monomials = generatePascalTriangle(4);
+ F.points = gmshGeneratePointsTriangle(4, false);
+ break;
+ case MSH_TRI_15I :
+ F.monomials = generatePascalSerendipityTriangle(5);
+ F.points = gmshGeneratePointsTriangle(5, true);
+ break;
+ case MSH_TRI_21 :
+ F.monomials = generatePascalTriangle(5);
+ F.points = gmshGeneratePointsTriangle(5, false);
+ break;
+ case MSH_TET_4 :
+ F.monomials = generatePascalTetrahedron(1);
+ F.points = gmshGeneratePointsTetrahedron(1, false);
+ break;
+ case MSH_TET_10 :
+ F.monomials = generatePascalTetrahedron(2);
+ F.points = gmshGeneratePointsTetrahedron(2, false);
+ break;
+ case MSH_TET_20 :
+ F.monomials = generatePascalTetrahedron(3);
+ F.points = gmshGeneratePointsTetrahedron(3, false);
+ break;
+ case MSH_TET_35 :
+ F.monomials = generatePascalTetrahedron(4);
+ F.points = gmshGeneratePointsTetrahedron(4, false);
+ break;
+ case MSH_TET_34 :
+ F.monomials = generatePascalSerendipityTetrahedron(4);
+ F.points = gmshGeneratePointsTetrahedron(4, true);
+ break;
+ case MSH_TET_52 :
+ F.monomials = generatePascalSerendipityTetrahedron(5);
+ F.points = gmshGeneratePointsTetrahedron(5, true);
+ break;
+ case MSH_TET_56 :
+ F.monomials = generatePascalTetrahedron(5);
+ F.points = gmshGeneratePointsTetrahedron(5, false);
+ break;
+ default :
+ Msg::Error("Unknown function space %d: reverting to TET_4", tag);
+ F.monomials = generatePascalTetrahedron(1);
+ F.points = gmshGeneratePointsTetrahedron(1, false);
+ break;
+ }
+ F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
+ fs.insert(std::make_pair(tag, F));
+ return fs[tag];
+}
+
+std::map<std::pair<int, int>, fullMatrix<double> > functionSpaces::injector;
+
+const fullMatrix<double> &functionSpaces::findInjector(int tag1, int tag2)
+{
+ std::pair<int,int> key(tag1,tag2);
+ std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
+ if (it != injector.end()) return it->second;
+
+ const functionSpace& fs1 = find(tag1);
+ const functionSpace& fs2 = find(tag2);
+
+ fullMatrix<double> inj(fs1.points.size1(),fs2.points.size1());
+
+ double sf[256];
+
+ for (int i = 0; i < fs1.points.size1(); i++) {
+ fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
+ for (int j = 0; j < fs2.points.size1(); j++) inj(i, j) = sf[j];
+ }
+
+ injector.insert(std::make_pair(key, inj));
+ return injector[key];
+}
diff --git a/Numeric/functionSpace.h b/Numeric/functionSpace.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ba14b7d84ca57c12faf1a762819c1b7c3863208
--- /dev/null
+++ b/Numeric/functionSpace.h
@@ -0,0 +1,107 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _FUNCTION_SPACE_H_
+#define _FUNCTION_SPACE_H_
+
+#include <math.h>
+#include <map>
+#include "fullMatrix.h"
+
+struct functionSpace
+{
+ fullMatrix<double> points;
+ fullMatrix<double> monomials;
+ fullMatrix<double> coefficients;
+ inline void evaluateMonomials(double u, double v, double w, double p[]) const
+ {
+ for (int j = 0; j < monomials.size1(); j++) {
+ p[j] = pow(u, (int)monomials(j, 0));
+ if (monomials.size2() > 1) p[j] *= pow(v, (int)monomials(j, 1));
+ if (monomials.size2() > 2) p[j] *= pow(w, (int)monomials(j, 2));
+ }
+ }
+ inline void f(double u, double v, double w, double *sf) const
+ {
+ double p[256];
+ evaluateMonomials(u, v, w, p);
+ for (int i = 0; i < coefficients.size1(); i++) {
+ sf[i] = 0;
+ for (int j = 0; j < coefficients.size2(); j++) {
+ sf[i] += coefficients(i, j) * p[j];
+ }
+ }
+ }
+ inline void df(double u, double v, double w, double grads[][3]) const
+ {
+ switch (monomials.size2()) {
+ case 1:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0);
+ }
+ }
+ break;
+ case 2:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1));
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1);
+ }
+ }
+ break;
+ case 3:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 2) > 0)
+ grads[i][2] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2) - 1) * (monomials)(j, 2);
+ }
+ }
+ break;
+ }
+ }
+};
+
+class functionSpaces
+{
+ private:
+ static std::map<int, functionSpace> fs;
+ static std::map<std::pair<int, int>, fullMatrix<double> > injector;
+ public :
+ static const functionSpace &find(int);
+ static const fullMatrix<double> &findInjector(int, int);
+};
+
+#endif
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 629ecf8f0264d598c906cfec2944550903d328af..c82b6eb39635ee862cdbf05ee9f209259c494b54 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -93,7 +93,7 @@ bool PViewDataGModel::finalize()
void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
{
int type = e->getType();
- const gmshFunctionSpace *fs = e->getFunctionSpace();
+ const functionSpace *fs = e->getFunctionSpace();
if(fs){
if(e->getPolynomialOrder() > 1)
setInterpolationMatrices(type, fs->coefficients, fs->monomials,