// (c) Yao Koutsawa, 25 February 2018
// face centered cubic (FCC) ellipsoids unit cell
SetFactory("OpenCASCADE");

Macro FCCUnitCell    
    // Spheres    
    coeffs_1 = {1,-1,1,-1, 1,-1,1,-1,  1,-1, 0, 0, 0, 0};
    coeffs_2 = {1,1,-1,-1, 1,1,-1,-1,  0, 0, 1,-1, 0, 0};
    coeffs_3 = {1,1,1,1, -1,-1,-1,-1,  0, 0, 0, 0, 1, -1};
    
    size = #coeffs_1() - 1;
    
    For i In {0:size}
        center = {0.5*coeffs_1[i], 0.5*coeffs_2[i], 0.5*coeffs_3[i]};
        v = newv; Sphere(v) = {center(0),center(1), center(2), r};
        iVolumes() += v;
    EndFor

    // Box
    v_box = newreg; Box(v_box) = {-0.5,-0.5,-0.5, 1,1,1};
    BooleanIntersection{ Volume{v_box}; Delete;}{ Volume{iVolumes()}; Delete; }
    
    // Box
    v_box = newreg; Box(v_box) = {-0.5,-0.5,-0.5, 1,1,1};
    BooleanDifference{ Volume{v_box}; Delete;}{ Volume{iVolumes()};}
    BooleanFragments { Volume{v_box,iVolumes()}; Delete;}{}
    mVolume() += v_box;
    
    Dilate { {0.,0.,0.}, {1/aspect_ratio_31, 1/aspect_ratio_32, 1} } { 
        Volume{mVolume(), iVolumes()};
    }
     
    Physical Volume(matrix) = {mVolume()};
    Physical Volume(fiber) =  {iVolumes()};
Return

// Options
Mesh.HighOrderOptimize = 1;
Mesh.Optimize = 1;
Mesh.OptimizeNetgen = 1;
Mesh.SecondOrderIncomplete = 1;
Mesh.Algorithm3D = 1;
Mesh.Algorithm = 2; 
Geometry.AutoCoherence = 1;
Geometry.CopyMeshingMethod = 1;
Geometry.HighlightOrphans = 1;

// Initialization
iVolumes() = {};
mVolume() = {};

// Volume fraction: vf
cl = 0.1;
vf = 0.15;
fiber = 1;
matrix = 2;
nparticles = 4.;
r = (3.*vf/(nparticles*4.*Pi))^(1.0/3.0);

Mesh.CharacteristicLengthMin = cl;
Mesh.CharacteristicLengthMax = cl;

// Aspect ratios
aspect_ratio_31 = 1/5; 
aspect_ratio_32 = 1;

Call FCCUnitCell;
// gmsh -3 ellipsoids_fcc_ruc.geo -order 2 -o ellipsoids_fcc_ruc.msh

