G3.Reset();

%START BUILDING MODEL 

%Create soil domain as a box
s1 = G3.CreateBox(G3.Axis(-X/2, -Y/2, 0, 0, 0, 1), X, Y, Z);

%Draw rectangle for defining shallow foundation.
s2 = G3.CreateRectangle(G3.Axis(-L/2, -B/2, Z, 0, 0, 1), L, B);
G3.ZoomAll();

%Draw two lines for defining point loads in center of domain. Can also
%define vertex
s3_1 = G3.CreateLine(XYZ(-L/2,-B/2,Z),XYZ(L/2,B/2,Z))
s3_2 = G3.CreateLine(XYZ(-L/2,B/2,Z),XYZ(L/2,-B/2,Z))

%partition soil domain and rectangle
G3.Partition(G3.NewShapeCollection().Added(s1).Added(s2).Added(s3_1).Added(s3_2));

%Create materials
%Mohr-Coulomb material
MC = G3.CreateMaterial(MaterialCategory.Solids,'MC Sand');
%set properties
ps = G3.NewPropertyCollection()...
    .Added('Type', 'Mohr-Coulomb')...
    .Added('c', G3.MakeValueParameter(0))...    
    .Added('Phi', G3.MakeValueParameter(phi))...
    .Added('GammaDry', G3.MakeValueParameter(gamma))...
    .Added('Color', G3.Color(222, 200, 133));
%Add properties, ps, to MC material
G3.SetMaterialProperties(MC, ps)

%Mohr-Coulomb interface material
MCIF = G3.CreateMaterial(MaterialCategory.Solids,'MC IF');
%set properties
ps = G3.NewPropertyCollection()...
    .Added('Type', 'Mohr-Coulomb')...
    .Added('c', G3.MakeValueParameter(0))...    
    .Added('Phi', G3.MakeValueParameter(phi_inter))...
    .Added('GammaDry', G3.MakeValueParameter(gamma))...
    .Added('Color', G3.Color(32, 223, 89));
%ad properties, ps, to MCIF material
G3.SetMaterialProperties(MCIF, ps)

%FASD material 
FASD = G3.CreateMaterial(MaterialCategory.Solids,'FASD sand');
%set properties
ps = G3.NewPropertyCollection()...
    .Added('Type', 'FASD Sand')...
    .Added('Phi_tc', phitc)...
    .Added('Phi_ps_tc',phips_phitc)...
    .Added('GammaDry', G3.MakeValueParameter(gamma))...
    .Added('Color', G3.Color(200, 200, 100));
%ad properties, ps, to FASD material
G3.SetMaterialProperties(FASD, ps)

%Shell material - rigid
shellmat = G3.CreateMaterial(MaterialCategory.Shells,'rigid_shell');
%set properties
ps = G3.NewPropertyCollection()...
    .Added('Type', 'Shell')...
    .Added('ContinuumModel',ContinuumModel.Rigid)...
    .Added('Color', G3.Color(121, 137, 163));
G3.SetMaterialProperties(shellmat, ps)

%Apply features
%Boundbox defined as (x,y,z,dx,dy,dz) where dx,dy and dz >= 0
selection = G3.Selection(G3.BoundBox(-X/2,-Y/2,0,X,Y,Z), ShapeKind.Solid);

%Initilize empty features container
features = G3.NewFeatureCollection();
%Add material to feature container
if length(soil) == 2 & soil == 'MC'
    features.Add(G3.MakeSolid(selection.At(0), MC));
elseif soil == 'FASD'
   features.Add(G3.MakeSolid(selection.At(0), FASD));
end;
%and then material to model
G3.CreateFeatures(features);

%Shell for foundation
%Initilize empty features container
features = G3.NewFeatureCollection();
%Select surfaces
selection = G3.Selection(G3.BoundBox(0,0,Z), ShapeKind.Face);

%Define MC interfaces 
MCinter = G3.MakeInterface(Rinter, MCIF); %Rinter (= [0,1]) , material

%find the correct surface
for i=0:selection.Count-1
    xmin = selection.At(i).Bound.Min.X;
    xmax = selection.At(i).Bound.Max.X;
    ymin = selection.At(i).Bound.Min.Y;
    ymax = selection.At(i).Bound.Max.Y;
    if (abs(xmin--L/2)<1e-5 & abs(xmax-L/2)<1e-5) | (abs(ymin--B/2)<1e-5 & abs(ymax-B/2)<1e-5) 
        features.Add(G3.MakeShell(selection.At(i), shellmat,[],MCinter));
    end
end
%and then material to model
G3.CreateFeatures(features);

%Create loads
%Initilize empty features container
features = G3.NewFeatureCollection();
%Vertical force at centre
selection = G3.Selection(G3.BoundBox(0,0,Z), ShapeKind.Vertex);
features.Add(G3.MakeLoad(selection.At(0), LoadOption.Multiplier, 0, 0, -1));
G3.CreateFeatures(features);

% standard fixities
G3.CreateStandardFixities();

% Stage settings
if Adaptivity == 'Y'
ps = G3.NewPropertyCollection()...
    .Added('AnalysisType', AnalysisType.LimitAnalysis)...
    .Added('Name', 'LA01')...
    .Added('TargetElements',NoEL)...
    .Added('StartElements',StartEL)...
    .Added('UseAdaptivity',YesNo.Yes)...
    .Added('AdaptivitySteps',AdaptSteps);
else
ps = G3.NewPropertyCollection()...
    .Added('AnalysisType', AnalysisType.LimitAnalysis)...
    .Added('Name', 'LA01')...
    .Added('TargetElements',NoEL)...
    .Added('StartElements',StartEL)...
    .Added('UseAdaptivity',YesNo.No)...
    .Added('AdaptivitySteps',AdaptSteps);    
end;
G3.SetStageProperties(ps);

%saving file
if 0
    file = 'D:\tmp\example_c-with-results.g3x';
    % saving
    G3.SaveProject(file);
end

%Run analysis
result = G3.RunAnalysis();

%extract load multiplier from result set
LM = result.StageResults.At(0).CriticalValue.At(3).Value;

%save results to file
if 0
    file = 'D:\tmp\example_c-with-results.g3x';
    %saving
    G3.SaveProject(file);
end

Vmax = LM;