Skip to content
Snippets Groups Projects
Commit c1f05e38 authored by Amaury Johnen's avatar Amaury Johnen
Browse files

revert monomials to be fullMatrix<double> + Add pow_int(double, double)

parent 820e1f9e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
FunctionSpace/BasisLagrange.h
+ 2
− 2
View file @ c1f05e38
...@@ -58,7 +58,7 @@ class BasisLagrange: public BasisLocal{ ...@@ -58,7 +58,7 @@ class BasisLagrange: public BasisLocal{
getPreEvaluatedDerivatives(unsigned int orientation) const; getPreEvaluatedDerivatives(unsigned int orientation) const;
const fullMatrix<double>& getCoefficient(void) const; const fullMatrix<double>& getCoefficient(void) const;
const fullMatrix<int>& getMonomial(void) const; const fullMatrix<double>& getMonomial(void) const;
std::vector<double> std::vector<double>
project(const MElement& element, project(const MElement& element,
...@@ -126,7 +126,7 @@ getCoefficient(void) const{ ...@@ -126,7 +126,7 @@ getCoefficient(void) const{
return lBasis->coefficients; return lBasis->coefficients;
} }
inline const fullMatrix<int>& BasisLagrange:: inline const fullMatrix<double>& BasisLagrange::
getMonomial(void) const{ getMonomial(void) const{
return lBasis->monomials; return lBasis->monomials;
} }
......
Numeric/nodalBasis.h
+ 0
− 2
View file @ c1f05e38
...@@ -14,8 +14,6 @@ class nodalBasis { ...@@ -14,8 +14,6 @@ class nodalBasis {
int type, parentType, order, dimension, numFaces; int type, parentType, order, dimension, numFaces;
bool serendip; bool serendip;
fullMatrix<double> points; fullMatrix<double> points;
fullMatrix<int> monomials;
fullMatrix<double> coefficients;
nodalBasis(int tag); nodalBasis(int tag);
virtual ~nodalBasis() {} virtual ~nodalBasis() {}
......
Numeric/polynomialBasis.cpp
+ 94
− 87
View file @ c1f05e38
...@@ -46,7 +46,7 @@ static void printClosure(polynomialBasis::clCont &fullClosure, ...@@ -46,7 +46,7 @@ static void printClosure(polynomialBasis::clCont &fullClosure,
static fullMatrix<double> generateLagrangeMonomialCoefficients static fullMatrix<double> generateLagrangeMonomialCoefficients
(const fullMatrix<int>& monomial, const fullMatrix<double>& point) (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{ {
if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){ if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d", Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
...@@ -75,34 +75,41 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients ...@@ -75,34 +75,41 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag) polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
{ {
fullMatrix<int> monom;
switch (parentType) { switch (parentType) {
case TYPE_PNT : case TYPE_PNT :
monomials = gmshGenerateMonomialsLine(0); monom = gmshGenerateMonomialsLine(0);
break; break;
case TYPE_LIN : case TYPE_LIN :
monomials = gmshGenerateMonomialsLine(order); monom = gmshGenerateMonomialsLine(order);
break; break;
case TYPE_TRI : case TYPE_TRI :
monomials = gmshGenerateMonomialsTriangle(order, serendip); monom = gmshGenerateMonomialsTriangle(order, serendip);
break; break;
case TYPE_QUA : case TYPE_QUA :
monomials = serendip ? gmshGenerateMonomialsQuadSerendipity(order) : monom = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
gmshGenerateMonomialsQuadrangle(order); gmshGenerateMonomialsQuadrangle(order);
break; break;
case TYPE_TET : case TYPE_TET :
monomials = gmshGenerateMonomialsTetrahedron(order, serendip); monom = gmshGenerateMonomialsTetrahedron(order, serendip);
break; break;
case TYPE_PRI : case TYPE_PRI :
monomials = serendip ? gmshGenerateMonomialsPrismSerendipity(order) : monom = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
gmshGenerateMonomialsPrism(order); gmshGenerateMonomialsPrism(order);
break; break;
case TYPE_HEX : case TYPE_HEX :
monomials = serendip ? gmshGenerateMonomialsHexaSerendipity(order) : monom = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
gmshGenerateMonomialsHexahedron(order); gmshGenerateMonomialsHexahedron(order);
break; break;
} }
monomials.resize(monom.size1(), monom.size2());
for (int i = 0; i < monom.size1(); ++i) {
for (int j = 0; j < monom.size2(); ++j) {
monomials(i, j) = static_cast<double>(monom(i, j));
}
}
coefficients = generateLagrangeMonomialCoefficients(monomials, points); coefficients = generateLagrangeMonomialCoefficients(monomials, points);
} }
...@@ -187,7 +194,7 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const ...@@ -187,7 +194,7 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0) if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) * grads[i][0] += coefficients(i, j) *
pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0); pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0);
} }
} }
break; break;
...@@ -199,12 +206,12 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const ...@@ -199,12 +206,12 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0) if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) * grads[i][0] += coefficients(i, j) *
pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) * pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0) *
pow_int(v, (int)monomials(j, 1)); pow_int(v, monomials(j, 1));
if (monomials(j, 1) > 0) if (monomials(j, 1) > 0)
grads[i][1] += coefficients(i, j) * grads[i][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1); pow_int(v, (monomials(j, 1) - 1)) * monomials(j, 1);
} }
} }
break; break;
...@@ -216,19 +223,19 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const ...@@ -216,19 +223,19 @@ void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0) if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) * grads[i][0] += coefficients(i, j) *
pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) * pow_int(u, (monomials(j, 0) - 1)) * monomials(j, 0) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if (monomials(j, 1) > 0) if (monomials(j, 1) > 0)
grads[i][1] += coefficients(i, j) * grads[i][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) * pow_int(v, (monomials(j, 1) - 1)) * monomials(j, 1) *
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if (monomials(j, 2) > 0) if (monomials(j, 2) > 0)
grads[i][2] += coefficients(i, j) * grads[i][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2); pow_int(w, (monomials(j, 2) - 1)) * monomials(j, 2);
} }
} }
break; break;
...@@ -248,7 +255,7 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con ...@@ -248,7 +255,7 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1) // second derivative !=0 if (monomials(j, 0) > 1) // second derivative !=0
hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) * hess[i][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 2)) *
monomials(j, 0) * (monomials(j, 0) - 1); monomials(j, 0) * (monomials(j, 0) - 1);
} }
} }
...@@ -260,16 +267,16 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con ...@@ -260,16 +267,16 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0; hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1) // second derivative !=0 if (monomials(j, 0) > 1) // second derivative !=0
hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) * hess[i][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 2)) *
monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1)); monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, monomials(j, 1));
if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0)) if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
hess[i][0][1] += coefficients(i, j) * hess[i][0][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) * pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1); pow_int(v, monomials(j, 1) - 1) * monomials(j, 1);
if (monomials(j, 1) > 1) if (monomials(j, 1) > 1)
hess[i][1][1] += coefficients(i, j) * hess[i][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1); pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
} }
hess[i][1][0] = hess[i][0][1]; hess[i][1][0] = hess[i][0][1];
} }
...@@ -282,35 +289,35 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con ...@@ -282,35 +289,35 @@ void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) con
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1) if (monomials(j, 0) > 1)
hess[i][0][0] += coefficients(i, j) * hess[i][0][0] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) * pow_int(u, monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)) if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
hess[i][0][1] += coefficients(i, j) * hess[i][0][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) * pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) * pow_int(v, monomials(j, 1) - 1) * monomials(j, 1) *
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0)) if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
hess[i][0][2] += coefficients(i, j) * hess[i][0][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) * pow_int(u, monomials(j, 0) - 1) * monomials(j, 0) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2); pow_int(w, monomials(j, 2) - 1) * monomials(j, 2);
if (monomials(j, 1) > 1) if (monomials(j, 1) > 1)
hess[i][1][1] += coefficients(i, j) * hess[i][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) * pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0)) if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
hess[i][1][2] += coefficients(i, j) * hess[i][1][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) * pow_int(v, monomials(j, 1) - 1) * monomials(j, 1) *
pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2); pow_int(w, monomials(j, 2) - 1) * monomials(j, 2);
if (monomials(j, 2) > 1) if (monomials(j, 2) > 1)
hess[i][2][2] += coefficients(i, j) * hess[i][2][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1); pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
} }
hess[i][1][0] = hess[i][0][1]; hess[i][1][0] = hess[i][0][1];
hess[i][2][0] = hess[i][0][2]; hess[i][2][0] = hess[i][0][2];
...@@ -334,7 +341,7 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3] ...@@ -334,7 +341,7 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 2) // third derivative !=0 if (monomials(j, 0) > 2) // third derivative !=0
third[i][0][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 3)) * third[i][0][0][0] += coefficients(i, j) * pow_int(u, (monomials(j, 0) - 3)) *
monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2); monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2);
} }
} }
...@@ -348,20 +355,20 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3] ...@@ -348,20 +355,20 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 2) // second derivative !=0 if (monomials(j, 0) > 2) // second derivative !=0
third[i][0][0][0] += coefficients(i, j) * third[i][0][0][0] += coefficients(i, j) *
pow_int(u, (int)(monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) * pow_int(u, (monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) *
pow_int(v, (int)monomials(j, 1)); pow_int(v, monomials(j, 1));
if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0)) if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
third[i][0][0][1] += coefficients(i, j) * third[i][0][0][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) * pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1); pow_int(v, monomials(j, 1) - 1) * monomials(j, 1);
if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1)) if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
third[i][0][1][1] += coefficients(i, j) * third[i][0][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)* pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1); pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
if (monomials(j, 1) > 2) if (monomials(j, 1) > 2)
third[i][1][1][1] += coefficients(i, j) * third[i][1][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2); pow_int(v, monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2);
} }
third[i][0][1][0] = third[i][0][0][1]; third[i][0][1][0] = third[i][0][0][1];
third[i][1][0][0] = third[i][0][0][1]; third[i][1][0][0] = third[i][0][0][1];
...@@ -378,62 +385,62 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3] ...@@ -378,62 +385,62 @@ void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]
for(int j = 0; j < coefficients.size2(); j++){ for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 2) // second derivative !=0 if (monomials(j, 0) > 2) // second derivative !=0
third[i][0][0][0] += coefficients(i, j) * third[i][0][0][0] += coefficients(i, j) *
pow_int(u, (int)(monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) * pow_int(u, (monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) *
pow_int(v, (int)monomials(j, 1))* pow_int(v, monomials(j, 1))*
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0)) if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
third[i][0][0][1] += coefficients(i, j) * third[i][0][0][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) * pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1)* pow_int(v, monomials(j, 1) - 1) * monomials(j, 1)*
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 0) > 1) && (monomials(j, 2) > 0)) if ((monomials(j, 0) > 1) && (monomials(j, 2) > 0))
third[i][0][0][2] += coefficients(i, j) * third[i][0][0][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) * pow_int(u, monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
pow_int(v, (int)monomials(j, 1)) * pow_int(v, monomials(j, 1)) *
pow_int(w, (int)monomials(j, 2) - 1)* monomials(j, 2); pow_int(w, monomials(j, 2) - 1)* monomials(j, 2);
if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1)) if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
third[i][0][1][1] += coefficients(i, j) * third[i][0][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)* pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)* pow_int(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)*
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)&& (monomials(j, 2) > 0)) if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)&& (monomials(j, 2) > 0))
third[i][0][1][2] += coefficients(i, j) * third[i][0][1][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)* pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
pow_int(v, (int)monomials(j, 1) - 1) *monomials(j, 1) * pow_int(v, monomials(j, 1) - 1) *monomials(j, 1) *
pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2); pow_int(w, monomials(j, 2) - 1) *monomials(j, 2);
if ((monomials(j, 0) > 0) && (monomials(j, 2) > 1)) if ((monomials(j, 0) > 0) && (monomials(j, 2) > 1))
third[i][0][2][2] += coefficients(i, j) * third[i][0][2][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)* pow_int(u, monomials(j, 0) - 1) *monomials(j, 0)*
pow_int(v, (int)monomials(j, 1))* pow_int(v, monomials(j, 1))*
pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1); pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
if ((monomials(j, 1) > 2)) if ((monomials(j, 1) > 2))
third[i][1][1][1] += coefficients(i, j) * third[i][1][1][1] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)* pow_int(v, monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)*
pow_int(w, (int)monomials(j, 2)); pow_int(w, monomials(j, 2));
if ((monomials(j, 1) > 1) && (monomials(j, 2) > 0)) if ((monomials(j, 1) > 1) && (monomials(j, 2) > 0))
third[i][1][1][2] += coefficients(i, j) * third[i][1][1][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)* pow_int(v, monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)*
pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2) ; pow_int(w, monomials(j, 2) - 1) *monomials(j, 2) ;
if ((monomials(j, 1) > 0) && (monomials(j, 2) > 1)) if ((monomials(j, 1) > 0) && (monomials(j, 2) > 1))
third[i][1][2][2] += coefficients(i, j) * third[i][1][2][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1)-1) *monomials(j, 1)* pow_int(v, monomials(j, 1)-1) *monomials(j, 1)*
pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ; pow_int(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ;
if ((monomials(j, 2) > 2)) if ((monomials(j, 2) > 2))
third[i][2][2][2] += coefficients(i, j) * third[i][2][2][2] += coefficients(i, j) *
pow_int(u, (int)monomials(j, 0)) * pow_int(u, monomials(j, 0)) *
pow_int(v, (int)monomials(j, 1))* pow_int(v, monomials(j, 1))*
pow_int(w, (int)monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2); pow_int(w, monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2);
} }
......
Numeric/polynomialBasis.h
+ 10
− 1
View file @ c1f05e38
...@@ -60,15 +60,24 @@ inline double pow_int(const double &a, const int &n) ...@@ -60,15 +60,24 @@ inline double pow_int(const double &a, const int &n)
return a4*a4*a2; return a4*a4*a2;
} }
default : default :
return pow_int(a,n-1)*a; return pow_int(a, n-9) * pow_int(a, 9);
} }
} }
inline double pow_int(const double &a, const double &d)
{
// Round double !
int n = static_cast<int>(d + .5);
return pow_int(a, n);
}
class polynomialBasis : public nodalBasis class polynomialBasis : public nodalBasis
{ {
public: public:
// for now the only implemented polynomial basis are nodal poly // for now the only implemented polynomial basis are nodal poly
// basis, we use the type of the corresponding gmsh element as type // basis, we use the type of the corresponding gmsh element as type
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
polynomialBasis(int tag); polynomialBasis(int tag);
~polynomialBasis(); ~polynomialBasis();
......
Numeric/pyramidalBasis.h
+ 1
− 0
View file @ c1f05e38
...@@ -18,6 +18,7 @@ class pyramidalBasis: public nodalBasis ...@@ -18,6 +18,7 @@ class pyramidalBasis: public nodalBasis
private: private:
// Orthogonal basis for the pyramid // Orthogonal basis for the pyramid
BergotBasis *bergot; BergotBasis *bergot;
fullMatrix<double> coefficients;
public: public:
pyramidalBasis(int tag); pyramidalBasis(int tag);
......
Post/PViewData.cpp
+ 0
− 64
View file @ c1f05e38
...@@ -122,23 +122,6 @@ void PViewData::setInterpolationMatrices(int type, ...@@ -122,23 +122,6 @@ void PViewData::setInterpolationMatrices(int type,
_interpolation[type].push_back(new fullMatrix<double>(expVal)); _interpolation[type].push_back(new fullMatrix<double>(expVal));
} }
void PViewData::setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<int> &expVal)
{
if(!type || _interpolation[type].size()) return;
fullMatrix<double> *dExpVal;
dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
for (int i = 0; i < expVal.size1(); ++i) {
for (int j = 0; i < expVal.size2(); ++j) {
(*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
}
}
_interpolation[type].push_back(new fullMatrix<double>(coefVal));
_interpolation[type].push_back(dExpVal);
}
void PViewData::setInterpolationMatrices(int type, void PViewData::setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal, const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal, const fullMatrix<double> &expVal,
...@@ -151,53 +134,6 @@ void PViewData::setInterpolationMatrices(int type, ...@@ -151,53 +134,6 @@ void PViewData::setInterpolationMatrices(int type,
_interpolation[type].push_back(new fullMatrix<double>(coefGeo)); _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
_interpolation[type].push_back(new fullMatrix<double>(expGeo)); _interpolation[type].push_back(new fullMatrix<double>(expGeo));
} }
void PViewData::setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<int> &expVal,
const fullMatrix<double> &coefGeo,
const fullMatrix<int> &expGeo)
{
if(!type || _interpolation[type].size()) return;
fullMatrix<double> *dExpVal, *dExpGeo;
dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
for (int i = 0; i < expVal.size1(); ++i) {
for (int j = 0; i < expVal.size2(); ++j) {
(*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
}
}
for (int i = 0; i < expGeo.size1(); ++i) {
for (int j = 0; i < expGeo.size2(); ++j) {
(*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
}
}
_interpolation[type].push_back(new fullMatrix<double>(coefVal));
_interpolation[type].push_back(dExpVal);
_interpolation[type].push_back(new fullMatrix<double>(coefGeo));
_interpolation[type].push_back(dExpGeo);
}
void PViewData::setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal,
const fullMatrix<double> &coefGeo,
const fullMatrix<int> &expGeo)
{
if(!type || _interpolation[type].size()) return;
fullMatrix<double> *dExpGeo;
dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
for (int i = 0; i < expGeo.size1(); ++i) {
for (int j = 0; i < expGeo.size2(); ++j) {
(*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
}
}
_interpolation[type].push_back(new fullMatrix<double>(coefVal));
_interpolation[type].push_back(new fullMatrix<double>(expVal));
_interpolation[type].push_back(new fullMatrix<double>(coefGeo));
_interpolation[type].push_back(dExpGeo);
}
int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p) int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p)
{ {
......
Post/PViewData.h
+ 0
− 13
View file @ c1f05e38
...@@ -205,24 +205,11 @@ class PViewData { ...@@ -205,24 +205,11 @@ class PViewData {
void setInterpolationMatrices(int type, void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal, const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal); const fullMatrix<double> &expVal);
void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<int> &expVal);
void setInterpolationMatrices(int type, void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal, const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal, const fullMatrix<double> &expVal,
const fullMatrix<double> &coefGeo, const fullMatrix<double> &coefGeo,
const fullMatrix<double> &expGeo); const fullMatrix<double> &expGeo);
void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<int> &expVal,
const fullMatrix<double> &coefGeo,
const fullMatrix<int> &expGeo);
void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal,
const fullMatrix<double> &coefGeo,
const fullMatrix<int> &expGeo);
int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p); int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
bool haveInterpolationMatrices(int type=0); bool haveInterpolationMatrices(int type=0);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment