
// BOX WITH SPHERES - 3D //
// --------------------- //
//
// Set global parameters

DOM_SIZE = 30. ;
lc       =  1.5 ;

// Lists of bubble center coordinates (x,y,z) and bubble radii

bubble_xs[] = { -4. ,  2. } ;
bubble_ys[] = {  7. , -6. } ;
bubble_zs[] = { -4. ,  7. } ;
bubble_rs[] = {  5. ,  5. } ;

// Definition of a macro generating a bubble

Macro GenerateBubble

	pp = newp ;
	Point( pp   ) = { xc     , yc     , zc     , lc } ;
	Point( pp+1 ) = { xc + rp, yc     , zc     , lc } ;
	Point( pp+2 ) = { xc     , yc + rp, zc     , lc } ;
	Point( pp+3 ) = { xc - rp, yc     , zc     , lc } ;
	Point( pp+4 ) = { xc     , yc - rp, zc     , lc } ;
	Point( pp+5 ) = { xc     , yc     , zc - rp, lc } ;
	Point( pp+6 ) = { xc     , yc     , zc + rp, lc } ;

	lp = newl ;
	Circle( lp    ) = { pp+1 , pp , pp+2} ;
	Circle( lp+1  ) = { pp+2 , pp , pp+3} ;
	Circle( lp+2  ) = { pp+3 , pp , pp+4} ;
	Circle( lp+3  ) = { pp+4 , pp , pp+1} ;
	Circle( lp+4  ) = { pp+2 , pp , pp+5} ;
	Circle( lp+5  ) = { pp+5 , pp , pp+4} ;
	Circle( lp+6  ) = { pp+4 , pp , pp+6} ;
	Circle( lp+7  ) = { pp+6 , pp , pp+2} ;
	Circle( lp+8  ) = { pp+1 , pp , pp+6} ;
	Circle( lp+9  ) = { pp+6 , pp , pp+3} ;
	Circle( lp+10 ) = { pp+3 , pp , pp+5} ;
	Circle( lp+11 ) = { pp+5 , pp , pp+1} ;

	
	llp = newll ;
	Line Loop( llp   ) = {   lp+1  ,   lp+7  , -(lp+9) } ;
	Line Loop( llp+1 ) = {   lp+9  ,   lp+2  ,   lp+6 } ;
	Line Loop( llp+2 ) = { -(lp+7) , -(lp+8) ,   lp   } ;
	Line Loop( llp+3 ) = { -(lp+10), -(lp+1) ,   lp+4 } ;
	Line Loop( llp+4 ) = { -(lp+4) , -(lp+11), -lp   } ;
	Line Loop( llp+5 ) = { -(lp+2) ,   lp+10 ,   lp+5 } ;
	Line Loop( llp+6 ) = { -(lp+6) ,   lp+3  ,   lp+8 } ;
	Line Loop( llp+7 ) = { -(lp+3) ,   lp+11 , -(lp+5) } ;

	
	sp = news ;
	Surface( sp   ) = { llp   } ;
	Surface( sp+1 ) = { llp+1 } ;
	Surface( sp+2 ) = { llp+2 } ;
	Surface( sp+3 ) = { llp+3 } ;
	Surface( sp+4 ) = { llp+4 } ;
	Surface( sp+5 ) = { llp+5 } ;
	Surface( sp+6 ) = { llp+6 } ;
	Surface( sp+7 ) = { llp+7 } ;

	slp = newsl ;
	Surface Loop( slp ) = {sp+7, sp+6, sp+1, sp, sp+3, sp+5, sp+4, sp+2} ;

	vp = newv ;
	Volume( vp ) = { slp } ;

	bubble_surf_loops[] += { slp } ;
	bubble_volumes[]    += { vp } ;

Return

// Domain construction

xd  = DOM_SIZE / 2.  ;
yd  = DOM_SIZE / 2.  ;
zd  = DOM_SIZE / 2.  ;
lcd = DOM_SIZE / 15. ;

pd = newp ;
Point( pd   ) = { -xd, -yd, -zd, lcd } ;
Point( pd+1 ) = {  xd, -yd, -zd, lcd } ;
Point( pd+2 ) = {  xd,  yd, -zd, lcd } ;
Point( pd+3 ) = { -xd,  yd, -zd, lcd } ;
Point( pd+4 ) = { -xd, -yd,  zd, lcd } ;
Point( pd+5 ) = {  xd, -yd,  zd, lcd } ;
Point( pd+6 ) = {  xd,  yd,  zd, lcd } ;
Point( pd+7 ) = { -xd,  yd,  zd, lcd } ;

ld = newl ;
Line( ld    ) = { pd+3 , pd+2 } ;
Line( ld+1  ) = { pd+2 , pd+1 } ;
Line( ld+2  ) = { pd+1 , pd   } ;
Line( ld+3  ) = { pd   , pd+3 } ;
Line( ld+4  ) = { pd+4 , pd+5 } ;
Line( ld+5  ) = { pd+5 , pd+6 } ;
Line( ld+6  ) = { pd+6 , pd+7 } ;
Line( ld+7  ) = { pd+7 , pd+4 } ;
Line( ld+8  ) = { pd   , pd+4 } ;
Line( ld+9  ) = { pd+3 , pd+7 } ;
Line( ld+10 ) = { pd+1 , pd+5 } ;
Line( ld+11 ) = { pd+2 , pd+6 } ;

lld = newll ;
Line Loop( lld   ) = { ld+2,  ld+3   ,  ld     ,  ld+1    } ;
Line Loop( lld+1 ) = { ld+4,  ld+5   ,  ld+6   ,  ld+7    } ;
Line Loop( lld+2 ) = { ld+8, -(ld+7) , -(ld+9) , -(ld+3)  } ;
Line Loop( lld+3 ) = { ld+6, -(ld+9) ,  ld     ,  ld+11   } ;
Line Loop( lld+4 ) = { ld+10, ld+5   , -(ld+11),  ld+1    } ;
Line Loop( lld+5 ) = { ld+4, -(ld+10),  ld+2   ,  ld+8    } ;

psd = news ;
Plane Surface( psd   ) = { lld   } ;
Plane Surface( psd+1 ) = { lld+1 } ;
Plane Surface( psd+2 ) = { lld+2 } ;
Plane Surface( psd+3 ) = { lld+3 } ;
Plane Surface( psd+4 ) = { lld+4 } ;
Plane Surface( psd+5 ) = { lld+5 } ;

sld = newsl ;
Surface Loop( sld ) = { psd+1, -(psd+5), -(psd+4), -(psd+3), psd+2, psd } ;

// Generate bubbles => ready storage lists

bubble_surf_loops[] = {} ;
bubble_volumes[]    = {} ;

// Generate bubbles => Loop on bubbles
For i_bubble In { 0:(#bubble_xs[] - 1) }

	// Generate bubbles => set bubble parameter
	xc = bubble_xs[ i_bubble ] ;
	yc = bubble_ys[ i_bubble ] ;
	zc = bubble_zs[ i_bubble ] ;
	rp = bubble_rs[ i_bubble ] ;

	// Generate bubbles => Call macro
	Call GenerateBubble ;

EndFor

// Finalize domain volume accounting for the holes created
// by the particles

vd = newv ;
Volume( vd ) = { sld, bubble_surf_loops[] } ;

Recombine Volume{ bubble_volumes[], vd } ;

// Coloring

Recursive Color Orange { Volume{ bubble_volumes[] } ; }
Recursive Color SkyBlue { Volume{ vd } ; }

// Physical tags

Physical Volume( "bubbles", 1000 ) = { bubble_volumes[] } ;
Physical Volume( "medium", 3000 ) = { vd } ;

// Mesh Configuration

Mesh.Algorithm = 8;
//Mesh.Algorithm3D = 1;
Mesh.RecombinationAlgorithm = 3;
Mesh.SubdivisionAlgorithm = 2;