From 9b9979f6134c3446243cb3cf7ac4fd7dbe970b40 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 24 May 2011 15:13:38 +0000
Subject: [PATCH] fix drawing and shape functions for order N
---
Geo/MHexahedron.cpp | 71 ++++++++++++++++-----------------------------
Geo/MHexahedron.h | 27 ++++++-----------
2 files changed, 34 insertions(+), 64 deletions(-)
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index c008a557db..38a392ac0f 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -102,8 +102,8 @@ void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3
double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
SPoint3 pnt1, pnt2;
- pnt(u1,v1,w1,pnt1);
- pnt(u2,v2,w2,pnt2);
+ pnt(u1, v1, w1, pnt1);
+ pnt(u2, v2, w2, pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -113,20 +113,19 @@ void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3
n[0] = n[1] = 1 ;
}
-int MHexahedronN::getNumEdgesRep(){
+int MHexahedronN::getNumEdgesRep()
+{
return 12 * CTX::instance()->mesh.numSubEdges;
}
-//int MHexaedronN::getNumFacesRep(){
-// return 8 * SQU(CTX::instance()->mesh.numSubEdges);
-//}
-
const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int nv = getNumVolumeVertices();
+ printf("nv = %d\n", nv);
+
if ((nv == 0) && (o == -1)) {
switch (order) {
case 0: return polynomialBases::find(MSH_HEX_1);
@@ -147,7 +146,7 @@ const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
case 0: return polynomialBases::find(MSH_HEX_1);
case 1: return polynomialBases::find(MSH_HEX_8);
case 2: return polynomialBases::find(MSH_HEX_27);
- case 3: return polynomialBases::find(MSH_HEX_64);
+ case 3: printf("BBBBBBBBBBB\n"); return polynomialBases::find(MSH_HEX_64);
case 4: return polynomialBases::find(MSH_HEX_125);
case 5: return polynomialBases::find(MSH_HEX_216);
case 6: return polynomialBases::find(MSH_HEX_343);
@@ -178,7 +177,6 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
SPoint3 pnt1, pnt2, pnt3;
double J1[3][3], J2[3][3], J3[3][3];
-
/*
0
0 1
@@ -252,12 +250,9 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
- hex->pnt(U1,V1,W1, pnt1);
- hex->pnt(U2,V2,W2, pnt2);
- hex->pnt(U3,V3,W3, pnt3);
- // hex->getJacobian(U1,V1,W1, J1);
- // hex->getJacobian(U2,V2,W2, J2);
- // hex->getJacobian(U3,V3,W3, J3);
+ hex->pnt(U1, V1, W1, pnt1);
+ hex->pnt(U2, V2, W2, pnt2);
+ hex->pnt(U3, V3, W3, pnt3);
}
else{
double U1 =
@@ -313,45 +308,29 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
- hex->pnt(U1,V1,W1, pnt1);
- hex->pnt(U2,V2,W2, pnt2);
- hex->pnt(U3,V3,W3, pnt3);
- // hex->getJacobian(U1,V1,W1, J1);
- // hex->getJacobian(U2,V2,W2, J2);
- // hex->getJacobian(U3,V3,W3, J3);
- }
- /*
- {
- SVector3 d1(J1[0][0], J1[0][1], J1[0][2]);
- SVector3 d2(J1[1][0], J1[1][1], J1[1][2]);
- n[0] = crossprod(d1, d2);
- n[0].normalize();
- }
- {
- SVector3 d1(J2[0][0], J2[0][1], J2[0][2]);
- SVector3 d2(J2[1][0], J2[1][1], J2[1][2]);
- n[1] = crossprod(d1, d2);
- n[1].normalize();
+ hex->pnt(U1, V1, W1, pnt1);
+ hex->pnt(U2, V2, W2, pnt2);
+ hex->pnt(U3, V3, W3, pnt3);
}
- {
- SVector3 d1(J3[0][0], J3[0][1], J3[0][2]);
- SVector3 d2(J3[1][0], J3[1][1], J3[1][2]);
- n[2] = crossprod(d1, d2);
- n[2].normalize();
- }
- */
- n[0] = 1;
- n[1] = 1;
- n[2] = 1;
+
x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
-}
+ SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
+ SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
+ n[0] = crossprod(d1, d2);
+ n[0].normalize();
+ n[1] = n[0];
+ n[2] = n[0];
+}
void MHexahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
-int MHexahedronN::getNumFacesRep(){ return 6 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2); }
+int MHexahedronN::getNumFacesRep()
+{
+ return 6 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index b90a83b29d..74e1fb0d0d 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -502,17 +502,16 @@ class MHexahedronN : public MHexahedron {
MHexahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
MVertex *v4, MVertex *v5, MVertex *v6, MVertex *v7,
const std::vector<MVertex*> &v, char order, int num=0, int part=0)
- : MHexahedron(v0,v1,v2,v3,v4,v5,v6,v7, num, part), _order(order), _vs(v)
- {
- for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
- }
+ : MHexahedron(v0, v1, v2, v3, v4, v5, v6, v7, num, part), _order(order), _vs(v)
+ {
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
MHexahedronN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
: MHexahedron(v[0], v[1], v[2], v[3],v[4], v[5], v[6], v[7], num, part), _order(order)
{
for(unsigned int i = 8; i < v.size(); i++) _vs.push_back(v[i]);
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
-
~MHexahedronN(){}
virtual int getPolynomialOrder() const { return (int)_order; }
virtual int getNumVertices() const { return 8 + _vs.size(); }
@@ -541,7 +540,7 @@ class MHexahedronN : public MHexahedron {
{
v.resize(_order+1);
MHexahedron::_getEdgeVertices(num, v);
- for (int i=0;i<_order-1;i++) v[2+i] = _vs[(_order-1)*num+i];
+ for (int i = 0; i < _order - 1; i++) v[2+i] = _vs[(_order-1)*num+i];
}
//virtual int getNumFacesRep();
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
@@ -557,12 +556,12 @@ class MHexahedronN : public MHexahedron {
{4, 5, 6, 7}
};
int count = 4;
- for (int i=0;i<4;i++){
- for (int j=0;j<_order-1;j++) v[count++] = _vs[(_order-1)*f[num][i]+j];
+ for (int i = 0; i < 4; i++){
+ for (int j = 0; j < _order - 1; j++) v[count++] = _vs[(_order-1)*f[num][i]+j];
}
- for (int i=0;i<(_order+1)*(_order+1);i++){
+ for (int i = 0; i < (_order + 1) * (_order + 1); i++){
int N = _order - 1;
- int start = 8 + 12 * N + num * (_order-1)*(_order-1);
+ int start = 8 + 12 * N + num * (_order - 1) * (_order - 1);
v[count++] = _vs[start + i];
}
}
@@ -585,16 +584,8 @@ class MHexahedronN : public MHexahedron {
if(_order == 9 && _vs.size() + 8 == 488) return MSH_HEX_488;
return 0;
}
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o){
- MElement::getShapeFunctions (u,v,w,s,o);
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
- {
- MElement::getGradShapeFunctions (u,v,w,s,o);
- }
virtual int getNumFacesRep();
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
-
virtual void revert()
{
Msg::Error("not done yet reverse hexN");
--
GitLab