diff --git a/Plugin/Gradient.cpp b/Plugin/Gradient.cpp
index 75655036dd74cba1d56d997c73088e4262d302cc..43fef86084b01bde08fe0ed0411467e1dfe29c84 100644
--- a/Plugin/Gradient.cpp
+++ b/Plugin/Gradient.cpp
@@ -1,4 +1,4 @@
-// $Id: Gradient.cpp,v 1.1 2004-11-26 10:31:54 remacle Exp $
+// $Id: Gradient.cpp,v 1.2 2004-11-26 16:16:39 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -30,10 +30,6 @@
#include <math.h>
#include <stdio.h>
-#if defined(HAVE_MATH_EVAL)
-#include "matheval.h"
-#endif
-
extern Context_T CTX;
StringXNumber GradientOptions_Number[] = {
@@ -48,11 +44,6 @@ extern "C"
}
}
-int GMSH_GradientPlugin::getNbOptionsStr() const
-{
- return 0;
-}
-
int GMSH_GradientPlugin::getNbOptions() const
{
return sizeof(GradientOptions_Number) / sizeof(StringXNumber);
@@ -63,13 +54,6 @@ StringXNumber *GMSH_GradientPlugin::getOption(int iopt)
return &GradientOptions_Number[iopt];
}
-StringXString *GMSH_GradientPlugin::getOptionStr(int iopt)
-{
- throw;
-}
-
-
-
GMSH_GradientPlugin::GMSH_GradientPlugin()
{
;
@@ -86,36 +70,25 @@ void GMSH_GradientPlugin::getInfos(char *author, char *copyright,
strcpy(author, "E. Marchandise");
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
- "Plugin(Gradient) computes\n"
- "the eigenvalues Lambda (1,2,3)\n"
- "of the tensor S_ik S_kj + Om_ik Om_kj,\n"
- "where S_ij = 0.5 (ui,j + uj,i)) and\n"
- "Om_ij = 0.5 (ui,j - uj,i) are respectively\n"
- "the symmetric and antisymmetric parts\n"
- "of the velocity gradient tensor.\n"
- "Vortices are well-represented by regions\n"
- "where Lambda (2) is negative.\n"
+ "Plugin(Gradient) computes XXX\n"
"If `iView' < 0, the plugin is run\n"
"on the current view.\n"
"\n"
- "Plugin(Lambda2) is executed in-place.\n");
+ "Plugin(Gradient) creates one new view.\n");
}
-
void GMSH_GradientPlugin::catchErrorMessage(char *errorMessage) const
{
strcpy(errorMessage, "Gradient failed...");
}
-//-----------------------------------------------------------------------------
-static void computeGradient( List_T *inList, int inNb,
- List_T *outListVector, int *outNbVector,
- int nbTime, int nbNod, int nbComp)
+static void computeGradient(List_T *inList, int inNb,
+ List_T *outListVector, int *outNbVector,
+ int nbTime, int nbNod, int nbComp)
{
-
- if(!inNb)
+ if(!inNb)
return;
-
+
int outNbComp, *outNb;
List_T *outList;
@@ -128,132 +101,127 @@ static void computeGradient( List_T *inList, int inNb,
int nb = List_Nbr(inList) / inNb;
- //on parcourt les elements
- //------------------------
- //------------------------
+ // on parcourt les elements
for(int i = 0; i < List_Nbr(inList); i += nb) {
-
+
double *yy = (double *)List_Pointer_Fast(inList, i + nbNod);
- if (yy[0] > 0)
- {
- for(int j = 0; j < 3 * nbNod; j++)
- List_Add(outList, List_Pointer_Fast(inList, i + j));
+ if(yy[0] > 0){
+ for(int j = 0; j < 3 * nbNod; j++)
+ List_Add(outList, List_Pointer_Fast(inList, i + j));
+
+ for(int t = 0; t < nbTime; t++){
+
+ //on lit vecteurs positions et vitesses
+ //-------------------------------------
+ double *x = (double *)List_Pointer_Fast(inList, i);
+ double *y = (double *)List_Pointer_Fast(inList, i + nbNod);
+ double *z = (double *)List_Pointer_Fast(inList, i + 2 * nbNod);
+ double val[MAX_COMP][MAX_NOD];
+ for(int k = 0; k < nbNod; k++){
+ double *v = (double *)List_Pointer_Fast(inList, i + 3 * nbNod +
+ nbNod * nbComp * t + nbComp * k);
+ for(int l = 0; l < nbComp; l++){
+ val[l][k] = v[l];
+ }
+ }
+
+ //on calcule le gradient des fonctions de forme
+ //---------------------------------------------
+ double GradPhi_x[4][3];
+ double GradPhi_ksi[4][3];
+ double dx_dksi[3][3];
+ double dksi_dx[3][3];
+ double det;
+
+ //triangles
+ //+++++++++
+ if (nbNod == 3){
+ double a[3], b[3], cross[3];
+ a[0]= x[1]-x[0]; a[1]= y[1]-y[0]; a[2]= z[1]-z[0];
+ b[0]= x[2]-x[0]; b[1]= y[2]-y[0]; b[2]= z[2]-z[0];
+ prodve(a, b, cross);
+ dx_dksi[0][0] = x[1] - x[0]; dx_dksi[0][1] = x[2]-x[0]; dx_dksi[0][2] = cross[0];
+ dx_dksi[1][0] = y[1] - y[0]; dx_dksi[1][1] = y[2]-y[0]; dx_dksi[1][2] = cross[1];
+ dx_dksi[2][0] = z[1] - z[0]; dx_dksi[2][1] = z[2]-z[0]; dx_dksi[2][2] = cross[2];
+ inv3x3(dx_dksi, dksi_dx, &det);
+
+ GradPhi_ksi[0][0]= -1; GradPhi_ksi[0][1]= -1; GradPhi_ksi[0][2]= 0;
+ GradPhi_ksi[1][0]= 1; GradPhi_ksi[1][1]= 0; GradPhi_ksi[1][2]= 0;
+ GradPhi_ksi[2][0]= 0; GradPhi_ksi[2][1]= 1; GradPhi_ksi[2][2]= 0;
+ }
- for(int t = 0; t < nbTime; t++){
-
- //on lit vecteurs positions et vitesses
- //-------------------------------------
- double *x = (double *)List_Pointer_Fast(inList, i);
- double *y = (double *)List_Pointer_Fast(inList, i + nbNod);
- double *z = (double *)List_Pointer_Fast(inList, i + 2 * nbNod);
- double val[MAX_COMP][MAX_NOD];
- for(int k = 0; k < nbNod; k++){
- double *v = (double *)List_Pointer_Fast(inList, i + 3 * nbNod + nbNod * nbComp * t + nbComp * k);
- for(int l = 0; l < nbComp; l++){
- val[l][k] = v[l];
- }
+ //tetraedre
+ //++++++++++
+ if (nbNod == 4){
+ dx_dksi[0][0] = x[1] - x[0]; dx_dksi[0][1] = x[2]-x[0]; dx_dksi[0][2] = x[3]-x[0];
+ dx_dksi[1][0] = y[1] - y[0]; dx_dksi[1][1] = y[2]-y[0]; dx_dksi[1][2] = y[3]-y[0];
+ dx_dksi[2][0] = z[1] - z[0]; dx_dksi[2][1] = z[2]-z[0]; dx_dksi[2][2] = z[3]-z[0];
+ inv3x3(dx_dksi, dksi_dx, &det);
+ GradPhi_ksi[0][0]= -1; GradPhi_ksi[0][1]= -1; GradPhi_ksi[0][2]= -1;
+ GradPhi_ksi[1][0]= 1; GradPhi_ksi[1][1]= 0; GradPhi_ksi[1][2]= 0;
+ GradPhi_ksi[2][0]= 0; GradPhi_ksi[2][1]= 1; GradPhi_ksi[2][2]= 0;
+ GradPhi_ksi[3][0]= 0; GradPhi_ksi[3][1]= 0; GradPhi_ksi[3][2]= 1;
+ }
+
+ for(int k = 0; k < nbNod ; k++){
+ for(int l = 0; l < 3; l++){
+ GradPhi_x[k][l] = 0.0;
+ for(int m = 0; m < 3; m++){
+ GradPhi_x[k][l] += GradPhi_ksi[k][m] * dksi_dx[l][m];
}
-
-
- //on calcule le gradient des fonctions de forme
- //---------------------------------------------
- double GradPhi_x[4][3];
- double GradPhi_ksi[4][3];
- double dx_dksi[3][3];
- double dksi_dx[3][3];
- double det;
-
- //triangles
- //+++++++++
- if (nbNod == 3){
- double a[3], b[3], cross[3];
- a[0]= x[1]-x[0]; a[1]= y[1]-y[0]; a[2]= z[1]-z[0];
- b[0]= x[2]-x[0]; b[1]= y[2]-y[0]; b[2]= z[2]-z[0];
- prodve(a, b, cross);
- dx_dksi[0][0] = x[1] - x[0]; dx_dksi[0][1] = x[2]-x[0]; dx_dksi[0][2] = cross[0];
- dx_dksi[1][0] = y[1] - y[0]; dx_dksi[1][1] = y[2]-y[0]; dx_dksi[1][2] = cross[1];
- dx_dksi[2][0] = z[1] - z[0]; dx_dksi[2][1] = z[2]-z[0]; dx_dksi[2][2] = cross[2];
- inv3x3(dx_dksi, dksi_dx, &det);
-
- GradPhi_ksi[0][0]= -1; GradPhi_ksi[0][1]= -1; GradPhi_ksi[0][2]= 0;
- GradPhi_ksi[1][0]= 1; GradPhi_ksi[1][1]= 0; GradPhi_ksi[1][2]= 0;
- GradPhi_ksi[2][0]= 0; GradPhi_ksi[2][1]= 1; GradPhi_ksi[2][2]= 0;
- }
-
- //tetraedre
- //++++++++++
- if (nbNod == 4){
- dx_dksi[0][0] = x[1] - x[0]; dx_dksi[0][1] = x[2]-x[0]; dx_dksi[0][2] = x[3]-x[0];
- dx_dksi[1][0] = y[1] - y[0]; dx_dksi[1][1] = y[2]-y[0]; dx_dksi[1][2] = y[3]-y[0];
- dx_dksi[2][0] = z[1] - z[0]; dx_dksi[2][1] = z[2]-z[0]; dx_dksi[2][2] = z[3]-z[0];
- inv3x3(dx_dksi, dksi_dx, &det);
- GradPhi_ksi[0][0]= -1; GradPhi_ksi[0][1]= -1; GradPhi_ksi[0][2]= -1;
- GradPhi_ksi[1][0]= 1; GradPhi_ksi[1][1]= 0; GradPhi_ksi[1][2]= 0;
- GradPhi_ksi[2][0]= 0; GradPhi_ksi[2][1]= 1; GradPhi_ksi[2][2]= 0;
- GradPhi_ksi[3][0]= 0; GradPhi_ksi[3][1]= 0; GradPhi_ksi[3][2]= 1;
- }
-
- for(int k = 0; k < nbNod ; k++){
- for(int l = 0; l < 3; l++){
- GradPhi_x[k][l] = 0.0;
- for(int m = 0; m < 3; m++){
- GradPhi_x[k][l] += GradPhi_ksi[k][m] * dksi_dx[l][m];
- }
- }
- }
-
- //on calcule le gradient des vitesses
- //----------------------------------
-
- //par noeud on a max 9 composantes
- double res, grad[3][3], sym[3][3], asym[3][3];
-// for(int l = outNbComp; l < 9; l++)
-// d[l] = 0.0;
- for(int k = 0; k < nbComp; k++){
- for(int l = 0; l < 3; l++){
- grad[k][l] = 0.0;
- for(int m = 0; m < nbNod; m++){
- grad[k][l] += val[k][m]* GradPhi_x[m][l];
+ }
+ }
+
+ //on calcule le gradient des vitesses
+ //----------------------------------
+
+ //par noeud on a max 9 composantes
+ double res, grad[3][3], sym[3][3], asym[3][3];
+ // for(int l = outNbComp; l < 9; l++)
+ // d[l] = 0.0;
+ for(int k = 0; k < nbComp; k++){
+ for(int l = 0; l < 3; l++){
+ grad[k][l] = 0.0;
+ for(int m = 0; m < nbNod; m++){
+ grad[k][l] += val[k][m]* GradPhi_x[m][l];
+ }
+ }
+ }
+
+ for(int k=0; k<nbComp; k++){
+ for(int l=0; l<3; l++){
+ sym[k][l] = 0.5*(grad[k][l]+grad[l][k]);
+ asym[k][l] = 0.5*(grad[k][l]-grad[l][k]);
+ }
}
+
+ double a[3][3];
+ for(int i=0; i<3; i++)
+ for(int j=0; j<3; j++){
+ a[i][j] = 0.0;
+ for(int k=0; k<3; k++)
+ a[i][j] += sym[i][k]*sym[k][j] + asym[i][k]*asym[k][j];
+ }
+
+ //calcul de lambda avec tri
+ double lambda[3];
+ eigenvalue(a, lambda);
+
+ for(int k = 0; k < nbNod; k++){
+ for(int l = 0; l < outNbComp; l++){
+ res = lambda[1];
+ List_Add(outList, &res);
+ }
+ }
}
- }
-
- for(int k=0; k<nbComp; k++){
- for(int l=0; l<3; l++){
- sym[k][l] = 0.5*(grad[k][l]+grad[l][k]);
- asym[k][l] = 0.5*(grad[k][l]-grad[l][k]);
- }
- }
-
- double a[3][3];
- for(int i=0; i<3; i++)
- for(int j=0; j<3; j++){
- a[i][j] = 0.0;
- for(int k=0; k<3; k++)
- a[i][j] += sym[i][k]*sym[k][j] + asym[i][k]*asym[k][j];
- }
-
- //calcul de lambda avec tri
- double lambda[3];
- eigenvalue(a, lambda);
-
- for(int k = 0; k < nbNod; k++){
- for(int l = 0; l < outNbComp; l++){
- res = lambda[1];
- List_Add(outList, &res);
- }
- }
- }
- (*outNb)++;
+ (*outNb)++;
}
}
}
-//-----------------------------------------------------------------------------
Post_View *GMSH_GradientPlugin::execute(Post_View * v)
{
-
int iView = (int)GradientOptions_Number[0].def;
if(iView < 0)
@@ -264,73 +232,60 @@ Post_View *GMSH_GradientPlugin::execute(Post_View * v)
return v;
}
- Post_View *v2 = BeginView(1);
-
- Post_View *v1 = (Post_View*)List_Pointer(CTX.post.list, iView);
+ Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
-// points
-//--------
-// computeGradient( v1->SP, v1->NbSP, v2->VP, &v2->NbVP, v1->NbTimeStep, 1, 1);
-// computeGradient( v1->VP, v1->NbVP, v2->TP, &v2->NbTP, v1->NbTimeStep, 1, 3);
+ Post_View *v2 = BeginView(1);
-// lines
-//-------
+ // points
+ //--------
+ // computeGradient( v1->SP, v1->NbSP, v2->VP, &v2->NbVP, v1->NbTimeStep, 1, 1);
+ // computeGradient( v1->VP, v1->NbVP, v2->TP, &v2->NbTP, v1->NbTimeStep, 1, 3);
+
+ // lines
+ //-------
// computeGradient( v1->SL, v1->NbSL, v2->VL, &v2->NbVL, v1->NbTimeStep, 2, 1);
// computeGradient( v1->VL, v1->NbVL, v2->TL, &v2->NbTL, v1->NbTimeStep, 2, 3);
-
-
-// triangles
-//------------
-// computeGradient( v1->ST, v1->NbST, v2->VT, &v2->NbVT, v1->NbTimeStep, 3, 1);
-// computeGradient( v1->VT, v1->NbVT, v2->TT, &v2->NbTT, v1->NbTimeStep, 3, 3);
-
-
-// tets
-//------
-// computeGradient( v1->SS, v1->NbSS, v2->VS, &v2->NbVS, v1->NbTimeStep, 4, 1);
+
+ // triangles
+ //------------
+ // computeGradient( v1->ST, v1->NbST, v2->VT, &v2->NbVT, v1->NbTimeStep, 3, 1);
+ // computeGradient( v1->VT, v1->NbVT, v2->TT, &v2->NbTT, v1->NbTimeStep, 3, 3);
+
+ // tets
+ //------
+ // computeGradient( v1->SS, v1->NbSS, v2->VS, &v2->NbVS, v1->NbTimeStep, 4, 1);
computeGradient( v1->VS, v1->NbVS, v2->SS, &v2->NbSS, v1->NbTimeStep, 4, 3);
-
-
-/*
-// quadrangles
-//-----------
-computeGradient( v1->SQ, v1->NbSQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 1);
-computeGradient( v1->VQ, v1->NbVQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 3);
-computeGradient( v1->TQ, v1->NbTQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 9);
-// hexas
-//--------
-computeGradient( v1->SH, v1->NbSH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 1);
-computeGradient( v1->VH, v1->NbVH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 3);
-computeGradient( v1->TH, v1->NbTH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 9);
-// prisms
-//--------
-computeGradient( v1->SI, v1->NbSI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 1);
-computeGradient( v1->VI, v1->NbVI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 3);
-computeGradient( v1->TI, v1->NbTI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 9);
-// pyramids
-//----------
-computeGradient( v1->SY, v1->NbSY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 1);
-computeGradient( v1->VY, v1->NbVY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 3);
-computeGradient( v1->TY, v1->NbTY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 9);
-*/
-
- if(v2->empty()) {
- RemoveViewByNumber(v2->Num);
- return v1;
- }
- else{
- // copy time data
- for(int i = 0; i < List_Nbr(v1->Time); i++)
- List_Add(v2->Time, List_Pointer(v1->Time, i));
- // finalize
- char name[1024], filename[1024];
- sprintf(name, "%s_Gradient", v1->Name);
- sprintf(filename, "%s_Gradient.pos", v1->Name);
- EndView(v2, 1, filename, name);
- return v2;
- }
+
+ /*
+ // quadrangles
+ //-----------
+ computeGradient( v1->SQ, v1->NbSQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 1);
+ computeGradient( v1->VQ, v1->NbVQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 3);
+ computeGradient( v1->TQ, v1->NbTQ, v2->VQ, &v2->NbVQ, v1->NbTimeStep, 4, 9);
+ // hexas
+ //--------
+ computeGradient( v1->SH, v1->NbSH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 1);
+ computeGradient( v1->VH, v1->NbVH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 3);
+ computeGradient( v1->TH, v1->NbTH, v2->VH, &v2->NbVH, v1->NbTimeStep, 8, 9);
+ // prisms
+ //--------
+ computeGradient( v1->SI, v1->NbSI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 1);
+ computeGradient( v1->VI, v1->NbVI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 3);
+ computeGradient( v1->TI, v1->NbTI, v2->VI, &v2->NbVI, v1->NbTimeStep, 6, 9);
+ // pyramids
+ //----------
+ computeGradient( v1->SY, v1->NbSY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 1);
+ computeGradient( v1->VY, v1->NbVY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 3);
+ computeGradient( v1->TY, v1->NbTY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 9);
+ */
+
+ // copy time data
+ for(int i = 0; i < List_Nbr(v1->Time); i++)
+ List_Add(v2->Time, List_Pointer(v1->Time, i));
+ // finalize
+ char name[1024], filename[1024];
+ sprintf(name, "%s_Gradient", v1->Name);
+ sprintf(filename, "%s_Gradient.pos", v1->Name);
+ EndView(v2, 1, filename, name);
+ return v2;
}
-
-
-
-//----------------------------------------------------------------------
diff --git a/Plugin/Gradient.h b/Plugin/Gradient.h
index 19ec8b91d71c5cc28145ce465ec244957f7c3733..516ef94946a69276f4dfe6ccecb859705fb8f163 100644
--- a/Plugin/Gradient.h
+++ b/Plugin/Gradient.h
@@ -37,8 +37,6 @@ public:
void catchErrorMessage(char *errorMessage) const;
int getNbOptions() const;
StringXNumber* getOption(int iopt);
- int getNbOptionsStr() const;
- StringXString* getOptionStr(int iopt);
Post_View *execute(Post_View *);
};
diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp
index 0062f1ae289b9b29dd137e80bdb1b9687593b3be..5300a519de8f302bb5c834c2cc9c41728c902fa3 100644
--- a/Plugin/Integrate.cpp
+++ b/Plugin/Integrate.cpp
@@ -1,4 +1,4 @@
-// $Id: Integrate.cpp,v 1.8 2004-11-26 15:06:50 remacle Exp $
+// $Id: Integrate.cpp,v 1.9 2004-11-26 16:16:39 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -30,235 +30,199 @@
extern Context_T CTX;
StringXNumber IntegrateOptions_Number[] = {
- {GMSH_FULLRC, "iView", NULL, -1.},
- {GMSH_FULLRC, "computeLevelsetPositive", NULL, 0.}
+ {GMSH_FULLRC, "iView", NULL, -1.},
+ {GMSH_FULLRC, "computeLevelsetPositive", NULL, 0.}
};
extern "C"
{
- GMSH_Plugin *GMSH_RegisterIntegratePlugin()
- {
- return new GMSH_IntegratePlugin();
- }
+ GMSH_Plugin *GMSH_RegisterIntegratePlugin()
+ {
+ return new GMSH_IntegratePlugin();
+ }
}
GMSH_IntegratePlugin::GMSH_IntegratePlugin()
{
- ;
+ ;
}
void GMSH_IntegratePlugin::getName(char *name) const
{
- strcpy(name, "Integrate");
+ strcpy(name, "Integrate");
}
void GMSH_IntegratePlugin::getInfos(char *author, char *copyright,
char *help_text) const
{
- strcpy(author, "C. Geuzaine (geuz@geuz.org)");
- strcpy(copyright, "DGR (www.multiphysics.com)");
- strcpy(help_text,
- "Plugin(Integrate) integrates scalar fields over\n"
- "all the elements in the view `iView', as well\n"
- "as the circulation/flux of vector fields over\n"
- "line/surface elements. If `iView' < 0, the\n"
- "plugin is run on the current view. If\n"
- "`computeLevelsetPositive' is set, the plugin\n"
- "computes the positive area (volume) of the map\n"
- "(for triangle/tetrahedron maps only).\n"
- "\n"
- "Plugin(Integrate) creates one new view.\n");
+ strcpy(author, "C. Geuzaine (geuz@geuz.org)");
+ strcpy(copyright, "DGR (www.multiphysics.com)");
+ strcpy(help_text,
+ "Plugin(Integrate) integrates scalar fields over\n"
+ "all the elements in the view `iView', as well\n"
+ "as the circulation/flux of vector fields over\n"
+ "line/surface elements. If `iView' < 0, the\n"
+ "plugin is run on the current view. If\n"
+ "`computeLevelsetPositive' is set, the plugin\n"
+ "computes the positive area (volume) of the map\n"
+ "(for triangle/tetrahedron maps only).\n"
+ "\n"
+ "Plugin(Integrate) creates one new view.\n");
}
int GMSH_IntegratePlugin::getNbOptions() const
{
- return sizeof(IntegrateOptions_Number) / sizeof(StringXNumber);
+ return sizeof(IntegrateOptions_Number) / sizeof(StringXNumber);
}
StringXNumber *GMSH_IntegratePlugin::getOption(int iopt)
{
- return &IntegrateOptions_Number[iopt];
+ return &IntegrateOptions_Number[iopt];
}
void GMSH_IntegratePlugin::catchErrorMessage(char *errorMessage) const
{
- strcpy(errorMessage, "Integrate failed...");
+ strcpy(errorMessage, "Integrate failed...");
}
static double integrate(int nbList, List_T *list, int dim,
int nbNod, int nbComp, int step)
{
- if(!nbList) return 0.;
+ if(!nbList) return 0.;
- const int levelsetPositive = (int)IntegrateOptions_Number[1].def;
-
- double res = 0.;
- int nb = List_Nbr(list) / nbList;
- for(int i = 0; i < List_Nbr(list); i += nb) {
- double *x = (double *)List_Pointer_Fast(list, i);
- double *y = (double *)List_Pointer_Fast(list, i + nbNod);
- double *z = (double *)List_Pointer_Fast(list, i + 2 * nbNod);
- double *v = (double *)List_Pointer_Fast(list, i + 3 * nbNod +
- nbNod * nbComp * step);
- if(dim == 0){
- if(nbNod == 1){
- point p(x, y, z);
- if(nbComp == 1) res += p.integrate(v);
- }
+ const int levelsetPositive = (int)IntegrateOptions_Number[1].def;
+
+ double res = 0.;
+ int nb = List_Nbr(list) / nbList;
+ for(int i = 0; i < List_Nbr(list); i += nb) {
+ double *x = (double *)List_Pointer_Fast(list, i);
+ double *y = (double *)List_Pointer_Fast(list, i + nbNod);
+ double *z = (double *)List_Pointer_Fast(list, i + 2 * nbNod);
+ double *v = (double *)List_Pointer_Fast(list, i + 3 * nbNod +
+ nbNod * nbComp * step);
+ if(dim == 0){
+ if(nbNod == 1){
+ point p(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += p.integrate(v);
+ else res += p.integrateLevelsetPositive(v);
+ }
+ }
+ }
+ else if(dim == 1){
+ if(nbNod == 2){
+ line l(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += l.integrate(v);
+ else res += l.integrateLevelsetPositive(v);
+ }
+ else if(nbComp == 3) res += l.integrateCirculation(v);
+ }
+ }
+ else if(dim == 2){
+ if(nbNod == 3){
+ triangle t(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += t.integrate(v);
+ else res += t.integrateLevelsetPositive(v);
+ }
+ else if(nbComp == 3) res += t.integrateFlux(v);
+ }
+ else if(nbNod == 4){
+ quadrangle q(x, y, z);
+ if(nbComp == 1) {
+ if(!levelsetPositive) res += q.integrate(v);
+ else res += q.integrateLevelsetPositive(v);
+ }
+ else if(nbComp == 3) res += q.integrateFlux(v);
+ }
+ }
+ else if(dim == 3){
+ if(nbNod == 4){
+ tetrahedron t(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += t.integrate(v);
+ else res += t.integrateLevelsetPositive(v);
}
- else if(dim == 1){
- if(nbNod == 2){
- line l(x, y, z);
- if(nbComp == 1) res += l.integrate(v);
- else if(nbComp == 3) res += l.integrateCirculation(v);
- }
+ }
+ else if(nbNod == 8){
+ hexahedron h(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += h.integrate(v);
+ else res += h.integrateLevelsetPositive(v);
}
- else if(dim == 2){
- if(nbNod == 3){
- triangle t(x, y, z);
- if(nbComp == 1){
- if(!levelsetPositive)
- res += t.integrate(v);
- else{
- double ONES[] = { 1.0 , 1.0 , 1.0 };
- const double area = t.integrate(ONES);
- const double SUM = v[0] + v[1] + v[2];
- const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]);
- if(SUMABS){
- const double XI = SUM / SUMABS;
- res += area * (1 - XI) * 0.5 ;
- }
- }
- }
- else if(nbComp == 3) res += t.integrateFlux(v);
- }
- else if(nbNod == 4){
- quadrangle q(x, y, z);
- if(nbComp == 1) {
- if(!levelsetPositive)
- {
- res += q.integrate(v);
- }
- else
- {
- double ONES[] = { 1.0 , 1.0 , 1.0 , 1.0};
- const double area = q.integrate(ONES);
- const double SUM = v[0] + v[1] + v[2] + v[3];
- const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]) + fabs (v[3]);
- if(SUMABS != 0.0){
- const double XI = SUM / SUMABS;
- res += area * (1 - XI) * 0.5 ;
- }
- }
- }
- else if(nbComp == 3) res += q.integrateFlux(v);
- }
+ }
+ else if(nbNod == 6){
+ prism p(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += p.integrate(v);
+ else res += p.integrateLevelsetPositive(v);
}
- else if(dim == 3)
- {
- if(nbNod == 4)
- {
- tetrahedron t(x, y, z);
- if(nbComp == 1)
- {
- if ( ! levelsetPositive )
- res += t.integrate(v);
- else
- {
- double ONES[] = { 1.0 , 1.0 , 1.0, 1.0 };
- const double area = fabs(t.integrate (ONES));
- const double SUM = v[0] + v[1] + v[2] + v[3];
- const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]) + fabs (v[3]);
- const double XI = SUM / SUMABS;
- res += area * (1 - XI) * 0.5 ;
- }
- }
- }
- else if(nbNod == 8){
- hexahedron h(x, y, z);
- if(nbComp == 1)
- {
- if ( ! levelsetPositive )
- res += h.integrate(v);
- else
- {
- double ONES[] = { 1.0 , 1.0 , 1.0, 1.0 ,1,1,1,1};
- const double area = fabs(h.integrate (ONES));
- const double SUM = v[0] + v[1] + v[2] + v[3]+v[4] + v[5] + v[6] + v[7];
- const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]) + fabs (v[3]) +
- fabs(v[4]) + fabs(v[5]) + fabs (v[6]) + fabs (v[7]);
- if (SUMABS)
- {
- const double XI = SUM / SUMABS;
- res += area * (1 - XI) * 0.5 ;
- }
- }
- }
- }
- else if(nbNod == 6){
- prism p(x, y, z);
- if(nbComp == 1) res += p.integrate(v);
- }
- else if(nbNod == 5){
- pyramid p(x, y, z);
- if(nbComp == 1) res += p.integrate(v);
- }
+ }
+ else if(nbNod == 5){
+ pyramid p(x, y, z);
+ if(nbComp == 1){
+ if(!levelsetPositive) res += p.integrate(v);
+ else res += p.integrateLevelsetPositive(v);
}
+ }
}
-// printf("integration res = %22.15E\n",res);
- return res;
+ }
+ // printf("integration res = %22.15E\n",res);
+ return res;
}
+
Post_View *GMSH_IntegratePlugin::execute(Post_View * v)
{
- int iView = (int)IntegrateOptions_Number[0].def;
-
- if(iView < 0)
- iView = v ? v->Index : 0;
-
- if(!List_Pointer_Test(CTX.post.list, iView)) {
- Msg(GERROR, "View[%d] does not exist", iView);
- return v;
- }
-
- Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
- Post_View *v2 = BeginView(1);
-
- double x = (v1->BBox[0]+v1->BBox[1])/2.;
- double y = (v1->BBox[2]+v1->BBox[3])/2.;
- double z = (v1->BBox[4]+v1->BBox[5])/2.;
- List_Add(v2->SP, &x);
- List_Add(v2->SP, &y);
- List_Add(v2->SP, &z);
- for(int ts = 0; ts < v1->NbTimeStep; ts++){
- double val = 0;
-// scalar fields
- val += integrate(v1->NbSP, v1->SP, 0, 1, 1, ts);
- val += integrate(v1->NbSL, v1->SL, 1, 2, 1, ts);
- val += integrate(v1->NbST, v1->ST, 2, 3, 1, ts);
- val += integrate(v1->NbSQ, v1->SQ, 2, 4, 1, ts);
- val += integrate(v1->NbSS, v1->SS, 3, 4, 1, ts);
- val += integrate(v1->NbSH, v1->SH, 3, 8, 1, ts);
- val += integrate(v1->NbSI, v1->SI, 3, 6, 1, ts);
- val += integrate(v1->NbSY, v1->SY, 3, 5, 1, ts);
-// circulations
- val += integrate(v1->NbVL, v1->VL, 1, 2, 3, ts);
-// fluxes
- val += integrate(v1->NbVT, v1->VT, 2, 3, 3, ts);
- val += integrate(v1->NbVQ, v1->VQ, 2, 4, 3, ts);
- Msg(INFO, "Step %d: integral = %.16g", ts, val);
- List_Add(v2->SP, &val);
- }
- v2->NbSP = 1;
- v2->IntervalsType = DRAW_POST_NUMERIC;
-
-// copy time data
- for(int i = 0; i < List_Nbr(v1->Time); i++)
- List_Add(v2->Time, List_Pointer(v1->Time, i));
-// finalize
- char name[1024], filename[1024];
- sprintf(name, "%s_Integrate", v1->Name);
- sprintf(filename, "%s_Integrate.pos", v1->Name);
- EndView(v2, 1, filename, name);
-
- return v2;
+ int iView = (int)IntegrateOptions_Number[0].def;
+
+ if(iView < 0)
+ iView = v ? v->Index : 0;
+
+ if(!List_Pointer_Test(CTX.post.list, iView)) {
+ Msg(GERROR, "View[%d] does not exist", iView);
+ return v;
+ }
+
+ Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
+ Post_View *v2 = BeginView(1);
+
+ double x = (v1->BBox[0]+v1->BBox[1])/2.;
+ double y = (v1->BBox[2]+v1->BBox[3])/2.;
+ double z = (v1->BBox[4]+v1->BBox[5])/2.;
+ List_Add(v2->SP, &x);
+ List_Add(v2->SP, &y);
+ List_Add(v2->SP, &z);
+ for(int ts = 0; ts < v1->NbTimeStep; ts++){
+ double val = 0;
+ // scalar fields
+ val += integrate(v1->NbSP, v1->SP, 0, 1, 1, ts);
+ val += integrate(v1->NbSL, v1->SL, 1, 2, 1, ts);
+ val += integrate(v1->NbST, v1->ST, 2, 3, 1, ts);
+ val += integrate(v1->NbSQ, v1->SQ, 2, 4, 1, ts);
+ val += integrate(v1->NbSS, v1->SS, 3, 4, 1, ts);
+ val += integrate(v1->NbSH, v1->SH, 3, 8, 1, ts);
+ val += integrate(v1->NbSI, v1->SI, 3, 6, 1, ts);
+ val += integrate(v1->NbSY, v1->SY, 3, 5, 1, ts);
+ // circulations
+ val += integrate(v1->NbVL, v1->VL, 1, 2, 3, ts);
+ // fluxes
+ val += integrate(v1->NbVT, v1->VT, 2, 3, 3, ts);
+ val += integrate(v1->NbVQ, v1->VQ, 2, 4, 3, ts);
+ Msg(INFO, "Step %d: integral = %.16g", ts, val);
+ List_Add(v2->SP, &val);
+ }
+ v2->NbSP = 1;
+ v2->IntervalsType = DRAW_POST_NUMERIC;
+
+ // copy time data
+ for(int i = 0; i < List_Nbr(v1->Time); i++)
+ List_Add(v2->Time, List_Pointer(v1->Time, i));
+ // finalize
+ char name[1024], filename[1024];
+ sprintf(name, "%s_Integrate", v1->Name);
+ sprintf(filename, "%s_Integrate.pos", v1->Name);
+ EndView(v2, 1, filename, name);
+
+ return v2;
}
diff --git a/Plugin/Lambda2.cpp b/Plugin/Lambda2.cpp
index 2be82101660e858572ee46e85edf20cb14431b2b..056c602634e537381300ed013035e2b660c18418 100644
--- a/Plugin/Lambda2.cpp
+++ b/Plugin/Lambda2.cpp
@@ -1,4 +1,4 @@
-// $Id: Lambda2.cpp,v 1.1 2004-11-26 10:31:54 remacle Exp $
+// $Id: Lambda2.cpp,v 1.2 2004-11-26 16:16:39 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -30,10 +30,6 @@
#include <math.h>
#include <stdio.h>
-#if defined(HAVE_MATH_EVAL)
-#include "matheval.h"
-#endif
-
extern Context_T CTX;
StringXNumber Lambda2Options_Number[] = {
@@ -41,7 +37,6 @@ StringXNumber Lambda2Options_Number[] = {
{GMSH_FULLRC, "iView", NULL, -1.}
};
-
extern "C"
{
GMSH_Plugin *GMSH_RegisterLambda2Plugin()
@@ -50,7 +45,6 @@ extern "C"
}
}
-
GMSH_Lambda2Plugin::GMSH_Lambda2Plugin()
{
;
@@ -79,7 +73,7 @@ void GMSH_Lambda2Plugin::getInfos(char *author, char *copyright,
"If `iView' < 0, the plugin is run\n"
"on the current view.\n"
"\n"
- "Plugin(Lambda2) is executed in-place.\n");
+ "Plugin(Lambda2) creates one noew view.\n");
}
int GMSH_Lambda2Plugin::getNbOptions() const
@@ -97,45 +91,44 @@ void GMSH_Lambda2Plugin::catchErrorMessage(char *errorMessage) const
strcpy(errorMessage, "Lambda2 failed...");
}
-//-----------------------------------------------------------------------------
static void eigen(List_T *inList, int inNb,
List_T *outListScalar, int *outNbScalar,
int nbTime, int nbNod, int nbComp, int lam)
{
-
if(!inNb)
return;
int outNbComp, *outNb;
List_T *outList;
-
+
outNbComp = 1;
outNb = outNbScalar;
outList = outListScalar;
-
+
int nb = List_Nbr(inList) / inNb;
-
+
//boucle sur les elements
//-----------------------
for(int i = 0; i < List_Nbr(inList); i += nb) {
-
+
for(int j = 0; j < 3 * nbNod; j++)
List_Add(outList, List_Pointer_Fast(inList, i + j));
-
+
for(int j = 0; j < nbTime; j++){
double *x = (double *)List_Pointer_Fast(inList, i);
double *y = (double *)List_Pointer_Fast(inList, i + nbNod);
double *z = (double *)List_Pointer_Fast(inList, i + 2 * nbNod);
+ double *v = (double *)List_Pointer_Fast(inList, i + 3 * nbNod +
+ nbNod * nbComp * j + nbComp * 0);
double GradVel[3][3];
- double *v = (double *)List_Pointer_Fast(inList, i + 3 * nbNod + nbNod * nbComp * j + nbComp * 0);
GradVel[0][0] = v[0]; GradVel[0][1] = v[1]; GradVel[0][2] = v[2];
GradVel[1][0] = v[3]; GradVel[1][1] = v[4]; GradVel[1][2] = v[5];
GradVel[2][0] = v[6]; GradVel[2][1] = v[7]; GradVel[2][2] = v[8];
-
+
//par element
//------------
-
+
//calcul partie sym et antisymetrique
double sym[3][3];
double asym[3][3];
@@ -156,7 +149,7 @@ static void eigen(List_T *inList, int inNb,
//calcul de lambda avec tri
double lambda[3];
eigenvalue(a, lambda);
-
+
double res;
for(int k = 0; k < nbNod; k++){
for(int l = 0; l < outNbComp; l++){
@@ -167,46 +160,34 @@ static void eigen(List_T *inList, int inNb,
}
(*outNb)++;
}
-
}
-//-----------------------------------------------------------------------------
Post_View *GMSH_Lambda2Plugin::execute(Post_View * v)
{
int ev = (int)Lambda2Options_Number[0].def;
int iView = (int)Lambda2Options_Number[1].def;
-
+
if(iView < 0)
iView = v ? v->Index : 0;
-
+
if(!List_Pointer_Test(CTX.post.list, iView)) {
Msg(GERROR, "View[%d] does not exist", iView);
return v;
}
+ Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
Post_View *v2 = BeginView(1);
- Post_View *v1 = (Post_View*)List_Pointer(CTX.post.list, iView);
-
eigen(v1->TT, v1->NbTT, v2->ST, &v2->NbSS, v1->NbTimeStep, 3, 9, ev);
-
eigen(v1->TS, v1->NbTS, v2->SS, &v2->NbSS, v1->NbTimeStep, 4, 9, ev);
-
- if(v2->empty()) {
- RemoveViewByNumber(v2->Num);
- return v1;
- }
- else{
- // copy time data
- for(int i = 0; i < List_Nbr(v1->Time); i++)
- List_Add(v2->Time, List_Pointer(v1->Time, i));
- // finalize
- char name[1024], filename[1024];
- sprintf(name, "%s_Extract", v1->Name);
- sprintf(filename, "%s_Extract.pos", v1->Name);
- EndView(v2, 1, filename, name);
- return v2;
-
- }
+
+ // copy time data
+ for(int i = 0; i < List_Nbr(v1->Time); i++)
+ List_Add(v2->Time, List_Pointer(v1->Time, i));
+ // finalize
+ char name[1024], filename[1024];
+ sprintf(name, "%s_Extract", v1->Name);
+ sprintf(filename, "%s_Extract.pos", v1->Name);
+ EndView(v2, 1, filename, name);
+ return v2;
}
-//------------------------------------------------------------------
diff --git a/Plugin/ShapeFunctions.h b/Plugin/ShapeFunctions.h
index 68b1d1a3dcc906769793bfd260138a2088a6e2c0..8e5c246ce79b32533cd873875ee27aa1a28120e9 100644
--- a/Plugin/ShapeFunctions.h
+++ b/Plugin/ShapeFunctions.h
@@ -93,6 +93,20 @@ public:
}
return sum;
}
+ double integrateLevelsetPositive(double val[])
+ {
+ double ones[8] = {1., 1., 1., 1., 1., 1., 1., 1.}; // FIXME: 8-node max
+ double area = integrate(ones);
+ double sum = 0, sumabs = 0.;
+ for(int i = 0; i < getNumNodes(); i++){
+ sum += val[i];
+ sumabs += fabs(val[i]);
+ }
+ double res = 0.;
+ if(sumabs)
+ res = area * (1 - sum/sumabs) * 0.5 ;
+ return res;
+ }
};
class point : public element{