// -----------------------------------------------------------------------------
//
//  Gmsh Packed-bed
//
// -----------------------------------------------------------------------------

SetFactory("OpenCASCADE"); //In order to make the box

// We first create two cubes:
Box(1) = {0,0,0, 1,1,1};


// We then create the five spheres:
x0 = 0.02 ; y0 = 0.02 ; z0 = 0.02 ; r = 0.16 ; dx = 0.15; n = 3;
x = x0; y = y0; z = z0;
For t In {1:n}
  x += dx*1.06;
  For u In {1:n}
    y += dx*1.06;
    For v In {1:n}
      z += dx*1.06;

      Sphere(1 + (t-1)*n*n + (u-1)*n + v) = {x,y,z,r};
      Physical Volume((t-1)*n*n + (u-1)*n + v) = {1 + (t-1)*n*n + (u-1)*n + v};     
 
      z += dx;
    EndFor
    y += dx;
    z = z0;
  EndFor
  x += dx;
  y = y0;
EndFor

// Deleting the spheres from the geometry
BooleanDifference(n*n*n+2) = { Volume{1}; Delete; }{ Volume{2 : n*n*n+1}; Delete;};

Physical Volume(n*n*n+3) = v(#v()-1);

// Creating entities using constructive solid geometry is very powerful, but can
// lead to practical issues for e.g. setting mesh sizes at points, or
// identifying boundaries.

// To identify points or other bounding entities you can take advantage of the
// `PointfsOf' (a special case of the more general `Boundary' command) and the
// `In BoundingBox' commands.
//lcar1 = .100;
//lcar2 = .020;
//lcar3 = .030;
//eps = 1e-3;

Mesh.CharacteristicLengthMin = 0.025;
Mesh.CharacteristicLengthMax = 0.045;

// Assign a mesh size to all the points of all the volumes:
//Characteristic Length{ PointsOf{ Volume{:}; } } = lcar1;

// Override this constraint on the points of the five spheres:
//Characteristic Length{ PointsOf{ Volume{2:126}; } } = lcar3;

// Select the corner point by searching for it geometrically:
//p() = Point In BoundingBox{0.5-eps, 0.5-eps, 0.5-eps,
//                           0.5+eps, 0.5+eps, 0.5+eps};
//Characteristic Length{ p() } = lcar2;

Mesh.SubdivisionAlgorithm = 2;
Mesh.ElementOrder = 2;
Mesh.HighOrderOptimize = 1;


//+
Physical Surface("Inlet") = {4};
//+
Physical Surface("Outlet") = {2};
//+
Physical Surface("Walls") = {3, 1, 5, 6};
//+
Physical Surface("spheres") = {28, 27, 30, 29, 7, 26, 10, 25, 21, 20, 22, 11, 33, 14, 8, 31, 9, 32, 16, 13, 17, 18, 15, 12, 24, 19, 23};
