diff --git a/COPYING b/COPYING
index b076a2ed1c7e9e86673290b6b85921cf64073621..07a2ee7721774b6751250061c1980ec9d615a09a 100644
--- a/COPYING
+++ b/COPYING
@@ -13,6 +13,7 @@ Copyright holders:
2010--2012 Christoph Grüninger
2006 Bernhard Haasdonk
2012 Olaf Ippisch
+2012 Arne Morten Kvarving
2009 Leonard Kern
2005--2007 Sreejith Pulloor Kuttanikkad
2003--2012 Robert Klöfkorn
diff --git a/dune/common/Makefile.am b/dune/common/Makefile.am
index 4cd44bc9916e804d3c637be4fee6314e4695a39d..3a71534cec123dd9b7d0edbd602c4eadf0298e85 100644
--- a/dune/common/Makefile.am
+++ b/dune/common/Makefile.am
@@ -8,6 +8,7 @@ noinst_LTLIBRARIES = libcommon.la
libcommon_la_SOURCES = \
debugallocator.cc \
fmatrixev.cc \
+ dynmatrixev.cc \
ios_state.cc \
parametertree.cc \
parametertreeparser.cc \
diff --git a/dune/common/dynmatrixev.cc b/dune/common/dynmatrixev.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1a6eed61ba765cf9548066b455dafbcc402dd203
--- /dev/null
+++ b/dune/common/dynmatrixev.cc
@@ -0,0 +1,154 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC
+#define DUNE_FMATRIXEIGENVALUES_EXT_CC
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#if HAVE_LAPACK
+
+// nonsymetric matrices
+#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
+
+// dsyev declaration (in liblapack)
+extern "C" {
+
+ /*
+ *
+ ** purpose
+ ** =======
+ **
+ ** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+ ** eigenvalues and, optionally, the left and/or right eigenvectors.
+ **
+ ** The right eigenvector v(j) of A satisfies
+ ** A * v(j) = lambda(j) * v(j)
+ ** where lambda(j) is its eigenvalue.
+ ** The left eigenvector u(j) of A satisfies
+ ** u(j)**T * A = lambda(j) * u(j)**T
+ ** where u(j)**T denotes the transpose of u(j).
+ **
+ ** The computed eigenvectors are normalized to have Euclidean norm
+ ** equal to 1 and largest component real.
+ **
+ ** arguments
+ ** =========
+ **
+ ** jobvl (input) char
+ ** = 'n': left eigenvectors of a are not computed;
+ ** = 'v': left eigenvectors of a are computed.
+ **
+ ** jobvr (input) char
+ ** = 'n': right eigenvectors of a are not computed;
+ ** = 'v': right eigenvectors of a are computed.
+ **
+ ** n (input) long int
+ ** the order of the matrix v. v >= 0.
+ **
+ ** a (input/output) BASE DATA TYPE array, dimension (lda,n)
+ ** on entry, the n-by-n matrix a.
+ ** on exit, a has been overwritten.
+ **
+ ** lda (input) long int
+ ** the leading dimension of the array a. lda >= max(1,n).
+ **
+ ** wr (output) BASE DATA TYPE array, dimension (n)
+ ** wi (output) BASE DATA TYPE array, dimension (n)
+ ** wr and wi contain the real and imaginary parts,
+ ** respectively, of the computed eigenvalues. complex
+ ** conjugate pairs of eigenvalues appear consecutively
+ ** with the eigenvalue having the positive imaginary part
+ ** first.
+ **
+ ** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+ ** if jobvl = 'v', the left eigenvectors u(j) are stored one
+ ** after another in the columns of vl, in the same order
+ ** as their eigenvalues.
+ ** if jobvl = 'n', vl is not referenced.
+ ** if the j-th eigenvalue is real, then u(j) = vl(:,j),
+ ** the j-th column of vl.
+ ** if the j-th and (j+1)-st eigenvalues form a complex
+ ** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+ ** u(j+1) = vl(:,j) - i*vl(:,j+1).
+ **
+ ** ldvl (input) long int
+ ** the leading dimension of the array vl. ldvl >= 1; if
+ ** jobvl = 'v', ldvl >= n.
+ **
+ ** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+ ** if jobvr = 'v', the right eigenvectors v(j) are stored one
+ ** after another in the columns of vr, in the same order
+ ** as their eigenvalues.
+ ** if jobvr = 'n', vr is not referenced.
+ ** if the j-th eigenvalue is real, then v(j) = vr(:,j),
+ ** the j-th column of vr.
+ ** if the j-th and (j+1)-st eigenvalues form a complex
+ ** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+ ** v(j+1) = vr(:,j) - i*vr(:,j+1).
+ **
+ ** ldvr (input) long int
+ ** the leading dimension of the array vr. ldvr >= 1; if
+ ** jobvr = 'v', ldvr >= n.
+ **
+ ** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+ ** on exit, if info = 0, work(1) returns the optimal lwork.
+ **
+ ** lwork (input) long int
+ ** the dimension of the array work. lwork >= max(1,3*n), and
+ ** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good
+ ** performance, lwork must generally be larger.
+ **
+ ** if lwork = -1, then a workspace query is assumed; the routine
+ ** only calculates the optimal size of the work array, returns
+ ** this value as the first entry of the work array, and no error
+ ** message related to lwork is issued by xerbla.
+ **
+ ** info (output) long int
+ ** = 0: successful exit
+ ** < 0: if info = -i, the i-th argument had an illegal value.
+ ** > 0: if info = i, the qr algorithm failed to compute all the
+ ** eigenvalues, and no eigenvectors have been computed;
+ ** elements i+1:n of wr and wi contain eigenvalues which
+ ** have converged.
+ **
+ **/
+
+ extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info);
+
+} // end extern C
+#endif
+
+namespace Dune {
+
+ namespace DynamicMatrixHelp {
+
+ void eigenValuesNonsymLapackCall(
+ const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info)
+ {
+#if HAVE_LAPACK
+ // call LAPACK dgeev
+ DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+ work, lwork, info);
+#else
+ DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+ }
+
+ } // end namespace FMatrixHelp
+
+} // end namespace Dune
+#endif
diff --git a/dune/common/dynmatrixev.hh b/dune/common/dynmatrixev.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5ad91d34712ce57e3dfd7741698cbd107ceda71c
--- /dev/null
+++ b/dune/common/dynmatrixev.hh
@@ -0,0 +1,70 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
+#define DUNE_DYNMATRIXEIGENVALUES_HH
+
+#include "dynmatrix.hh"
+
+namespace Dune {
+
+ namespace DynamicMatrixHelp {
+
+ // defined in fmatrixev_ext.cpp
+ extern void eigenValuesNonsymLapackCall(
+ const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info);
+
+ template <typename K, class C>
+ static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
+ DynamicVector<C>& eigenValues)
+ {
+ {
+ const long int N = matrix.rows();
+ const char jobvl = 'n';
+ const char jobvr = 'n';
+
+
+ // matrix to put into dgeev
+ double matrixVector[N * N];
+
+ // copy matrix
+ int row = 0;
+ for(int i=0; i<N; ++i)
+ {
+ for(int j=0; j<N; ++j, ++row)
+ {
+ matrixVector[ row ] = matrix[ i ][ j ];
+ }
+ }
+
+ // working memory
+ double eigenR[N];
+ double eigenI[N];
+ double work[3*N];
+
+ // return value information
+ long int info = 0;
+ long int lwork = 3*N;
+
+ // call LAPACK routine (see fmatrixev_ext.cc)
+ eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
+ &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+ &lwork, &info);
+
+ if( info != 0 )
+ {
+ std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+ DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+ }
+ for (int i=0; i<N; ++i)
+ eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
+ }
+ }
+
+ }
+
+}
+
+#endif
diff --git a/dune/common/fmatrixev.cc b/dune/common/fmatrixev.cc
index feba39854b3db2c43d4f909bd5bdc020f3499de1..d720d509caa8a6b43247e6bb30c7a47e2d8903a8 100644
--- a/dune/common/fmatrixev.cc
+++ b/dune/common/fmatrixev.cc
@@ -78,6 +78,111 @@ extern "C" {
int* n, double* a, const long int* lda, double* w,
double* work, const long int* lwork, long int* info);
+ /*
+ *
+ ** purpose
+ ** =======
+ **
+ ** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+ ** eigenvalues and, optionally, the left and/or right eigenvectors.
+ **
+ ** The right eigenvector v(j) of A satisfies
+ ** A * v(j) = lambda(j) * v(j)
+ ** where lambda(j) is its eigenvalue.
+ ** The left eigenvector u(j) of A satisfies
+ ** u(j)**T * A = lambda(j) * u(j)**T
+ ** where u(j)**T denotes the transpose of u(j).
+ **
+ ** The computed eigenvectors are normalized to have Euclidean norm
+ ** equal to 1 and largest component real.
+ **
+ ** arguments
+ ** =========
+ **
+ ** jobvl (input) char
+ ** = 'n': left eigenvectors of a are not computed;
+ ** = 'v': left eigenvectors of a are computed.
+ **
+ ** jobvr (input) char
+ ** = 'n': right eigenvectors of a are not computed;
+ ** = 'v': right eigenvectors of a are computed.
+ **
+ ** n (input) long int
+ ** the order of the matrix v. v >= 0.
+ **
+ ** a (input/output) BASE DATA TYPE array, dimension (lda,n)
+ ** on entry, the n-by-n matrix a.
+ ** on exit, a has been overwritten.
+ **
+ ** lda (input) long int
+ ** the leading dimension of the array a. lda >= max(1,n).
+ **
+ ** wr (output) BASE DATA TYPE array, dimension (n)
+ ** wi (output) BASE DATA TYPE array, dimension (n)
+ ** wr and wi contain the real and imaginary parts,
+ ** respectively, of the computed eigenvalues. complex
+ ** conjugate pairs of eigenvalues appear consecutively
+ ** with the eigenvalue having the positive imaginary part
+ ** first.
+ **
+ ** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+ ** if jobvl = 'v', the left eigenvectors u(j) are stored one
+ ** after another in the columns of vl, in the same order
+ ** as their eigenvalues.
+ ** if jobvl = 'n', vl is not referenced.
+ ** if the j-th eigenvalue is real, then u(j) = vl(:,j),
+ ** the j-th column of vl.
+ ** if the j-th and (j+1)-st eigenvalues form a complex
+ ** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+ ** u(j+1) = vl(:,j) - i*vl(:,j+1).
+ **
+ ** ldvl (input) long int
+ ** the leading dimension of the array vl. ldvl >= 1; if
+ ** jobvl = 'v', ldvl >= n.
+ **
+ ** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+ ** if jobvr = 'v', the right eigenvectors v(j) are stored one
+ ** after another in the columns of vr, in the same order
+ ** as their eigenvalues.
+ ** if jobvr = 'n', vr is not referenced.
+ ** if the j-th eigenvalue is real, then v(j) = vr(:,j),
+ ** the j-th column of vr.
+ ** if the j-th and (j+1)-st eigenvalues form a complex
+ ** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+ ** v(j+1) = vr(:,j) - i*vr(:,j+1).
+ **
+ ** ldvr (input) long int
+ ** the leading dimension of the array vr. ldvr >= 1; if
+ ** jobvr = 'v', ldvr >= n.
+ **
+ ** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+ ** on exit, if info = 0, work(1) returns the optimal lwork.
+ **
+ ** lwork (input) long int
+ ** the dimension of the array work. lwork >= max(1,3*n), and
+ ** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good
+ ** performance, lwork must generally be larger.
+ **
+ ** if lwork = -1, then a workspace query is assumed; the routine
+ ** only calculates the optimal size of the work array, returns
+ ** this value as the first entry of the work array, and no error
+ ** message related to lwork is issued by xerbla.
+ **
+ ** info (output) long int
+ ** = 0: successful exit
+ ** < 0: if info = -i, the i-th argument had an illegal value.
+ ** > 0: if info = i, the qr algorithm failed to compute all the
+ ** eigenvalues, and no eigenvectors have been computed;
+ ** elements i+1:n of wr and wi contain eigenvalues which
+ ** have converged.
+ **
+ **/
+
+ extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info);
+
} // end extern C
#endif
@@ -98,6 +203,21 @@ namespace Dune {
#endif
}
+ void eigenValuesNonsymLapackCall(
+ const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info)
+ {
+#if HAVE_LAPACK
+ // call LAPACK dgeev
+ DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+ work, lwork, info);
+#else
+ DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+ }
+
} // end namespace FMatrixHelp
} // end namespace Dune
diff --git a/dune/common/fmatrixev.hh b/dune/common/fmatrixev.hh
index 2f244eae4980bd3184f10af89d1b5f05cbc7c223..cee1aa190ae14820d68b436d73a4d33a45a98f68 100644
--- a/dune/common/fmatrixev.hh
+++ b/dune/common/fmatrixev.hh
@@ -30,6 +30,12 @@ namespace Dune {
int* n, double* a, const long int* lda, double* w,
double* work, const long int* lwork, long int* info);
+ extern void eigenValuesNonsymLapackCall(
+ const char* jobvl, const char* jobvr, const long
+ int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+ const long int* ldvl, double* vr, const long int* ldvr, double* work,
+ const long int* lwork, const long int* info);
+
/** \brief calculates the eigenvalues of a symetric field matrix
\param[in] matrix matrix eigenvalues are calculated for
\param[out] eigenvalues FieldVector that contains eigenvalues in
@@ -122,6 +128,58 @@ namespace Dune {
}
}
+ template <int dim, typename K, class C>
+ static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
+ FieldVector<C, dim>& eigenValues)
+ {
+ {
+ const long int N = dim ;
+ const char jobvl = 'n';
+ const char jobvr = 'n';
+
+ // length of matrix vector
+ const long int w = N * N ;
+
+ // matrix to put into dgeev
+ double matrixVector[dim * dim];
+
+ // copy matrix
+ int row = 0;
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j, ++row)
+ {
+ matrixVector[ row ] = matrix[ i ][ j ];
+ }
+ }
+
+ // working memory
+ double eigenR[dim];
+ double eigenI[dim];
+ double work[3*dim];
+
+ // return value information
+ long int info = 0;
+ long int lwork = 3*dim;
+
+ // call LAPACK routine (see fmatrixev_ext.cc)
+ eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
+ &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+ &lwork, &info);
+
+ if( info != 0 )
+ {
+ std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+ DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+ }
+ for (int i=0; i<N; ++i) {
+ eigenValues[i].real = eigenR[i];
+ eigenValues[i].imag = eigenI[i];
+ }
+ }
+
+ }
+
} // end namespace FMatrixHelp
/** @} end documentation */
diff --git a/dune/common/test/Makefile.am b/dune/common/test/Makefile.am
index 9b676024c87b486f8b0091777c2a88d3144c34d1..f6a7e38f1d1f73fb43f8062616996e75a5c69d03 100644
--- a/dune/common/test/Makefile.am
+++ b/dune/common/test/Makefile.am
@@ -11,6 +11,7 @@ TESTPROGS = \
diagonalmatrixtest \
dynmatrixtest \
dynvectortest \
+ eigenvaluestest \
enumsettest \
fassigntest \
fmatrixtest \
@@ -161,6 +162,9 @@ dynmatrixtest_SOURCES = dynmatrixtest.cc
dynvectortest_SOURCES = dynvectortest.cc
+eigenvaluestest_SOURCES = eigenvaluestest.cc
+eigenvaluestest_LDADD = $(LAPACK_LIBS) $(LDADD) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+
fmatrixtest_SOURCES = fmatrixtest.cc
fmatrixtest_LDADD = $(LAPACK_LIBS) $(LDADD) $(BLAS_LIBS) $(LIBS) $(FLIBS)
diff --git a/dune/common/test/eigenvaluestest.cc b/dune/common/test/eigenvaluestest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..dea8906f6395eb8a0d695c5e195818d0bb96066e
--- /dev/null
+++ b/dune/common/test/eigenvaluestest.cc
@@ -0,0 +1,259 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+//==============================================================================
+//!
+//! \file shapefunctions.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Classes for shape functions. Loosely based on code in dune-grid-howto
+//!
+//==============================================================================
+
+#include <dune/common/fvector.hh>
+#include <dune/common/dynmatrixev.hh>
+
+#include <complex>
+
+
+//! \brief Represents a cardinal function on a line
+template<class ctype, class rtype>
+class LagrangeCardinalFunction
+{
+public:
+ //! \brief Empty default constructor
+ LagrangeCardinalFunction() {}
+
+ //! \brief Construct a cardinal function with the given nodes
+ //! \param[in] nodes_ The nodes
+ //! \param[in] i The node this function is associated with
+ LagrangeCardinalFunction(const std::vector<rtype>& nodes_,
+ size_t i)
+ : nodes(nodes_), node(i) {}
+
+ //! \brief Evaluate the shape function
+ //! \param[in] local The local coordinates
+ rtype evaluateFunction(const ctype& local) const
+ {
+ rtype result = 1;
+ for (size_t i=0; i < nodes.size(); ++i) {
+ if (i != node)
+ result *= (local-nodes[i])/(nodes[node]-nodes[i]);
+ }
+
+ return result;
+ }
+
+ //! \brief Evaluate the derivative of the cardinal function
+ //! \param[in] local The local coordinates
+ rtype evaluateGradient(const ctype& local) const
+ {
+ rtype result = 0;
+ for (size_t i=0; i < nodes.size(); ++i) {
+ rtype f = 1;
+ for (int j=0; j < nodes.size(); ++j) {
+ if (i != j && j != node)
+ f *= (local-nodes[j])/(nodes[node]-nodes[j]);
+ }
+ result += f/(nodes[node]-nodes[i]);
+ }
+
+ return result;
+ }
+
+private:
+ //! \brief The nodes
+ std::vector<rtype> nodes;
+
+ size_t node;
+};
+
+//! \brief Represents a tensor-product of 1D functions
+template<class rtype, class ctype, class ftype, int dim>
+class TensorProductFunction
+{
+public:
+ //! \brief The dimension of the function
+ enum { dimension = dim };
+
+ //! \brief Empty default constructor
+ TensorProductFunction() {}
+
+ //! \brief Construct a tensor-product function
+ //! \param[in] funcs_ The functions
+ TensorProductFunction(const Dune::FieldVector<ftype, dim>& funcs_)
+ : funcs(funcs_) {}
+
+ //! \brief Evaluate the function
+ //! \param[in] local The local coordinates
+ rtype evaluateFunction(const Dune::FieldVector<ctype,dim>& local) const
+ {
+ rtype result = 1;
+ for (int i=0; i < dim; ++i)
+ result *= funcs[i].evaluateFunction(local[i]);
+
+ return result;
+ }
+
+ Dune::FieldVector<rtype, dim>
+ evaluateGradient(const Dune::FieldVector<ctype,dim>& local) const
+ {
+ Dune::FieldVector<rtype, dim> result;
+ for (int i=0; i < dim; ++i)
+ result[i] = funcs[i].evaluateGradient(local[i]);
+ }
+private:
+ Dune::FieldVector<ftype, dim> funcs;
+};
+
+template<int dim>
+class PNShapeFunctionSet
+{
+public:
+ typedef LagrangeCardinalFunction<double, double> CardinalFunction;
+
+ typedef TensorProductFunction<double, double, CardinalFunction, dim>
+ ShapeFunction;
+
+ PNShapeFunctionSet(int n1, int n2, int n3=0)
+ {
+ int dims[3] = {n1, n2, n3};
+ cfuncs.resize(dim);
+ for (int i=0; i < dim; ++i) {
+ std::vector<double> grid;
+ grid = gaussLobattoLegendreGrid(dims[i]);
+ for (int j=0; j<dims[i]; ++j)
+ cfuncs[i].push_back(CardinalFunction(grid,j));
+ }
+ int l=0;
+ Dune::FieldVector<CardinalFunction,dim> fs;
+ if (dim == 3) {
+ f.resize(n1*n2*n3);
+ for (int k=0; k < n3; ++k) {
+ for (int j=0; j < n2; ++j)
+ for (int i=0; i< n1; ++i) {
+ fs[0] = cfuncs[0][i];
+ fs[1] = cfuncs[1][j];
+ fs[2] = cfuncs[2][k];
+ f[l++] = ShapeFunction(fs);
+ }
+ }
+ } else {
+ f.resize(n1*n2);
+ for (int j=0; j < n2; ++j) {
+ for (int i=0; i< n1; ++i) {
+ fs[0] = cfuncs[0][i];
+ fs[1] = cfuncs[1][j];
+ f[l++] = ShapeFunction(fs);
+ }
+ }
+ }
+ }
+
+ //! \brief Obtain a given shape function
+ //! \param[in] i The requested shape function
+ const ShapeFunction& operator[](int i) const
+ {
+ return f[i];
+ }
+
+ int size()
+ {
+ return f.size();
+ }
+protected:
+ std::vector< std::vector<CardinalFunction> > cfuncs;
+ std::vector<ShapeFunction> f;
+
+ double legendre(double x, int n)
+ {
+ std::vector<double> Ln;
+ Ln.resize(n+1);
+ Ln[0] = 1.f;
+ Ln[1] = x;
+ if( n > 1 ) {
+ for( int i=1; i<n; i++ )
+ Ln[i+1] = (2*i+1.0)/(i+1.0)*x*Ln[i]-i/(i+1.0)*Ln[i-1];
+ }
+
+ return Ln[n];
+ }
+
+ double legendreDerivative(double x, int n)
+ {
+ std::vector<double> Ln;
+ Ln.resize(n+1);
+
+ Ln[0] = 1.0; Ln[1] = x;
+
+ if( (x == 1.0) || (x == -1.0) )
+ return( pow(x,n-1)*n*(n+1.0)/2.0 );
+ else {
+ for( int i=1; i<n; i++ )
+ Ln[i+1] = (2.0*i+1.0)/(i+1.0)*x*Ln[i]-(double)i/(i+1.0)*Ln[i-1];
+ return( (double)n/(1.0-x*x)*Ln[n-1]-n*x/(1-x*x)*Ln[n] );
+ }
+ }
+
+
+ std::vector<double> gaussLegendreGrid(int n)
+ {
+ Dune::DynamicMatrix<double> A(n,n,0.0);
+
+ A[0][1] = 1.f;
+ for (int i=1; i<n-1; ++i) {
+ A[i][i-1] = (double)i/(2.0*(i+1.0)-1.0);
+ A[i][i+1] = (double)(i+1.0)/(2*(i+1.0)-1.0);
+ }
+ A[n-1][n-2] = (n-1.0)/(2.0*n-1.0);
+
+ Dune::DynamicVector<std::complex<double> > eigenValues(n);
+ Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues);
+
+ std::vector<double> result(n);
+ for (int i=0; i < n; ++i)
+ result[i] = std::real(eigenValues[i]);
+ std::sort(result.begin(),result.begin()+n);
+
+ return result;
+ }
+
+ std::vector<double> gaussLobattoLegendreGrid(int n)
+ {
+ assert(n > 1);
+ const double tolerance = 1.e-15;
+
+ std::vector<double> result(n);
+ result[0] = 0.0;
+ result[n-1] = 1.0;
+ if (n == 3)
+ result[1] = 0.5;
+
+ if (n < 4)
+ return result;
+
+ std::vector<double> glgrid = gaussLegendreGrid(n-1);
+ for (int i=1; i<n-1; ++i) {
+ result[i] = (glgrid[i-1]+glgrid[i])/2.0;
+ double old = 0.0;
+ while (std::abs(old-result[i]) > tolerance) {
+ old = result[i];
+ double L = legendre(old, n-1);
+ double Ld = legendreDerivative(old, n-1);
+ result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L);
+ }
+ result[i] = (result[i]+1.0)/2.0;
+ }
+
+ return result;
+ }
+};
+
+int main()
+{
+ // Not really a test but better than nothing.
+ PNShapeFunctionSet<2> lbasis(2, 3);
+ return 0;
+}