Generalization of the levelset plugin, take II.

One can now compute cuts/maps of abritrary views (any dimension, any kind of
element, any kind of field) + one can plot any other view on the resulting
cut/map.

Plugin/DecomposeInSimplex.cpp

-// $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
 }

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

benchmarks/misc/levelset/levelsetTest.geo

+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};

benchmarks/misc/levelset/levelsetTest.opt

+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;
+

benchmarks/misc/levelset/levelsetTest.pro

+// 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"  ] ;
+}

doc/VERSIONS

-$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