From c36677f3feb0b6f2fce150f87a137e7219f9b444 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 26 Nov 2004 10:31:54 +0000
Subject: [PATCH] *** empty log message ***
---
Plugin/Gradient.cpp | 336 +++++++++++++++++++++++++++++++++++++++++
Plugin/Gradient.h | 45 ++++++
Plugin/Integrate.cpp | 347 +++++++++++++++++++++++--------------------
Plugin/Lambda2.cpp | 212 ++++++++++++++++++++++++++
Plugin/Lambda2.h | 52 +++++++
5 files changed, 835 insertions(+), 157 deletions(-)
create mode 100644 Plugin/Gradient.cpp
create mode 100644 Plugin/Gradient.h
create mode 100644 Plugin/Lambda2.cpp
create mode 100644 Plugin/Lambda2.h
diff --git a/Plugin/Gradient.cpp b/Plugin/Gradient.cpp
new file mode 100644
index 0000000000..75655036dd
--- /dev/null
+++ b/Plugin/Gradient.cpp
@@ -0,0 +1,336 @@
+// $Id: Gradient.cpp,v 1.1 2004-11-26 10:31:54 remacle Exp $
+//
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "Gradient.h"
+#include "List.h"
+#include "Views.h"
+#include "Context.h"
+#include "Numeric.h"
+#include <algorithm>
+#include <functional>
+#include <math.h>
+#include <stdio.h>
+
+#if defined(HAVE_MATH_EVAL)
+#include "matheval.h"
+#endif
+
+extern Context_T CTX;
+
+StringXNumber GradientOptions_Number[] = {
+ {GMSH_FULLRC, "iView", NULL, -1.}
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterGradientPlugin()
+ {
+ return new GMSH_GradientPlugin();
+ }
+}
+
+int GMSH_GradientPlugin::getNbOptionsStr() const
+{
+ return 0;
+}
+
+int GMSH_GradientPlugin::getNbOptions() const
+{
+ return sizeof(GradientOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_GradientPlugin::getOption(int iopt)
+{
+ return &GradientOptions_Number[iopt];
+}
+
+StringXString *GMSH_GradientPlugin::getOptionStr(int iopt)
+{
+ throw;
+}
+
+
+
+GMSH_GradientPlugin::GMSH_GradientPlugin()
+{
+ ;
+}
+
+void GMSH_GradientPlugin::getName(char *name) const
+{
+ strcpy(name, "Gradient");
+}
+
+void GMSH_GradientPlugin::getInfos(char *author, char *copyright,
+ char *help_text) const
+{
+ 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"
+ "If `iView' < 0, the plugin is run\n"
+ "on the current view.\n"
+ "\n"
+ "Plugin(Lambda2) is executed in-place.\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)
+{
+
+ if(!inNb)
+ return;
+
+ int outNbComp, *outNb;
+ List_T *outList;
+
+ outNbComp = 1;
+ outNb = outNbVector;
+ outList = outListVector;
+
+ int MAX_NOD = 4;
+ int MAX_COMP= 3;
+
+ int nb = List_Nbr(inList) / inNb;
+
+ //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));
+
+ 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;
+ }
+
+ //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];
+ }
+ }
+ }
+
+ 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)++;
+ }
+ }
+}
+
+//-----------------------------------------------------------------------------
+Post_View *GMSH_GradientPlugin::execute(Post_View * v)
+{
+
+ int iView = (int)GradientOptions_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 *v2 = BeginView(1);
+
+ 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);
+
+// 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);
+ 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;
+ }
+}
+
+
+
+//----------------------------------------------------------------------
diff --git a/Plugin/Gradient.h b/Plugin/Gradient.h
new file mode 100644
index 0000000000..19ec8b91d7
--- /dev/null
+++ b/Plugin/Gradient.h
@@ -0,0 +1,45 @@
+#ifndef _GRADIENT_H_
+#define _GRADIENT_H
+
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "Numeric.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterGradientPlugin();
+}
+
+class GMSH_GradientPlugin : public GMSH_Post_Plugin
+{
+public:
+ GMSH_GradientPlugin();
+ void getName(char *name) const;
+ void getInfos(char *author, char *copyright, char *helpText) const;
+ void catchErrorMessage(char *errorMessage) const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+ int getNbOptionsStr() const;
+ StringXString* getOptionStr(int iopt);
+ Post_View *execute(Post_View *);
+};
+
+#endif
diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp
index 2784c9258d..8601528c2b 100644
--- a/Plugin/Integrate.cpp
+++ b/Plugin/Integrate.cpp
@@ -1,4 +1,4 @@
-// $Id: Integrate.cpp,v 1.6 2004-11-25 22:07:53 geuzaine Exp $
+// $Id: Integrate.cpp,v 1.7 2004-11-26 10:29:20 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -30,203 +30,236 @@
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
+ 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;
+ const int levelsetPositive = (int)IntegrateOptions_Number[1].def;
- double res = 0.;
- double res2 = 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);
- }
- }
- 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(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 ;
+ 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);
}
- }
}
- else if(nbComp == 3) res += t.integrateFlux(v);
- }
- else if(nbNod == 4){
- quadrangle q(x, y, z);
- if(nbComp == 1) res += q.integrate(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{
- 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]);
- if(SUMABS){
- const double XI = SUM / SUMABS;
- res += area * (1 - XI) * 0.5 ;
- res2 += area * (1 + XI) * 0.5 ;
+ 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) res += h.integrate(v);
- }
- 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(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){
+ const double XI = SUM / SUMABS;
+ res += area * (1 - XI) * 0.5 ;
+ }
+ }
+ }
+ 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
+ {
+ 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);
+ }
+ }
+
+// 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
new file mode 100644
index 0000000000..2be8210166
--- /dev/null
+++ b/Plugin/Lambda2.cpp
@@ -0,0 +1,212 @@
+// $Id: Lambda2.cpp,v 1.1 2004-11-26 10:31:54 remacle Exp $
+//
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "Lambda2.h"
+#include "List.h"
+#include "Views.h"
+#include "Context.h"
+#include "Numeric.h"
+#include <algorithm>
+#include <functional>
+#include <math.h>
+#include <stdio.h>
+
+#if defined(HAVE_MATH_EVAL)
+#include "matheval.h"
+#endif
+
+extern Context_T CTX;
+
+StringXNumber Lambda2Options_Number[] = {
+ {GMSH_FULLRC, "Eigenvalue", NULL, 2.},
+ {GMSH_FULLRC, "iView", NULL, -1.}
+};
+
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterLambda2Plugin()
+ {
+ return new GMSH_Lambda2Plugin();
+ }
+}
+
+
+GMSH_Lambda2Plugin::GMSH_Lambda2Plugin()
+{
+ ;
+}
+
+void GMSH_Lambda2Plugin::getName(char *name) const
+{
+ strcpy(name, "Lambda2");
+}
+
+void GMSH_Lambda2Plugin::getInfos(char *author, char *copyright,
+ char *help_text) const
+{
+ strcpy(author, "E. Marchandise");
+ strcpy(copyright, "DGR (www.multiphysics.com)");
+ strcpy(help_text,
+ "Plugin(Lambda2) 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"
+ "If `iView' < 0, the plugin is run\n"
+ "on the current view.\n"
+ "\n"
+ "Plugin(Lambda2) is executed in-place.\n");
+}
+
+int GMSH_Lambda2Plugin::getNbOptions() const
+{
+ return sizeof(Lambda2Options_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_Lambda2Plugin::getOption(int iopt)
+{
+ return &Lambda2Options_Number[iopt];
+}
+
+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 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];
+ for(int i=0; i<3; i++){
+ for(int j=0; j<3; j++){
+ sym[i][j] = 0.5*(GradVel[i][j]+GradVel[j][i]);
+ asym[i][j] = 0.5*(GradVel[i][j]-GradVel[j][i]);
+ }
+ }
+ 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);
+
+ double res;
+ for(int k = 0; k < nbNod; k++){
+ for(int l = 0; l < outNbComp; l++){
+ res = lambda[lam-1];
+ List_Add(outList, &res);
+ }
+ }
+ }
+ (*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 *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;
+
+ }
+}
+//------------------------------------------------------------------
diff --git a/Plugin/Lambda2.h b/Plugin/Lambda2.h
new file mode 100644
index 0000000000..f7d9f93174
--- /dev/null
+++ b/Plugin/Lambda2.h
@@ -0,0 +1,52 @@
+#ifndef _LAMBDA2_H_
+#define _LAMBDA2_H
+
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "List.h"
+#include "Numeric.h"
+#include <algorithm>
+#include <functional>
+#include <math.h>
+#include <stdio.h>
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterLambda2Plugin();
+}
+
+class GMSH_Lambda2Plugin : public GMSH_Post_Plugin
+{
+public:
+ GMSH_Lambda2Plugin();
+ void getName(char *name) const;
+ void getInfos(char *author, char *copyright, char *helpText) const;
+ void catchErrorMessage(char *errorMessage) const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+
+ Post_View *execute(Post_View *);
+
+
+
+//---------------------------------------------------------
+};
+#endif
--
GitLab