diff --git a/Plugin/DecomposeInSimplex.cpp b/Plugin/DecomposeInSimplex.cpp index bab2995e08bd3826094a580da08628cca0e6e43b..f92c30d61b2a766a408866a050f0a29a0e6de0ba 100644 --- a/Plugin/DecomposeInSimplex.cpp +++ b/Plugin/DecomposeInSimplex.cpp @@ -1,4 +1,4 @@ -// $Id: DecomposeInSimplex.cpp,v 1.3 2003-11-21 07:56:32 geuzaine Exp $ +// $Id: DecomposeInSimplex.cpp,v 1.4 2003-11-22 01:45:48 geuzaine Exp $ // // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle // @@ -111,17 +111,19 @@ Post_View *GMSH_DecomposeInSimplexPlugin::execute(Post_View * v) } // Utility class -int DecomposeInSimplex::num() + +int DecomposeInSimplex::numSimplices() { - switch(_numNodes){ - case 4 : return 2; // quad -> 2 tris - case 5 : return 2; // pyramid -> 2 tets - case 6 : return 3; // prism -> 3 tets - case 8 : return 6; // hexa -> 6 tets + switch(_numNodes) { + case 4: return 2; // quad -> 2 tris + case 5: return 2; // pyramid -> 2 tets + case 6: return 3; // prism -> 3 tets + case 8: return 6; // hexa -> 6 tets } + return 0; } -int DecomposeInSimplex::numNodes() +int DecomposeInSimplex::numSimplexNodes() { if(_numNodes == 4) return 3; // quad -> tris @@ -129,31 +131,37 @@ int DecomposeInSimplex::numNodes() return 4; // all others -> tets } -void DecomposeInSimplex::decompose() +void DecomposeInSimplex::reorder(int map[4], int n, + double *x, double *y, double *z, double *val, + double *xn, double *yn, double *zn, double *valn) { -#if 0 - switch(_numNodes){ - case 4: // quad - 0 1 2 - 0 2 3 - break ; - - case 8: // hexa - 0 1 2 5 - 0 2 5 6 - 0 4 5 6 - 0 2 3 6 - 0 4 6 7 - 0 3 6 7 - - case 6: // prism - 0 1 2 4 - 0 2 4 5 - 0 3 4 5 - - case 5: // pyramid - 0 1 3 4 - 1 2 3 4 + for(int i = 0; i < n; i++) { + xn[i] = x[map[i]]; + yn[i] = y[map[i]]; + zn[i] = z[map[i]]; + for(int j = 0; j < _numComponents; j++) + valn[i*_numComponents+j] = val[map[i]*_numComponents+j]; + } +} + +void DecomposeInSimplex::decompose(int num, + double *x, double *y, double *z, double *val, + double *xn, double *yn, double *zn, double *valn) +{ + int quadTri[2][4] = {{0,1,2,-1}, {0,2,3,-1}}; + int hexaTet[6][4] = {{0,1,2,5}, {0,2,5,6}, {0,4,5,6}, {0,2,3,6}, {0,4,6,7}, {0,3,6,7}}; + int prisTet[3][4] = {{0,1,2,4}, {0,2,4,5}, {0,3,4,5}}; + int pyraTet[2][4] = {{0,1,3,4}, {1,2,3,4}}; + + if(num < 0 || num > numSimplices()-1) { + Msg(GERROR, "Invalid decomposition"); + num = 0; + } + + switch(_numNodes) { + case 4: reorder(quadTri[num], 3, x, y, z, val, xn, yn, zn, valn); break ; + case 8: reorder(hexaTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ; + case 6: reorder(prisTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ; + case 5: reorder(pyraTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ; } -#endif } diff --git a/Plugin/DecomposeInSimplex.h b/Plugin/DecomposeInSimplex.h index 70dcac93ec0ec263a7aba22091bd8a0184a60fd5..cb469f2a56d0e4caa0e5f9b34685ed44bd5b1afa 100644 --- a/Plugin/DecomposeInSimplex.h +++ b/Plugin/DecomposeInSimplex.h @@ -45,19 +45,22 @@ class DecomposeInSimplex{ int _numNodes; // how many field components int _numComponents; - // how many time steps - int _numTimeSteps; + // create a simplex + void reorder(int map[4], int n, + double *x, double *y, double *z, double *val, + double *xn, double *yn, double *zn, double *valn); public: // default constructor - DecomposeInSimplex(int numNodes, int numComponents, int numTimeSteps) - : _numNodes(numNodes), _numComponents(numComponents), - _numTimeSteps(numTimeSteps) { ; } + DecomposeInSimplex(int numNodes, int numComponents) + : _numNodes(numNodes), _numComponents(numComponents) { ; } // the number of simplices into which the element is decomposed - int num(); + int numSimplices(); // the number of nodes of the simplex - int numNodes(); + int numSimplexNodes(); // returns the i-th simplex in the decomposition - void decompose(); + void decompose(int num, + double *x, double *y, double *z, double *val, + double *xn, double *yn, double *zn, double *valn); }; #endif diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp index 0f08e0545d9f0dd6baf267bc44a1a35b71190439..ef463c4f2b62ee54bef24bf025d4bb5cef47e374 100644 --- a/Plugin/Levelset.cpp +++ b/Plugin/Levelset.cpp @@ -1,4 +1,4 @@ -// $Id: Levelset.cpp,v 1.1 2003-11-21 07:56:32 geuzaine Exp $ +// $Id: Levelset.cpp,v 1.2 2003-11-22 01:45:48 geuzaine Exp $ // // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle // @@ -40,6 +40,17 @@ GMSH_LevelsetPlugin::GMSH_LevelsetPlugin() _orientation = GMSH_LevelsetPlugin::NONE; } +static void affect(double *xpi, double *ypi, double *zpi, double valpi[12][9], int i, + double *xp, double *yp, double *zp, double valp[12][9], int j, + int nb) +{ + xpi[i] = xp[j]; + ypi[i] = yp[j]; + zpi[i] = zp[j]; + for(int k = 0; k < nb; k++) + valpi[i][k] = valp[j][k]; +} + int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, int nbVert, int nbEdg, int exn[12][2], double *x, double *y, double *z, @@ -50,10 +61,9 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, double levels[8], scalarVal[8]; // compute the value of the levelset function at each node - if(_valueIndependent){ - for(int k = 0; k < nbVert; k++) { + if(_valueIndependent) { + for(int k = 0; k < nbVert; k++) levels[k] = levelset(x[k], y[k], z[k], 0.0); - } } else{ for(int k = 0; k < nbVert; k++) { @@ -78,14 +88,13 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, int np = 0; for(int k = 0; k < nbEdg; k++) { if(levels[exn[k][0]] * levels[exn[k][1]] <= 0.0) { - if(iVal && dVal){ + if(iVal && dVal) { double coef = InterpolateIso(x, y, z, levels, 0.0, exn[k][0], exn[k][1], &xp[np], &yp[np], &zp[np]); double *val1 = &dVal[dNbComp * exn[k][0]]; double *val2 = &dVal[dNbComp * exn[k][1]]; - for(int l = 0; l < dNbComp; l++) { + for(int l = 0; l < dNbComp; l++) valp[np][l] = coef * (val2[l] - val1[l]) + val1[l]; - } } np++; } @@ -94,28 +103,48 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, if(!iVal || !dVal) return np; + double xpi[12], ypi[12], zpi[12], valpi[12][9]; + + // Remove identical nodes (can happen only if an edge actually + // belongs to the zero levelset, i.e., if levels[] * levels[] == + // 0). Actuallym, we should be doing this even for np < 4, but it + // takes some time... And we don't really care if some nodes in + // a postprocessing element are identical. + if(np > 4){ + int npi; + affect(xpi, ypi, zpi, valpi, 0, xp, yp, zp, valp, 0, dNbComp); + npi = 1; + for(int j = 1; j < np; j++){ + for(int i = 0; i < npi; i++){ + if(fabs(xp[j] - xpi[i]) < 1.e-12 && + fabs(yp[j] - ypi[i]) < 1.e-12 && + fabs(zp[j] - zpi[i]) < 1.e-12){ + break; + } + if(i == npi-1){ + affect(xpi, ypi, zpi, valpi, i+1, xp, yp, zp, valp, j, dNbComp); + npi++; + } + } + } + for(int i = 0; i < npi; i++) + affect(xp, yp, zp, valp, i, xpi, ypi, zpi, valpi, i, dNbComp); + np = npi; + } + + // can't deal with this--just return... if(np < 1 || np > 4) return 0; // avoid ``butterflies'' - if(np == 4){ - double xx = xp[3], yy = yp[3], zz = zp[3], vv[9]; - for(int l = 0; l < dNbComp; l++) - vv[l] = valp[3][l]; - xp[3] = xp[2]; - yp[3] = yp[2]; - zp[3] = zp[2]; - for(int l = 0; l < dNbComp; l++) - valp[3][l] = valp[2][l]; - xp[2] = xx; - yp[2] = yy; - zp[2] = zz; - for(int l = 0; l < dNbComp; l++) - valp[2][l] = vv[l]; + if(np == 4) { + affect(xpi, ypi, zpi, valpi, 0, xp, yp, zp, valp, 3, dNbComp); + affect(xp, yp, zp, valp, 3, xp, yp, zp, valp, 2, dNbComp); + affect(xp, yp, zp, valp, 2, xpi, ypi, zpi, valpi, 0, dNbComp); } // orient the triangles and the quads to get the normals right - if(np == 3 || np == 4){ + if(np == 3 || np == 4) { if(!timeStep || !_valueIndependent) { // test this only once for spatially-fixed views double v1[3] = { xp[2] - xp[0], yp[2] - yp[0], zp[2] - zp[0] }; @@ -142,20 +171,10 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, } if(_invert > 0.) { double xpi[12], ypi[12], zpi[12], valpi[12][9]; - for(int k = 0; k < np; k++) { - xpi[k] = xp[k]; - ypi[k] = yp[k]; - zpi[k] = zp[k]; - for(int l = 0; l < dNbComp; l++) - valpi[k][l] = valp[k][l]; - } - for(int k = 0; k < np; k++) { - xp[k] = xpi[np - k - 1]; - yp[k] = ypi[np - k - 1]; - zp[k] = zpi[np - k - 1]; - for(int l = 0; l < dNbComp; l++) - valp[k][l] = valpi[np - k - 1][l]; - } + for(int k = 0; k < np; k++) + affect(xpi, ypi, zpi, valpi, k, xp, yp, zp, valp, k, dNbComp); + for(int k = 0; k < np; k++) + affect(xp, yp, zp, valp, k, xpi, ypi, zpi, valpi, np-k-1, dNbComp); } } @@ -163,32 +182,35 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, Post_View *view = _valueIndependent ? out[0] : out[timeStep]; List_T *list; int *nbPtr; - if(np == 1){ - if(dNbComp == 1) { list = view->SP; nbPtr = &view->NbSP; } - else if(dNbComp == 3){ list = view->VP; nbPtr = &view->NbVP; } - else { list = view->TP; nbPtr = &view->NbTP; } + if(np == 1) { + if(dNbComp == 1) { list = view->SP; nbPtr = &view->NbSP; } + else if(dNbComp == 3) { list = view->VP; nbPtr = &view->NbVP; } + else { list = view->TP; nbPtr = &view->NbTP; } } - else if(np == 2){ - if(dNbComp == 1) { list = view->SL; nbPtr = &view->NbSL; } - else if(dNbComp == 3){ list = view->VL; nbPtr = &view->NbVL; } - else { list = view->TL; nbPtr = &view->NbTL; } + else if(np == 2) { + if(dNbComp == 1) { list = view->SL; nbPtr = &view->NbSL; } + else if(dNbComp == 3) { list = view->VL; nbPtr = &view->NbVL; } + else { list = view->TL; nbPtr = &view->NbTL; } } - else if(np == 3){ - if(dNbComp == 1) { list = view->ST; nbPtr = &view->NbST; } - else if(dNbComp == 3){ list = view->VT; nbPtr = &view->NbVT; } - else { list = view->TT; nbPtr = &view->NbTT; } + else if(np == 3) { + if(dNbComp == 1) { list = view->ST; nbPtr = &view->NbST; } + else if(dNbComp == 3) { list = view->VT; nbPtr = &view->NbVT; } + else { list = view->TT; nbPtr = &view->NbTT; } } else{ - if(dNbComp == 1) { list = view->SQ; nbPtr = &view->NbSQ; } - else if(dNbComp == 3){ list = view->VQ; nbPtr = &view->NbVQ; } - else { list = view->TQ; nbPtr = &view->NbTQ; } + if(dNbComp == 1) { list = view->SQ; nbPtr = &view->NbSQ; } + else if(dNbComp == 3) { list = view->VQ; nbPtr = &view->NbVQ; } + else { list = view->TQ; nbPtr = &view->NbTQ; } } // copy the elements in the output view if(!timeStep || !_valueIndependent) { - for(int k = 0; k < np; k++) List_Add(list, &xp[k]); - for(int k = 0; k < np; k++) List_Add(list, &yp[k]); - for(int k = 0; k < np; k++) List_Add(list, &zp[k]); + for(int k = 0; k < np; k++) + List_Add(list, &xp[k]); + for(int k = 0; k < np; k++) + List_Add(list, &yp[k]); + for(int k = 0; k < np; k++) + List_Add(list, &zp[k]); (*nbPtr)++; } for(int k = 0; k < np; k++) @@ -200,10 +222,10 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, static int getTimeStep(int timeStep, int ts, Post_View *dView) { - if(timeStep < 0){ + if(timeStep < 0) { return ts; } - if(timeStep >= dView->NbTimeStep){ + if(timeStep >= dView->NbTimeStep) { Msg(WARNING, "Wrong TimeStep %d in View[%d]: reverting to TimeStep 0", timeStep, dView->Index); return 0; @@ -235,7 +257,7 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList, double *y = (double *)List_Pointer_Fast(iList, i + nbVert); double *z = (double *)List_Pointer_Fast(iList, i + 2 * nbVert); - if(nbVert == 2 || nbVert == 3 || (nbVert == 4 && nbEdg == 6)){ + if(nbVert == 2 || nbVert == 3 || (nbVert == 4 && nbEdg == 6)) { // easy for simplices: at most one element is created per time step for(int iTS = 0; iTS < iView->NbTimeStep; iTS++) { int dTS = getTimeStep(_timeStep, iTS, dView); @@ -248,38 +270,36 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList, } } else{ - Msg(GERROR, "Arghh! Non-tet levelset not finished!"); -#if 0 // we decompose the element into simplices - DecomposeInSimplex iDec(nbVert, iNbComp, iView->NbTimeStep); - DecomposeInSimplex dDec(nbVert, dNbComp, dView->NbTimeStep); + DecomposeInSimplex iDec(nbVert, iNbComp); + DecomposeInSimplex dDec(nbVert, dNbComp); - int nbVert_new = iDec.numNodes(); - int nbEdg_new = (nbVert_new==4)?6:3; + int nbVertNew = iDec.numSimplexNodes(); + int nbEdgNew = (nbVertNew == 4) ? 6 : 3; int exnTri[12][2] = {{0,1},{0,2},{1,2}}; int exnTet[12][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; - double x_new[4], y_new[4], z_new[4]; - double *iVal_new = new double[4 * iNbComp * iView->NbTimeStep]; - double *dVal_new = new double[4 * dNbComp * dView->NbTimeStep]; + double xNew[4], yNew[4], zNew[4]; + double *iValNew = new double[4 * iNbComp]; + double *dValNew = new double[4 * dNbComp]; - if(_valueIndependent){ + if(_valueIndependent) { // since the intersection is value-independent, we can loop on // the time steps so that only one element gets created per // time step. This allows us to still generate a *single* // multi-timestep view. - if(zeroLevelset(0, nbVert, nbEdg, exn, x, y, z, NULL, 0, NULL, 0, out)){ - for(int k = 0; k < iDec.num(); k++){ + if(zeroLevelset(0, nbVert, nbEdg, exn, x, y, z, NULL, 0, + NULL, 0, out)) { + for(int k = 0; k < iDec.numSimplices(); k++) { for(int iTS = 0; iTS < iView->NbTimeStep; iTS++) { int dTS = getTimeStep(_timeStep, iTS, dView); double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbVert + iNbComp * nbVert * iTS); double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbVert + dNbComp * nbVert * dTS); - iDec.decompose(k, x, y, z, iVal, x_new, y_new, z_new, iVal_new); - dDec.decompose(k, x, y, z, dVal, x_new, y_new, z_new, dVal_new); - zeroLevelset(iTS, nbVert_new, nbEdg_new, (nbVert_new==4)?exnTet:exnTri, - x_new, y_new, z_new, iVal_new, iNbComp, dVal_new, dNbComp, - out); + iDec.decompose(k, x, y, z, iVal, xNew, yNew, zNew, iValNew); + dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew); + zeroLevelset(iTS, nbVertNew, nbEdgNew, (nbVertNew == 4) ? exnTet : exnTri, + xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out); } } } @@ -288,27 +308,25 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList, // since we generate one view for each time step, we can // generate multiple elements per time step without problem. for(int iTS = 0; iTS < iView->NbTimeStep; iTS++) { - int dTS = getTimeStep(_timeStep, iTS, dView); double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbVert + iNbComp * nbVert * iTS); - double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbVert + - dNbComp * nbVert * dTS); - if(zeroLevelset(1, iTS, nbVert, nbEdg, exn, x, y, z, - iVal, iNbComp, dVal, dNbComp, out)){ - for(int k = 0; k < iDec.num(); k++){ - iDec.decompose(k, x, y, z, iVal, x_new, y_new, z_new, iVal_new); - dDec.decompose(k, x, y, z, dVal, x_new, y_new, z_new, dVal_new); - zeroLevelset(iTS, nbVert_new, nbEdg_new, (nbVert_new==4)?exnTet:exnTri, - x_new, y_new, z_new, iVal_new, iNbComp, dVal_new, dNbComp, - out); + if(zeroLevelset(iTS, nbVert, nbEdg, exn, x, y, z, iVal, iNbComp, + NULL, 0, out)) { + int dTS = getTimeStep(_timeStep, iTS, dView); + double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbVert + + dNbComp * nbVert * dTS); + for(int k = 0; k < iDec.numSimplices(); k++) { + iDec.decompose(k, x, y, z, iVal, xNew, yNew, zNew, iValNew); + dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew); + zeroLevelset(iTS, nbVertNew, nbEdgNew, (nbVertNew == 4) ? exnTet : exnTri, + xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out); } } } } - delete [] iVal_new; - delete [] dVal_new; -#endif + delete [] iValNew; + delete [] dValNew; } } @@ -319,7 +337,7 @@ Post_View *GMSH_LevelsetPlugin::execute(Post_View * v) Post_View *w; vector<Post_View *> out; - if(_targetView < 0){ + if(_targetView < 0) { w = v; } else if(!(w = (Post_View *)List_Pointer_Test(CTX.post.list, _targetView))) { @@ -336,60 +354,109 @@ Post_View *GMSH_LevelsetPlugin::execute(Post_View * v) out.push_back(BeginView(1)); } + // We should definitely recode the View interface in C++ (and define + // iterators on the different kinds of elements...) + // lines int exnLin[12][2] = {{0,1}}; executeList(v, v->SL, v->NbSL, 1, w, w->SL, w->NbSL, 1, 2, 1, exnLin, out); - executeList(v, v->SL, v->NbSL, 1, w, w->SL, w->NbVL, 3, 2, 1, exnLin, out); - executeList(v, v->SL, v->NbSL, 1, w, w->SL, w->NbTL, 9, 2, 1, exnLin, out); + executeList(v, v->SL, v->NbSL, 1, w, w->VL, w->NbVL, 3, 2, 1, exnLin, out); + executeList(v, v->SL, v->NbSL, 1, w, w->TL, w->NbTL, 9, 2, 1, exnLin, out); executeList(v, v->VL, v->NbVL, 3, w, w->SL, w->NbSL, 1, 2, 1, exnLin, out); - executeList(v, v->VL, v->NbVL, 3, w, w->SL, w->NbVL, 3, 2, 1, exnLin, out); - executeList(v, v->VL, v->NbVL, 3, w, w->SL, w->NbTL, 9, 2, 1, exnLin, out); + executeList(v, v->VL, v->NbVL, 3, w, w->VL, w->NbVL, 3, 2, 1, exnLin, out); + executeList(v, v->VL, v->NbVL, 3, w, w->TL, w->NbTL, 9, 2, 1, exnLin, out); executeList(v, v->TL, v->NbTL, 9, w, w->SL, w->NbSL, 1, 2, 1, exnLin, out); - executeList(v, v->TL, v->NbTL, 9, w, w->SL, w->NbVL, 3, 2, 1, exnLin, out); - executeList(v, v->TL, v->NbTL, 9, w, w->SL, w->NbTL, 9, 2, 1, exnLin, out); + executeList(v, v->TL, v->NbTL, 9, w, w->VL, w->NbVL, 3, 2, 1, exnLin, out); + executeList(v, v->TL, v->NbTL, 9, w, w->TL, w->NbTL, 9, 2, 1, exnLin, out); // triangles - int exnTri[12][2] = {{0,1},{0,2},{1,2}}; + int exnTri[12][2] = {{0,1}, {0,2}, {1,2}}; executeList(v, v->ST, v->NbST, 1, w, w->ST, w->NbST, 1, 3, 3, exnTri, out); - executeList(v, v->ST, v->NbST, 1, w, w->ST, w->NbVT, 3, 3, 3, exnTri, out); - executeList(v, v->ST, v->NbST, 1, w, w->ST, w->NbTT, 9, 3, 3, exnTri, out); + executeList(v, v->ST, v->NbST, 1, w, w->VT, w->NbVT, 3, 3, 3, exnTri, out); + executeList(v, v->ST, v->NbST, 1, w, w->TT, w->NbTT, 9, 3, 3, exnTri, out); executeList(v, v->VT, v->NbVT, 3, w, w->ST, w->NbST, 1, 3, 3, exnTri, out); - executeList(v, v->VT, v->NbVT, 3, w, w->ST, w->NbVT, 3, 3, 3, exnTri, out); - executeList(v, v->VT, v->NbVT, 3, w, w->ST, w->NbTT, 9, 3, 3, exnTri, out); + executeList(v, v->VT, v->NbVT, 3, w, w->VT, w->NbVT, 3, 3, 3, exnTri, out); + executeList(v, v->VT, v->NbVT, 3, w, w->TT, w->NbTT, 9, 3, 3, exnTri, out); executeList(v, v->TT, v->NbTT, 9, w, w->ST, w->NbST, 1, 3, 3, exnTri, out); - executeList(v, v->TT, v->NbTT, 9, w, w->ST, w->NbVT, 3, 3, 3, exnTri, out); - executeList(v, v->TT, v->NbTT, 9, w, w->ST, w->NbTT, 9, 3, 3, exnTri, out); + executeList(v, v->TT, v->NbTT, 9, w, w->VT, w->NbVT, 3, 3, 3, exnTri, out); + executeList(v, v->TT, v->NbTT, 9, w, w->TT, w->NbTT, 9, 3, 3, exnTri, out); // tets int exnTet[12][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; executeList(v, v->SS, v->NbSS, 1, w, w->SS, w->NbSS, 1, 4, 6, exnTet, out); - executeList(v, v->SS, v->NbSS, 1, w, w->SS, w->NbVS, 3, 4, 6, exnTet, out); - executeList(v, v->SS, v->NbSS, 1, w, w->SS, w->NbTS, 9, 4, 6, exnTet, out); + executeList(v, v->SS, v->NbSS, 1, w, w->VS, w->NbVS, 3, 4, 6, exnTet, out); + executeList(v, v->SS, v->NbSS, 1, w, w->TS, w->NbTS, 9, 4, 6, exnTet, out); executeList(v, v->VS, v->NbVS, 3, w, w->SS, w->NbSS, 1, 4, 6, exnTet, out); - executeList(v, v->VS, v->NbVS, 3, w, w->SS, w->NbVS, 3, 4, 6, exnTet, out); - executeList(v, v->VS, v->NbVS, 3, w, w->SS, w->NbTS, 9, 4, 6, exnTet, out); + executeList(v, v->VS, v->NbVS, 3, w, w->VS, w->NbVS, 3, 4, 6, exnTet, out); + executeList(v, v->VS, v->NbVS, 3, w, w->TS, w->NbTS, 9, 4, 6, exnTet, out); executeList(v, v->TS, v->NbTS, 9, w, w->SS, w->NbSS, 1, 4, 6, exnTet, out); - executeList(v, v->TS, v->NbTS, 9, w, w->SS, w->NbVS, 3, 4, 6, exnTet, out); - executeList(v, v->TS, v->NbTS, 9, w, w->SS, w->NbTS, 9, 4, 6, exnTet, out); + executeList(v, v->TS, v->NbTS, 9, w, w->VS, w->NbVS, 3, 4, 6, exnTet, out); + executeList(v, v->TS, v->NbTS, 9, w, w->TS, w->NbTS, 9, 4, 6, exnTet, out); // quads - //int exnQua[12][2] = + int exnQua[12][2] = {{0,1}, {0,3}, {1,2}, {2,3}}; + executeList(v, v->SQ, v->NbSQ, 1, w, w->SQ, w->NbSQ, 1, 4, 4, exnQua, out); + executeList(v, v->SQ, v->NbSQ, 1, w, w->VQ, w->NbVQ, 3, 4, 4, exnQua, out); + executeList(v, v->SQ, v->NbSQ, 1, w, w->TQ, w->NbTQ, 9, 4, 4, exnQua, out); + + executeList(v, v->VQ, v->NbVQ, 3, w, w->SQ, w->NbSQ, 1, 4, 4, exnQua, out); + executeList(v, v->VQ, v->NbVQ, 3, w, w->VQ, w->NbVQ, 3, 4, 4, exnQua, out); + executeList(v, v->VQ, v->NbVQ, 3, w, w->TQ, w->NbTQ, 9, 4, 4, exnQua, out); + + executeList(v, v->TQ, v->NbTQ, 9, w, w->SQ, w->NbSQ, 1, 4, 4, exnQua, out); + executeList(v, v->TQ, v->NbTQ, 9, w, w->VQ, w->NbVQ, 3, 4, 4, exnQua, out); + executeList(v, v->TQ, v->NbTQ, 9, w, w->TQ, w->NbTQ, 9, 4, 4, exnQua, out); // hexes - //int exnHex[12][2] = + int exnHex[12][2] = {{0,1}, {0,3}, {0,4}, {1,2}, {1,5}, {2,3}, + {2,6}, {3,7}, {4,5}, {4,7}, {5,6}, {6,7}}; + executeList(v, v->SH, v->NbSH, 1, w, w->SH, w->NbSH, 1, 8, 12, exnHex, out); + executeList(v, v->SH, v->NbSH, 1, w, w->VH, w->NbVH, 3, 8, 12, exnHex, out); + executeList(v, v->SH, v->NbSH, 1, w, w->TH, w->NbTH, 9, 8, 12, exnHex, out); + + executeList(v, v->VH, v->NbVH, 3, w, w->SH, w->NbSH, 1, 8, 12, exnHex, out); + executeList(v, v->VH, v->NbVH, 3, w, w->VH, w->NbVH, 3, 8, 12, exnHex, out); + executeList(v, v->VH, v->NbVH, 3, w, w->TH, w->NbTH, 9, 8, 12, exnHex, out); + + executeList(v, v->TH, v->NbTH, 9, w, w->SH, w->NbSH, 1, 8, 12, exnHex, out); + executeList(v, v->TH, v->NbTH, 9, w, w->VH, w->NbVH, 3, 8, 12, exnHex, out); + executeList(v, v->TH, v->NbTH, 9, w, w->TH, w->NbTH, 9, 8, 12, exnHex, out); // prisms - //int exnPri[12][2] = + int exnPri[12][2] = {{0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5}, + {3,4}, {3,5}, {4,5}}; + executeList(v, v->SI, v->NbSI, 1, w, w->SI, w->NbSI, 1, 6, 9, exnPri, out); + executeList(v, v->SI, v->NbSI, 1, w, w->VI, w->NbVI, 3, 6, 9, exnPri, out); + executeList(v, v->SI, v->NbSI, 1, w, w->TI, w->NbTI, 9, 6, 9, exnPri, out); + + executeList(v, v->VI, v->NbVI, 3, w, w->SI, w->NbSI, 1, 6, 9, exnPri, out); + executeList(v, v->VI, v->NbVI, 3, w, w->VI, w->NbVI, 3, 6, 9, exnPri, out); + executeList(v, v->VI, v->NbVI, 3, w, w->TI, w->NbTI, 9, 6, 9, exnPri, out); + + executeList(v, v->TI, v->NbTI, 9, w, w->SI, w->NbSI, 1, 6, 9, exnPri, out); + executeList(v, v->TI, v->NbTI, 9, w, w->VI, w->NbVI, 3, 6, 9, exnPri, out); + executeList(v, v->TI, v->NbTI, 9, w, w->TI, w->NbTI, 9, 6, 9, exnPri, out); // pyramids - //int exnPyr[12][2] = - + int exnPyr[12][2] = {{0,1}, {0,3}, {0,4}, {1,2}, {1,4}, {2,3}, {2,4}, {3,4}}; + executeList(v, v->SY, v->NbSY, 1, w, w->SY, w->NbSY, 1, 5, 8, exnPyr, out); + executeList(v, v->SY, v->NbSY, 1, w, w->VY, w->NbVY, 3, 5, 8, exnPyr, out); + executeList(v, v->SY, v->NbSY, 1, w, w->TY, w->NbTY, 9, 5, 8, exnPyr, out); + + executeList(v, v->VY, v->NbVY, 3, w, w->SY, w->NbSY, 1, 5, 8, exnPyr, out); + executeList(v, v->VY, v->NbVY, 3, w, w->VY, w->NbVY, 3, 5, 8, exnPyr, out); + executeList(v, v->VY, v->NbVY, 3, w, w->TY, w->NbTY, 9, 5, 8, exnPyr, out); + + executeList(v, v->TY, v->NbTY, 9, w, w->SY, w->NbSY, 1, 5, 8, exnPyr, out); + executeList(v, v->TY, v->NbTY, 9, w, w->VY, w->NbVY, 3, 5, 8, exnPyr, out); + executeList(v, v->TY, v->NbTY, 9, w, w->TY, w->NbTY, 9, 5, 8, exnPyr, out); + for(unsigned int i = 0; i < out.size(); i++) { char name[1024], filename[1024]; sprintf(name, "cut-%s-%d", v->Name, i); @@ -399,7 +466,7 @@ Post_View *GMSH_LevelsetPlugin::execute(Post_View * v) // remove empty views (this is a bit ugly because, due to the // dynamic GUI events, this should actually be locked... - for(int i = List_Nbr(CTX.post.list) - 1; i >= 0; --i){ + for(int i = List_Nbr(CTX.post.list) - 1; i >= 0; --i) { w = (Post_View*) List_Pointer_Test(CTX.post.list, i); if(w && w->empty()) RemoveViewByIndex(i); diff --git a/benchmarks/misc/levelset/levelsetTest.geo b/benchmarks/misc/levelset/levelsetTest.geo new file mode 100644 index 0000000000000000000000000000000000000000..3ee8155d7fad98b11ee54a1ff99dc4ed59bcb04b --- /dev/null +++ b/benchmarks/misc/levelset/levelsetTest.geo @@ -0,0 +1,39 @@ +lc = 0.1; + +Point(1) = {0,0,0,lc}; +Point(2) = {0,1,0,lc}; +Point(3) = {1,1,0,lc}; +Point(4) = {1,0,0,lc}; +Line(1) = {1,4}; +Line(2) = {4,3}; +Line(3) = {3,2}; +Line(4) = {2,1}; +Line Loop(5) = {3,4,1,2}; + +Plane Surface(6) = {5}; + +Transfinite Surface{6} = {1,2,3,4}; +Recombine Surface{6}; + +Translate {2,0,0} { + Duplicata{ Surface{6}; } +} +Translate {4,0,0} { + Duplicata{ Surface{6}; } +} + +Extrude Surface {6, {-1,1,0}, {2.5,2.5,1}, Pi/4}{ + Layers { {2/lc}, {100}, {1} }; Recombine; +}; + +Extrude Surface {7, {0,0,1}}{ + Layers { {1/lc}, {200}, {1} }; Recombine; +}; + +Extrude Surface {12, {-1,-1,0}, {2.5,2.5,1}, Pi/4}{ + Layers { {1/lc}, {300}, {1} }; +}; + +Physical Volume(1) = {100,200,300}; +Physical Surface(2) = {1:100}; +Physical Line(3) = {1:100}; diff --git a/benchmarks/misc/levelset/levelsetTest.opt b/benchmarks/misc/levelset/levelsetTest.opt new file mode 100644 index 0000000000000000000000000000000000000000..e9062d6db0af1d2f1504153f5338792a2d3e18ed --- /dev/null +++ b/benchmarks/misc/levelset/levelsetTest.opt @@ -0,0 +1,18 @@ +Plugin(CutPlane).A = 1; +Plugin(CutPlane).B = 0; +Plugin(CutPlane).C = 0; +Plugin(CutPlane).D = 0; +Plugin(CutPlane).Run; + +Plugin(CutPlane).D = -2.5; +Plugin(CutPlane).Run; + +Plugin(CutPlane).D = -4.1; +Plugin(CutPlane).Run; + +Plugin(CutMap).A = 1; +Plugin(CutMap).Run; + +Plugin(CutMap).A = 22; +Plugin(CutMap).Run; + diff --git a/benchmarks/misc/levelset/levelsetTest.pro b/benchmarks/misc/levelset/levelsetTest.pro new file mode 100644 index 0000000000000000000000000000000000000000..f0ef74f62a9f69b81f31537b5a36d34e02b75aa2 --- /dev/null +++ b/benchmarks/misc/levelset/levelsetTest.pro @@ -0,0 +1,79 @@ +// Moronic getdp input file to create some test post-proessing maps + +Group { + Omega = Region[ {1,2,3,4} ]; +} + +Function { + Fct[] = Complex[ X[]^2-Sin[2*Y[]]+Sqrt[Z[]] , Sin[2*Z[]] + X[]^2 ]; +} + +Jacobian { + { Name JVol ; Case { { Region All ; Jacobian Vol ; } } } + { Name JSur ; Case { { Region All ; Jacobian Sur ; } } } +} + +Integration { + { Name I1 ; + Case { {Type Gauss ; + Case { { GeoElement Line ; NumberOfPoints 4 ; } + { GeoElement Triangle ; NumberOfPoints 12 ; } + { GeoElement Quadrangle ; NumberOfPoints 4 ; } + { GeoElement Tetrahedron ; NumberOfPoints 15 ; } + { GeoElement Hexahedron ; NumberOfPoints 34 ; } + { GeoElement Prism ; NumberOfPoints 9 ; } + { GeoElement Pyramid ; NumberOfPoints 8 ; } } + } + } + } +} + +FunctionSpace { + { Name Proj ; Type Form0 ; + BasisFunction { + { Name sn ; NameOfCoef phin ; Function BF_Node ; + Support Omega ; Entity NodesOf[ All ] ; } + } + } +} + +Formulation { + { Name Proj ; Type FemEquation ; + Quantity { + { Name v ; Type Local ; NameOfSpace Proj ; } + } + Equation { + Galerkin { [ Dof{v} , {v} ] ; In Omega ; Jacobian JVol ; Integration I1 ; } + Galerkin { [ -Fct[] , {v} ] ; In Omega ; Jacobian JVol ; Integration I1 ; } + } + } + +} + +Resolution { + { Name Proj ; + System { + { Name Proj ; NameOfFormulation Proj ; Type Complex; Frequency 1; } + } + Operation { + Generate Proj ; Solve Proj ; SaveSolution Proj ; + } + } + +} + +PostProcessing { + { Name Proj ; NameOfFormulation Proj ; + Quantity { + { Name v ; Value { Term { [ {v} ] ; In Omega ; } } } + { Name v3 ; Value { Term { [ Vector[0,{v},0] ] ; In Omega ; } } } + } + } + +} + +PostOperation v UsingPost Proj { + Print[ v , OnElementsOf Omega , File "levelsetTest1.pos" ] ; + Print[ v , OnElementsOf Omega , File "levelsetTest2.pos" , DecomposeInSimplex ] ; + Print[ v3 , OnElementsOf Omega , File "levelsetTest3.pos" ] ; +} diff --git a/doc/VERSIONS b/doc/VERSIONS index 90aa0c764887b50b319b9b751aa1c3c4e210980a..81f8937fde3938421c7b50595b5198e170ff38bc 100644 --- a/doc/VERSIONS +++ b/doc/VERSIONS @@ -1,8 +1,10 @@ -$Id: VERSIONS,v 1.160 2003-11-21 08:00:04 geuzaine Exp $ +$Id: VERSIONS,v 1.161 2003-11-22 01:45:48 geuzaine Exp $ New since 1.47: new DisplacementRaise plugin to plot arbitrary fields -on a deformed mesh; genearlized CutMap, CutPlane and CutSphere -plugins; small bug fixes (configure script); +on deformed meshes; generalized CutMap, CutPlane and CutSphere +plugins; small bug fixes (configure tests for libpng, handling of +erroneous options, multi timestep scalar prism drawings, copy of +surface mesh attributes); New in 1.47: fixed extrusion of surfaces defined by only two curves; new syntax to retrieve point coordinates and indices of entities