// ================================================================
// to build the mesh in command line :
// gmsh script_Domaine.geo  -3 -format med -o myMesh.med
// ================================================================

// This defines the CAD Kernel - since we operate on brep, it *must* be this one
SetFactory("OpenCASCADE");

// First we import the brep
Merge "Domaine.brep";

// Define some basic groups
Physical Volume("Rocher") = {1};
Physical Volume("Pliocene") = {3};
Physical Volume("Alluvion") = {2};
Physical Surface("Base") = {2};

psrc() = {};

// Create new points and embed them in the mesh
c0[] = Point{4}; // c0 contains the coords in indexes 0,1,2
c1[] = Point{3};
p=newp; // newp is command to generate a new point number - we store it in p
Point(p)={0,0, (c1[2]-c0[2])/2};
psrc() += p;
Physical Point("PGHM") = {p};
//
c0[] = Point{3};
c1[] = Point{16};
p=newp;
Point(p)={(c1[0]-c0[0])/2,0, -(c1[2]-c0[2])/2};
psrc() += p;
Physical Point("PMM") = {p};
//
c0[] = Point{3};
c1[] = Point{11};
p=newp;
Point(p)={(c1[0]-c0[0])/2,0, c0[2]};
psrc() += p;
Physical Point("PMT") = {p};
//
c0[] = Point{4};
c1[] = Point{16};
p=newp;
Point(p)={(c1[0]-c0[0])/2,0, 0};
psrc() += p;
Physical Point("PMB") = {p};

BooleanFragments{ Volume{:}; Delete; }{ Point{psrc()}; Delete; }

// ================================
// Rocher
// ================================
// Get the bounding box
bbox() = BoundingBox Volume{1};
xmin = bbox(0);
rch_ymin = bbox(1);
zmin = bbox(2);
xmax = bbox(3);
rch_ymax = bbox(4);
zmax = bbox(5);

// Elements size according to linear law - V=800->2500 (*3.1)
rch_min = 80; // This is the only parameter to set the mesh size
rch_max = 3.1*rch_min;
Field[1] = MathEval;
Field[1].F = Sprintf("(%g) * y/(%g) + (%g) - (%g)*y/(%g)", rch_max, rch_ymin, rch_min, rch_min, rch_ymin);

// Assign sizes to considered volume, faces and edges
Field[2] = Restrict;
Field[2].IField = 1;
Field[2].RegionsList = {1};
fac() = Abs(Boundary { Volume{1}; });
Field[2].FacesList = {fac()};
edg() = Abs(Boundary { Surface{fac()}; });
Field[2].EdgesList = {edg()};

// ================================
// Pliocene
// ================================
bbox() = BoundingBox Volume{3};
xmin = bbox(0);
plc_ymin = bbox(1);
zmin = bbox(2);
xmax = bbox(3);
plc_ymax = bbox(4);
zmax = bbox(5);

// Elements size according to linear law - V=350->800 (*2.3)
plc_min = rch_min/2.3;
plc_max = rch_min;
Field[3] = MathEval;
Field[3].F = Sprintf("(%g) * y/(%g) + (%g) - (%g)*y/(%g)", plc_max, plc_ymin, plc_min, plc_min, plc_ymin);

Field[4] = Restrict;
Field[4].IField = 3;
Field[4].RegionsList = {3};
fac() = Abs(Boundary { Volume{3}; });
Field[4].FacesList = {fac()};
edg() = Abs(Boundary { Surface{fac()}; });
Field[4].EdgesList = {edg()};

// ================================
// Alluvion
// ================================
bbox() = BoundingBox Volume{2};
xmin = bbox(0);
alv_ymin = bbox(1);
zmin = bbox(2);
xmax = bbox(3);
alv_ymax = bbox(4);
zmax = bbox(5);

// Elements size according to linear law - V=200->350 (*1.8)
alv_min = plc_min/1.8;
alv_max = plc_min;
Field[5] = MathEval;
Field[5].F = Sprintf("(%g) * y/(%g) + (%g) - (%g)*y/(%g)", alv_max, alv_ymin, alv_min, alv_min, alv_ymin);

Field[6] = Restrict;
Field[6].IField = 5;
Field[6].RegionsList = {2};
fac() = Abs(Boundary { Volume{2}; });
Field[6].FacesList = {fac()};
edg() = Abs(Boundary { Surface{fac()}; });
Field[6].EdgesList = {edg()};

// ================================
// Synthesis
// ================================
Field[7] = Min;
Field[7].FieldsList = {2, 4, 6};
Background Field = 7;

Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.ElementOrder = 1;
//Mesh.OptimizeNetgen = 1;
Mesh.Algorithm3D = 10;
//Mesh.SurfaceFaces = 1;
