Select Git revision
TetNodeBasis.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Nicolas Marsic authored-- Add Point Mapping from ABC space to XYZ space -- Handle Elements with no Face ** LineBasis: -- Works again ** Jacobian and GroupOfJacobian: -- Computes Good Jacobian from Basis ReferenceSpace ** Formulation: -- Addapted for good Jacobian ** TermProjection: -- Addapted to compute function to project at good point (see ReferenceSpace ABCtoXYZ mapping) ** Laplace, Poisson, SteadyWave: -- Passed with all basis ** WARNING: -- Tet Basis NOT Retested (but seems to work) -- FunctionSpace::getDof doesn't work if wrong element type (because of basis implication) -- SystemAbstract::fixCoef doesn't work anymore (see above) -- Problem with cureved Element (getTypeForMSH gives HO element for FunctionSpace::getDof)Nicolas Marsic authored-- Add Point Mapping from ABC space to XYZ space -- Handle Elements with no Face ** LineBasis: -- Works again ** Jacobian and GroupOfJacobian: -- Computes Good Jacobian from Basis ReferenceSpace ** Formulation: -- Addapted for good Jacobian ** TermProjection: -- Addapted to compute function to project at good point (see ReferenceSpace ABCtoXYZ mapping) ** Laplace, Poisson, SteadyWave: -- Passed with all basis ** WARNING: -- Tet Basis NOT Retested (but seems to work) -- FunctionSpace::getDof doesn't work if wrong element type (because of basis implication) -- SystemAbstract::fixCoef doesn't work anymore (see above) -- Problem with cureved Element (getTypeForMSH gives HO element for FunctionSpace::getDof)
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