I feel like you have helped me every step of the way so I just want to say thank you 
Setting “Soil as Solid” as default was my initial ide, does not work for me (changed in both net6.0 folder / FEM-Design API and when iterating). Have you tested it?
When I test the same node number is chosen every iteration.
What I find interesting is that any of the calculation tabs Analysis → Composite generate a mesh before calculation with node number 558 at my result point.
But when using Finite elements “Generate” button the result points gets a different node number (1808).
My code so far:
using FemDesign;
using FemDesign.Calculate;
using FemDesign.Materials;
using FemDesign.Results;
using FemDesign.Shells;
using FemDesign.Soil;
using StruSoft.Interop.StruXml.Data;
namespace FemDesign.SoilIter //ITERATE SOIL STRENGTH
{
internal class Program
{
static void Main()
{
// READ AND SAVE MODEL NAME
string struxmlPath = "GEO.struxml";
Model model = Model.DeserializeFromFilePath(struxmlPath);
string struxmlPath_saved = "GEO_saved.str";
// 1. READ SOIL MATERIAL TO ANALYSE
var soil = model.Materials.Material.Where(m => m.Name == "Soil_ver").ToList()[0];
var drained = (Material_typeStratumBehaviourDrained)soil.Stratum.Behaviour.Item;
var C_k_org = drained.C_k;
var Phi_k_org = drained.Phi_k;
var Phi_cvk_org = drained.Phi_cvk;
int i_failure = 1; // SAVEING .STRUXML FILE AT FAILURE (ADJUST AFTER FAILIURE ITERATION IS KNOW)
List<double> MsfList_SDL = new List<double> { 0.0, 0.4, 0.5, 0.55, 0.56, 0.57, 0.58, 0.6, 0.66 };
//List<double> MsfList_SUR = new List<double> { 0.0, 0.1, 0.2, 0.23, 0.25, 0.27, 0.30, 0.35 };
// ANALYSIS SETTINGS
var combItem = new Calculate.CombItem
{
ImpfRqd = 0,
StabRqd = 0,
NLE = false,
PL = false,
NLS = true,
Cr = false,
f2nd = false,
Im = 0,
Waterlevel = 0,
};
var combItems = new List<Calculate.CombItem>
{
combItem,
};
var comb = new Calculate.Comb
{
NLEmaxiter = 30,
PLdefloadstep = 20,
PLminloadstep = 2,
PlKeepLoadStep = true,
PlTolerance = 1,
PLmaxeqiter = 50,
PlShellLayers = 10,
NLSMohr = true,
NLSinitloadstep = 10,
NLSminloadstep = 10,
NLSactiveelemratio = 5,
NLSplasticelemratio = 5,
CRloadstep = 20,
CRmaxiter = 30,
CRstifferror = 2,
CombItem = combItems,
};
Analysis analysis = new Analysis(comb: comb, calcCase: true, calcComb: true);
double myPointX = 4;
double myPointY = 1.5;
double myPointZ = 0;
string LCName = "LC1";
// ITERATION & ANALYSIS PROCESS
using (var femDesign = new FemDesignConnection(keepOpen: true))
{
for (int i = 0; i < MsfList_SDL.Count; i++)
//for (int i = 0; i < MsfList_SUR.Count; i++)
{
Console.WriteLine($"Iteration {i}");
// Setting material safty factor and reduce soil strength
double Msf = 1 + MsfList_SDL[i];
// double Msf = 1 + MsfList_SUR[i];
drained.C_k = C_k_org / Msf;
drained.Phi_k = Math.Atan(Math.Tan(Phi_k_org) / Msf);
drained.Phi_cvk = Math.Atan(Math.Tan(Phi_cvk_org) / Msf);
// Print of friction angle and material safty factor
Console.WriteLine($"Msf: {Math.Round(Msf, 2)}");
Console.WriteLine($"C_k: {Math.Round(drained.C_k, 2)}");
Console.WriteLine($"Phi_k: {Math.Round(180 / Math.PI * drained.Phi_k, 2)}");
Console.WriteLine($"Phi_cvk: {Math.Round(180 / Math.PI * drained.Phi_cvk, 2)}");
// Run analysis
femDesign.Open(model);
// --- Run calculation with soil as solid elements (NOT IMPLEMENTED YET) ---
femDesign.RunAnalysis(analysis); // ADD BREAKPOINT AND MANUALLY SETT "CALCULATE SOIL AS SOILD ELEMENT" (TEMPERORY FIX)
// Find node number at myPoint
var feaNode = femDesign.GetFeaNodes();
var myPoint = new Geometry.Point3d(myPointX, myPointY, myPointZ);
var myNodeId = feaNode.Where(n => (new Geometry.Point3d(n.X, n.Y, n.Z) - myPoint).Length() < 0.001).FirstOrDefault().NodeId;
Console.WriteLine($"Node Id: {myNodeId}");
var disp = femDesign.GetResults<NodalDisplacement>();
var myRes = disp.Where(r => r.CaseIdentifier == LCName).Where(x => x.NodeId == myNodeId).ToList();
// Print node displacement
foreach (var myRe in myRes)
{
var dispID = myRe.Id;
var dispX = Math.Round(myRe.Ex,3);
var dispY = Math.Round(myRe.Ey,3);
var dispZ = Math.Round(myRe.Ez,3);
var dispAbs = Math.Round(Math.Sqrt(Math.Pow(dispX, 2) + Math.Pow(dispY, 2) + Math.Pow(dispZ, 2)), 3);
var dispCId = myRe.CaseIdentifier;
Console.WriteLine($"{dispID,5} | ex={dispX,5}, ey={dispY,5}, ez={dispZ,5}, eabs={dispAbs,5} | {dispCId,5}");
}
// Save failure iteration
if (i == i_failure)
{
Console.WriteLine($"Saving model at i = {i}");
femDesign.Save(struxmlPath_saved);
}
Console.WriteLine("\n");
}
Console.WriteLine("FEM-Design is closing.");
}
Console.WriteLine("Press any key to exit.");
Console.ReadKey();
}
}
}
Comparing manual calculation with the iteration loops is correct for Msf 1 / 1.5 but differ for 1.4 / 1.5 but maybe this is just because of the large displacement 
Comparing it to the example in Geotechnical module there is a bit of difference that I will need think about further …
