From 3b95283bd3a2c898797395a09be733c946b252ef Mon Sep 17 00:00:00 2001
From: Koen Hillewaert <koen.hillewaert@cenaero.be>
Date: Mon, 25 Aug 2008 13:30:56 +0000
Subject: [PATCH] integrated higher order mappings for triangles corrected bug
for higher order mappings in post views (supports q up to 4)
---
Geo/MElement.cpp | 67 ++++++++++++++++++++++++++++++++++++++++++++-
Geo/MElement.h | 41 +++++++++++++++------------
Parser/Gmsh.tab.cpp | 2 +-
Parser/Gmsh.y | 4 +--
4 files changed, 92 insertions(+), 22 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 125af1ae6d..e500de751c 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -719,7 +719,72 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
if(name) *name = "Unknown";
return 0;
}
-}
+}
+
+void MTriangle::getShapeFunction(int num,double uu, double vv,double ww,double& s) {
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+
+ int nf = getNumFaceVertices();
+ double fv[256];
+
+ if (!nf){
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, fv); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, fv); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, fv); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, fv); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, fv); break;
+ default: Msg::Error("Order %d triangle shape function not implemented", getPolynomialOrder()); break;
+ }
+ }
+ else{
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, fv); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, fv); break;
+ 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;
+ }
+ }
+ s = fv[num];
+#endif
+
+}
+
+void MTriangle::getGradShapeFunction(int num,double uu,double vv,double ww,double s[3]) {
+
+
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double grads[256][2];
+ int nf = getNumFaceVertices();
+
+ if (!nf){
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
+ 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;
+ }
+ }
+ else{
+ switch(getPolynomialOrder()){
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
+ 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;
+ }
+ }
+
+ for (int i=0;i<3;i++) s[i] = grads[num][i];
+
+#endif
+}
void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double ww,
double j[2][3])
diff --git a/Geo/MElement.h b/Geo/MElement.h
index fc96c63628..806cba53e7 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -566,24 +566,28 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s)
- {
- switch(num){
- case 0 : s = 1. - u - v; break;
- case 1 : s = u ; break;
- case 2 : s = v; break;
- default : s = 0.; break;
- }
- }
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
- {
- switch(num) {
- case 0 : s[0] = -1.; s[1] = -1.; s[2] = 0.; break;
- case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
- case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
- default : s[0] = s[1] = s[2] = 0.; break;
- }
- }
+
+ 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 getShapeFunction(int num, double u, double v, double w, double &s) */
+/* { */
+/* switch(num){ */
+/* case 0 : s = 1. - u - v; break; */
+/* case 1 : s = u ; break; */
+/* case 2 : s = v; break; */
+/* default : s = 0.; break; */
+/* } */
+/* } */
+/* virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]) */
+/* { */
+/* switch(num) { */
+/* case 0 : s[0] = -1.; s[1] = -1.; s[2] = 0.; break; */
+/* case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break; */
+/* case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break; */
+/* default : s[0] = s[1] = s[2] = 0.; break; */
+/* } */
+/* } */
virtual bool isInside(double u, double v, double w, double tol=1.e-8)
{
if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
@@ -704,6 +708,7 @@ class MTriangle6 : public MTriangle {
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
+
virtual void jac(double u, double v, double w, double j[2][3])
{
MTriangle::jac(2, _vs, u, v, w, j);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 3d2a661fe4..9a22a17e6b 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -376,7 +376,7 @@ static PViewDataList *ViewData;
static ExtrudeParams extr;
static gmshSurface *myGmshSurface = 0;
static List_T *ViewValueList = 0;
-static double ViewCoord[100];
+static double ViewCoord[105];
static int *ViewNumList = 0, ViewCoordIdx = 0;
#define MAX_RECUR_LOOPS 100
static int ImbricatedLoop = 0;
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index cb013f861f..3806432de1 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -51,7 +51,7 @@ static PViewDataList *ViewData;
static ExtrudeParams extr;
static gmshSurface *myGmshSurface = 0;
static List_T *ViewValueList = 0;
-static double ViewCoord[100];
+static double ViewCoord[105]; // KH: support up to order 4 mappings
static int *ViewNumList = 0, ViewCoordIdx = 0;
#define MAX_RECUR_LOOPS 100
static int ImbricatedLoop = 0;
@@ -453,7 +453,7 @@ Element :
#if !defined(HAVE_NO_POST)
if(ViewValueList){
for(int i = 0; i < 3; i++)
- for(int j = 0; j < ViewCoordIdx / 3; j++)
+ for(int j = 0; j < ViewCoordIdx / 3; j++)
List_Add(ViewValueList, &ViewCoord[3 * j + i]);
}
#endif
--
GitLab