Select Git revision
functionSpace.h
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
functionSpace.h 3.29 KiB
// 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