12#ifndef SHAPEFUNCTIONS_HPP_
13#define SHAPEFUNCTIONS_HPP_
15#include <dune/common/fvector.hh>
16#include <dune/common/dynmatrixev.hh>
25 template<
class ctype,
class rtype,
int dim>
30 enum { dimension = dim };
39 const Dune::FieldVector<rtype,dim>& coeff1_)
40 : coeff0(coeff0_), coeff1(coeff1_) {}
45 void setCoeff(
const Dune::FieldVector<rtype,dim>& coeff0_,
const Dune::FieldVector<rtype,dim>& coeff1_)
56 for (
int i = 0; i < dim; ++i)
57 result *= (coeff0[i]+coeff1[i] * local[i]);
63 Dune::FieldVector<rtype,dim>
66 Dune::FieldVector<rtype, dim> result;
67 for (
int i=0;i<dim;++i) {
69 for (
int j=0;j<dim;++j) {
71 result[i] *= coeff1[j];
73 result[i] *= (coeff0[j]+coeff1[j]*local[j]);
82 Dune::FieldVector<rtype,dim> coeff0;
85 Dune::FieldVector<rtype,dim> coeff1;
89 template<
class ctype,
class rtype>
101 : nodes(nodes_), node(i) {}
108 for (
size_t i=0; i < nodes.size(); ++i) {
110 result *= (local-nodes[i])/(nodes[node]-nodes[i]);
121 for (
size_t i=0; i < nodes.size(); ++i) {
123 for (
int j=0; j < nodes.size(); ++j) {
124 if (i != j && j != node)
125 f *= (local-nodes[j])/(nodes[node]-nodes[j]);
127 result += f/(nodes[node]-nodes[i]);
135 std::vector<rtype> nodes;
141 template<
class rtype,
class ctype,
class ftype,
int dim>
146 enum { dimension = dim };
161 for (
int i=0; i < dim; ++i)
167 Dune::FieldVector<rtype, dim>
168 evaluateGradient(
const Dune::FieldVector<ctype,dim>& local)
const
170 Dune::FieldVector<rtype, dim> result;
171 for (
int i=0; i < dim; ++i)
172 result[i] = funcs[i].evaluateGradient(local[i]);
175 Dune::FieldVector<ftype, dim> funcs;
179 template<
class ctype,
class rtype,
int dim>
184 enum { n = dim==2?4:8 };
210 static rtype coeffs11[] = {0,
212 static rtype coeffs12[] = {1,
215 static rtype coeffs21[] = { 1, 1,
219 static rtype coeffs22[] = {-1,-1,
224 static rtype coeffs31[] = { 1, 1, 1,
232 static rtype coeffs32[] = {-1,-1,-1,
256 Dune::FieldVector<rtype,dim> c1;
257 Dune::FieldVector<rtype,dim> c2;
258 for (
int i=0;i<n;++i) {
259 for (
int j=0;j<dim;++j) {
260 c1[j] = coeffs1[i*dim+j];
261 c2[j] = coeffs2[i*dim+j];
263 f[i].setCoeff(c1,c2);
268 std::array<ShapeFunction, n> f;
282 const int dims[3] = {n1, n2, n3};
284 for (
int i=0; i < dim; ++i) {
285 std::vector<double> grid;
286 grid = gaussLobattoLegendreGrid(dims[i]);
287 for (
int j=0;j<dims[i];++j)
291 Dune::FieldVector<CardinalFunction,dim> fs;
294 for (
int k=0; k < n3; ++k) {
295 for (
int j=0; j < n2; ++j)
296 for (
int i=0; i< n1; ++i) {
297 fs[0] = cfuncs[0][i];
298 fs[1] = cfuncs[1][j];
299 fs[2] = cfuncs[2][k];
305 for (
int j=0; j < n2; ++j) {
306 for (
int i=0; i< n1; ++i) {
307 fs[0] = cfuncs[0][i];
308 fs[1] = cfuncs[1][j];
327 std::vector< std::vector<CardinalFunction> > cfuncs;
328 std::vector<ShapeFunction> f;
330 double legendre(
double x,
int n)
332 std::vector<double> Ln;
337 for(
int i=1;i<n;i++ )
338 Ln[i+1] = (2*i+1.0)/(i+1.0)*x*Ln[i]-i/(i+1.0)*Ln[i-1];
344 double legendreDerivative(
double x,
int n)
346 std::vector<double> Ln;
349 Ln[0] = 1.0; Ln[1] = x;
351 if( (x == 1.0) || (x == -1.0) )
352 return( pow(x,n-1)*n*(n+1.0)/2.0 );
354 for(
int i=1;i<n;i++ )
355 Ln[i+1] = (2.0*i+1.0)/(i+1.0)*x*Ln[i]-(
double)i/(i+1.0)*Ln[i-1];
356 return( (
double)n/(1.0-x*x)*Ln[n-1]-n*x/(1-x*x)*Ln[n] );
360 std::vector<double> gaussLegendreGrid(
int n)
362 Dune::DynamicMatrix<double> A(n,n,0.0);
365 for (
int i=1;i<n-1;++i) {
366 A[i][i-1] = (double)i/(2.0*(i+1.0)-1.0);
367 A[i][i+1] = (double)(i+1.0)/(2*(i+1.0)-1.0);
369 A[n-1][n-2] = (n-1.0)/(2.0*n-1.0);
371 Dune::DynamicVector<std::complex<double> > eigenValues(n);
372 Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues);
374 std::vector<double> result(n);
375 for (
int i=0; i < n; ++i)
376 result[i] = std::real(eigenValues[i]);
377 std::sort(result.begin(),result.begin()+n);
382 std::vector<double> gaussLobattoLegendreGrid(
int n)
385 const double tolerance = 1.e-15;
387 std::vector<double> result(n);
396 std::vector<double> glgrid = gaussLegendreGrid(n-1);
397 for (
int i=1;i<n-1;++i) {
398 result[i] = (glgrid[i-1]+glgrid[i])/2.0;
400 while (std::abs(old-result[i]) > tolerance) {
402 double L = legendre(old, n-1);
403 double Ld = legendreDerivative(old, n-1);
404 result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L);
406 result[i] = (result[i]+1.0)/2.0;
Represents a cardinal function on a line.
Definition shapefunctions.hpp:91
rtype evaluateGradient(const ctype &local) const
Evaluate the derivative of the cardinal function.
Definition shapefunctions.hpp:118
rtype evaluateFunction(const ctype &local) const
Evaluate the shape function.
Definition shapefunctions.hpp:105
LagrangeCardinalFunction()
Empty default constructor.
Definition shapefunctions.hpp:94
LagrangeCardinalFunction(const std::vector< rtype > &nodes_, size_t i)
Construct a cardinal function with the given nodes.
Definition shapefunctions.hpp:99
Represents a linear shape function on a Q4/Q8 element.
Definition shapefunctions.hpp:27
void setCoeff(const Dune::FieldVector< rtype, dim > &coeff0_, const Dune::FieldVector< rtype, dim > &coeff1_)
Set the given conefficients.
Definition shapefunctions.hpp:45
LinearShapeFunction()
Default constructor.
Definition shapefunctions.hpp:33
Dune::FieldVector< rtype, dim > evaluateGradient(const Dune::FieldVector< ctype, dim > &local) const
Evaluate the gradient of the shape function.
Definition shapefunctions.hpp:64
rtype evaluateFunction(const Dune::FieldVector< ctype, dim > &local) const
Evaluate the shape function.
Definition shapefunctions.hpp:53
LinearShapeFunction(const Dune::FieldVector< rtype, dim > &coeff0_, const Dune::FieldVector< rtype, dim > &coeff1_)
Construct a shape function with the given coefficients.
Definition shapefunctions.hpp:38
Singleton handler for the set of LinearShapeFunction.
Definition shapefunctions.hpp:181
rtype resulttype
The type of the return value from a shape function.
Definition shapefunctions.hpp:190
LinearShapeFunction< ctype, rtype, dim > ShapeFunction
A single shape function.
Definition shapefunctions.hpp:187
static const P1ShapeFunctionSet & instance()
Get the only instance of this class.
Definition shapefunctions.hpp:193
const ShapeFunction & operator[](int i) const
Obtain a given shape function.
Definition shapefunctions.hpp:201
Definition shapefunctions.hpp:273
const ShapeFunction & operator[](int i) const
Obtain a given shape function.
Definition shapefunctions.hpp:317
Represents a tensor-product of 1D functions.
Definition shapefunctions.hpp:143
TensorProductFunction()
Empty default constructor.
Definition shapefunctions.hpp:149
rtype evaluateFunction(const Dune::FieldVector< ctype, dim > &local) const
Evaluate the function.
Definition shapefunctions.hpp:158
TensorProductFunction(const Dune::FieldVector< ftype, dim > &funcs_)
Construct a tensor-product function.
Definition shapefunctions.hpp:153
Inverting small matrices.
Definition ImplicitAssembly.hpp:43