From 9ecb53d51f7d1dfb05c1d2486de5354c2e3b3e29 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Thu, 15 Nov 2012 18:39:51 +0000
Subject: [PATCH] Fixed bug (wrong basis type) in OptHOM. Fixed bug in
computation of Bergot basis gradient, still wrong.
---
Numeric/BergotBasis.cpp | 14 +++++++-------
Numeric/polynomialBasis.cpp | 2 --
Numeric/pyramidalBasis.cpp | 12 ++++++------
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 4 ++--
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 2 +-
5 files changed, 16 insertions(+), 18 deletions(-)
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index 2638b5bd13..0c018405e6 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -27,7 +27,7 @@ BergotBasis::BergotBasis(int p): order(p)
for (int i=0;i<=order;i++) {
for (int j=0;j<=order;j++) {
int mIJ = std::max(i,j);
- for (int k=0;k<= (int) (order - mIJ);k++,index++) {
+ for (int k=0; k<=order-mIJ; k++,index++) {
iOrder[index] = i;
jOrder[index] = j;
@@ -56,13 +56,13 @@ BergotBasis::~BergotBasis()
void BergotBasis::f(double u, double v, double w, double* val) const
{
- double uhat = (w == 1.) ? 1 : u/(1.-w);
+ double uhat = (w == 1.) ? 1. : u/(1.-w);
std::vector<double> uFcts(order+1);
//double uFcts[order+1];
legendre.f(uhat,&(uFcts[0]));
- double vhat = (w == 1.) ? 1 : v/(1.-w);
+ double vhat = (w == 1.) ? 1. : v/(1.-w);
std::vector<double> vFcts(order+1);
legendre.f(vhat,&(vFcts[0]));
@@ -83,9 +83,9 @@ void BergotBasis::f(double u, double v, double w, double* val) const
for (int i=0;i<=order;i++) {
for (int j=0;j<=order;j++) {
int mIJ = std::max(i,j);
- double fact = pow_int(1-w,(int) mIJ);
+ double fact = pow(1-w, mIJ);
double* wf = (wFcts.find(mIJ))->second;
- for (int k=0;k<=(int) (order - mIJ);k++,index++) {
+ for (int k=0; k<=order-mIJ; k++,index++) {
val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
}
}
@@ -140,11 +140,11 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const
double* wf = (wFcts .find(mIJ))->second;
double* wg = (wGrads.find(mIJ))->second;
- double wPowM2 = pow_int(1.-w,((int) mIJ) - 2);
+ double wPowM2 = pow(1.-w, mIJ-2);
double wPowM1 = wPowM2*(1.-w);
double wPowM0 = wPowM1*(1.-w);
- for (int k=0;k<= (int) (order - mIJ);k++,index++) {
+ for (int k=0; k<=order-mIJ; k++,index++) {
grads[index][0] = uGrads[i] * vFcts[j] * wf[k] * wPowM1;
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index c32c1adec6..bda418798d 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -422,7 +422,6 @@ void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) con
{
double dfv[1256][3];
dfm.resize (coefficients.size1(), coord.size1() * 3, false);
- int ii = 0;
int dimCoord = coord.size2();
for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
df(coord(iPoint,0), dimCoord > 1 ? coord(iPoint, 1) : 0., dimCoord > 2 ? coord(iPoint, 2) : 0., dfv);
@@ -430,7 +429,6 @@ void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) con
dfm(i, iPoint * 3 + 0) = dfv[i][0];
dfm(i, iPoint * 3 + 1) = dfv[i][1];
dfm(i, iPoint * 3 + 2) = dfv[i][2];
- ++ii;
}
}
}
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 24a382d29a..84f12f22e2 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -43,7 +43,7 @@ pyramidalBasis::~pyramidalBasis()
-inline void pyramidalBasis::f(double u, double v, double w, double *val) const
+void pyramidalBasis::f(double u, double v, double w, double *val) const
{
const int N = bergot->size();
@@ -60,7 +60,7 @@ inline void pyramidalBasis::f(double u, double v, double w, double *val) const
-inline void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
+void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
{
const int N = bergot->size(), NPts = coord.size1();
@@ -81,7 +81,7 @@ inline void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
-inline void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
+void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
{
const int N = bergot->size();
@@ -102,18 +102,18 @@ inline void pyramidalBasis::df(double u, double v, double w, double grads[][3])
-inline void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
+void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
{
const int N = bergot->size(), NPts = coord.size1();
double dfv[N][3];
- dfm.resize (NPts, 3*N, false);
+ dfm.resize (N, 3*NPts, false);
for (int iPt=0; iPt<NPts; iPt++) {
df(coord(iPt,0), coord(iPt,1), coord(iPt,2), dfv);
for (int i = 0; i < N; i++) {
- dfm(i, 3*iPt+0) = dfv[i][0];
+ dfm(i, 3*iPt) = dfv[i][0];
dfm(i, 3*iPt+1) = dfv[i][1];
dfm(i, 3*iPt+2) = dfv[i][2];
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index a37d265549..a74db6c0ea 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -13,7 +13,7 @@ std::map<int, fullMatrix<double> > Mesh::_lag2Bez;
-fullMatrix<double> Mesh::computeGSF(const polynomialBasis *lagrange, const bezierBasis *bezier)
+fullMatrix<double> Mesh::computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier)
{
// bezier points are defined in the [0,1] x [0,1] quad
fullMatrix<double> bezierPoints = bezier->points;
@@ -70,7 +70,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
for(std::set<MElement*>::const_iterator it = els.begin(); it != els.end(); ++it, ++iEl) {
MElement *el = *it;
_el[iEl] = el;
- const polynomialBasis *lagrange = (polynomialBasis*) el->getFunctionSpace();
+ const nodalBasis *lagrange = el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_lag2Bez.find(lagrange->type) == _lag2Bez.end()) {
_gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 61c89a72ca..87bdd18ffd 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -85,7 +85,7 @@ private:
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
SVector3 getNormalEl(int iEl);
- static fullMatrix<double> computeGSF(const polynomialBasis *lagrange, const bezierBasis *bezier);
+ static fullMatrix<double> computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier);
static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }
--
GitLab