From a29b84d6057eda129638cadf19aeb05ba184e1b3 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 2 Oct 2009 15:12:54 +0000
Subject: [PATCH] High order quads
---
Common/GmshDefines.h | 7 ++
Geo/GModelIO_Mesh.cpp | 8 +-
Geo/MElement.cpp | 11 ++-
Geo/MQuadrangle.cpp | 145 +++++++++++++++++++++++++++++
Geo/MQuadrangle.h | 187 +++++++++++++++++++++++++++-----------
Mesh/HighOrder.cpp | 176 +++++++++++++++--------------------
Numeric/functionSpace.cpp | 156 +++++++++++++++++++++++++++++++
7 files changed, 529 insertions(+), 161 deletions(-)
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index d4b66d6be4..bbd732ae8b 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -90,6 +90,13 @@
#define MSH_TET_52 33
#define MSH_POLYG_ 34
#define MSH_POLYH_ 35
+#define MSH_QUA_16 36
+#define MSH_QUA_25 37
+#define MSH_QUA_36 38
+#define MSH_QUA_12 39
+#define MSH_QUA_16I 40
+#define MSH_QUA_20 41
+
#define MSH_NUM_TYPE 35
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index b7745f466a..78969fc777 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -477,7 +477,7 @@ int GModel::readMSH(const std::string &name)
static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elements,
int t1, int t2=0, int t3=0, int t4=0,
- int t5=0, int t6=0, int t7=0, int t8=0)
+ int t5=0, int t6=0, int t7=0, int t8=0, int t9=0)
{
if(!binary) return;
@@ -515,6 +515,10 @@ static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elem
data[0] = t8; data[1] = elements[t8]; data[2] = numTags;
fwrite(data, sizeof(int), 3, fp);
}
+ else if(t9 && elements.count(t9)){
+ data[0] = t9; data[1] = elements[t9]; data[2] = numTags;
+ fwrite(data, sizeof(int), 3, fp);
+ }
}
template<class T>
@@ -704,7 +708,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
- writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8);
+ writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8, MSH_QUA_16, MSH_QUA_25, MSH_QUA_36,MSH_QUA_12,MSH_QUA_16I,MSH_QUA_20);
for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
if(itP->first->getType() == TYPE_QUA) {
writeElementsMSH(fp, itP->first, saveAll, version, binary, num,
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 5231fff5e4..14c715ed05 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -684,7 +684,13 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_TRI_21 : if(name) *name = "Triangle 21"; return 3 + 12 + 6;
case MSH_QUA_4 : if(name) *name = "Quadrilateral 4"; return 4;
case MSH_QUA_8 : if(name) *name = "Quadrilateral 8"; return 4 + 4;
- case MSH_QUA_9 : if(name) *name = "Quadrilateral 9"; return 4 + 4 + 1;
+ case MSH_QUA_9 : if(name) *name = "Quadrilateral 9"; return 9;
+ case MSH_QUA_16 : if(name) *name = "Quadrilateral 16"; return 16;
+ case MSH_QUA_25 : if(name) *name = "Quadrilateral 25"; return 25;
+ case MSH_QUA_36 : if(name) *name = "Quadrilateral 36"; return 36;
+ case MSH_QUA_12 : if(name) *name = "Quadrilateral 12"; return 12;
+ case MSH_QUA_16I : if(name) *name = "Quadrilateral 16I"; return 16;
+ case MSH_QUA_20 : if(name) *name = "Quadrilateral 20"; return 20;
case MSH_POLYG_ : if(name) *name = "Polygon"; return 0;
case MSH_TET_4 : if(name) *name = "Tetrahedron 4"; return 4;
case MSH_TET_10 : if(name) *name = "Tetrahedron 10"; return 4 + 6;
@@ -731,6 +737,9 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_QUA_4: return new MQuadrangle(v, num, part);
case MSH_QUA_8: return new MQuadrangle8(v, num, part);
case MSH_QUA_9: return new MQuadrangle9(v, num, part);
+ case MSH_QUA_16: return new MQuadrangleN(v, 3, num, part);
+ case MSH_QUA_25: return new MQuadrangleN(v, 4, num, part);
+ case MSH_QUA_36: return new MQuadrangleN(v, 5, num, part);
case MSH_POLYG_: return new MPolygon(v, num, part);
case MSH_TET_4: return new MTetrahedron(v, num, part);
case MSH_TET_10: return new MTetrahedron10(v, num, part);
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 9c689220bc..6880925b40 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -7,6 +7,151 @@
#include "GaussLegendre1D.h"
#include "Context.h"
#include "qualityMeasures.h"
+#define SQU(a) ((a)*(a))
+const functionSpace* MQuadrangle::getFunctionSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nf = getNumFaceVertices();
+
+ if ((nf == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return &functionSpaces::find(MSH_QUA_4);
+ case 2: return &functionSpaces::find(MSH_QUA_8);
+ case 3: return &functionSpaces::find(MSH_QUA_12);
+ case 4: return &functionSpaces::find(MSH_QUA_16I);
+ case 5: return &functionSpaces::find(MSH_QUA_20);
+ }
+ }
+ switch (order) {
+ case 1: return &functionSpaces::find(MSH_QUA_4);
+ case 2: return &functionSpaces::find(MSH_QUA_9);
+ case 3: return &functionSpaces::find(MSH_QUA_16);
+ case 4: return &functionSpaces::find(MSH_QUA_25);
+ case 5: return &functionSpaces::find(MSH_QUA_36);
+ default: Msg::Error("Order %d triangle function space not implemented", order);
+ }
+ return 0;
+}
+
+int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
+int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
+int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
+
+static void _myGetEdgeRep(MQuadrangle *q, int num, double *x, double *y, double *z,
+ SVector3 *n, int numSubEdges)
+{
+ n[0] = n[1] = n[2] = n[3] = q->getFace(0).normal();
+ int ie = num/numSubEdges;
+ int isub = num%numSubEdges;
+ double xi1 = -1. + (2.*isub)/numSubEdges;
+ double xi2 = -1. + (2.*(isub+1))/numSubEdges;
+ SPoint3 pnt1, pnt2;
+ switch(ie){
+ case 0:
+ q->pnt( xi1, -1., 0., pnt1);
+ q->pnt( xi2, -1., 0., pnt2);
+ break;
+ case 1:
+ q->pnt( 1., xi1, 0., pnt1);
+ q->pnt( 1., xi2, 0., pnt2);
+ break;
+ case 2:
+ q->pnt( xi1, 1., 0., pnt1);
+ q->pnt( xi2, 1., 0., pnt2);
+ break;
+ case 3:
+ q->pnt(-1., xi1, 0., pnt1);
+ q->pnt(-1., xi2, 0., pnt2);
+ break;
+ }
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
+}
+
+void MQuadrangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+void MQuadrangle8::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+void MQuadrangle9::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+
+int MQuadrangleN::getNumFacesRep(){ return 2*SQU(CTX::instance()->mesh.numSubEdges); }
+int MQuadrangle8::getNumFacesRep(){ return 2*SQU(CTX::instance()->mesh.numSubEdges); }
+int MQuadrangle9::getNumFacesRep(){ return 2*SQU(CTX::instance()->mesh.numSubEdges); }
+
+static void _myGetFaceRep(MQuadrangle *t, int num, double *x, double *y, double *z,
+ SVector3 *n, int numSubEdges)
+{
+ int io = num%2;
+ int ix = (num/2)/numSubEdges;
+ int iy = (num/2)%numSubEdges;
+
+ const double d = 2. / numSubEdges;
+ const double ox = -1. + d*ix;
+ const double oy = -1. + d*iy;
+
+ SPoint3 pnt1, pnt2, pnt3;
+ double J1[3][3], J2[3][3], J3[3][3];
+ if (io == 0){
+ t->pnt(ox, oy, 0, pnt1);
+ t->pnt(ox + d, oy, 0, pnt2);
+ t->pnt(ox + d, oy + d, 0, pnt3);
+ t->getJacobian(ox, oy, 0, J1);
+ t->getJacobian(ox + d, oy, 0, J2);
+ t->getJacobian(ox + d, oy + d, 0, J3);
+ } else{
+ t->pnt(ox, oy, 0, pnt1);
+ t->pnt(ox + d, oy + d, 0, pnt2);
+ t->pnt(ox, oy + d, 0, pnt3);
+ t->getJacobian(ox, oy, 0, J1);
+ t->getJacobian(ox + d, oy + d, 0, J2);
+ t->getJacobian(ox, oy + d, 0, 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();
+ }
+
+ 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();
+}
+
+void MQuadrangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+void MQuadrangle8::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+void MQuadrangle9::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 376265d3fb..62603daf10 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -115,24 +115,11 @@ class MQuadrangle : public MElement {
virtual const char *getStringForPOS() const { return "SQ"; }
virtual const char *getStringForBDF() const { return "CQUAD4"; }
virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
+ virtual const functionSpace* getFunctionSpace(int o=-1) const;
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
}
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
- {
- s[0] = (1. - u) * (1. - v) * 0.25;
- s[1] = (1. + u) * (1. - v) * 0.25;
- s[2] = (1. + u) * (1. + v) * 0.25;
- s[3] = (1. - u) * (1. + v) * 0.25;
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
- {
- s[0][0] = -0.25 * (1. - v); s[0][1] = -0.25 * (1. - u); s[0][2] = 0.;
- s[1][0] = 0.25 * (1. - v); s[1][1] = -0.25 * (1. + u); s[1][2] = 0.;
- s[2][0] = 0.25 * (1. + v); s[2][1] = 0.25 * (1. + u); s[2][2] = 0.;
- s[3][0] = -0.25 * (1. + v); s[3][1] = 0.25 * (1. - u); s[3][2] = 0.;
- }
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
@@ -203,32 +190,16 @@ class MQuadrangle8 : public MQuadrangle {
return getVertex(map[num]);
}
virtual int getNumEdgeVertices() const { return 4; }
- virtual int getNumEdgesRep(){ return 8; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[8][2] = {
- {0, 4}, {4, 1},
- {1, 5}, {5, 2},
- {2, 6}, {6, 3},
- {3, 7}, {7, 0}
- };
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MQuadrangle::_getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 6; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[6][3] = {
- {0, 4, 7}, {1, 5, 4}, {2, 6, 5}, {3, 7, 6}, {4, 5, 6}, {4, 6, 7}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+ virtual int getNumFacesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(8);
@@ -292,33 +263,16 @@ class MQuadrangle9 : public MQuadrangle {
}
virtual int getNumEdgeVertices() const { return 4; }
virtual int getNumFaceVertices() const { return 1; }
- virtual int getNumEdgesRep(){ return 8; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[8][2] = {
- {0, 4}, {4, 1},
- {1, 5}, {5, 2},
- {2, 6}, {6, 3},
- {3, 7}, {7, 0}
- };
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MQuadrangle::getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 8; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[8][4] = {
- {0, 4, 8}, {0, 8, 7}, {1, 5, 8}, {1, 8, 4},
- {2, 6, 8}, {2, 8, 5}, {3, 7, 8}, {3, 8, 6}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+ virtual int getNumFacesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(9);
@@ -341,4 +295,127 @@ class MQuadrangle9 : public MQuadrangle {
}
};
+/*
+ * MQuadrangle
+ *
+ * 3--3+3E-...--4+2E--2
+ * | | E = order - 1;
+ * | | N = total number of vertices
+ * 4+3E 3+2E
+ * | | Interior vertex numbers
+ * ... 4+4e to N-1 ... for edge 0 <= i <= 3: 4+i*E to 3+(i+1)*E
+ * | | in volume : 4+4*E to N-1
+ * 3+4E 4+E
+ * | |
+ * | |
+ * 0---4--...---3+E---1
+ *
+ */
+class MQuadrangleN : public MQuadrangle {
+ protected:
+ std::vector<MVertex *> _vs;
+ const char _order;
+ public:
+ MQuadrangleN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
+ std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MQuadrangle(v0, v1, v2, v3, num, part), _vs(v), _order(order)
+ {
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ MQuadrangleN(std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MQuadrangle(v[0], v[1], v[2], v[3], num, part), _order(order)
+ {
+ for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ ~MQuadrangleN(){}
+ virtual int getPolynomialOrder() const { return _order; }
+ virtual int getNumVertices() const {return 4 + _vs.size(); }
+ virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual int getNumFaceVertices() const
+ {
+ if( _order>1 && _vs.size()+4==(_order+1)*(_order+1))
+ return (_order-1)*(_order-1);
+ else
+ return 0;
+ }
+ virtual int getNumEdgeVertices() const { return 4 * (_order - 1); }
+ virtual int getNumEdgesRep();
+ virtual int getNumFacesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize(_order + 1);
+ MQuadrangle::_getEdgeVertices(num, v);
+ int j = 2;
+ const int ie = (num + 1) * (_order - 1);
+ for(int i = num * (_order-1); i != ie; ++i)
+ v[j++] = _vs[i];
+ }
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize(4 + _vs.size());
+ MQuadrangle::_getFaceVertices(v);
+ for(unsigned int i = 0; i != _vs.size(); ++i) v[i + 4] = _vs[i];
+ }
+ virtual int getTypeForMSH() const
+ {
+ if(_order==2 && _vs.size()+4==8) return MSH_QUA_8;
+ if(_order==3 && _vs.size()+4==12) return MSH_QUA_12;
+ if(_order==3 && _vs.size()+4==16) return MSH_QUA_16;
+ if(_order==4 && _vs.size()+4==16) return MSH_QUA_16I;
+ if(_order==4 && _vs.size()+4==25) return MSH_QUA_25;
+ if(_order==5 && _vs.size()+4==20) return MSH_QUA_20;
+ if(_order==5 && _vs.size()+4==36) return MSH_QUA_36;
+ return 0;
+ }
+ virtual void revert()
+ {
+ MVertex *tmp;
+ tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
+ std::vector<MVertex*> inv;
+ inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
+ _vs = inv;
+ }
+};
+
+template <class T>
+void inline sort2(T &a, T &b){
+ if(b<a){
+ T t=b;
+ b=a;
+ a=t;
+ }
+}
+
+template <class T>
+void sort4(T *t[3])
+{
+ sort2<T*>(t[0],t[1]);
+ sort2<T*>(t[2],t[3]);
+ sort2<T*>(t[0],t[2]);
+ sort2<T*>(t[1],t[3]);
+ sort2<T*>(t[1],t[2]);
+}
+
+struct compareMQuadrangleLexicographic
+{
+ bool operator () (MQuadrangle *t1, MQuadrangle *t2) const
+ {
+ MVertex *_v1[] = {t1->getVertex(0), t1->getVertex(1), t1->getVertex(2), t1->getVertex(3)};
+ MVertex *_v2[] = {t2->getVertex(0), t2->getVertex(1), t2->getVertex(2), t2->getVertex(3)};
+ sort4(_v1);
+ sort4(_v2);
+ if(_v1[0] < _v2[0]) return true;
+ if(_v1[0] > _v2[0]) return false;
+ if(_v1[1] < _v2[1]) return true;
+ if(_v1[1] > _v2[1]) return false;
+ if(_v1[2] < _v2[2]) return true;
+ if(_v1[2] > _v2[2]) return false;
+ if(_v1[3] < _v1[3]) return true;
+ return false;
+ }
+};
+
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index e5c10cbadd..cc2af1b450 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -408,28 +408,6 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
if(gf->geomType() == GEntity::DiscreteSurface ||
gf->geomType() == GEntity::BoundaryLayerSurface)
linear = true;
-
- fullMatrix<double> points;
- int start = 0;
-
- switch (nPts){
- case 2 :
- points = functionSpaces::find(MSH_TRI_10).points;
- start = 9;
- break;
- case 3 :
- points = functionSpaces::find(MSH_TRI_15).points;
- start = 12;
- break;
- case 4 :
- points = functionSpaces::find(MSH_TRI_21).points;
- start = 15;
- break;
- default :
- // do nothing (e.g. for 2nd order tri faces or for quad faces)
- break;
- }
-
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
faceContainer::iterator fIter = faceVertices.find(face);
@@ -444,81 +422,54 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
for(int k = 0; k < incomplete->getNumVertices(); k++)
reparamOK &= reparamMeshVertexOnFace(incomplete->getVertex(k), gf, pts[k]);
}
- if(face.getNumVertices() == 3 && nPts > 1){ // tri face
- for(int k = start; k < points.size1(); k++){
- MVertex *v;
- const double t1 = points(k, 0);
- const double t2 = points(k, 1);
- if(linear){
- SPoint3 pc = face.interpolate(t1, t2);
- v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
- }
- else{
- double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
- double sf[256];
- incomplete->getShapeFunctions(t1, t2, 0, sf);
- for (int j = 0; j < incomplete->getNumVertices(); j++){
- MVertex *vt = incomplete->getVertex(j);
- X += sf[j] * vt->x();
- Y += sf[j] * vt->y();
- Z += sf[j] * vt->z();
- if (reparamOK){
- GUESS[0] += sf[j] * pts[j][0];
- GUESS[1] += sf[j] * pts[j][1];
- }
+ int start = face.getNumVertices()*(nPts+1);
+ const fullMatrix<double> &points = ele->getFunctionSpace(nPts+1)->points;
+ for(int k = start; k < points.size1(); k++){
+ MVertex *v;
+ const double t1 = points(k, 0);
+ const double t2 = points(k, 1);
+ if(linear){
+ SPoint3 pc = face.interpolate(t1, t2);
+ v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
+ }
+ else{
+ double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
+ double sf[256];
+ incomplete->getShapeFunctions(-1., -1., 0, sf);
+ for (int j = 0; j < incomplete->getNumVertices(); j++)
+ incomplete->getShapeFunctions(t1, t2, 0, sf);
+ for (int j = 0; j < incomplete->getNumVertices(); j++){
+ MVertex *vt = incomplete->getVertex(j);
+ X += sf[j] * vt->x();
+ Y += sf[j] * vt->y();
+ Z += sf[j] * vt->z();
+ if (reparamOK){
+ GUESS[0] += sf[j] * pts[j][0];
+ GUESS[1] += sf[j] * pts[j][1];
}
- if(reparamOK){
- GPoint gp = gf->closestPoint(SPoint3(X, Y, Z), GUESS);
- if (gp.g()){
- v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
- }
- else{
- v = new MVertex(X, Y, Z, gf);
- }
+ }
+ if(reparamOK){
+ GPoint gp = gf->closestPoint(SPoint3(X, Y, Z), GUESS);
+ if (gp.g()){
+ v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
}
else{
v = new MVertex(X, Y, Z, gf);
}
- if(displ3D || displ2D){
- SPoint3 pc2 = face.interpolate(t1, t2);
- if(displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
- if(displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
- }
}
- // should be expensive -> induces a new search each time
- vtcs.push_back(v);
- gf->mesh_vertices.push_back(v);
- vf.push_back(v);
- }
- }
- else if(face.getNumVertices() == 4){ // quad face
- for(int j = 0; j < nPts; j++){
- for(int k = 0; k < nPts; k++){
- MVertex *v;
- // parameters are between -1 and 1
- double t1 = 2. * (double)(j + 1) / (nPts + 1) - 1.;
- double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
- if(linear || !reparamOK){
- SPoint3 pc = face.interpolate(t1, t2);
- v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
- }
- else{
- double uc = 0.25 * ((1 - t1) * (1 - t2) * pts[0][0] +
- (1 + t1) * (1 - t2) * pts[1][0] +
- (1 + t1) * (1 + t2) * pts[2][0] +
- (1 - t1) * (1 + t2) * pts[3][0]);
- double vc = 0.25 * ((1 - t1) * (1 - t2) * pts[0][1] +
- (1 + t1) * (1 - t2) * pts[1][1] +
- (1 + t1) * (1 + t2) * pts[2][1] +
- (1 - t1) * (1 + t2) * pts[3][1]);
- GPoint pc = gf->point(uc, vc);
- v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
- }
- vtcs.push_back(v);
- gf->mesh_vertices.push_back(v);
- vf.push_back(v);
+ else{
+ v = new MVertex(X, Y, Z, gf);
}
+ if(displ3D || displ2D){
+ SPoint3 pc2 = face.interpolate(t1, t2);
+ if(displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if(displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ }
}
+ // should be expensive -> induces a new search each time
+ vtcs.push_back(v);
+ gf->mesh_vertices.push_back(v);
+ vf.push_back(v);
}
}
}
@@ -737,6 +688,37 @@ MTriangle* setHighOrder(MTriangle *t, GFace *gf,
}
}
}
+MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts,
+ highOrderSmoother *displ2D,
+ highOrderSmoother *displ3D)
+{
+ std::vector<MVertex*> ve, vf;
+ getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+ if(incomplete){
+ if(nPts==1){
+ return new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2),q->getVertex(3), ve[0],ve[1],ve[2],ve[3]);
+ }else{
+ return new MQuadrangleN(q->getVertex(0), q->getVertex(1), q->getVertex(2), q->getVertex(3), ve, nPts + 1);
+ }
+ } else {
+ if (displ2D && gf->geomType() == GEntity::Plane){
+ MQuadrangle incpl(q->getVertex(0), q->getVertex(1), q->getVertex(2), q->getVertex(3));
+ getFaceVertices(gf, &incpl, q, vf, faceVertices, linear, nPts, displ2D, displ3D);
+ }else{
+ MQuadrangleN incpl(q->getVertex(0), q->getVertex(1), q->getVertex(2), q->getVertex(3), ve, nPts + 1);
+ getFaceVertices(gf, &incpl, q, vf, faceVertices, linear, nPts, displ2D, displ3D);
+ }
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ if(nPts==1){
+ return new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),q->getVertex(3), ve[0], ve[1], ve[2], ve[3], vf[0]);
+ }else{
+ return new MQuadrangleN(q->getVertex(0), q->getVertex(1), q->getVertex(2), q->getVertex(3), ve, nPts + 1);
+ }
+ }
+}
static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
faceContainer &faceVertices, bool linear, bool incomplete,
@@ -752,23 +734,11 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
delete t;
}
gf->triangles = triangles2;
-
std::vector<MQuadrangle*> quadrangles2;
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
- std::vector<MVertex*> ve, vf;
- getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts, displ2D, displ3D);
- if(incomplete){
- quadrangles2.push_back
- (new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2),
- q->getVertex(3), ve[0], ve[1], ve[2], ve[3]));
- }
- else{
- getFaceVertices(gf, q, q, vf, faceVertices, linear, nPts);
- quadrangles2.push_back
- (new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),
- q->getVertex(3), ve[0], ve[1], ve[2], ve[3], vf[0]));
- }
+ MQuadrangle *qNew = setHighOrder(q,gf,edgeVertices,faceVertices,linear,incomplete,nPts,displ2D,displ3D);
+ quadrangles2.push_back(qNew);
delete q;
}
gf->quadrangles = quadrangles2;
diff --git a/Numeric/functionSpace.cpp b/Numeric/functionSpace.cpp
index fe26349a3d..d54c96c461 100644
--- a/Numeric/functionSpace.cpp
+++ b/Numeric/functionSpace.cpp
@@ -72,6 +72,83 @@ static fullMatrix<double> generatePascalSerendipityTriangle(int order)
return monomials;
}
+// generate all monomials xi^m * eta^n with n and m <= order
+static fullMatrix<double> generatePascalQuad(int order)
+{
+
+ fullMatrix<double> monomials( (order+1)*(order+1), 2);
+ int index = 0;
+ for (int p = 0; p <= order; p++) {
+ for(int i = 0; i < p; i++, index++) {
+ monomials(index, 0) = p;
+ monomials(index, 1) = i;
+ }
+ for(int i = 0; i <= p; i++, index++) {
+ monomials(index, 0) = p-i;
+ monomials(index, 1) = p;
+ }
+ }
+ return monomials;
+}
+/*
+00 10 20 30 40 ⋯
+01 11 21 31 41 ⋯
+02 12
+03 13
+04 14
+â‹® â‹®
+*/
+static fullMatrix<double> generatePascalQuadSerendip(int order)
+{
+ fullMatrix<double> monomials( (order)*4, 2);
+ monomials(0,0)=0.;
+ monomials(0,1)=0.;
+ monomials(1,0)=1.;
+ monomials(1,1)=0.;
+ monomials(2,0)=0.;
+ monomials(2,1)=1.;
+ monomials(3,0)=1.;
+ monomials(3,1)=1.;
+ int index = 4;
+ for (int p = 2; p <= order; p++) {
+ monomials(index, 0) = p;
+ monomials(index, 1) = 0;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = p;
+ index++;
+ monomials(index, 0) = p;
+ monomials(index, 1) = 1;
+ index++;
+ monomials(index, 0) = 1;
+ monomials(index, 1) = p;
+ index++;
+ }
+ return monomials;
+}
+
+/*static fullMatrix<double> generatePascalQuadSerendip(int order)
+{
+
+ fullMatrix<double> monomials( order*4, 2);
+ int index = 0;
+ for (int p = 0; p < order; p++) {
+ monomials(p, 0) = p;
+ monomials(p, 1) = 0;
+
+ monomials(p+order, 0) = order;
+ monomials(p+order, 1) = p;
+
+ monomials(p+3*order, 0) = order-p;
+ monomials(p+3*order, 1) = order;
+
+ monomials(p+2*order, 0) = 0;
+ monomials(p+2*order, 1) = order-p;
+ }
+ monomials.print();
+ return monomials;
+}*/
+
// generate the monomials subspace of all monomials of order exactly == p
static fullMatrix<double> generateMonomialSubspace(int dim, int p)
@@ -153,6 +230,8 @@ static fullMatrix<double> generatePascalTetrahedron(int order)
return monomials;
}
+
+
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
//static int nbdoftriangleserendip(int order) { return 3 * order; }
@@ -517,6 +596,47 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
return point;
}
+static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
+{
+ int nbPoints = serendip ? order*4 : (order+1)*(order+1);
+ fullMatrix<double> point(nbPoints, 2);
+
+ double dd = 1. / order;
+
+ if (order > 0) {
+ point(0, 0) = -1;
+ point(0, 1) = -1;
+ point(1, 0) = 1;
+ point(1, 1) = -1;
+ point(2, 0) = 1;
+ point(2, 1) = 1;
+ point(3, 0) = -1;
+ point(3, 1) = 1;
+
+ if (order > 1) {
+ int index = 4;
+ const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
+ for (int iedge=0; iedge<4; iedge++) {
+ int p0 = edges[iedge][0];
+ int p1 = edges[iedge][1];
+ for (int i = 1; i < order; i++, index++) {
+ point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
+ point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
+ }
+ }
+ if (order > 2 && !serendip) {
+ fullMatrix<double> inner = gmshGeneratePointsQuad(order - 2, false);
+ inner.scale(1. - 2./order);
+ point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ } else {
+ point(0, 0) = 0;
+ point(0, 1) = 0;
+ }
+ return point;
+}
+
static fullMatrix<double> generateLagrangeMonomialCoefficients
(const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{
@@ -747,6 +867,42 @@ const functionSpace &functionSpaces::find(int tag)
F.points = gmshGeneratePointsTetrahedron(5, false);
generate3dFaceClosure (F.faceClosure,5);
break;
+ case MSH_QUA_4 :
+ F.monomials = generatePascalQuad(1);
+ F.points = gmshGeneratePointsQuad(1,false);
+ break;
+ case MSH_QUA_9 :
+ F.monomials = generatePascalQuad(2);
+ F.points = gmshGeneratePointsQuad(2,false);
+ break;
+ case MSH_QUA_16 :
+ F.monomials = generatePascalQuad(3);
+ F.points = gmshGeneratePointsQuad(3,false);
+ break;
+ case MSH_QUA_25 :
+ F.monomials = generatePascalQuad(4);
+ F.points = gmshGeneratePointsQuad(4,false);
+ break;
+ case MSH_QUA_36 :
+ F.monomials = generatePascalQuad(5);
+ F.points = gmshGeneratePointsQuad(5,false);
+ break;
+ case MSH_QUA_8 :
+ F.monomials = generatePascalQuadSerendip(2);
+ F.points = gmshGeneratePointsQuad(2,true);
+ break;
+ case MSH_QUA_12 :
+ F.monomials = generatePascalQuadSerendip(3);
+ F.points = gmshGeneratePointsQuad(3,true);
+ break;
+ case MSH_QUA_16I :
+ F.monomials = generatePascalQuadSerendip(4);
+ F.points = gmshGeneratePointsQuad(4,true);
+ break;
+ case MSH_QUA_20 :
+ F.monomials = generatePascalQuadSerendip(5);
+ F.points = gmshGeneratePointsQuad(5,true);
+ break;
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.monomials = generatePascalTetrahedron(1);
--
GitLab