diff --git a/Plugin/CutGrid.cpp b/Plugin/CutGrid.cpp
index c600a59460d2a63668828815b245e572cb2ccf08..f1b60402c21751e1552732fe0ab7d7b2a8c44a8c 100644
--- a/Plugin/CutGrid.cpp
+++ b/Plugin/CutGrid.cpp
@@ -1,4 +1,4 @@
-// $Id: CutGrid.cpp,v 1.15 2005-01-09 02:18:59 geuzaine Exp $
+// $Id: CutGrid.cpp,v 1.16 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -192,17 +192,16 @@ void GMSH_CutGridPlugin::getInfos(char *author, char *copyright,
strcpy(author, "J.-F. Remacle (remacle@scorec.rpi.edu)");
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
- "Plugin(CutGrid) cuts a triangle/tetrahedron view\n"
- "with a rectangular grid defined by the 3 points\n"
+ "Plugin(CutGrid) cuts the view `iView' with a\n"
+ "rectangular grid defined by the 3 points\n"
"(`X0',`Y0',`Z0') (origin), (`X1',`Y1',`Z1') (axis of U)\n"
"and (`X2',`Y2',`Z2') (axis of V). The number of points\n"
"along U and V is set with the options `nPointsU'\n"
"and `nPointsV'. If `ConnectPoints' is zero, the\n"
- "plugin creates scalar or vector points; otherwise,\n"
- "the plugin generates scalar or vector quadrangles,\n"
- "lines or points depending on the values of `nPointsU'\n"
- "and `nPointsV'. If `iView' < 0, the plugin is run on\n"
- "the current view.\n"
+ "plugin creates points; otherwise, the plugin\n"
+ "generates quadrangles, lines or points depending\n"
+ "on the values of `nPointsU' and `nPointsV'. If\n"
+ "`iView' < 0, the plugin is run on the current view.\n"
"\n"
"Plugin(CutGrid) creates one new view.\n");
}
@@ -247,7 +246,85 @@ void GMSH_CutGridPlugin::getPoint(int iU, int iV, double *X)
v * (CutGridOptions_Number[8].def-CutGridOptions_Number[2].def) ;
}
-Post_View * GMSH_CutGridPlugin::GenerateView(Post_View * v, int connect) const
+void GMSH_CutGridPlugin::addInView(Post_View *v, int connect, int nbcomp,
+ double ***pnts, double ***vals,
+ List_T *P, int *nP,
+ List_T *L, int *nL,
+ List_T *Q, int *nQ)
+{
+ if(!connect || (getNbU() == 1 && getNbV() == 1)){ // generate points
+
+ for(int i = 0; i < getNbU(); ++i){
+ for(int j = 0; j < getNbV(); ++j){
+ List_Add(P, &pnts[i][j][0]);
+ List_Add(P, &pnts[i][j][1]);
+ List_Add(P, &pnts[i][j][2]);
+ (*nP)++;
+ for(int k = 0; k < v->NbTimeStep; ++k){
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(P, &vals[i][j][nbcomp*k+l]);
+ }
+ }
+ }
+
+ }
+ else{ // generate lines or quads
+
+ if(getNbU() == 1){
+ for(int i = 0; i < getNbV()-1; ++i){
+ List_Add(L, &pnts[0][i][0]); List_Add(L, &pnts[0][i+1][0]);
+ List_Add(L, &pnts[0][i][1]); List_Add(L, &pnts[0][i+1][1]);
+ List_Add(L, &pnts[0][i][2]); List_Add(L, &pnts[0][i+1][2]);
+ (*nL)++;
+ for(int k = 0; k < v->NbTimeStep; ++k){
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &vals[0][i ][nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &vals[0][i+1][nbcomp*k+l]);
+ }
+ }
+ }
+ else if(getNbV() == 1){
+ for(int i = 0; i < getNbU()-1; ++i){
+ List_Add(L, &pnts[i][0][0]); List_Add(L, &pnts[i+1][0][0]);
+ List_Add(L, &pnts[i][0][1]); List_Add(L, &pnts[i+1][0][1]);
+ List_Add(L, &pnts[i][0][2]); List_Add(L, &pnts[i+1][0][2]);
+ (*nL)++;
+ for(int k = 0; k < v->NbTimeStep; ++k){
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &vals[i ][0][nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &vals[i+1][0][nbcomp*k+l]);
+ }
+ }
+ }
+ else{
+ for(int i = 0; i < getNbU()-1; ++i){
+ for(int j = 0; j < getNbV()-1; ++j){
+ List_Add(Q, &pnts[i ][j ][0]); List_Add(Q, &pnts[i+1][j ][0]);
+ List_Add(Q, &pnts[i+1][j+1][0]); List_Add(Q, &pnts[i ][j+1][0]);
+ List_Add(Q, &pnts[i ][j ][1]); List_Add(Q, &pnts[i+1][j ][1]);
+ List_Add(Q, &pnts[i+1][j+1][1]); List_Add(Q, &pnts[i ][j+1][1]);
+ List_Add(Q, &pnts[i ][j ][2]); List_Add(Q, &pnts[i+1][j ][2]);
+ List_Add(Q, &pnts[i+1][j+1][2]); List_Add(Q, &pnts[i ][j+1][2]);
+ (*nQ)++;
+ for(int k = 0; k < v->NbTimeStep; ++k){
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(Q, &vals[i ][j ][nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(Q, &vals[i+1][j ][nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(Q, &vals[i+1][j+1][nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(Q, &vals[i ][j+1][nbcomp*k+l]);
+ }
+ }
+ }
+ }
+ }
+}
+
+Post_View * GMSH_CutGridPlugin::GenerateView(Post_View * v, int connect)
{
if(getNbU() <= 0 || getNbV() <= 0)
return v;
@@ -256,6 +333,17 @@ Post_View * GMSH_CutGridPlugin::GenerateView(Post_View * v, int connect) const
OctreePost o(v);
+ int nbs = 0, nbv = 0, nbt = 0;
+
+ if(v->NbST || v->NbSQ || v->NbSS || v->NbSH || v->NbSI || v->NbSY)
+ nbs = 1;
+ if(v->NbVT || v->NbVQ || v->NbVS || v->NbVH || v->NbVI || v->NbVY)
+ nbv = 1;
+ if(v->NbTT || v->NbTQ || v->NbTS || v->NbTH || v->NbTI || v->NbTY)
+ nbt = 1;
+
+ int maxcomp = nbt ? 9 : (nbv ? 3 : 1);
+
double ***pnts = new double** [getNbU()];
double ***vals = new double** [getNbU()];
for(int i = 0; i < getNbU(); i++){
@@ -263,190 +351,35 @@ Post_View * GMSH_CutGridPlugin::GenerateView(Post_View * v, int connect) const
vals[i] = new double* [getNbV()];
for(int j = 0; j < getNbV(); j++){
pnts[i][j] = new double[3];
- vals[i][j] = new double[3 * v->NbTimeStep]; // change to 9 for tensors
+ vals[i][j] = new double[maxcomp * v->NbTimeStep];
getPoint(i, j, pnts[i][j]);
}
}
- if(v->NbST || v->NbSS){
-
+ if(nbs){
for(int i = 0; i < getNbU(); i++)
for(int j = 0; j < getNbV(); j++)
- o.searchScalar(pnts[i][j][0], pnts[i][j][1], pnts[i][j][2],
- vals[i][j]);
-
- if(!connect){ // generate points
-
- for(int i = 0; i < getNbU(); ++i){
- for(int j = 0; j < getNbV(); ++j){
- List_Add(View->SP, &pnts[i][j][0]);
- List_Add(View->SP, &pnts[i][j][1]);
- List_Add(View->SP, &pnts[i][j][2]);
- View->NbSP ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->SP, &vals[i][j][k]);
- }
- }
- }
-
- }
- else{
-
- if(getNbU() == 1 && getNbV() == 1){
- List_Add(View->SP, &pnts[0][0][0]);
- List_Add(View->SP, &pnts[0][0][1]);
- List_Add(View->SP, &pnts[0][0][2]);
- View->NbSP ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->SP, &vals[0][0][k]);
- }
- }
- else if(getNbU() == 1){
- for(int i = 0; i < getNbV()-1; ++i){
- List_Add(View->SL, &pnts[0][i][0]); List_Add(View->SL, &pnts[0][i+1][0]);
- List_Add(View->SL, &pnts[0][i][1]); List_Add(View->SL, &pnts[0][i+1][1]);
- List_Add(View->SL, &pnts[0][i][2]); List_Add(View->SL, &pnts[0][i+1][2]);
- View->NbSL ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->SL, &vals[0][i ][k]);
- List_Add(View->SL, &vals[0][i+1][k]);
- }
- }
- }
- else if(getNbV() == 1){
- for(int i = 0; i < getNbU()-1; ++i){
- List_Add(View->SL, &pnts[i][0][0]); List_Add(View->SL, &pnts[i+1][0][0]);
- List_Add(View->SL, &pnts[i][0][1]); List_Add(View->SL, &pnts[i+1][0][1]);
- List_Add(View->SL, &pnts[i][0][2]); List_Add(View->SL, &pnts[i+1][0][2]);
- View->NbSL ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->SL, &vals[i ][0][k]);
- List_Add(View->SL, &vals[i+1][0][k]);
- }
- }
- }
- else{
- for(int i = 0; i < getNbU()-1; ++i){
- for(int j = 0; j < getNbV()-1; ++j){
- List_Add(View->SQ, &pnts[i ][j ][0]); List_Add(View->SQ, &pnts[i+1][j ][0]);
- List_Add(View->SQ, &pnts[i+1][j+1][0]); List_Add(View->SQ, &pnts[i ][j+1][0]);
- List_Add(View->SQ, &pnts[i ][j ][1]); List_Add(View->SQ, &pnts[i+1][j ][1]);
- List_Add(View->SQ, &pnts[i+1][j+1][1]); List_Add(View->SQ, &pnts[i ][j+1][1]);
- List_Add(View->SQ, &pnts[i ][j ][2]); List_Add(View->SQ, &pnts[i+1][j ][2]);
- List_Add(View->SQ, &pnts[i+1][j+1][2]); List_Add(View->SQ, &pnts[i ][j+1][2]);
- View->NbSQ ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->SQ, &vals[i ][j ][k]);
- List_Add(View->SQ, &vals[i+1][j ][k]);
- List_Add(View->SQ, &vals[i+1][j+1][k]);
- List_Add(View->SQ, &vals[i ][j+1][k]);
- }
- }
- }
- }
- }
-
+ o.searchScalar(pnts[i][j][0], pnts[i][j][1], pnts[i][j][2], vals[i][j]);
+ addInView(v, connect, 1, pnts, vals,
+ View->SP, &View->NbSP, View->SL, &View->NbSL, View->SQ, &View->NbSQ);
}
- if(v->NbVT || v->NbVS){
-
+ if(nbv){
double sizeElem;
for(int i = 0; i < getNbU(); i++)
for(int j = 0; j < getNbV(); j++)
- o.searchVector(pnts[i][j][0], pnts[i][j][1], pnts[i][j][2],
- vals[i][j], &sizeElem);
-
- if(!connect){ // generate points
-
- for(int i = 0; i < getNbU(); ++i){
- for(int j = 0; j < getNbV(); ++j){
- List_Add(View->VP, &pnts[i][j][0]);
- List_Add(View->VP, &pnts[i][j][1]);
- List_Add(View->VP, &pnts[i][j][2]);
- View->NbVP ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->VP, &vals[i][j][3*k]);
- List_Add(View->VP, &vals[i][j][3*k+1]);
- List_Add(View->VP, &vals[i][j][3*k+2]);
- }
- }
- }
-
- }
- else{
-
- if(getNbU() == 1 && getNbV() == 1){
- List_Add(View->VP, &pnts[0][0][0]);
- List_Add(View->VP, &pnts[0][0][1]);
- List_Add(View->VP, &pnts[0][0][2]);
- View->NbVP ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->VP, &vals[0][0][3*k]);
- List_Add(View->VP, &vals[0][0][3*k+1]);
- List_Add(View->VP, &vals[0][0][3*k+2]);
- }
- }
- else if(getNbU() == 1){
- for(int i = 0; i < getNbV()-1; ++i){
- List_Add(View->VL, &pnts[0][i][0]); List_Add(View->VL, &pnts[0][i+1][0]);
- List_Add(View->VL, &pnts[0][i][1]); List_Add(View->VL, &pnts[0][i+1][1]);
- List_Add(View->VL, &pnts[0][i][2]); List_Add(View->VL, &pnts[0][i+1][2]);
- View->NbVL ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->VL, &vals[0][i ][3*k]);
- List_Add(View->VL, &vals[0][i ][3*k+1]);
- List_Add(View->VL, &vals[0][i ][3*k+2]);
- List_Add(View->VL, &vals[0][i+1][3*k]);
- List_Add(View->VL, &vals[0][i+1][3*k+1]);
- List_Add(View->VL, &vals[0][i+1][3*k+2]);
- }
- }
- }
- else if(getNbV() == 1){
- for(int i = 0; i < getNbU()-1; ++i){
- List_Add(View->VL, &pnts[i][0][0]); List_Add(View->VL, &pnts[i+1][0][0]);
- List_Add(View->VL, &pnts[i][0][1]); List_Add(View->VL, &pnts[i+1][0][1]);
- List_Add(View->VL, &pnts[i][0][2]); List_Add(View->VL, &pnts[i+1][0][2]);
- View->NbVL ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->VL, &vals[i ][0][3*k]);
- List_Add(View->VL, &vals[i ][0][3*k+1]);
- List_Add(View->VL, &vals[i ][0][3*k+2]);
- List_Add(View->VL, &vals[i+1][0][3*k]);
- List_Add(View->VL, &vals[i+1][0][3*k+1]);
- List_Add(View->VL, &vals[i+1][0][3*k+2]);
- }
- }
- }
- else{
- for(int i = 0; i < getNbU()-1; ++i){
- for(int j = 0; j < getNbV()-1; ++j){
- List_Add(View->VQ, &pnts[i ][j ][0]); List_Add(View->VQ, &pnts[i+1][j ][0]);
- List_Add(View->VQ, &pnts[i+1][j+1][0]); List_Add(View->VQ, &pnts[i ][j+1][0]);
- List_Add(View->VQ, &pnts[i ][j ][1]); List_Add(View->VQ, &pnts[i+1][j ][1]);
- List_Add(View->VQ, &pnts[i+1][j+1][1]); List_Add(View->VQ, &pnts[i ][j+1][1]);
- List_Add(View->VQ, &pnts[i ][j ][2]); List_Add(View->VQ, &pnts[i+1][j ][2]);
- List_Add(View->VQ, &pnts[i+1][j+1][2]); List_Add(View->VQ, &pnts[i ][j+1][2]);
- View->NbVQ ++;
- for(int k = 0; k < v->NbTimeStep; ++k){
- List_Add(View->VQ, &vals[i ][j ][3*k]);
- List_Add(View->VQ, &vals[i ][j ][3*k+1]);
- List_Add(View->VQ, &vals[i ][j ][3*k+2]);
- List_Add(View->VQ, &vals[i+1][j ][3*k]);
- List_Add(View->VQ, &vals[i+1][j ][3*k+1]);
- List_Add(View->VQ, &vals[i+1][j ][3*k+2]);
- List_Add(View->VQ, &vals[i+1][j+1][3*k]);
- List_Add(View->VQ, &vals[i+1][j+1][3*k+1]);
- List_Add(View->VQ, &vals[i+1][j+1][3*k+2]);
- List_Add(View->VQ, &vals[i ][j+1][3*k]);
- List_Add(View->VQ, &vals[i ][j+1][3*k+1]);
- List_Add(View->VQ, &vals[i ][j+1][3*k+2]);
- }
- }
- }
- }
- }
+ o.searchVector(pnts[i][j][0], pnts[i][j][1], pnts[i][j][2], vals[i][j],
+ &sizeElem);
+ addInView(v, connect, 3, pnts, vals,
+ View->VP, &View->NbVP, View->VL, &View->NbVL, View->VQ, &View->NbVQ);
+ }
+ if(nbt){
+ for(int i = 0; i < getNbU(); i++)
+ for(int j = 0; j < getNbV(); j++)
+ o.searchTensor(pnts[i][j][0], pnts[i][j][1], pnts[i][j][2], vals[i][j]);
+ addInView(v, connect, 9, pnts, vals,
+ View->TP, &View->NbTP, View->TL, &View->NbTL, View->TQ, &View->NbTQ);
}
for(int i = 0; i < getNbU(); i++){
diff --git a/Plugin/CutGrid.h b/Plugin/CutGrid.h
index 2f9b90721391c3dab21e0c645a720d8e53b19112..880a5f55bfe24d515d46a6ef2b6e5892ef387071 100644
--- a/Plugin/CutGrid.h
+++ b/Plugin/CutGrid.h
@@ -31,6 +31,12 @@ class GMSH_CutGridPlugin : public GMSH_Post_Plugin
{
static double callback(int num, int action, double value, double *opt,
double step, double min, double max);
+ void addInView(Post_View *v, int connect, int nbcomp,
+ double ***pnts, double ***vals,
+ List_T *P, int *nP,
+ List_T *L, int *nL,
+ List_T *Q, int *nQ);
+ Post_View * GenerateView (Post_View * v, int connectPoints);
public:
GMSH_CutGridPlugin();
void getName (char *name) const;
@@ -41,7 +47,6 @@ public:
int getNbOptions() const;
StringXNumber *getOption (int iopt);
Post_View *execute (Post_View *);
- virtual Post_View * GenerateView (Post_View * v, int connectPoints) const ;
static int getNbU ();
static int getNbV ();
diff --git a/Plugin/CutParametric.cpp b/Plugin/CutParametric.cpp
index dfd618da718a9b6c36b883a44c78bff8a5e0cfc2..c0f9935339da98514feec1c0330a7d1b61068cb9 100644
--- a/Plugin/CutParametric.cpp
+++ b/Plugin/CutParametric.cpp
@@ -1,4 +1,4 @@
-// $Id: CutParametric.cpp,v 1.10 2005-01-09 02:18:59 geuzaine Exp $
+// $Id: CutParametric.cpp,v 1.11 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -76,14 +76,13 @@ void GMSH_CutParametricPlugin::getInfos(char *author, char *copyright,
strcpy(author, "C. Geuzaine (geuzaine@acm.caltech.edu)");
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
- "Plugin(CutParametric) cuts a scalar triangle/\n"
- "tetrahedron view `iView' with the parametric function\n"
- "(`X'(u), `Y'(u), `Z'(u)), using `nPointsU' values of\n"
- "the parameter u in [`MinU', `MaxU']. If\n"
- "`ConnectPoints' is set, the plugin creates scalar\n"
- "line elements; otherwise, the plugin generates\n"
- "scalar points. If `iView' < 0, the plugin is run on\n"
- "the current view.\n"
+ "Plugin(CutParametric) cuts the scalar view `iView'\n"
+ "with the parametric function (`X'(u), `Y'(u), `Z'(u)),\n"
+ "using `nPointsU' values of the parameter u in\n"
+ "[`MinU', `MaxU']. If `ConnectPoints' is set, the\n"
+ "plugin creates scalar line elements; otherwise,\n"
+ "the plugin generates scalar points. If `iView' < 0,\n"
+ "the plugin is run on the current view.\n"
"\n"
"Plugin(CutParametric) creates one new view.\n");
}
@@ -113,6 +112,36 @@ void GMSH_CutParametricPlugin::catchErrorMessage(char *errorMessage) const
strcpy(errorMessage, "CutParametric failed...");
}
+static void addInView(int connect, int i, int nbcomp, int nbtime,
+ double x0, double y0, double z0, double *res0,
+ double x, double y, double z, double *res,
+ List_T *P, int *nP, List_T *L, int *nL)
+{
+ if(connect){
+ if(i){
+ List_Add(L, &x0); List_Add(L, &x);
+ List_Add(L, &y0); List_Add(L, &y);
+ List_Add(L, &z0); List_Add(L, &z);
+ for(int k = 0; k < nbtime; ++k){
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &res0[nbcomp*k+l]);
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(L, &res[nbcomp*k+l]);
+ }
+ (*nL)++;
+ }
+ }
+ else{
+ List_Add(P, &x);
+ List_Add(P, &y);
+ List_Add(P, &z);
+ for(int k = 0; k < nbtime; ++k)
+ for(int l = 0; l < nbcomp; ++l)
+ List_Add(P, &res[nbcomp*k+l]);
+ (*nP)++;
+ }
+}
+
Post_View *GMSH_CutParametricPlugin::execute(Post_View * v)
{
int iView = (int)CutParametricOptions_Number[4].def;
@@ -167,45 +196,46 @@ Post_View *GMSH_CutParametricPlugin::execute(Post_View * v)
OctreePost o(v1);
Post_View *v2 = BeginView(1);
- double *res0 = new double[v1->NbTimeStep];
- double *res = new double[v1->NbTimeStep];
- double x, y, z, x0, y0, z0;
+ double *res0 = new double[9*v1->NbTimeStep];
+ double *res = new double[9*v1->NbTimeStep];
+ double x = 0., y = 0., z = 0., x0 = 0., y0 = 0., z0 = 0.;
+
+ for(int k = 0; k < 9*v1->NbTimeStep; ++k) res0[k] = res[k] = 0.;
+
for(int i = 0; i < nbU; ++i){
if(i && connect){
x0 = x;
y0 = y;
z0 = z;
- for(int k = 0; k < v1->NbTimeStep; ++k) res0[k] = res[k];
+ for(int k = 0; k < 9*v1->NbTimeStep; ++k) res0[k] = res[k];
}
- double u = minU + (double)(i)/(double)(nbU-1) * (maxU - minU);
+
+ double u;
+ if(nbU == 1 || maxU == minU)
+ u = minU;
+ else
+ u = minU + (double)(i)/(double)(nbU-1) * (maxU - minU);
char *names[] = { "u" };
double values[] = { u };
x = evaluator_evaluate(fx, sizeof(names)/sizeof(names[0]), names, values);
y = evaluator_evaluate(fy, sizeof(names)/sizeof(names[0]), names, values);
z = evaluator_evaluate(fz, sizeof(names)/sizeof(names[0]), names, values);
- o.searchScalar(x, y, z, res);
- if(connect){
- if(i){
- List_Add(v2->SL, &x0);
- List_Add(v2->SL, &x);
- List_Add(v2->SL, &y0);
- List_Add(v2->SL, &y);
- List_Add(v2->SL, &z0);
- List_Add(v2->SL, &z);
- for(int k = 0; k < v1->NbTimeStep; ++k){
- List_Add(v2->SL, &res0[k]);
- List_Add(v2->SL, &res[k]);
- }
- v2->NbSL++;
- }
+
+ if(v->NbST || v->NbSQ || v->NbSS || v->NbSH || v->NbSI || v->NbSY){
+ o.searchScalar(x, y, z, res);
+ addInView(connect, i, 1, v1->NbTimeStep, x0, y0, z0, res0, x, y, z, res,
+ v2->SP, &v2->NbSP, v2->SL, &v2->NbSL);
+ }
+ if(v->NbVT || v->NbVQ || v->NbVS || v->NbVH || v->NbVI || v->NbVY){
+ double size;
+ o.searchVector(x, y, z, res, &size);
+ addInView(connect, i, 3, v1->NbTimeStep, x0, y0, z0, res0, x, y, z, res,
+ v2->VP, &v2->NbVP, v2->VL, &v2->NbVL);
}
- else{
- List_Add(v2->SP, &x);
- List_Add(v2->SP, &y);
- List_Add(v2->SP, &z);
- for(int k = 0; k < v1->NbTimeStep; ++k)
- List_Add(v2->SP, &res[k]);
- v2->NbSP++;
+ if(v->NbTT || v->NbTQ || v->NbTS || v->NbTH || v->NbTI || v->NbTY){
+ o.searchTensor(x, y, z, res);
+ addInView(connect, i, 9, v1->NbTimeStep, x0, y0, z0, res0, x, y, z, res,
+ v2->TP, &v2->NbTP, v2->TL, &v2->NbTL);
}
}
diff --git a/Plugin/Evaluate.cpp b/Plugin/Evaluate.cpp
index 81e43046efff5f8c71bcaa00a911b9ff1501aa47..fcd08ebdbe1bb643fcb1596d65fcd74f69163b95 100644
--- a/Plugin/Evaluate.cpp
+++ b/Plugin/Evaluate.cpp
@@ -1,4 +1,4 @@
-// $Id: Evaluate.cpp,v 1.16 2005-01-12 19:38:01 geuzaine Exp $
+// $Id: Evaluate.cpp,v 1.17 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -77,8 +77,8 @@ void GMSH_EvaluatePlugin::getInfos(char *author, char *copyright,
"- the usual mathematical functions (Log, Sqrt\n"
"Sin, Cos, Fabs, ...) and operators (+, -, *, /, ^);\n"
"\n"
- "- the symbols x, y and z, to retrieve the coordinates\n"
- "of the current node;\n"
+ "- the symbols x, y and z, to retrieve the\n"
+ "coordinates of the current node;\n"
"\n"
"- the symbols Time and TimeStep, to retrieve the\n"
"current time and time step values;\n"
@@ -93,11 +93,11 @@ void GMSH_EvaluatePlugin::getInfos(char *author, char *copyright,
"\n"
"- the symbol w, to retrieve the `Component'-th\n"
"component of the field in `ExternalView' at the\n"
- "`ExternalTimeStep'-th time step. `ExternalView' and\n"
- "`iView' must be of the same type (scalar, vector\n"
- "or tensor); if `ExternalView' and `iView' are not\n"
- "based on the same spatial grid, `ExternalView' is\n"
- "interpolated onto `iView';\n"
+ "`ExternalTimeStep'-th time step. `ExternalView'\n"
+ "and `iView' must be of the same type (scalar,\n"
+ "vector or tensor); if `ExternalView' and `iView'\n"
+ "are not based on the same spatial grid,\n"
+ "`ExternalView' is interpolated onto `iView';\n"
"\n"
"- the symbols w0, w1, w2, ..., w8, to retrieve each\n"
"component of the field in `ExternalView' at the\n"
@@ -185,8 +185,10 @@ void GMSH_EvaluatePlugin::evaluate(Post_View *v1, List_T *list1, int nbElm1,
val2 = tmp;
if(nbComp == 1)
_octree->searchScalar(x[j], y[j], z[j], val2, timeStep2);
- else
+ else if(nbComp == 3)
_octree->searchVector(x[j], y[j], z[j], val2, &sizeElm, timeStep2);
+ else
+ _octree->searchTensor(x[j], y[j], z[j], val2, timeStep2);
}
else
val2 = (double *)List_Pointer_Fast(list2,
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 2b3a4804eb7af48e83cc7186f49e26c3061823d7..f7df8d5847242455a3b0291866597de480c2ca5a 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.79 2005-01-12 19:10:41 geuzaine Exp $
+# $Id: Makefile,v 1.80 2005-03-02 07:49:41 geuzaine Exp $
#
# Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
#
@@ -147,7 +147,7 @@ OctreePost.o: OctreePost.cpp Octree.h OctreeInternals.h OctreePost.h \
../DataStr/List.h ../Common/Views.h ../Common/ColorTable.h \
../Common/VertexArray.h ../Common/SmoothNormals.h \
../Common/GmshMatrix.h ../Common/AdaptiveViews.h ../Numeric/Numeric.h \
- ../Common/Message.h
+ ../Common/Message.h ShapeFunctions.h
StreamLines.o: StreamLines.cpp OctreePost.h Octree.h OctreeInternals.h \
StreamLines.h Plugin.h ../Common/Options.h ../Common/Message.h \
../Common/Views.h ../Common/ColorTable.h ../DataStr/List.h \
diff --git a/Plugin/OctreePost.cpp b/Plugin/OctreePost.cpp
index 06f567b3924fe5fad7ba0428865049f5390ecefc..1022ee03ad247351f96a051539b693959dd7a9fb 100644
--- a/Plugin/OctreePost.cpp
+++ b/Plugin/OctreePost.cpp
@@ -1,4 +1,4 @@
-// $Id: OctreePost.cpp,v 1.13 2005-01-24 17:39:28 geuzaine Exp $
+// $Id: OctreePost.cpp,v 1.14 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -25,51 +25,10 @@
#include "Views.h"
#include "Numeric.h"
#include "Message.h"
+#include "ShapeFunctions.h"
-static double computeBarycentricTriangle(double *X, double *Y, double *Z,
- double *P, double *U)
-{
- // FIXME: only valid for triangles in the XY-plane
- double mat[2][2], b[2];
- U[2] = 0.0;
- mat[0][0] = X[1]-X[0];
- mat[0][1] = X[2]-X[0];
- mat[1][0] = Y[1]-Y[0];
- mat[1][1] = Y[2]-Y[0];
- b[0] = P[0] -X[0];
- b[1] = P[1] -Y[0];
-
- sys2x2(mat, b, U);
- return 0.5 * ( mat[0][0] * mat[1][1] - mat[1][0] *mat[0][1]);
-}
-
-static void computeBarycentricSimplex(double *X, double *Y, double *Z,
- double *P, double *U)
-{
- double mat[3][3], b[3];
-
- mat[0][0] = X[1]-X[0];
- mat[0][1] = X[2]-X[0];
- mat[0][2] = X[3]-X[0];
- mat[1][0] = Y[1]-Y[0];
- mat[1][1] = Y[2]-Y[0];
- mat[1][2] = Y[3]-Y[0];
- mat[2][0] = Z[1]-Z[0];
- mat[2][1] = Z[2]-Z[0];
- mat[2][2] = Z[3]-Z[0];
- b[0] = P[0] - X[0];
- b[1] = P[1] - Y[0];
- b[2] = P[2] - Z[0];
- double det;
- sys3x3(mat, b, U, &det);
-}
-
-static void minmax(int n,
- double *X,
- double *Y,
- double *Z,
- double *min,
- double *max)
+static void minmax(int n, double *X, double *Y, double *Z,
+ double *min, double *max)
{
min[0] = X[0];
min[1] = Y[0];
@@ -78,27 +37,23 @@ static void minmax(int n,
max[1] = Y[0];
max[2] = Z[0];
- for (int i = 1; i < n; ++i) {
- min[0] = (X[i] < min[0]) ? X[i]:min[0];
- min[1] = (Y[i] < min[1]) ? Y[i]:min[1];
- min[2] = (Z[i] < min[2]) ? Z[i]:min[2];
- max[0] = (X[i] > max[0]) ? X[i]:max[0];
- max[1] = (Y[i] > max[1]) ? Y[i]:max[1];
- max[2] = (Z[i] > max[2]) ? Z[i]:max[2];
+ for(int i = 1; i < n; i++) {
+ min[0] = (X[i] < min[0]) ? X[i] : min[0];
+ min[1] = (Y[i] < min[1]) ? Y[i] : min[1];
+ min[2] = (Z[i] < min[2]) ? Z[i] : min[2];
+ max[0] = (X[i] > max[0]) ? X[i] : max[0];
+ max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
+ max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
}
}
-static void centroid(int n,
- double *X,
- double *Y,
- double *Z,
- double *c)
+static void centroid(int n, double *X, double *Y, double *Z, double *c)
{
const double oc = 1./(double)n;
c[0] = X[0];
c[1] = Y[0];
c[2] = Z[0];
- for (int i = 1; i < n; ++i) {
+ for(int i = 1; i < n; i++) {
c[0] += X[i];
c[1] += Y[i];
c[2] += Z[i];
@@ -108,79 +63,127 @@ static void centroid(int n,
c[2] *= oc;
}
-void PostTriangleBB(void *a, double *min, double *max)
+void triBB(void *a, double *min, double *max)
{
- double *X = (double*) a;
- double *Y = &X[3];
- double *Z = &X[6];
-
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6];
minmax(3, X, Y, Z, min, max);
+}
+
+void quaBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ minmax(4, X, Y, Z, min, max);
+}
+
+void tetBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ minmax(4, X, Y, Z, min, max);
+}
+
+void hexBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16];
+ minmax(8, X, Y, Z, min, max);
+}
+
+void priBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12];
+ minmax(6, X, Y, Z, min, max);
+}
- // FIXME: only valid for triangles in the XY-plane
- min[2] = -1;
- max[2] = 1;
+void pyrBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10];
+ minmax(5, X, Y, Z, min, max);
+}
+
+int triInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6], uvw[3];
+ triangle tri(X, Y, Z);
+ tri.xyz2uvw(x, uvw);
+ return tri.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int quaInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
+ quadrangle qua(X, Y, Z);
+ qua.xyz2uvw(x, uvw);
+ return qua.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int tetInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
+ tetrahedron tet(X, Y, Z);
+ tet.xyz2uvw(x, uvw);
+ return tet.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int hexInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16], uvw[3];
+ hexahedron hex(X, Y, Z);
+ hex.xyz2uvw(x, uvw);
+ return hex.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int priInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12], uvw[3];
+ prism pri(X, Y, Z);
+ pri.xyz2uvw(x, uvw);
+ return pri.isInside(uvw[0], uvw[1], uvw[2]);
}
-int PostTriangleInEle(void *a, double *x)
+int pyrInEle(void *a, double *x)
{
- const double eps = 1.e-3;
- double U[3];
- double *X = (double*) a;
- double *Y = &X[3];
- double *Z = &X[6];
- computeBarycentricTriangle(X, Y, Z, x, U);
- double W = 1.-U[0]-U[1];
- if (U[0] < -eps || U[0] > 1+eps ||
- U[1] < -eps || U[1] > 1+eps ||
- W < -eps || W > 1+eps) return 0;
- return 1;
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10], uvw[3];
+ pyramid pyr(X, Y, Z);
+ pyr.xyz2uvw(x, uvw);
+ return pyr.isInside(uvw[0], uvw[1], uvw[2]);
}
-void PostTriangleCentroid(void *a, double *x)
+void triCentroid(void *a, double *x)
{
- double *X = (double*) a;
- double *Y = &X[3];
- double *Z = &X[6];
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6];
centroid(3, X, Y, Z, x);
}
-void PostSimplexBB(void *a, double *min, double *max)
+void quaCentroid(void *a, double *x)
{
- double *X = (double*) a;
- double *Y = &X[4];
- double *Z = &X[8];
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ centroid(4, X, Y, Z, x);
+}
- minmax(4, X, Y, Z, min, max);
+void tetCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ centroid(4, X, Y, Z, x);
+}
+
+void hexCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16];
+ centroid(8, X, Y, Z, x);
}
-int PostSimplexInEle(void *a, double *x)
+void priCentroid(void *a, double *x)
{
- const double eps = 1.e-5;
- double U[3];
- double *X = (double*) a;
- double *Y = &X[4];
- double *Z = &X[8];
- computeBarycentricSimplex(X, Y, Z, x, U);
- double W = 1.-U[0]-U[1]-U[2];
-
- if(U[0] < -eps || U[0] > 1+eps ||
- U[1] < -eps || U[1] > 1+eps ||
- U[2] < -eps || U[2] > 1+eps ||
- W < -eps || W > 1+eps) return 0;
- return 1;
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12];
+ centroid(6, X, Y, Z, x);
}
-void PostSimplexCentroid(void *a, double *x)
+void pyrCentroid(void *a, double *x)
{
- double *X = (double*) a;
- double *Y = &X[4];
- double *Z = &X[8];
- centroid (4, X, Y, Z, x);
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10];
+ centroid(5, X, Y, Z, x);
}
-static void addListOfStuff(Octree *o,
- List_T *l,
- int nbelm)
+static void addListOfStuff(Octree *o, List_T *l, int nbelm)
{
if(!l) return;
@@ -192,251 +195,178 @@ static void addListOfStuff(Octree *o,
OctreePost::~OctreePost()
{
- Octree_Delete(ST);
- Octree_Delete(VT);
- Octree_Delete(TT);
- Octree_Delete(SS);
- Octree_Delete(VS);
- Octree_Delete(TS);
+ Octree_Delete(ST); Octree_Delete(VT); Octree_Delete(TT);
+ Octree_Delete(SQ); Octree_Delete(VQ); Octree_Delete(TQ);
+ Octree_Delete(SS); Octree_Delete(VS); Octree_Delete(TS);
+ Octree_Delete(SH); Octree_Delete(VH); Octree_Delete(TH);
+ Octree_Delete(SI); Octree_Delete(VI); Octree_Delete(TI);
+ Octree_Delete(SY); Octree_Delete(VY); Octree_Delete(TY);
}
OctreePost::OctreePost(Post_View *v)
: theView(v)
{
- double min [3] = {v->BBox[0],v->BBox[2],v->BBox[4]};
+ double min [3] = {v->BBox[0], v->BBox[2], v->BBox[4]};
- double size[3] = {v->BBox[1]-v->BBox[0],
- v->BBox[3]-v->BBox[2],
- v->BBox[5]-v->BBox[4]};
+ double size[3] = {v->BBox[1] - v->BBox[0],
+ v->BBox[3] - v->BBox[2],
+ v->BBox[5] - v->BBox[4]};
const int maxElePerBucket = 100; // memory vs. speed trade-off
- ST = Octree_Create(maxElePerBucket, min, size,
- PostTriangleBB,
- PostTriangleCentroid,
- PostTriangleInEle);
+ ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
addListOfStuff(ST, v->ST, 9 + 3 * v->NbTimeStep);
-
- VT = Octree_Create(maxElePerBucket, min, size,
- PostTriangleBB,
- PostTriangleCentroid,
- PostTriangleInEle);
+ Octree_Arrange(ST);
+ VT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
addListOfStuff(VT, v->VT, 9 + 9 * v->NbTimeStep);
-
- TT = Octree_Create(maxElePerBucket, min, size,
- PostTriangleBB,
- PostTriangleCentroid,
- PostTriangleInEle);
+ Octree_Arrange(VT);
+ TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
addListOfStuff(TT, v->TT, 9 + 27 * v->NbTimeStep);
+ Octree_Arrange(TT);
- SS = Octree_Create(maxElePerBucket, min, size,
- PostSimplexBB,
- PostSimplexCentroid,
- PostSimplexInEle);
+ SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(SQ, v->SQ, 12 + 4 * v->NbTimeStep);
+ Octree_Arrange(SQ);
+ VQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(VQ, v->VQ, 12 + 12 * v->NbTimeStep);
+ Octree_Arrange(VQ);
+ TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(TQ, v->TQ, 12 + 36 * v->NbTimeStep);
+ Octree_Arrange(TQ);
+
+ SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
addListOfStuff(SS, v->SS, 12 + 4 * v->NbTimeStep);
-
- VS = Octree_Create(maxElePerBucket, min, size,
- PostSimplexBB,
- PostSimplexCentroid,
- PostSimplexInEle);
- addListOfStuff(VS, v->VS, 12 + 12 * v->NbTimeStep);
-
- TS = Octree_Create(maxElePerBucket, min, size,
- PostSimplexBB,
- PostSimplexCentroid,
- PostSimplexInEle);
- addListOfStuff(TS, v->TS, 12 + 36 * v->NbTimeStep);
-
- Octree_Arrange(ST);
- Octree_Arrange(VT);
- Octree_Arrange(TT);
Octree_Arrange(SS);
+ VS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(VS, v->VS, 12 + 12 * v->NbTimeStep);
Octree_Arrange(VS);
+ TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(TS, v->TS, 12 + 36 * v->NbTimeStep);
Octree_Arrange(TS);
+
+ SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(SH, v->SH, 24 + 8 * v->NbTimeStep);
+ Octree_Arrange(SH);
+ VH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(VH, v->VH, 24 + 24 * v->NbTimeStep);
+ Octree_Arrange(VH);
+ TH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(TH, v->TH, 24 + 72 * v->NbTimeStep);
+ Octree_Arrange(TH);
+
+ SI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(SI, v->SI, 18 + 6 * v->NbTimeStep);
+ Octree_Arrange(SI);
+ VI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(VI, v->VI, 18 + 18 * v->NbTimeStep);
+ Octree_Arrange(VI);
+ TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(TI, v->TI, 18 + 54 * v->NbTimeStep);
+ Octree_Arrange(TI);
+
+ SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(SY, v->SY, 15 + 5 * v->NbTimeStep);
+ Octree_Arrange(SY);
+ VY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(VY, v->VY, 15 + 15 * v->NbTimeStep);
+ Octree_Arrange(VY);
+ TY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(TY, v->TY, 15 + 45 * v->NbTimeStep);
+ Octree_Arrange(TY);
+}
+
+bool OctreePost::getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int timestep, double *values)
+{
+ if(!in) return false;
+
+ double *X = (double*)in, *Y = &X[nbNod], *Z = &X[2*nbNod], *V = &X[3*nbNod], U[3];
+
+ elementFactory factory;
+ element *e = factory.create(nbNod, dim, X, Y, Z);
+ if(!e) return false;
+
+ e->xyz2uvw(P, U);
+ if(timestep < 0){
+ for(int i = 0; i < theView->NbTimeStep; i++)
+ for(int j = 0; j < nbComp; j++)
+ values[nbComp * i + j] = e->interpolate(&V[nbNod * nbComp * i + j],
+ U[0], U[1], U[2], nbComp);
+ }
+ else{
+ for(int j = 0; j < nbComp; j++)
+ values[j] = e->interpolate(&V[nbNod * nbComp * timestep + j],
+ U[0], U[1], U[2], nbComp);
+ }
+
+ delete e;
+ return true;
+}
+
+bool OctreePost::searchScalar(double x, double y, double z,
+ double * values, int timestep)
+{
+ double P[3] = {x, y, z};
+
+ if(timestep < 0)
+ for(int i = 0; i < theView->NbTimeStep; i++)
+ values[i] = 0.0;
+ else
+ values[0] = 0.0;
+
+ if(getValue(Octree_Search(P, SS), 3, 4, 1, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, SH), 3, 8, 1, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, SI), 3, 6, 1, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, SY), 3, 5, 1, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, ST), 2, 3, 1, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, SQ), 2, 4, 1, P, timestep, values)) return true;
+
+ return false;
}
-bool OctreePost::searchVector(double x,
- double y,
- double z,
- double * values,
- double * size_elem,
- int timestep)
+bool OctreePost::searchVector(double x, double y, double z,
+ double * values, double * size_elem, int timestep)
{
- double P[3] = {x,y,z};
+ double P[3] = {x, y, z};
if(timestep < 0)
- for (int i = 0; i < 3*theView->NbTimeStep; ++i)
+ for(int i = 0; i < 3 * theView->NbTimeStep; i++)
values[i] = 0.0;
else
- values[0] = values[1] = values[2] = 0.0;
-
- void * inVT = Octree_Search(P, VT);
-
- // values[0] = -0.5*y;
- // values[1] = 0.5*x;
- // values[2] = 1;
- // return true;
-
- if(inVT){
- double U[3];
- double *X = (double*) inVT;
- double *Y = &X[3];
- double *Z = &X[6];
- double *V = &X[9];
- *size_elem = fabs(computeBarycentricTriangle(X, Y, Z, P, U));
- // bof
- *size_elem = sqrt(*size_elem);
- if(timestep < 0){
- for (int i = 0; i < theView->NbTimeStep; ++i){
- values[3*i] =
- V[9*i+3] * U[0] +
- V[9*i+6] * U[1] +
- V[9*i+0] * (1-U[0]-U[1]);
- values[3*i+1] =
- V[9*i+4] * U[0] +
- V[9*i+7] * U[1] +
- V[9*i+1] * (1-U[0]-U[1]);
- values[3*i+2] =
- V[9*i+5] * U[0] +
- V[9*i+8] * U[1] +
- V[9*i+2] * (1-U[0]-U[1]);
- }
- }
- else{
- values[0] =
- V[9*timestep+3] * U[0] +
- V[9*timestep+6] * U[1] +
- V[9*timestep+0] * (1-U[0]-U[1]);
- values[1] =
- V[9*timestep+4] * U[0] +
- V[9*timestep+7] * U[1] +
- V[9*timestep+1] * (1-U[0]-U[1]);
- values[2] =
- V[9*timestep+5] * U[0] +
- V[9*timestep+8] * U[1] +
- V[9*timestep+2] * (1-U[0]-U[1]);
- }
- return true;
- }
-
- void * inVS = Octree_Search(P, VS);
-
- if(inVS){
- double U[3];
- double *X = (double*) inVS;
- double *Y = &X[4];
- double *Z = &X[8];
- double *V = &X[12];
- computeBarycentricSimplex(X, Y, Z, P, U);
- // bof
- *size_elem = .1;
- if(timestep < 0){
- for (int i=0;i<theView->NbTimeStep;++i){
- values[3*i] =
- V[12*i+3] * U[0] +
- V[12*i+6] * U[1] +
- V[12*i+9] * U[2] +
- V[12*i+0] * (1-U[0]-U[1]-U[2]);
- values[3*i+1] =
- V[12*i+4] * U[0] +
- V[12*i+9] * U[1] +
- V[12*i+10] * U[2] +
- V[12*i+1] * (1-U[0]-U[1]-U[2]);
- values[3*i+2] =
- V[12*i+5] * U[0] +
- V[12*i+7] * U[1] +
- V[12*i+11] * U[2] +
- V[12*i+2] * (1-U[0]-U[1]-U[2]);
- }
- }
- else{
- values[0] =
- V[12*timestep+3] * U[0] +
- V[12*timestep+6] * U[1] +
- V[12*timestep+9] * U[2] +
- V[12*timestep+0] * (1-U[0]-U[1]-U[2]);
- values[1] =
- V[12*timestep+4] * U[0] +
- V[12*timestep+7] * U[1] +
- V[12*timestep+10] * U[2] +
- V[12*timestep+1 ] * (1-U[0]-U[1]-U[2]);
- values[2] =
- V[12*timestep+5 ] * U[0] +
- V[12*timestep+8 ] * U[1] +
- V[12*timestep+11] * U[2] +
- V[12*timestep+2 ] * (1-U[0]-U[1]-U[2]);
- }
- return true;
- }
+ for(int i = 0; i < 3; i++)
+ values[i] = 0.0;
+
+ // FIXME: compute this!
+ *size_elem = 1.;
+
+ if(getValue(Octree_Search(P, VS), 3, 4, 3, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, VH), 3, 8, 3, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, VI), 3, 6, 3, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, VY), 3, 5, 3, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, VT), 2, 3, 3, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, VQ), 2, 4, 3, P, timestep, values)) return true;
+
return false;
}
-bool OctreePost::searchScalar(double x,
- double y,
- double z,
- double * values,
- int timestep)
+bool OctreePost::searchTensor(double x, double y, double z,
+ double * values, int timestep)
{
- double P[3] = {x,y,z};
- void * inST = Octree_Search(P, ST);
+ double P[3] = {x, y, z};
if(timestep < 0)
- for(int i = 0; i <theView->NbTimeStep; ++i)
+ for(int i = 0; i < 9 * theView->NbTimeStep; i++)
values[i] = 0.0;
else
- values[0] = 0.0;
+ for(int i = 0; i < 9; i++)
+ values[i] = 0.0;
+
+ if(getValue(Octree_Search(P, TS), 3, 4, 9, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, TH), 3, 8, 9, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, TI), 3, 6, 9, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, TY), 3, 5, 9, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, TT), 2, 3, 9, P, timestep, values)) return true;
+ if(getValue(Octree_Search(P, TQ), 2, 4, 9, P, timestep, values)) return true;
- if(inST){
- double U[3];
- double *X = (double*) inST;
- double *Y = &X[3];
- double *Z = &X[6];
- double *V = &X[9];
- computeBarycentricTriangle(X, Y, Z, P, U);
- if(timestep < 0){
- for (int i = 0; i < theView->NbTimeStep; ++i){
- values[i] =
- V[3*i+1] * U[0] +
- V[3*i+2] * U[1] +
- V[3*i+0] * (1-U[0]-U[1]);
- }
- }
- else{
- values[0] =
- V[3*timestep+1] * U[0] +
- V[3*timestep+2] * U[1] +
- V[3*timestep+0] * (1-U[0]-U[1]);
- }
- return true;
- }
-
- void * inSS = Octree_Search(P, SS);
-
- if(inSS){
- double U[3];
- double *X = (double*) inSS;
- double *Y = &X[4];
- double *Z = &X[8];
- double *V = &X[12];
- computeBarycentricSimplex(X, Y, Z, P, U);
- if(timestep < 0){
- for (int i = 0; i < theView->NbTimeStep; ++i){
- values[i] =
- V[4*i+1] * U[0] +
- V[4*i+2] * U[1] +
- V[4*i+3] * U[2] +
- V[4*i+0] * (1.-U[0]-U[1]-U[2]);
- }
- }
- else{
- values[0] =
- V[4*timestep+1] * U[0] +
- V[4*timestep+2] * U[1] +
- V[4*timestep+3] * U[2] +
- V[4*timestep ] * (1-U[0]-U[1]-U[2]);
- }
- return true;
- }
-
return false;
}
diff --git a/Plugin/OctreePost.h b/Plugin/OctreePost.h
index 6432639d24938f90e8c74b6c53e225702f2123ad..fcbab6b0d237e6c84688b910a3a6f3179a5a8eed 100644
--- a/Plugin/OctreePost.h
+++ b/Plugin/OctreePost.h
@@ -26,28 +26,29 @@ class Post_View;
class OctreePost
{
- Octree *ST,*VT,*TT;
- Octree *SS,*VS,*TS;
+ Octree *ST, *VT, *TT;
+ Octree *SQ, *VQ, *TQ;
+ Octree *SS, *VS, *TS;
+ Octree *SH, *VH, *TH;
+ Octree *SI, *VI, *TI;
+ Octree *SY, *VY, *TY;
Post_View *theView;
+ bool getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int timestep, double *values);
public :
- OctreePost ( Post_View * );
- ~OctreePost ();
+ OctreePost(Post_View *);
+ ~OctreePost();
// search for the value of the View at point
// x, y, z. Values is interpolated linearly in
// the post element. If several time steps
// are present, they are all interpolated unless
// time step is set to a different value than -1.
- bool searchScalar ( double x ,
- double y ,
- double z,
- double * values ,
- int timestep = -1);
- bool searchVector ( double x ,
- double y ,
- double z,
- double * values ,
- double * size_elem ,
- int timestep = -1);
+ bool searchScalar(double x, double y, double z,
+ double * values, int timestep = -1);
+ bool searchVector(double x, double y, double z,
+ double * values, double * size_elem, int timestep = -1);
+ bool searchTensor(double x, double y, double z,
+ double * values, int timestep = -1);
};
#endif
diff --git a/Plugin/Probe.cpp b/Plugin/Probe.cpp
index c11ae9a637e05c169fce881b0ed67b06063d6b25..cc45e2db39fe1307e5d3e66ef1fcee1c02ebadc9 100644
--- a/Plugin/Probe.cpp
+++ b/Plugin/Probe.cpp
@@ -1,4 +1,4 @@
-// $Id: Probe.cpp,v 1.7 2005-01-09 02:19:00 geuzaine Exp $
+// $Id: Probe.cpp,v 1.8 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -131,9 +131,9 @@ void GMSH_ProbePlugin::getInfos(char *author, char *copyright,
strcpy(author, "C. Geuzaine (geuzaine@acm.caltech.edu)");
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
- "Plugin(Probe) gets the value of the simplectic view\n"
- "`iView' at the point (`X',`Y',`Z'). If `iView' < 0,\n"
- "the plugin is run on the current view.\n"
+ "Plugin(Probe) gets the value of the view `iView' at\n"
+ "the point (`X',`Y',`Z'). If `iView' < 0, the plugin is\n"
+ "run on the current view.\n"
"\n"
"Plugin(Probe) creates one new view.\n");
}
@@ -189,14 +189,22 @@ Post_View *GMSH_ProbePlugin::execute(Post_View * v)
List_Add(v2->VP, &y);
List_Add(v2->VP, &z);
for(int i = 0; i < v1->NbTimeStep; i++){
- List_Add(v2->VP, &val[3*i]);
- List_Add(v2->VP, &val[3*i+1]);
- List_Add(v2->VP, &val[3*i+2]);
+ for(int j = 0; j < 3; j++)
+ List_Add(v2->VP, &val[3*i+j]);
}
v2->NbVP++;
}
- // FIXME: do the tensor stuff
+ if(o.searchTensor(x, y, z, val)){
+ List_Add(v2->TP, &x);
+ List_Add(v2->TP, &y);
+ List_Add(v2->TP, &z);
+ for(int i = 0; i < v1->NbTimeStep; i++){
+ for(int j = 0; j < 9; j++)
+ List_Add(v2->TP, &val[9*i+j]);
+ }
+ v2->NbTP++;
+ }
delete [] val;
diff --git a/Plugin/ShapeFunctions.h b/Plugin/ShapeFunctions.h
index 3b4585b72ad443ce62f64e9740c9e0545d8982fd..ccc11d8ce0ecc23622a2ae28cc9820125c3deab8 100644
--- a/Plugin/ShapeFunctions.h
+++ b/Plugin/ShapeFunctions.h
@@ -22,6 +22,9 @@
#include "Numeric.h"
+#define ONE (1. + 1.e-6)
+#define ZERO (-1.e-6)
+
class element{
protected:
double *_x, *_y, *_z;
@@ -186,6 +189,47 @@ public:
Msg(GERROR, "integrateFlux not available for this element");
return 0.;
}
+ virtual void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ // general newton routine for the nonlinear case (more efficient
+ // routines are implemented for simplices, where the basis
+ // functions are linear)
+ uvw[0] = uvw[1] = uvw[2] = 0.0;
+
+ int iter = 1, maxiter = 20;
+ double error = 1., tol = 1.e-6;
+
+ while (error > tol && iter < maxiter){
+ double jac[3][3];
+ if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
+
+ double xn = 0., yn = 0., zn = 0.;
+ for (int i = 0; i < getNumNodes(); i++) {
+ double s;
+ getShapeFunction(i, uvw[0], uvw[1], uvw[2], s);
+ xn += _x[i] * s;
+ yn += _y[i] * s;
+ zn += _z[i] * s;
+ }
+ double inv[3][3];
+ inv3x3(jac, inv);
+
+ double un = uvw[0] +
+ inv[0][0] * (xyz[0] - xn) + inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
+ double vn = uvw[1] +
+ inv[0][1] * (xyz[0] - xn) + inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn) ;
+ double wn = uvw[2] +
+ inv[0][2] * (xyz[0] - xn) + inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn) ;
+
+ error = sqrt(SQR(un - uvw[0]) + SQR(vn - uvw[1]) + SQR(wn - uvw[2]));
+ uvw[0] = un;
+ uvw[1] = vn;
+ uvw[2] = wn;
+ iter++ ;
+ }
+ //if(error > tol) Msg(WARNING, "Newton did not converge in xyz2uvw") ;
+ }
+ virtual int isInside(double u, double v, double w) = 0;
};
class point : public element{
@@ -214,6 +258,14 @@ public:
{
s[0] = s[1] = s[2] = 0.;
}
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ uvw[0] = uvw[1] = uvw[2] = 0.;
+ }
+ int isInside(double u, double v, double w)
+ {
+ return 1;
+ }
};
class line : public element{
@@ -264,6 +316,12 @@ public:
prosca(t, v, &d);
return d;
}
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || u > ONE)
+ return 0;
+ return 1;
+ }
};
class triangle : public element{
@@ -325,6 +383,26 @@ public:
prosca(n, v, &d);
return d;
}
+#if 0 // faster, but only valid for triangles in the z=0 plane
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ double mat[2][2], b[2];
+ mat[0][0] = _x[1] - _x[0];
+ mat[0][1] = _x[2] - _x[0];
+ mat[1][0] = _y[1] - _y[0];
+ mat[1][1] = _y[2] - _y[0];
+ b[0] = xyz[0] - _x[0];
+ b[1] = xyz[1] - _y[0];
+ sys2x2(mat, b, uvw);
+ uvw[2] = 0.;
+ }
+#endif
+ int isInside(double u, double v, double w)
+ {
+ if(u < ZERO || v < ZERO || u > (ONE - v))
+ return 0;
+ return 1;
+ }
};
class quadrangle : public element{
@@ -389,6 +467,12 @@ public:
prosca(n, v, &d);
return d;
}
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || v < -ONE || u > ONE || v > ONE)
+ return 0;
+ return 1;
+ }
};
class tetrahedron : public element{
@@ -439,6 +523,30 @@ public:
default : s[0] = s[1] = s[2] = 0.; break;
}
}
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ double mat[3][3], b[3];
+ mat[0][0] = _x[1] - _x[0];
+ mat[0][1] = _x[2] - _x[0];
+ mat[0][2] = _x[3] - _x[0];
+ mat[1][0] = _y[1] - _y[0];
+ mat[1][1] = _y[2] - _y[0];
+ mat[1][2] = _y[3] - _y[0];
+ mat[2][0] = _z[1] - _z[0];
+ mat[2][1] = _z[2] - _z[0];
+ mat[2][2] = _z[3] - _z[0];
+ b[0] = xyz[0] - _x[0];
+ b[1] = xyz[1] - _y[0];
+ b[2] = xyz[2] - _z[0];
+ double det;
+ sys3x3(mat, b, uvw, &det);
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < ZERO || v < ZERO || w < ZERO || u > (ONE - v - w))
+ return 0;
+ return 1;
+ }
};
class hexahedron : public element{
@@ -521,6 +629,12 @@ public:
default : s[0] = s[1] = s[2] = 0.; break;
}
}
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || v < -ONE || w < -ONE || u > ONE || v > ONE || w > ONE)
+ return 0;
+ return 1;
+ }
};
class prism : public element{
@@ -593,6 +707,12 @@ public:
default : s[0] = s[1] = s[2] = 0.; break;
}
}
+ int isInside(double u, double v, double w)
+ {
+ if(w > ONE || w < -ONE || u < ZERO || v < ZERO || u > (ONE - v))
+ return 0;
+ return 1;
+ }
};
class pyramid : public element{
@@ -678,6 +798,13 @@ public:
}
}
}
+ int isInside(double u, double v, double w)
+ {
+ if(u < (w - ONE) || u > (ONE - w) || v < (w - ONE) || v > (ONE - w) ||
+ w < ZERO || w > ONE)
+ return 0;
+ return 1;
+ }
};
class elementFactory{
@@ -701,4 +828,7 @@ class elementFactory{
}
};
+#undef ONE
+#undef ZERO
+
#endif
diff --git a/Plugin/StreamLines.cpp b/Plugin/StreamLines.cpp
index 1f6d67fff91e8f54845aff434453c8a0c7844677..e89465956e8cd52408b6f11dae8658f14bed0d89 100644
--- a/Plugin/StreamLines.cpp
+++ b/Plugin/StreamLines.cpp
@@ -1,4 +1,4 @@
-// $Id: StreamLines.cpp,v 1.20 2005-01-09 02:19:00 geuzaine Exp $
+// $Id: StreamLines.cpp,v 1.21 2005-03-02 07:49:41 geuzaine Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -171,9 +171,9 @@ void GMSH_StreamLinesPlugin::getInfos(char *author, char *copyright,
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
"Plugin(StreamLines) computes stream lines\n"
- "from a triangle/tetrahedron vector view `iView'\n"
- "and optionally interpolates the scalar view `dView'\n"
- "on the resulting stream lines. It takes as input a\n"
+ "from a vector view `iView' and optionally\n"
+ "interpolates the scalar view `dView' on the\n"
+ "resulting stream lines. It takes as input a\n"
"grid defined by the 3 points (`X0',`Y0',`Z0')\n"
"(origin), (`X1',`Y1',`Z1') (axis of U) and\n"
"(`X2',`Y2',`Z2') (axis of V). The number of points\n"
diff --git a/doc/VERSIONS b/doc/VERSIONS
index bf01705104eb1ff987a60a98ff8b5de2d4b83ced..45c7d5b4f850d1698d7a26f39758a373e887125e 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,7 +1,8 @@
-$Id: VERSIONS,v 1.310 2005-02-20 06:36:58 geuzaine Exp $
+$Id: VERSIONS,v 1.311 2005-03-02 07:49:41 geuzaine Exp $
New since 1.59: added support for discrete curves; new Window menu on
-Mac OS X; fixed small bugs.
+Mac OS X; generalized octree-based plugins (CutGrid, StreamLines,
+Probe, etc.) to all element types; fixed small bugs.
New in 1.59: added support for discrete (triangulated) surfaces,
either in STL format or with the new "Discrete Surface" command; added
diff --git a/doc/texinfo/opt_plugin.texi b/doc/texinfo/opt_plugin.texi
index f81c5d757b68c1e7d4394f7cb97b0178b754e89a..cd6ffcb7c6eaf60c944f703c6cd2d0c84a82ed93 100644
--- a/doc/texinfo/opt_plugin.texi
+++ b/doc/texinfo/opt_plugin.texi
@@ -51,17 +51,16 @@ Default value: @code{-1}
@end table
@item Plugin(CutGrid)
-Plugin(CutGrid) cuts a triangle/tetrahedron view
-with a rectangular grid defined by the 3 points
+Plugin(CutGrid) cuts the view `iView' with a
+rectangular grid defined by the 3 points
(`X0',`Y0',`Z0') (origin), (`X1',`Y1',`Z1') (axis of U)
and (`X2',`Y2',`Z2') (axis of V). The number of points
along U and V is set with the options `nPointsU'
and `nPointsV'. If `ConnectPoints' is zero, the
-plugin creates scalar or vector points; otherwise,
-the plugin generates scalar or vector quadrangles,
-lines or points depending on the values of `nPointsU'
-and `nPointsV'. If `iView' < 0, the plugin is run on
-the current view.
+plugin creates points; otherwise, the plugin
+generates quadrangles, lines or points depending on
+the values of `nPointsU' and `nPointsV'. If `iView' < 0,
+the plugin is run on the current view.
Plugin(CutGrid) creates one new view.
@@ -131,14 +130,13 @@ Default value: @code{-1}
@end table
@item Plugin(CutParametric)
-Plugin(CutParametric) cuts a scalar triangle/
-tetrahedron view `iView' with the parametric function
-(`X'(u), `Y'(u), `Z'(u)), using `nPointsU' values of
-the parameter u in [`MinU', `MaxU']. If
-`ConnectPoints' is set, the plugin creates scalar
-line elements; otherwise, the plugin generates
-scalar points. If `iView' < 0, the plugin is run on
-the current view.
+Plugin(CutParametric) cuts the scalar view `iView'
+with the parametric function (`X'(u), `Y'(u), `Z'(u)),
+using `nPointsU' values of the parameter u in
+[`MinU', `MaxU']. If `ConnectPoints' is set, the
+plugin creates scalar line elements; otherwise, the
+plugin generates scalar points. If `iView' < 0, the
+plugin is run on the current view.
Plugin(CutParametric) creates one new view.
@@ -493,9 +491,9 @@ Default value: @code{-1}
@end table
@item Plugin(Probe)
-Plugin(Probe) gets the value of the simplectic view
-`iView' at the point (`X',`Y',`Z'). If `iView' < 0,
-the plugin is run on the current view.
+Plugin(Probe) gets the value of the view `iView' at
+the point (`X',`Y',`Z'). If `iView' < 0, the plugin
+is run on the current view.
Plugin(Probe) creates one new view.
@@ -608,9 +606,9 @@ Default value: @code{-1}
@item Plugin(StreamLines)
Plugin(StreamLines) computes stream lines
-from a triangle/tetrahedron vector view `iView'
-and optionally interpolates the scalar view `dView'
-on the resulting stream lines. It takes as input a
+from a vector view `iView' and optionally
+interpolates the scalar view `dView' on the
+resulting stream lines. It takes as input a
grid defined by the 3 points (`X0',`Y0',`Z0')
(origin), (`X1',`Y1',`Z1') (axis of U) and
(`X2',`Y2',`Z2') (axis of V). The number of points
diff --git a/doc/texinfo/opt_view.texi b/doc/texinfo/opt_view.texi
index 7574abf4637112178643d123f51021285a431255..db490edcad9a55a11be9e5b8ed6ba154d2610c3e 100644
--- a/doc/texinfo/opt_view.texi
+++ b/doc/texinfo/opt_view.texi
@@ -6,7 +6,7 @@ Saved in: @code{General.OptionsFileName}
@item View.AbscissaFormat
Abscissa number format for 2D graphs (in standard C form)@*
-Default value: @code{"%.3e"}@*
+Default value: @code{"%.3g"}@*
Saved in: @code{General.OptionsFileName}
@item View.FileName
@@ -16,7 +16,7 @@ Saved in: @code{-}
@item View.Format
Number format (in standard C form)@*
-Default value: @code{"%.3e"}@*
+Default value: @code{"%.3g"}@*
Saved in: @code{General.OptionsFileName}
@item View.GeneralizedRaiseX