From 6000d7f2fcfd8e576a1f092ae79a3b7fd94b0e36 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 18 Apr 2006 03:13:30 +0000
Subject: [PATCH] *** empty log message ***
---
benchmarks/3d/percolation.geo | 229 +++++++++++++++++-----------------
1 file changed, 113 insertions(+), 116 deletions(-)
diff --git a/benchmarks/3d/percolation.geo b/benchmarks/3d/percolation.geo
index 6f2053c173..836624cde8 100644
--- a/benchmarks/3d/percolation.geo
+++ b/benchmarks/3d/percolation.geo
@@ -1,138 +1,135 @@
-lc=0.25;
-//dimensions
-xdim=5;
-ydim=5;
-zdim=5;
-//electrodes
-a=newp; Point(a) = {-1, -1, -Cos(45),lc};
-a=newp; Point(a) = {xdim, -1, -Cos(45),lc};
-a=newp; Point(a) = {xdim, ydim, -Cos(45),lc};
-a=newp; Point(a) = {-1, ydim, -Cos(45),lc};
-b1=newl; Line (b1) = {a-3,a-2};
-b2=newl; Line (b2) = {a-2,a-1};
-b3=newl; Line (b3) = {a-1,a};
-b4=newl; Line (b4) = {a,a-3};
-c=newll; Line Loop(c) = {b1,b2,b3,b4};
-ss=news; Plane Surface(ss) = {c};
+sphere_dist = 0.05; // spheres don't touch
+// sphere_dist = 0; // spheres touch
-Extrude {0,0,-1} { Surface{ss}; } //lower electrode
-aa[] = Extrude {0,0,zdim-1+2*Cos(45)} { Surface{ss}; } ;
-Extrude {0,0,1} { Surface{aa[0]}; } // top electrode
+lc = 0.5;
-Delete{ Volume{aa[1]}; } // this volume is wrong: we will be carving holes in it
+xdim = 5;
+ydim = 5;
+zdim = 5;
-//functions
-Function newsphere
-rr = Cos(45)-0.05;
-a=newp; Point(a) = {i,j,k,lc};
-a=newp; Point(a) = {i,j+ rr,k,lc};
-a=newp; Point(a) = {i,j- rr,k,lc};
-a=newp; Point(a) = {i+rr,j,k,lc};
-a=newp; Point(a) = {i-rr,j,k,lc};
-a=newp; Point(a) = {i,j,k+rr,lc};
-a=newp; Point(a) = {i,j,k-rr,lc};
-b=newl; Circle(b) = {a-1,a-6,a-3};
-b=newl; Circle(b) = {a-3,a-6,a};
-b=newl; Circle(b) = {a,a-6,a-2};
-b=newl; Circle(b) = {a-2,a-6,a-1};
-b=newl; Circle(b) = {a-5,a-6,a};
-b=newl; Circle(b) = {a,a-6,a-4};
-b=newl; Circle(b) = {a-4,a-6,a-1};
-b=newl; Circle(b) = {a-1,a-6,a-5};
-b=newl; Circle(b) = {a-5,a-6,a-3};
-b=newl; Circle(b) = {a-3,a-6,a-4};
-b=newl; Circle(b) = {a-4,a-6,a-2};
-b=newl; Circle(b) = {a-2,a-6,a-5};
-c=newll; Line Loop(c) = {b-6,-(b-2),b-10};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-2,b-5,(b-11)};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-5,-(b-8),-(b-1)};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-1,-(b-9),b-6};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b,b-7,(b-9)};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-7,-(b-10),-(b-3)};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-3,-(b-11),b-4};
-d=news; Ruled Surface(d) = {c};
-c=newll; Line Loop(c) = {b-4,-b,b-8};
-d=news; Ruled Surface(d) = {c};
-e=newsl; Surface Loop(e) = {(d-14),-(d-8),-(d-10),(d-12),(d-2),(d-4),-(d-6),-d};
-vv=newv; Volume(vv)={e};
+Point(1) = {-1,-1,-Cos(45),lc};
+Point(2) = {xdim,-1,-Cos(45),lc};
+Point(3) = {xdim,ydim,-Cos(45),lc};
+Point(4) = {-1,ydim,-Cos(45),lc};
+Line (1) = {1,2};
+Line (2) = {2,3};
+Line (3) = {3,4};
+Line (4) = {4,1};
+Line Loop(1) = {1,2,3,4};
+Plane Surface(1) = {1};
-holes[]+=e;
+Extrude {0, 0, -1} { Surface{1}; } // lower electrode
+aa[] = Extrude {0, 0, zdim - 1 + 2 * Cos(45)} { Surface{1}; } ;
+Extrude {0, 0, 1} { Surface{aa[0]}; } // top electrode
+
+Delete{ Volume{aa[1]}; } // this volume is wrong: we will be carving holes in it
+Function newsphere
+ rr = 0.5-sphere_dist;
+ a = newp;
+ b = newl;
+ c = newll;
+ d = news;
+ e = newsl;
+ f = newv;
+ Point(a) = {i,j,k,lc};
+ Point(a+1) = {i,j+rr,k,lc};
+ Point(a+2) = {i,j-rr,k,lc};
+ Point(a+3) = {i+rr,j,k,lc};
+ Point(a+4) = {i-rr,j,k,lc};
+ Point(a+5) = {i,j,k+rr,lc};
+ Point(a+6) = {i,j,k-rr,lc};
+ Circle(b) = {a+5,a,a+3};
+ Circle(b+1) = {a+3,a,a+6};
+ Circle(b+2) = {a+6,a,a+4};
+ Circle(b+3) = {a+4,a,a+5};
+ Circle(b+4) = {a+1,a,a+6};
+ Circle(b+5) = {a+6,a,a+2};
+ Circle(b+6) = {a+2,a,a+5};
+ Circle(b+7) = {a+5,a,a+1};
+ Circle(b+8) = {a+1,a,a+3};
+ Circle(b+9) = {a+3,a,a+2};
+ Circle(b+10) = {a+2,a,a+4};
+ Circle(b+11) = {a+4,a,a+1};
+ Line Loop(c) = {b+5,-(b+9),b+1}; Ruled Surface(d) = {c};
+ Line Loop(c+1) = {b+9,b+6,b}; Ruled Surface(d+1) = {c+1};
+ Line Loop(c+2) = {b+6,-(b+3),-(b+10)}; Ruled Surface(d+2) = {c+2};
+ Line Loop(c+3) = {b+10,-(b+2),b+5}; Ruled Surface(d+3) = {c+3};
+ Line Loop(c+4) = {b+11,b+4,b+2}; Ruled Surface(d+4) = {c+4};
+ Line Loop(c+5) = {b+4,-(b+1),-(b+8)}; Ruled Surface(d+5) = {c+5};
+ Line Loop(c+6) = {b+8,-b,b+7}; Ruled Surface(d+6) = {c+6};
+ Line Loop(c+7) = {b+7,-(b+11),b+3}; Ruled Surface(d+7) = {c+7};
+ Surface Loop(e) = {d+1,-(d+4),-(d+3),d+2,d,(d+6),-(d+5),-(d+7)};
+ Volume(f) = e;
+ holes[] += e;
Return
//random distribution
For i In {0:xdim-1:1}
-For j In {0:ydim-1:1}
-For k In {0:zdim-1:1}
-num=Fmod(Rand(1), 1);
-If (num>=0.5) num=1; EndIf
-If (num<0.5) num=0; EndIf
-array[i*ydim*zdim + j*zdim + k] =num;
-EndFor
-EndFor
+ For j In {0:ydim-1:1}
+ For k In {0:zdim-1:1}
+ num = Fmod(Rand(1), 1);
+ If(num >= 0.5) num = 1; EndIf
+ If(num < 0.5) num = 0; EndIf
+ array[i*ydim*zdim + j*zdim + k] = num;
+ EndFor
+ EndFor
EndFor
//remove isolated particles
For i In {1:xdim-2:1}
-For j In {1:ydim-2:1}
-For k In {1:zdim-2:1}
-If (array[i*ydim*zdim + j*zdim + k]==1)
- If (
- array[(i+1)*ydim*zdim + j*zdim + k]==0 &&
- array[(i+1)*ydim*zdim + (j+1)*zdim + k]==0 &&
- array[(i+1)*ydim*zdim + (j-1)*zdim + k]==0 &&
- array[(i-1)*ydim*zdim + j*zdim + k]==0 &&
- array[(i-1)*ydim*zdim + (j+1)*zdim + k]==0 &&
- array[(i-1)*ydim*zdim + (j-1)*zdim + k]==0 &&
- array[(i)*ydim*zdim + (j+1)*zdim + k]==0 &&
- array[(i)*ydim*zdim + (j-1)*zdim + k]==0 &&
-
- array[(i+1)*ydim*zdim + j*zdim + (k-1)]==0 &&
- array[(i+1)*ydim*zdim + (j+1)*zdim + (k-1)]==0 &&
- array[(i+1)*ydim*zdim + (j-1)*zdim + (k-1)]==0 &&
- array[(i-1)*ydim*zdim + j*zdim + (k-1)]==0 &&
- array[(i-1)*ydim*zdim + (j+1)*zdim + (k-1)]==0 &&
- array[(i-1)*ydim*zdim + (j-1)*zdim + (k-1)]==0 &&
- array[(i)*ydim*zdim + (j+1)*zdim + (k-1)]==0 &&
- array[(i)*ydim*zdim + (j-1)*zdim + (k-1)]==0 &&
- array[i*ydim*zdim + j*zdim + (k-1)] ==0 &&
-
- array[(i+1)*ydim*zdim + j*zdim + (k+1)]==0 &&
- array[(i+1)*ydim*zdim + (j+1)*zdim + (k+1)]==0 &&
- array[(i+1)*ydim*zdim + (j-1)*zdim + (k+1)]==0 &&
- array[(i-1)*ydim*zdim + j*zdim + (k+1)]==0 &&
- array[(i-1)*ydim*zdim + (j+1)*zdim + (k+1)]==0 &&
- array[(i-1)*ydim*zdim + (j-1)*zdim + (k+1)]==0 &&
- array[(i)*ydim*zdim + (j+1)*zdim + (k+1)]==0 &&
- array[(i)*ydim*zdim + (j-1)*zdim + (k+1)]==0 &&
- array[i*ydim*zdim + j*zdim + (k+1)] ==0)
- array[i*ydim*zdim + j*zdim + k]=0;
-EndIf
-EndIf
-
-EndFor
-EndFor
+ For j In {1:ydim-2:1}
+ For k In {1:zdim-2:1}
+ If (array[i*ydim*zdim + j*zdim + k] == 1)
+ If (array[(i+1)*ydim*zdim + j*zdim + k] == 0 &&
+ array[(i+1)*ydim*zdim + (j+1)*zdim + k] == 0 &&
+ array[(i+1)*ydim*zdim + (j-1)*zdim + k] == 0 &&
+ array[(i-1)*ydim*zdim + j*zdim + k] == 0 &&
+ array[(i-1)*ydim*zdim + (j+1)*zdim + k] == 0 &&
+ array[(i-1)*ydim*zdim + (j-1)*zdim + k] == 0 &&
+ array[(i)*ydim*zdim + (j+1)*zdim + k] == 0 &&
+ array[(i)*ydim*zdim + (j-1)*zdim + k] == 0 &&
+
+ array[(i+1)*ydim*zdim + j*zdim + (k-1)] == 0 &&
+ array[(i+1)*ydim*zdim + (j+1)*zdim + (k-1)] == 0 &&
+ array[(i+1)*ydim*zdim + (j-1)*zdim + (k-1)] == 0 &&
+ array[(i-1)*ydim*zdim + j*zdim + (k-1)] == 0 &&
+ array[(i-1)*ydim*zdim + (j+1)*zdim + (k-1)] == 0 &&
+ array[(i-1)*ydim*zdim + (j-1)*zdim + (k-1)] == 0 &&
+ array[(i)*ydim*zdim + (j+1)*zdim + (k-1)] == 0 &&
+ array[(i)*ydim*zdim + (j-1)*zdim + (k-1)] == 0 &&
+ array[i*ydim*zdim + j*zdim + (k-1)] == 0 &&
+
+ array[(i+1)*ydim*zdim + j*zdim + (k+1)] == 0 &&
+ array[(i+1)*ydim*zdim + (j+1)*zdim + (k+1)] == 0 &&
+ array[(i+1)*ydim*zdim + (j-1)*zdim + (k+1)] == 0 &&
+ array[(i-1)*ydim*zdim + j*zdim + (k+1)] == 0 &&
+ array[(i-1)*ydim*zdim + (j+1)*zdim + (k+1)] == 0 &&
+ array[(i-1)*ydim*zdim + (j-1)*zdim + (k+1)] == 0 &&
+ array[(i)*ydim*zdim + (j+1)*zdim + (k+1)] == 0 &&
+ array[(i)*ydim*zdim + (j-1)*zdim + (k+1)] == 0 &&
+ array[i*ydim*zdim + j*zdim + (k+1)] == 0)
+ array[i*ydim*zdim + j*zdim + k] = 0;
+ EndIf
+ EndIf
+ EndFor
+ EndFor
EndFor
-//check
For i In {0:xdim-1:1}
-For j In {0:ydim-1:1}
-For k In {0:zdim-1:1}
-//Printf("array[%g][%g][%g] = %g", i, j, k,array[i*ydim*zdim + j*zdim + k]);
-If (array[i*ydim*zdim + j*zdim + k]==1)
+ For j In {0:ydim-1:1}
+ For k In {0:zdim-1:1}
+ //Printf("array[%g][%g][%g] = %g", i, j, k,array[i*ydim*zdim + j*zdim + k]);
+ If (array[i*ydim*zdim + j*zdim + k] == 1)
Call newsphere;
-EndIf
-EndFor
-EndFor
+ EndIf
+ EndFor
+ EndFor
EndFor
-c = newsl; Surface Loop(c) = {45,-6,37,41,50,49};
+// volume between electrodes with spheres carved out
+c = newsl; Surface Loop(c) = {43,-1,35,39,48,47};
v = newv; Volume(v) = {c, holes[]};
+Coherence; // make sure duplicate nodes are removed when sphere_dist=0
--
GitLab