diff --git a/Plugin/CutMap.cpp b/Plugin/CutMap.cpp index 2361fee279680ba408bdfd5609e69d9fe0620ec0..83d9bf971efe06d98f12b7491dc9052be10877aa 100644 --- a/Plugin/CutMap.cpp +++ b/Plugin/CutMap.cpp @@ -1,4 +1,4 @@ -// $Id: CutMap.cpp,v 1.21 2001-08-09 20:53:23 geuzaine Exp $ +// $Id: CutMap.cpp,v 1.22 2001-08-11 23:25:50 geuzaine Exp $ #include "CutMap.h" #include "List.h" @@ -68,6 +68,7 @@ Post_View *GMSH_CutMapPlugin::execute (Post_View *v) int iView = (int)CutMapOptions_Number[1].def; _ith_field_to_draw_on_the_iso = (int)CutMapOptions_Number[2].def; + _orientation = ORIENT_MAP; if(v && iView < 0) vv = v; diff --git a/Plugin/CutPlane.cpp b/Plugin/CutPlane.cpp index 825bc61209fbbb8a33e6f3a54cc8db415c476302..468231196ca777a5c132eb9380cb38b25903d367 100644 --- a/Plugin/CutPlane.cpp +++ b/Plugin/CutPlane.cpp @@ -1,4 +1,4 @@ -// $Id: CutPlane.cpp,v 1.17 2001-08-09 20:53:23 geuzaine Exp $ +// $Id: CutPlane.cpp,v 1.18 2001-08-11 23:25:50 geuzaine Exp $ #include "CutPlane.h" #include "List.h" @@ -69,6 +69,10 @@ Post_View *GMSH_CutPlanePlugin::execute (Post_View *v) Post_View *vv; int iView = (int)CutPlaneOptions_Number[4].def; + _orientation = ORIENT_PLANE; + _ref[0] = CutPlaneOptions_Number[0].def; + _ref[1] = CutPlaneOptions_Number[1].def; + _ref[2] = CutPlaneOptions_Number[2].def; if(v && iView < 0) vv = v; diff --git a/Plugin/CutSphere.cpp b/Plugin/CutSphere.cpp index bd60ccd59ccd06f129aada9fc26dff57dc4c8577..493b3c186a969f6ad28f39678964000b45d13113 100644 --- a/Plugin/CutSphere.cpp +++ b/Plugin/CutSphere.cpp @@ -1,4 +1,4 @@ -// $Id: CutSphere.cpp,v 1.16 2001-08-09 20:53:23 geuzaine Exp $ +// $Id: CutSphere.cpp,v 1.17 2001-08-11 23:25:50 geuzaine Exp $ #include <string.h> #include "CutSphere.h" @@ -69,7 +69,12 @@ extern List_T *Post_ViewList; Post_View *GMSH_CutSpherePlugin::execute (Post_View *v) { Post_View *vv; + int iView = (int)CutSphereOptions_Number[4].def; + _orientation = ORIENT_SPHERE; + _ref[0] = CutSphereOptions_Number[0].def; + _ref[1] = CutSphereOptions_Number[1].def; + _ref[2] = CutSphereOptions_Number[2].def; if(v && iView < 0) vv = v; diff --git a/Plugin/LevelsetPlugin.cpp b/Plugin/LevelsetPlugin.cpp index 5d50855246310d6802f95b8a70327778c60bfaf0..b4cabc8c0535c9d25bf7f3e7eb8365367c8c9f8d 100644 --- a/Plugin/LevelsetPlugin.cpp +++ b/Plugin/LevelsetPlugin.cpp @@ -1,20 +1,17 @@ -// $Id: LevelsetPlugin.cpp,v 1.15 2001-08-09 22:05:51 remacle Exp $ +// $Id: LevelsetPlugin.cpp,v 1.16 2001-08-11 23:25:50 geuzaine Exp $ #include "LevelsetPlugin.h" #include "List.h" #include "Views.h" #include "Iso.h" - -/// Includes are basdly designed, i prefer forward decls. -void prodve (double a[3], double b[3], double c[3]); -void prosca (double a[3], double b[3], double *c); -void norme (double a[3]); -int sys3x3 (double mat[3][3], double b[3], double res[3], double *det); +#include "Numeric.h" GMSH_LevelsetPlugin::GMSH_LevelsetPlugin() { processed = 0; _ith_field_to_draw_on_the_iso = 0; + _orientation = ORIENT_NONE; + _ref[0] = _ref[1] = _ref[2] = 0.; strcpy (OutputFileName,"levelset.pos"); } @@ -39,122 +36,130 @@ Post_View *GMSH_LevelsetPlugin::execute (Post_View *v) This plugin creates a new view which is the result of a cut of the actual view with a levelset. */ - int k,i,nb,edtet[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; - double Xpi[6],Ypi[6],Zpi[6],myValsi[6]; + int k,i,j,nb,edtet[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; + double *X, *Y, *Z, *Vals, levels[6], coef; + double Xp[6], Yp[6], Zp[6], myVals[6]; + double Xpi[6], Ypi[6], Zpi[6], myValsi[6]; Post_View *View; - // for all scalar simplices - - if(v->NbSS) - { - View = BeginView(1); - nb = List_Nbr(v->SS) / v->NbSS ; - for(i = 0 ; i < List_Nbr(v->SS) ; i+=nb) - { - double levels[6],Xp[6],Yp[6],Zp[6],myVals[6]; - double *X = (double*)List_Pointer_Fast(v->SS,i); - double *Y = (double*)List_Pointer_Fast(v->SS,i+4); - double *Z = (double*)List_Pointer_Fast(v->SS,i+8); - double *VAL = (double*)List_Pointer_Fast(v->SS,i+12); - for(int j=0;j<4;j++)levels[j] = levelset(X[j],Y[j],Z[j],VAL[j]); - int nx = 0; - for(k=0;k<6;k++) - { - if(levels[edtet[k][0]] * levels[edtet[k][1]] <= 0.0) - { - double coef = InterpolateIso(X,Y,Z,levels,0.0, - edtet[k][0],edtet[k][1], - &Xp[nx],&Yp[nx],&Zp[nx]); - myVals[nx] = what_to_draw (Xp[nx],Yp[nx],Zp[nx],edtet[k][0],edtet[k][1],coef,VAL); - nx++; - } - } - - if(nx == 4) - { - double xx = Xp[3]; - double yy = Yp[3]; - double zz = Zp[3]; - double vv = myVals[3]; - Xp[3] = Xp[2]; - Yp[3] = Yp[2]; - Zp[3] = Zp[2]; - myVals[3] = myVals[2]; - Xp[2] = xx; - Yp[2] = yy; - Zp[2] = zz; - myVals[2] = vv; - } - - double v1[3] = {Xp[2]-Xp[0],Yp[2]-Yp[0],Zp[2]-Zp[0]}; - double v2[3] = {Xp[1]-Xp[0],Yp[1]-Yp[0],Zp[1]-Zp[0]}; - double gr[3]; - double n[3],xx; - prodve(v1,v2,n); - norme(n); - gradSimplex(X,Y,Z,VAL,gr); - prosca(gr,n,&xx); - - if(n[0] + n[1] + n[2] > 0){ - for(k=0;k<nx;k++){ - Xpi[k] = Xp[k]; - Ypi[k] = Yp[k]; - Zpi[k] = Zp[k]; - myValsi[k] = myVals[k]; - } - for(k=0;k<nx;k++){ - Xp[k] = Xpi[nx-k-1]; - Yp[k] = Ypi[nx-k-1]; - Zp[k] = Zpi[nx-k-1]; - myVals[k] = myValsi[nx-k-1]; - } - } - - if(nx == 3 || nx == 4) - { - for(k=0;k<3;k++)List_Add(View->ST, &Xp[k]); - for(k=0;k<3;k++)List_Add(View->ST, &Yp[k]); - for(k=0;k<3;k++)List_Add(View->ST, &Zp[k]); - for(k=0;k<3;k++)List_Add(View->ST, &myVals[k]); - View->NbST++; - } - if(nx == 4) - { - for(k=2;k<5;k++)List_Add(View->ST, &Xp[k %4]); - for(k=2;k<5;k++)List_Add(View->ST, &Yp[k % 4]); - for(k=2;k<5;k++)List_Add(View->ST, &Zp[k % 4]); - for(k=2;k<5;k++)List_Add(View->ST, &myVals[k %4]); - View->NbST++; - } + // for all scalar tets + + if(v->NbSS){ + View = BeginView(1); + nb = List_Nbr(v->SS) / v->NbSS ; + for(i=0 ; i<List_Nbr(v->SS) ; i+=nb){ + X = (double*)List_Pointer_Fast(v->SS,i); + Y = (double*)List_Pointer_Fast(v->SS,i+4); + Z = (double*)List_Pointer_Fast(v->SS,i+8); + Vals = (double*)List_Pointer_Fast(v->SS,i+12); + for(j=0 ; j<4 ; j++) + levels[j] = levelset(X[j],Y[j],Z[j],Vals[j]); + int nx = 0; + for(k=0 ; k<6 ; k++){ + if(levels[edtet[k][0]] * levels[edtet[k][1]] <= 0.0){ + coef = InterpolateIso(X,Y,Z,levels,0.0, + edtet[k][0],edtet[k][1], + &Xp[nx],&Yp[nx],&Zp[nx]); + myVals[nx] = what_to_draw (Xp[nx],Yp[nx],Zp[nx], + edtet[k][0],edtet[k][1],coef,Vals); + nx++; } - char name[1024],filename[1024]; - - sprintf(name,"cut-%s",v->Name); - sprintf(filename,"cut-%s",v->FileName); - EndView(1, filename, name); + } + + if(nx == 4){ + double xx = Xp[3]; + double yy = Yp[3]; + double zz = Zp[3]; + double vv = myVals[3]; + Xp[3] = Xp[2]; + Yp[3] = Yp[2]; + Zp[3] = Zp[2]; + myVals[3] = myVals[2]; + Xp[2] = xx; + Yp[2] = yy; + Zp[2] = zz; + myVals[2] = vv; + } + + double v1[3] = {Xp[2]-Xp[0],Yp[2]-Yp[0],Zp[2]-Zp[0]}; + double v2[3] = {Xp[1]-Xp[0],Yp[1]-Yp[0],Zp[1]-Zp[0]}; + double gr[3]; + double n[3],test; + prodve(v1,v2,n); + + switch(_orientation){ + case ORIENT_MAP: + gradSimplex(X,Y,Z,Vals,gr); + prosca(gr,n,&test); + break; + case ORIENT_PLANE: + prosca(n,_ref,&test); + break; + case ORIENT_SPHERE: + gr[0] = _ref[0]-Xp[0]; + gr[1] = _ref[1]-Yp[0]; + gr[2] = _ref[2]-Zp[0]; + prosca(gr,n,&test); + break; + default: + test = 0.; + break; + } - Msg(INFO, "Created view '%s' (%d triangles)", name, View->NbST); - processed = View; - return View; + if(test<0){ + for(k=0;k<nx;k++){ + Xpi[k] = Xp[k]; + Ypi[k] = Yp[k]; + Zpi[k] = Zp[k]; + myValsi[k] = myVals[k]; + } + for(k=0;k<nx;k++){ + Xp[k] = Xpi[nx-k-1]; + Yp[k] = Ypi[nx-k-1]; + Zp[k] = Zpi[nx-k-1]; + myVals[k] = myValsi[nx-k-1]; + } + } + + if(nx == 3 || nx == 4){ + for(k=0 ; k<3 ; k++) List_Add(View->ST, &Xp[k]); + for(k=0 ; k<3 ; k++) List_Add(View->ST, &Yp[k]); + for(k=0 ; k<3 ; k++) List_Add(View->ST, &Zp[k]); + for(k=0 ; k<3 ; k++) List_Add(View->ST, &myVals[k]); + View->NbST++; + } + if(nx == 4){ + for(k=2 ; k<5 ; k++) List_Add(View->ST, &Xp[k % 4]); + for(k=2 ; k<5 ; k++) List_Add(View->ST, &Yp[k % 4]); + for(k=2 ; k<5 ; k++) List_Add(View->ST, &Zp[k % 4]); + for(k=2 ; k<5 ; k++) List_Add(View->ST, &myVals[k % 4]); + View->NbST++; + } } + char name[1024],filename[1024]; + sprintf(name,"cut-%s",v->Name); + sprintf(filename,"cut-%s",v->FileName); + EndView(View, 1, filename, name); + + Msg(INFO, "Created view '%s' (%d triangles)", name, View->NbST); + processed = View; + return View; + } + return 0; } -double GMSH_LevelsetPlugin::what_to_draw (double x, - double y, - double z, - int p1, - int p2, - double coef, - double *VAL) const +double GMSH_LevelsetPlugin::what_to_draw (double x, double y, double z, + int p1, int p2, + double coef, double *Vals) const { int offset = _ith_field_to_draw_on_the_iso * 4; // TEST JF, this would draw y coord on the iso // return y; p2 += offset; p1 += offset; - return coef * (VAL[p2] - VAL[p1]) + VAL[p1]; + return coef * (Vals[p2] - Vals[p1]) + Vals[p1]; } diff --git a/Plugin/LevelsetPlugin.h b/Plugin/LevelsetPlugin.h index f25ec33104c1fd52415d713af806e54a67e6c8d0..b2f5f16322f9eae6ddc5d93f68ff344f37642535 100644 --- a/Plugin/LevelsetPlugin.h +++ b/Plugin/LevelsetPlugin.h @@ -2,10 +2,17 @@ #define _LEVELSETPLUGIN_H_ #include "Plugin.h" +#define ORIENT_NONE 0 +#define ORIENT_MAP 1 +#define ORIENT_PLANE 2 +#define ORIENT_SPHERE 3 + class GMSH_LevelsetPlugin : public GMSH_Post_Plugin { protected: int _ith_field_to_draw_on_the_iso; + int _orientation; + double _ref[3]; private: virtual double levelset (double x, double y, double z, double val) const = 0; virtual double what_to_draw (double x, double y, double z, diff --git a/Plugin/Skin.cpp b/Plugin/Skin.cpp index 32e2e2ec2966a79f421e5fcaa68523167ce118b8..045def54c6c03abf5b5619b61c3930ba777fe305 100644 --- a/Plugin/Skin.cpp +++ b/Plugin/Skin.cpp @@ -1,4 +1,4 @@ -// $Id: Skin.cpp,v 1.7 2001-08-06 12:26:26 geuzaine Exp $ +// $Id: Skin.cpp,v 1.8 2001-08-11 23:25:50 geuzaine Exp $ #include "Plugin.h" #include "Skin.h" @@ -6,6 +6,7 @@ #include "Tree.h" #include "Views.h" #include "Context.h" +#include "Malloc.h" extern Context_T CTX; @@ -21,7 +22,6 @@ extern "C" } } - GMSH_SkinPlugin::GMSH_SkinPlugin() { } @@ -36,8 +36,9 @@ void GMSH_SkinPlugin::getInfos(char *author, char *copyright, char *help_text) c strcpy(author, "C. Geuzaine (geuz@geuz.org)"); strcpy(copyright, "DGR (www.multiphysics.com)"); strcpy(help_text, - "Gets the skin of a 3D view (eliminates all interior drawing).\n" - "Script name: Plugin(Skin)."); + "Gets the skin (i.e. the boundary) of a view,\n" + "eliminating all interior drawing).\n" + "Script name: Plugin(Skin).\n"); } int GMSH_SkinPlugin::getNbOptions() const @@ -55,63 +56,89 @@ void GMSH_SkinPlugin::CatchErrorMessage (char *errorMessage) const strcpy(errorMessage,"Skin failed..."); } -extern List_T *Post_ViewList; -struct elm{ - int nbnod; - double coord[9]; - double val[3]; -}; +static List_T * List; +static int * NbList, NbNod, NbComp, NbTime; + +typedef struct{ + double Coord[9]; + double *Val; +} Elm; -int fcmp_elm(const void *a, const void *b){ - struct elm *e1, *e2 ; - double s1, s2, TOL=CTX.lc*1.e-6 ; +static int fcmpElm(const void *a, const void *b){ + Elm *e1=(Elm*)a, *e2=(Elm*)b; + double s1, s2, TOL=CTX.lc*1.e-6; int i; - e1 = (struct elm*)a; e2 = (struct elm*)b; - s1 = s2 = 0.0 ; - for(i=0;i<e1->nbnod;i++){ s1 += e1->coord[i]; s2 += e2->coord[i]; } + for(i=0;i<NbNod-1;i++){ s1 += e1->Coord[i]; s2 += e2->Coord[i]; } if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1; - s1 = s2 = 0.0 ; - for(i=0;i<e1->nbnod;i++){ s1 += e1->coord[e1->nbnod+i]; s2 += e2->coord[e1->nbnod+i]; } + for(i=0;i<NbNod-1;i++){ s1 += e1->Coord[NbNod-1+i]; s2 += e2->Coord[NbNod-1+i]; } if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1; - s1 = s2 = 0.0 ; - for(i=0;i<e1->nbnod;i++){ s1 += e1->coord[2*e1->nbnod+i]; s2 += e2->coord[2*e1->nbnod+i]; } + for(i=0;i<NbNod-1;i++){ s1 += e1->Coord[2*(NbNod-1)+i]; s2 += e2->Coord[2*(NbNod-1)+i]; } if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1; return 0; } -void get_face(int *nod, int nbnod, int nbcomp, - double *C, double *V, double Cp[9], double Vp[3]){ - int i, j; - for(i=0; i<nbnod; i++) Cp[i] = C[nod[i]]; //x - for(i=0; i<nbnod; i++) Cp[nbnod+i] = C[(nbnod+1) + nod[i]]; //y - for(i=0; i<nbnod; i++) Cp[2*nbnod+i] = C[2*(nbnod+1) + nod[i]]; //z - for(i=0; i<nbnod; i++) - for(j=0; j<nbcomp; j++) Vp[nbcomp*i+j] = V[nbcomp*nod[i]+j]; //vals +static void getElm(int *Nod, double *Coord, double *Val, Elm *Elm){ + int i, j, k; + Elm->Val = (double*)Malloc((NbNod-1)*NbComp*NbTime*sizeof(double)); + for(i=0; i<NbNod-1; i++) Elm->Coord[i] = Coord[Nod[i]]; //x + for(i=0; i<NbNod-1; i++) Elm->Coord[NbNod-1+i] = Coord[NbNod + Nod[i]]; //y + for(i=0; i<NbNod-1; i++) Elm->Coord[2*(NbNod-1)+i] = Coord[2*NbNod + Nod[i]]; //z + for(i=0; i<NbTime; i++) + for(j=0; j<NbNod-1; j++) + for(k=0; k<NbComp; k++) + Elm->Val[(NbNod-1)*NbComp*i+NbComp*j+k] = + Val[NbNod*NbComp*i+NbComp*Nod[j]+k]; } -static Post_View *View; +static void addInView(void *a, void *b){ + int i, k; + Elm *e = (Elm*)a; + for(i=0; i<3*(NbNod-1); i++) List_Add(List, &e->Coord[i]); + for(i=0; i<NbTime; i++) + for(k=0;k<(NbNod-1)*NbComp;k++) + List_Add(List, &e->Val[(NbNod-1)*NbComp*i+k]); + Free(e->Val); + (*NbList)++; +} -void addSTinView(void *a, void *b){ - int k; - struct elm *elm = (struct elm*)a; - for(k=0;k<9;k++)List_Add(View->ST, &elm->coord[k]); - for(k=0;k<3;k++)List_Add(View->ST, &elm->val[k]); - View->NbST++; +static void skinSimplex(List_T *Simp, int NbSimp){ + double *Coords, *Vals; + int i, j; + int FacesTet[4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}}; + int EdgesTri[3][2] = {{0,1},{1,2},{2,0}}; + Elm e, *pe; + + Tree_T * Skin = Tree_Create(sizeof(Elm), fcmpElm); + for(i = 0 ; i < List_Nbr(Simp) ; i+=NbSimp){ + Coords = (double*)List_Pointer_Fast(Simp,i); + Vals = (double*)List_Pointer_Fast(Simp,i+3*NbNod); + for(j=0 ; j<NbNod ; j++){ + getElm(NbNod == 4 ? FacesTet[j] : EdgesTri[j],Coords,Vals,&e); + if(!(pe=(Elm*)Tree_PQuery(Skin, &e))) + Tree_Add(Skin, &e); + else{ + Free(pe->Val); + Free(e.Val); + Tree_Suppress(Skin, &e); + } + } + } + Tree_Action(Skin, addInView); + Tree_Delete(Skin); } +extern List_T * Post_ViewList; + Post_View *GMSH_SkinPlugin::execute (Post_View *v) { - Post_View *vv; - double *C, *V; - int faces_tet[4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}}; - int i, j, nb; - struct elm elm; + Post_View *vv, *View; + int iView = (int)SkinOptions_Number[0].def; if(v && iView < 0) @@ -124,34 +151,35 @@ Post_View *GMSH_SkinPlugin::execute (Post_View *v) } } - if(vv->NbSS){ + if(vv->NbSS || vv->NbVS || vv->NbST || vv->NbVT){ View = BeginView(1); - Tree_T * skin = Tree_Create(sizeof(struct elm), fcmp_elm); - nb = List_Nbr(vv->SS) / vv->NbSS ; - - for(i = 0 ; i < List_Nbr(vv->SS) ; i+=nb){ - C = (double*)List_Pointer_Fast(vv->SS,i); - V = (double*)List_Pointer_Fast(vv->SS,i+12); - - for(j=0 ; j<4 ; j++){//for each face - get_face(faces_tet[j],3,1,C,V,elm.coord,elm.val); - elm.nbnod = 3; - if(!Tree_PQuery(skin, &elm)) - Tree_Add(skin, &elm); - else - Tree_Suppress(skin, &elm); - } + NbTime = vv->NbTimeStep; + if(vv->NbSS){ + List = View->ST; NbList = &View->NbST; NbNod = 4; NbComp = 1; + skinSimplex(vv->SS, List_Nbr(vv->SS) / vv->NbSS); } - - Tree_Action(skin, addSTinView); - Tree_Delete(skin); - - char name[1024], filename[1024]; - sprintf(name,"skin-%s",vv->Name); - sprintf(filename,"skin-%s",vv->FileName); - EndView(1, filename, name); - Msg(INFO, "Created view '%s' (%d triangles)", name, View->NbST); - return View; + if(vv->NbVS){ + List = View->VT; NbList = &View->NbVT; NbNod = 4; NbComp = 3; + skinSimplex(vv->VS, List_Nbr(vv->VS) / vv->NbVS) ; + } + if(vv->NbST){ + List = View->SL; NbList = &View->NbSL; NbNod = 3; NbComp = 1; + skinSimplex(vv->ST, List_Nbr(vv->ST) / vv->NbST); + } + if(vv->NbVT){ + List = View->VL; NbList = &View->NbVL; NbNod = 3; NbComp = 3; + skinSimplex(vv->VT, List_Nbr(vv->VT) / vv->NbVT) ; + } + if(View->NbST || View->NbVT || View->NbSL || View->NbVL){ + char name[1024], filename[1024]; + sprintf(name,"skin-%s",vv->Name); + sprintf(filename,"skin-%s",vv->FileName); + EndView(View, 1, filename, name); + Msg(INFO, "Created view '%s'", name); + return View; + } + else + FreeView(View->Index); } return 0;