From f8a3ffbe6096f57008d28b06d3eebf355538f7bc Mon Sep 17 00:00:00 2001
From: Koen Hillewaert <koen.hillewaert@cenaero.be>
Date: Thu, 28 Aug 2008 09:47:07 +0000
Subject: [PATCH] bug fix: completed MElement interface (pnt, shapeFunction,
gradShapeFunction) bug fix: orientation of face vertices of 4th order
tetrahedra
---
Geo/MElement.cpp | 220 ++++++++++++++++++++++++++++++++++++--
Geo/MElement.h | 31 +++++-
Mesh/HighOrder.cpp | 69 +++++++-----
Numeric/FunctionSpace.cpp | 39 +++++++
4 files changed, 321 insertions(+), 38 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 7b0eb4a721..89c97adfbd 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -423,6 +423,26 @@ double MElement::interpolateDiv(double val[], double u, double v, double w, int
interpolateGrad(&val[2], u, v, w, fz, stride, inv);
return fx[0] + fy[1] + fz[2];
}
+
+
+// Expensive, generic version
+
+void MElement::pnt(double u,double v,double w,SPoint3& p) {
+
+ double x(0),y(0),z(0);
+
+ for (int i=0;i<getNumVertices();i++) {
+ double s;
+ getShapeFunction(i,u,v,w,s);
+
+ MVertex* v = getVertex(i);
+
+ x += v->x() * s;
+ y += v->y() * s;
+ z += v->z() * s;
+ }
+ p = SPoint3(x,y,z);
+}
void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
int elementary, int physical)
@@ -722,6 +742,118 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
}
}
+void MLine::pnt(double u,double v,double w,SPoint3& p) {
+ double f1,f2;
+ getShapeFunction(0,u,v,w,f1);
+ getShapeFunction(1,u,v,w,f2);
+
+ double x = f1 * _v[0]->x() + f2 * _v[1]->x();
+ double y = f1 * _v[0]->y() + f2 * _v[1]->y();
+ double z = f1 * _v[0]->z() + f2 * _v[1]->z();
+
+ p = SPoint3(x,y,z);
+}
+
+void MLine3::pnt(double u, double v, double w, SPoint3 &p) {
+
+ double sf[3];
+ gmshFunctionSpaces::find(MSH_LIN_3).f(u, v, w, sf);
+ double x = sf[0] * _v[0]->x() + sf[1] * _v[1]->x() + sf[2] * _vs[0]->x();
+ double y = sf[0] * _v[0]->y() + sf[1] * _v[1]->y() + sf[2] * _vs[0]->y();
+ double z = sf[0] * _v[0]->z() + sf[1] * _v[1]->z() + sf[2] * _vs[0]->z();
+
+ p = SPoint3(x,y,z);
+}
+
+void MLine3::getShapeFunction(int num,double uu, double vv,double ww,double& s) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+
+ if (num > 2) s = 0.;
+ double sf[3];
+ gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf);
+ s = sf[num];
+
+#endif
+
+}
+
+void MLine3::getGradShapeFunction(int num,double uu,double vv,double ww,double s[3]) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+
+ if (num > 2) for (int i=0;i<3;i++) s[i] = 0.;
+ double sf[3][3];
+ gmshFunctionSpaces::find(MSH_LIN_3).df(uu, vv, ww, sf);
+ for (int i=0;i<3;i++) s[i] = sf[num][i];
+
+#endif
+
+}
+
+
+void MLineN::pnt(double uu, double vv, double ww, SPoint3 &p) {
+
+ double sf[6];
+
+ switch (getPolynomialOrder()) {
+ case 1: gmshFunctionSpaces::find(MSH_LIN_2).f(uu, vv, ww, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_LIN_4).f(uu, vv, ww, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_LIN_5).f(uu, vv, ww, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_LIN_6).f(uu, vv, ww, sf); break;
+ default: Msg::Error("Order %d line point interpolation not implemented", getPolynomialOrder()); break;
+ }
+
+ double x = 0; for (int i=0;i<2;i++) x += sf[i] * _v[i]->x();
+ double y = 0; for (int i=0;i<2;i++) y += sf[i] * _v[i]->y();
+ double z = 0; for (int i=0;i<2;i++) z += sf[i] * _v[i]->z();
+
+ for (int i=0;i<_vs.size();i++) x += sf[i+2] * _vs[i]->x();
+ for (int i=0;i<_vs.size();i++) y += sf[i+2] * _vs[i]->y();
+ for (int i=0;i<_vs.size();i++) z += sf[i+2] * _vs[i]->z();
+
+ p = SPoint3(x,y,z);
+
+}
+
+void MLineN::getShapeFunction(int num,double uu, double vv,double ww,double& s) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+
+ double sf[6];
+ switch (getPolynomialOrder()) {
+ case 1: gmshFunctionSpaces::find(MSH_LIN_2).f(uu, vv, ww, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_LIN_4).f(uu, vv, ww, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_LIN_5).f(uu, vv, ww, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_LIN_6).f(uu, vv, ww, sf); break;
+ default: Msg::Error("Order %d line shape function not provided", getPolynomialOrder()); break;
+ }
+ s = sf[num];
+
+#endif
+
+}
+
+void MLineN::getGradShapeFunction(int num,double uu,double vv,double ww,double s[3]) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+
+ double sf[6][3];
+ switch (getPolynomialOrder()) {
+ case 1: gmshFunctionSpaces::find(MSH_LIN_2).df(uu, vv, ww, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_LIN_3).df(uu, vv, ww, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_LIN_4).df(uu, vv, ww, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_LIN_5).df(uu, vv, ww, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_LIN_6).df(uu, vv, ww, sf); break;
+ default: Msg::Error("Order %d line shape function gradient not implemented", getPolynomialOrder()); break;
+ }
+ for (int i=0;i<3;i++) s[i] = sf[num][i];
+#endif
+
+}
+
void MTriangle::getShapeFunction(int num,double uu, double vv,double ww,double& s) {
#if !defined(HAVE_GMSH_EMBEDDED)
@@ -746,7 +878,7 @@ void MTriangle::getShapeFunction(int num,double uu, double vv,double ww,double&
case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, fv); break;
case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, fv); break;
case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, fv); break;
- default: Msg::Error("Order %d triangle jac not implemented", getPolynomialOrder()); break;
+ default: Msg::Error("Order %d triangle shape function not implemented", getPolynomialOrder()); break;
}
}
s = fv[num];
@@ -768,7 +900,7 @@ void MTriangle::getGradShapeFunction(int num,double uu,double vv,double ww,doubl
case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, ww, grads); break;
case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, ww, grads); break;
case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, ww, grads); break;
- default: Msg::Error("Order %d triangle jac not implemented", getPolynomialOrder()); break;
+ default: Msg::Error("Order %d triangle shape function gradient not implemented", getPolynomialOrder()); break;
}
}
else{
@@ -778,7 +910,7 @@ void MTriangle::getGradShapeFunction(int num,double uu,double vv,double ww,doubl
case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, ww, grads); break;
case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, ww, grads); break;
case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, ww, grads); break;
- default: Msg::Error("Order %d triangle jac not implemented", getPolynomialOrder()); break;
+ default: Msg::Error("Order %d triangle shape function gradient not implemented", getPolynomialOrder()); break;
}
}
@@ -946,8 +1078,10 @@ void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,
double y = 0 ; for(int i = 0; i < 4; i++) y += sf[i] * _v[i]->y();
double z = 0 ; for(int i = 0; i < 4; i++) z += sf[i] * _v[i]->z();
- const int N = (ord+1)*(ord+2)*(ord+3)/6;
+ int nbComplete = (ord+1)*(ord+2)*(ord+3)/6;
+ int nbInterior = (ord-3)*(ord-2)*(ord-1)/6;
+ const int N = nv ? nbComplete : nbComplete - nbInterior;
for(int i = 4; i < N; i++) x += sf[i] * vs[i - 4]->x();
for(int i = 4; i < N; i++) y += sf[i] * vs[i - 4]->y();
@@ -1058,6 +1192,81 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
}
}
+void MTetrahedron10::getShapeFunction(int num, double uu, double vv, double ww, double &s) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double sf[256];
+ gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww, sf);
+ s = sf[num];
+#endif
+}
+
+void MTetrahedron10::getGradShapeFunction(int num, double uu, double vv, double ww, double gs[3]) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double gsf[256][3];
+ gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, gsf);
+ for (int i=0;i<3;i++) gs[i] = gsf[num][i];
+#endif
+}
+
+
+void MTetrahedronN::getShapeFunction(int num, double uu, double vv, double ww, double &s)
+{
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double sf[256];
+
+ int nv = getNumVolumeVertices();
+
+ if (!nv){
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_34).f(uu, vv, ww, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_52).f(uu, vv, ww, sf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder()); break;
+ }
+ }
+ else{
+ switch(getPolynomialOrder()){
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww, sf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder()); break;
+ }
+ }
+
+ s = sf[num];
+#endif
+}
+
+void MTetrahedronN::getGradShapeFunction(int num, double uu, double vv, double ww, double gs[3])
+{
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double gsf[256][3];
+
+ int nv = getNumVolumeVertices();
+
+ if (!nv){
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TET_4) .df(uu, vv, ww, gsf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, gsf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, gsf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_34).df(uu, vv, ww, gsf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_52).df(uu, vv, ww, gsf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder()); break;
+ }
+ }
+ else{
+ switch(getPolynomialOrder()){
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, gsf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, gsf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder()); break;
+ }
+ }
+ for (int i=0;i<3;i++) gs[i] = gsf[num][i];
+#endif
+}
int MTetrahedronN::getNumEdgesRep(){ return 6 * numSubEdges; }
@@ -1194,9 +1403,6 @@ void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector
}
}
-
-
-
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 8d66a06ad4..060c4189f3 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -169,11 +169,13 @@ class MElement
// coordinates
virtual void getGradShapeFunction(int num, double u, double v, double w,
double s[3]) = 0;
-
+
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
double getJacobian(double u, double v, double w, double jac[3][3]);
+ virtual void pnt(double u, double v, double w, SPoint3 &p);
+
// invert the parametrisation
virtual void xyz2uvw(double xyz[3], double uvw[3]);
@@ -340,6 +342,10 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
+
+
+ virtual void pnt(double u, double v, double w, SPoint3 &p);
+
virtual void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
@@ -405,12 +411,19 @@ class MLine3 : public MLine {
};
_getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n);
}
+
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MLine::_getEdgeVertices(v);
v[2] = _vs[0];
}
+
+ virtual void pnt(double u, double v, double w, SPoint3 &p);
+
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s);
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
+
virtual int getTypeForMSH(){ return MSH_LIN_3; }
virtual int getTypeForUNV(){ return 24; } // parabolic beam
virtual int getTypeForVTK(){ return 21; }
@@ -459,6 +472,12 @@ class MLineN : public MLine {
MLine::_getEdgeVertices(v);
for(unsigned int i = 0; i != _vs.size(); ++i) v[i+2] = _vs[i];
}
+
+ virtual void pnt(double u, double v, double w, SPoint3 &p) ;
+
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s);
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
+
virtual int getTypeForMSH()
{
if(_vs.size() == 2) return MSH_LIN_4;
@@ -1435,6 +1454,11 @@ class MTetrahedron10 : public MTetrahedron {
tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
tmp = _vs[5]; _vs[5] = _vs[3]; _vs[3] = tmp;
}
+
+
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s);
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
+
virtual void jac(double u,double v,double w, double jac[3][3])
{
MTetrahedron::jac(2,_vs,u,v,w,jac);
@@ -1607,7 +1631,7 @@ class MTetrahedronN : public MTetrahedron {
case MSH_TET_34:
{
std::vector<MVertex*> inv(30);
- for (int i=0;i<31;i++) inv[i] = _vs[reverseTet34[i+4]-4];
+ for (int i=0;i<30;i++) inv[i] = _vs[reverseTet34[i+4]-4];
_vs = inv;
break;
}
@@ -1620,6 +1644,9 @@ class MTetrahedronN : public MTetrahedron {
}
}
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s);
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
+
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumEdgesRep();
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index adc8fea423..272eb91a8e 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -668,7 +668,13 @@ void reorientTrianglePoints(std::vector<MVertex*>& vtcs,int orientation,bool swa
if (swap) {
// --- interior "principal vertices"
- for (int i=0;i<2;i++) for (int j=0;j<2-i;j++) vtcs[i*2+j] = tmp[i+j*2];
+
+
+ vtcs[orientation] = tmp[orientation];
+ vtcs[(orientation + 1) % 3] = tmp[(orientation+2)%3];
+ vtcs[(orientation + 2) % 3] = tmp[(orientation+1)%3];
+
+ // for (int i=0;i<2;i++) for (int j=0;j<2-i;j++) vtcs[i*2+j] = tmp[i+j*2];
}
// no swap
@@ -711,12 +717,6 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- // std::vector<MVertex*> p;
-// face.getOrderedVertices(p);
-// if(faceVertices.count(p)){
-// // checking orientation not possible ???
-// vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
-// }
faceContainer::iterator fIter = faceVertices.find(face);
@@ -756,7 +756,10 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
else hoEdgeNodes.insert(hoEdgeNodes.end(),eIter->second.rbegin(),eIter->second.rend());
}
- MTriangleN incomplete(face.getVertex(0),face.getVertex(1),face.getVertex(2),hoEdgeNodes,nPts+1);
+ MTriangleN incomplete(face.getVertex(0),
+ face.getVertex(1),
+ face.getVertex(2),
+ hoEdgeNodes,nPts+1);
double X(0),Y(0),Z(0);
@@ -765,18 +768,22 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
double t1 = points(k,0);
double t2 = points(k,1);
-
- for (int j=0; j<incomplete.getNumVertices(); j++){
+ SPoint3 pos;
+ incomplete.pnt(t1,t2,0,pos);
+ MVertex* v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
+
+
+ // for (int j=0; j<incomplete.getNumVertices(); j++){
- double sf ; incomplete.getShapeFunction(j,t1,t2,0,sf);
- MVertex *vt = incomplete.getVertex(j);
+// double sf ; incomplete.getShapeFunction(j,t1,t2,0,sf);
+// MVertex *vt = incomplete.getVertex(j);
- X += sf * vt->x();
- Y += sf * vt->y();
- Z += sf * vt->z();
- }
+// X += sf * vt->x();
+// Y += sf * vt->y();
+// Z += sf * vt->z();
+// }
- MVertex* v = new MVertex(X,Y,Z,gr);
+// MVertex* v = new MVertex(X,Y,Z,gr);
// SPoint3 pc = face.interpolate(t1, t2);
@@ -861,6 +868,10 @@ void getRegionVertices(GRegion *gr,
const double t2 = points(k, 1);
const double t3 = points(k, 2);
+ SPoint3 pos;
+ incomplete->pnt(t1,t2,t3,pos);
+ v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
+
// FIXME: KOEN - I had to comment this out (MElement does not have
// pnt() member) -- CG
@@ -868,16 +879,16 @@ void getRegionVertices(GRegion *gr,
// incomplete->pnt(t1,t2,t3,pos);
// v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
- double X(0),Y(0),Z(0);
- for (int j=0; j<incomplete->getNumVertices(); j++){
- double sf ; incomplete->getShapeFunction(j,t1,t2,t3,sf);
- MVertex *vt = incomplete->getVertex(j);
- X += sf * vt->x();
- Y += sf * vt->y();
- Z += sf * vt->z();
- }
- v = new MVertex(X,Y,Z, gr);
-
+ // double X(0),Y(0),Z(0);
+ // for (int j=0; j<incomplete->getNumVertices(); j++){
+ // double sf ; incomplete->getShapeFunction(j,t1,t2,t3,sf);
+ // MVertex *vt = incomplete->getVertex(j);
+ // X += sf * vt->x();
+ // Y += sf * vt->y();
+ // Z += sf * vt->z();
+ // }
+ // v = new MVertex(X,Y,Z, gr);
+
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
@@ -982,8 +993,8 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve.insert(ve.end(), vf.begin(), vf.end());
MTetrahedronN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
ve, nPts + 1);
- // getRegionVertices(gr, &incpl, t, vr, linear, nPts);
- // ve.insert(ve.end(), vr.begin(), vr.end());
+ getRegionVertices(gr, &incpl, t, vr, linear, nPts);
+ ve.insert(ve.end(), vr.begin(), vr.end());
tetrahedra2.push_back
(new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
ve, nPts + 1));
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 40c6ce22d8..76640414c3 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -7,6 +7,25 @@
#include "GmshDefines.h"
#include "Message.h"
+Double_Matrix generate1DMonomials(int order) {
+
+ Double_Matrix monomials(order+1,1);
+ for (int i=0;i<order+1;i++) monomials(i,0) = i;
+
+}
+
+Double_Matrix generate1DPoints(int order) {
+
+ Double_Matrix line(order+1,1);
+
+ line(0,0) = -1.;
+ line(1,0) = 1.;
+
+ double dd = 2./order;
+
+ for (int i=2;i<order;i++) line(i,0) = -1. + dd * (i-1);
+}
+
Double_Matrix generatePascalTriangle(int order)
{
Double_Matrix monomials((order + 1) * (order + 2) / 2, 2);
@@ -567,6 +586,26 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
gmshFunctionSpace F;
switch (tag){
+ case MSH_LIN_2 :
+ F.monomials = generate1DMonomials(1);
+ F.points = generate1DPoints(1);
+ break;
+ case MSH_LIN_3 :
+ F.monomials = generate1DMonomials(2);
+ F.points = generate1DPoints(2);
+ break;
+ case MSH_LIN_4:
+ F.monomials = generate1DMonomials(3);
+ F.points = generate1DPoints(3);
+ break;
+ case MSH_LIN_5:
+ F.monomials = generate1DMonomials(4);
+ F.points = generate1DPoints(4);
+ break;
+ case MSH_LIN_6:
+ F.monomials = generate1DMonomials(5);
+ F.points = generate1DPoints(5);
+ break;
case MSH_TRI_3 :
F.monomials = generatePascalTriangle(1);
F.points = gmshGeneratePointsTriangle(1, false);
--
GitLab