Skip to content
Snippets Groups Projects
Commit 9b9979f6 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix drawing and shape functions for order N

parent 714b5a65
No related branches found
No related tags found
No related merge requests found
...@@ -113,20 +113,19 @@ void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 ...@@ -113,20 +113,19 @@ void MHexahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3
n[0] = n[1] = 1 ; n[0] = n[1] = 1 ;
} }
int MHexahedronN::getNumEdgesRep(){ int MHexahedronN::getNumEdgesRep()
{
return 12 * CTX::instance()->mesh.numSubEdges; return 12 * CTX::instance()->mesh.numSubEdges;
} }
//int MHexaedronN::getNumFacesRep(){
// return 8 * SQU(CTX::instance()->mesh.numSubEdges);
//}
const polynomialBasis* MHexahedron::getFunctionSpace(int o) const const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
{ {
int order = (o == -1) ? getPolynomialOrder() : o; int order = (o == -1) ? getPolynomialOrder() : o;
int nv = getNumVolumeVertices(); int nv = getNumVolumeVertices();
printf("nv = %d\n", nv);
if ((nv == 0) && (o == -1)) { if ((nv == 0) && (o == -1)) {
switch (order) { switch (order) {
case 0: return polynomialBases::find(MSH_HEX_1); case 0: return polynomialBases::find(MSH_HEX_1);
...@@ -147,7 +146,7 @@ const polynomialBasis* MHexahedron::getFunctionSpace(int o) const ...@@ -147,7 +146,7 @@ const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
case 0: return polynomialBases::find(MSH_HEX_1); case 0: return polynomialBases::find(MSH_HEX_1);
case 1: return polynomialBases::find(MSH_HEX_8); case 1: return polynomialBases::find(MSH_HEX_8);
case 2: return polynomialBases::find(MSH_HEX_27); 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 4: return polynomialBases::find(MSH_HEX_125);
case 5: return polynomialBases::find(MSH_HEX_216); case 5: return polynomialBases::find(MSH_HEX_216);
case 6: return polynomialBases::find(MSH_HEX_343); case 6: return polynomialBases::find(MSH_HEX_343);
...@@ -178,7 +177,6 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl ...@@ -178,7 +177,6 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
SPoint3 pnt1, pnt2, pnt3; SPoint3 pnt1, pnt2, pnt3;
double J1[3][3], J2[3][3], J3[3][3]; double J1[3][3], J2[3][3], J3[3][3];
/* /*
0 0
0 1 0 1
...@@ -255,9 +253,6 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl ...@@ -255,9 +253,6 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
hex->pnt(U1, V1, W1, pnt1); hex->pnt(U1, V1, W1, pnt1);
hex->pnt(U2, V2, W2, pnt2); hex->pnt(U2, V2, W2, pnt2);
hex->pnt(U3, V3, W3, pnt3); hex->pnt(U3, V3, W3, pnt3);
// hex->getJacobian(U1,V1,W1, J1);
// hex->getJacobian(U2,V2,W2, J2);
// hex->getJacobian(U3,V3,W3, J3);
} }
else{ else{
double U1 = double U1 =
...@@ -316,42 +311,26 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl ...@@ -316,42 +311,26 @@ static void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, doubl
hex->pnt(U1, V1, W1, pnt1); hex->pnt(U1, V1, W1, pnt1);
hex->pnt(U2, V2, W2, pnt2); hex->pnt(U2, V2, W2, pnt2);
hex->pnt(U3, V3, W3, pnt3); 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();
} }
{
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(); x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y(); y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z(); 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) void MHexahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{ {
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges); _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);
}
...@@ -512,7 +512,6 @@ class MHexahedronN : public MHexahedron { ...@@ -512,7 +512,6 @@ class MHexahedronN : public MHexahedron {
for(unsigned int i = 8; i < v.size(); i++) _vs.push_back(v[i]); 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); for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
} }
~MHexahedronN(){} ~MHexahedronN(){}
virtual int getPolynomialOrder() const { return (int)_order; } virtual int getPolynomialOrder() const { return (int)_order; }
virtual int getNumVertices() const { return 8 + _vs.size(); } virtual int getNumVertices() const { return 8 + _vs.size(); }
...@@ -585,16 +584,8 @@ class MHexahedronN : public MHexahedron { ...@@ -585,16 +584,8 @@ class MHexahedronN : public MHexahedron {
if(_order == 9 && _vs.size() + 8 == 488) return MSH_HEX_488; if(_order == 9 && _vs.size() + 8 == 488) return MSH_HEX_488;
return 0; 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 int getNumFacesRep();
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n); virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void revert() virtual void revert()
{ {
Msg::Error("not done yet reverse hexN"); Msg::Error("not done yet reverse hexN");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment