Skip to content
Snippets Groups Projects
Commit 0801d0f9 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent 9b8fa83e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/GmshMatrix.cpp → Numeric/fullMatrix.cpp
+ 20
− 20
View file @ 0801d0f9
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
#include <complex> #include <complex>
#include <string.h> #include <string.h>
#include "GmshConfig.h" #include "GmshConfig.h"
#include "GmshMatrix.h" #include "fullMatrix.h"
#include "GmshMessage.h" #include "GmshMessage.h"
//#if defined(_MSC_VER) //#if defined(_MSC_VER)
...@@ -39,7 +39,7 @@ extern "C" { ...@@ -39,7 +39,7 @@ extern "C" {
} }
template<> 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 M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1(); int LDA = _r, LDB = b.size1(), LDC = c.size1();
...@@ -49,8 +49,8 @@ void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c ...@@ -49,8 +49,8 @@ void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c
} }
template<> template<>
void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<double> > &b, void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<double> > &b,
gmshMatrix<std::complex<double> > &c) fullMatrix<std::complex<double> > &c)
{ {
int M = c.size1(), N = c.size2(), K = _c; int M = c.size1(), N = c.size2(), K = _c;
int LDA = _r, LDB = b.size1(), LDC = c.size1(); int LDA = _r, LDB = b.size1(), LDC = c.size1();
...@@ -60,7 +60,7 @@ void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<doubl ...@@ -60,7 +60,7 @@ void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<doubl
} }
template<> 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) double alpha, double beta)
{ {
int M = size1(), N = size2(), K = a.size2(); int M = size1(), N = size2(), K = a.size2();
...@@ -70,8 +70,8 @@ void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b, ...@@ -70,8 +70,8 @@ void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b,
} }
template<> template<>
void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &a, void fullMatrix<std::complex<double> >::gemm(fullMatrix<std::complex<double> > &a,
gmshMatrix<std::complex<double> > &b, fullMatrix<std::complex<double> > &b,
std::complex<double> alpha, std::complex<double> alpha,
std::complex<double> beta) std::complex<double> beta)
{ {
...@@ -82,7 +82,7 @@ void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > & ...@@ -82,7 +82,7 @@ void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &
} }
template<> 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; int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
double alpha = 1., beta = 0.; double alpha = 1., beta = 0.;
...@@ -91,8 +91,8 @@ void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y ...@@ -91,8 +91,8 @@ void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y
} }
template<> template<>
void gmshMatrix<std::complex<double> >::mult(const gmshVector<std::complex<double> > &x, void fullMatrix<std::complex<double> >::mult(const fullVector<std::complex<double> > &x,
gmshVector<std::complex<double> > &y) fullVector<std::complex<double> > &y)
{ {
int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1; int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
std::complex<double> alpha = 1., beta = 0.; std::complex<double> alpha = 1., beta = 0.;
...@@ -123,7 +123,7 @@ extern "C" { ...@@ -123,7 +123,7 @@ extern "C" {
} }
template<> template<>
bool gmshMatrix<double>::invertInPlace() bool fullMatrix<double>::invertInPlace()
{ {
int N = size1(), nrhs = N, lda = N, ldb = N, info; int N = size1(), nrhs = N, lda = N, ldb = N, info;
int *ipiv = new int[N]; int *ipiv = new int[N];
...@@ -148,8 +148,8 @@ bool gmshMatrix<double>::invertInPlace() ...@@ -148,8 +148,8 @@ bool gmshMatrix<double>::invertInPlace()
template<> template<>
bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR, bool fullMatrix<double>::eig(fullMatrix<double> &VL, fullVector<double> &DR,
gmshVector<double> &DI, gmshMatrix<double> &VR, fullVector<double> &DI, fullMatrix<double> &VR,
bool sortRealPart) bool sortRealPart)
{ {
int N = size1(), info; int N = size1(), info;
...@@ -193,7 +193,7 @@ bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR, ...@@ -193,7 +193,7 @@ bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, gmshVector<double> &DR,
} }
template<> 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 N = size1(), nrhs = 1, lda = N, ldb = N, info;
int *ipiv = new int[N]; int *ipiv = new int[N];
...@@ -209,7 +209,7 @@ bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<doub ...@@ -209,7 +209,7 @@ bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<doub
} }
template<> template<>
bool gmshMatrix<double>::invert(gmshMatrix<double> &result) bool fullMatrix<double>::invert(fullMatrix<double> &result)
{ {
int M = size1(), N = size2(), lda = size1(), info; int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)]; int *ipiv = new int[std::min(M, N)];
...@@ -231,9 +231,9 @@ bool gmshMatrix<double>::invert(gmshMatrix<double> &result) ...@@ -231,9 +231,9 @@ bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
} }
template<> 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 M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)]; int *ipiv = new int[std::min(M, N)];
F77NAME(dgetrf)(&M, &N, tmp._data, &lda, ipiv, &info); F77NAME(dgetrf)(&M, &N, tmp._data, &lda, ipiv, &info);
...@@ -253,12 +253,12 @@ double gmshMatrix<double>::determinant() const ...@@ -253,12 +253,12 @@ double gmshMatrix<double>::determinant() const
} }
template<> 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 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)); 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, F77NAME(dgesvd)("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
VT._data, &LDVT, WORK._data, &LWORK, &info); VT._data, &LDVT, WORK._data, &LWORK, &info);
V = VT.transpose(); V = VT.transpose();
......
Numeric/GmshMatrix.h → Numeric/fullMatrix.h
+ 31
− 31
View file @ 0801d0f9
...@@ -11,27 +11,27 @@ ...@@ -11,27 +11,27 @@
#include "GmshConfig.h" #include "GmshConfig.h"
#include "GmshMessage.h" #include "GmshMessage.h"
template <class scalar> class gmshMatrix; template <class scalar> class fullMatrix;
template <class scalar> template <class scalar>
class gmshVector class fullVector
{ {
private: private:
int _r; int _r;
scalar *_data; scalar *_data;
friend class gmshMatrix<scalar>; friend class fullMatrix<scalar>;
public: public:
gmshVector(int r) : _r(r) fullVector(int r) : _r(r)
{ {
_data = new scalar[_r]; _data = new scalar[_r];
scale(0.); scale(0.);
} }
gmshVector(const gmshVector<scalar> &other) : _r(other._r) fullVector(const fullVector<scalar> &other) : _r(other._r)
{ {
_data = new scalar[_r]; _data = new scalar[_r];
for(int i = 0; i < _r; ++i) _data[i] = other._data[i]; 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 inline scalar operator () (int i) const
{ {
return _data[i]; return _data[i];
...@@ -55,7 +55,7 @@ class gmshVector ...@@ -55,7 +55,7 @@ class gmshVector
for(int i = 0; i < _r; ++i) _data[i] *= s; 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; scalar s = 0.0;
for(int i = 0; i < _r; ++i) s += _data[i]*other._data[i]; for(int i = 0; i < _r; ++i) s += _data[i]*other._data[i];
...@@ -73,27 +73,27 @@ class gmshVector ...@@ -73,27 +73,27 @@ class gmshVector
}; };
template <class scalar> template <class scalar>
class gmshMatrix class fullMatrix
{ {
private: private:
int _r, _c; int _r, _c;
scalar *_data; scalar *_data;
public: public:
gmshMatrix(int r, int c) : _r(r), _c(c) fullMatrix(int r, int c) : _r(r), _c(c)
{ {
_data = new scalar[_r * _c]; _data = new scalar[_r * _c];
scale(0.); 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]; _data = new scalar[_r * _c];
gM_memcpy(other); gM_memcpy(other);
} }
gmshMatrix() : _r(0), _c(0), _data(0) {} fullMatrix() : _r(0), _c(0), _data(0) {}
~gmshMatrix() { if(_data) delete [] _data; } ~fullMatrix() { if(_data) delete [] _data; }
inline int size1() const { return _r; } inline int size1() const { return _r; }
inline int size2() const { return _c; } inline int size2() const { return _c; }
gmshMatrix<scalar> & operator = (const gmshMatrix<scalar> &other) fullMatrix<scalar> & operator = (const fullMatrix<scalar> &other)
{ {
if(this != &other){ if(this != &other){
_r = other._r; _r = other._r;
...@@ -103,7 +103,7 @@ class gmshMatrix ...@@ -103,7 +103,7 @@ class gmshMatrix
} }
return *this; 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]; for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
} }
...@@ -115,14 +115,14 @@ class gmshMatrix ...@@ -115,14 +115,14 @@ class gmshMatrix
{ {
return _data[i + _r * j]; 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) int desti0, int destj0)
{ {
for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++) for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++) for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
(*this)(desti, destj) = a(i, j); (*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) #if !defined(HAVE_BLAS)
{ {
c.scale(0.); c.scale(0.);
...@@ -133,11 +133,11 @@ class gmshMatrix ...@@ -133,11 +133,11 @@ class gmshMatrix
} }
#endif #endif
; ;
void gemm(gmshMatrix<scalar> &a, gmshMatrix<scalar> &b, void gemm(fullMatrix<scalar> &a, fullMatrix<scalar> &b,
scalar alpha=1., scalar beta=1.) scalar alpha=1., scalar beta=1.)
#if !defined(HAVE_BLAS) #if !defined(HAVE_BLAS)
{ {
gmshMatrix<scalar> temp(a.size1(), b.size2()); fullMatrix<scalar> temp(a.size1(), b.size2());
a.mult(b, temp); a.mult(b, temp);
temp.scale(alpha); temp.scale(alpha);
scale(beta); scale(beta);
...@@ -160,13 +160,13 @@ class gmshMatrix ...@@ -160,13 +160,13 @@ class gmshMatrix
{ {
for(int i = 0; i < _r * _c; ++i) _data[i] += a; 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 i = 0; i < size1(); i++)
for(int j = 0; j < size2(); j++) for(int j = 0; j < size2(); j++)
(*this)(i, j) += m(i, 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) #if !defined(HAVE_BLAS)
{ {
y.scale(0.); y.scale(0.);
...@@ -176,9 +176,9 @@ class gmshMatrix ...@@ -176,9 +176,9 @@ class gmshMatrix
} }
#endif #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 i = 0; i < size1(); i++)
for(int j = 0; j < size2(); j++) for(int j = 0; j < size2(); j++)
T(j, i) = (*this)(i, j); T(j, i) = (*this)(i, j);
...@@ -197,7 +197,7 @@ class gmshMatrix ...@@ -197,7 +197,7 @@ class gmshMatrix
_data[j + _r * i] = t; _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) #if !defined(HAVE_LAPACK)
{ {
Msg::Error("LU factorization requires LAPACK"); Msg::Error("LU factorization requires LAPACK");
...@@ -213,10 +213,10 @@ class gmshMatrix ...@@ -213,10 +213,10 @@ class gmshMatrix
} }
#endif #endif
; ;
bool eig(gmshMatrix<scalar> &VL, // left eigenvectors bool eig(fullMatrix<scalar> &VL, // left eigenvectors
gmshVector<double> &DR, // Real part of eigenvalues fullVector<double> &DR, // Real part of eigenvalues
gmshVector<double> &DI, // Im part of eigen fullVector<double> &DI, // Im part of eigen
gmshMatrix<scalar> &VR, fullMatrix<scalar> &VR,
bool sortRealPart=false) // if true: sorted from min 'DR' to max 'DR' bool sortRealPart=false) // if true: sorted from min 'DR' to max 'DR'
#if !defined(HAVE_LAPACK) #if !defined(HAVE_LAPACK)
{ {
...@@ -225,7 +225,7 @@ class gmshMatrix ...@@ -225,7 +225,7 @@ class gmshMatrix
} }
#endif #endif
; ;
bool invert(gmshMatrix<scalar> &result) bool invert(fullMatrix<scalar> &result)
#if !defined(HAVE_LAPACK) #if !defined(HAVE_LAPACK)
{ {
Msg::Error("LU factorization requires LAPACK"); Msg::Error("LU factorization requires LAPACK");
...@@ -233,11 +233,11 @@ class gmshMatrix ...@@ -233,11 +233,11 @@ class gmshMatrix
} }
#endif #endif
; ;
gmshMatrix<scalar> cofactor(int i, int j) const fullMatrix<scalar> cofactor(int i, int j) const
{ {
int ni = size1(); int ni = size1();
int nj = size2(); 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 I = 0; I < ni; I++){
for(int J = 0; J < nj; J++){ for(int J = 0; J < nj; J++){
if(J != j && I != i) if(J != j && I != i)
...@@ -254,7 +254,7 @@ class gmshMatrix ...@@ -254,7 +254,7 @@ class gmshMatrix
} }
#endif #endif
; ;
bool svd(gmshMatrix<scalar> &V, gmshVector<scalar> &S) bool svd(fullMatrix<scalar> &V, fullVector<scalar> &S)
#if !defined(HAVE_LAPACK) #if !defined(HAVE_LAPACK)
{ {
Msg::Error("Singular value decomposition requires LAPACK"); Msg::Error("Singular value decomposition requires LAPACK");
......
Numeric/functionSpace.cpp 0 → 100644
+ 672
− 0
View file @ 0801d0f9
// 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];
}
Numeric/functionSpace.h 0 → 100644
+ 107
− 0
View file @ 0801d0f9
// 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
Post/PViewDataGModel.cpp
+ 1
− 1
View file @ 0801d0f9
...@@ -93,7 +93,7 @@ bool PViewDataGModel::finalize() ...@@ -93,7 +93,7 @@ bool PViewDataGModel::finalize()
void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e) void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
{ {
int type = e->getType(); int type = e->getType();
const gmshFunctionSpace *fs = e->getFunctionSpace(); const functionSpace *fs = e->getFunctionSpace();
if(fs){ if(fs){
if(e->getPolynomialOrder() > 1) if(e->getPolynomialOrder() > 1)
setInterpolationMatrices(type, fs->coefficients, fs->monomials, setInterpolationMatrices(type, fs->coefficients, fs->monomials,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment