Skip to content
Snippets Groups Projects
Commit cd1204ce authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

- fixed CTX.post.list error in new plugins

- cleanup integrateLevelsetPositive in Integrate.cpp

- re-indent
parent 40ec0d66
No related branches found
No related tags found
No related merge requests found
// $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,33 +70,22 @@ 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)
{
if(!inNb)
return;
......@@ -129,14 +102,11 @@ static void computeGradient( List_T *inList, int inNb,
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)
{
if(yy[0] > 0){
for(int j = 0; j < 3 * nbNod; j++)
List_Add(outList, List_Pointer_Fast(inList, i + j));
......@@ -149,13 +119,13 @@ static void computeGradient( List_T *inList, int inNb,
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);
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];
......@@ -250,10 +220,8 @@ static void computeGradient( List_T *inList, int inNb,
}
}
//-----------------------------------------------------------------------------
Post_View *GMSH_GradientPlugin::execute(Post_View * v)
{
int iView = (int)GradientOptions_Number[0].def;
if(iView < 0)
......@@ -264,9 +232,9 @@ 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);
Post_View *v2 = BeginView(1);
// points
//--------
......@@ -278,19 +246,16 @@ Post_View *GMSH_GradientPlugin::execute(Post_View * v)
// 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
//-----------
......@@ -314,11 +279,6 @@ 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));
......@@ -329,8 +289,3 @@ computeGradient( v1->TY, v1->NbTY, v2->VY, &v2->NbVY, v1->NbTimeStep, 5, 9);
EndView(v2, 1, filename, name);
return v2;
}
}
//----------------------------------------------------------------------
......@@ -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 *);
};
......
// $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
//
......@@ -103,13 +103,19 @@ static double integrate(int nbList, List_T *list, int dim,
if(dim == 0){
if(nbNod == 1){
point p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
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) res += l.integrate(v);
if(nbComp == 1){
if(!levelsetPositive) res += l.integrate(v);
else res += l.integrateLevelsetPositive(v);
}
else if(nbComp == 3) res += l.integrateCirculation(v);
}
}
......@@ -117,97 +123,55 @@ static double integrate(int nbList, List_T *list, int dim,
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 ;
}
}
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
{
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 ;
}
}
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)
{
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 ;
}
if(nbComp == 1){
if(!levelsetPositive) res += t.integrate(v);
else res += t.integrateLevelsetPositive(v);
}
}
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 ;
}
}
if(nbComp == 1){
if(!levelsetPositive) res += h.integrate(v);
else res += h.integrateLevelsetPositive(v);
}
}
else if(nbNod == 6){
prism p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
if(nbComp == 1){
if(!levelsetPositive) res += p.integrate(v);
else res += p.integrateLevelsetPositive(v);
}
}
else if(nbNod == 5){
pyramid p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
if(nbComp == 1){
if(!levelsetPositive) res += p.integrate(v);
else res += p.integrateLevelsetPositive(v);
}
}
}
}
// printf("integration res = %22.15E\n",res);
return res;
}
Post_View *GMSH_IntegratePlugin::execute(Post_View * v)
{
int iView = (int)IntegrateOptions_Number[0].def;
......
// $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,12 +91,10 @@ 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;
......@@ -127,8 +119,9 @@ static void eigen(List_T *inList, int inNb,
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];
......@@ -167,10 +160,8 @@ static void eigen(List_T *inList, int inNb,
}
(*outNb)++;
}
}
//-----------------------------------------------------------------------------
Post_View *GMSH_Lambda2Plugin::execute(Post_View * v)
{
int ev = (int)Lambda2Options_Number[0].def;
......@@ -184,19 +175,12 @@ Post_View *GMSH_Lambda2Plugin::execute(Post_View * v)
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));
......@@ -206,7 +190,4 @@ Post_View *GMSH_Lambda2Plugin::execute(Post_View * v)
sprintf(filename, "%s_Extract.pos", v1->Name);
EndView(v2, 1, filename, name);
return v2;
}
}
//------------------------------------------------------------------
......@@ -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{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment