How to change boundary conditions during AceFEM analysis?

Question

I'm trying to prestress a sample by a compressive strain of 10 percent and then stretch it by 30 percent. The way I intend to do this is apply the prestress in one step and then have the code run and strain the sample by 30 percent as usual. What I'm having issues with is that in every iteration my code both applies the prestress and the 30 percent strain, whereas I only want the prestress applied once in the begining of the code and then let the code apply the 30 percent stress using the increments outlined in the code below:

CompressionStrainX = -0.3;
CompressionStrainY = 0;
Prestress = -0.1;
AnalysisSwitch = 1; (* UNIFORM COMPRESSION Case *)

SolverSwitch = 1; (* 1 -in core; 2 -out-of-core *)
SolverMatrixType = -2;
For[AnalysisIndex = 1, AnalysisIndex <=  Length[AnalysisBatch], 
 AnalysisIndex++,

 Analysis name;
 AnalysisName = AnalysisBatch[[AnalysisIndex, 1]];
 E0StratumCorneum = AnalysisBatch[[AnalysisIndex, 2]];
 AnalysisSwitch = AnalysisBatch[[AnalysisIndex, 3]];

 Start analysis;
 SMTInputData[];

 SMTAddDomain[{
   {"Stratum-Corneum", EFormulationStratumCorneum,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], 
     "cSC*" -> cSC, "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, 
     "\[Rho] *" -> \[Rho] , "\[Beta] *" -> \[Beta] , 
     "\[Delta] *" -> \[Delta] , "Eg" -> Eg}},
   {"Living-Epidermis", EFormulationLivingEpidermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , 
     "Eg" -> Eg}},
   {"Papillary-Dermis", EFormulationPapillaryDermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , 
     "Eg" -> Eg}},
   {"Reticular-Dermis", EFormulationReticularDermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , "Eg" -> Eg}}
   }];
 (*{"Living-Epidermis",EFormulationLivingEpidermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q *"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}},
 {"Papillary-Dermis",EFormulationPapillaryDermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q *"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}},
 {"Reticular-Dermis",EFormulationReticularDermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q*"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}}*)

 SMTMesh["Stratum-Corneum", "H1", {Xnd, Ynd, ZndStratumCorneum}, 
  StratumCorneumMastermesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Living-Epidermis", "H1", {Xnd, Ynd, ZndLivingEpidermis}, 
  LivingEpidermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Papillary-Dermis", 
  "H1-rf", {Xnd/3, Ynd/3, ZndPapillaryDermis}, 
  PapillaryDermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Reticular-Dermis", 
  "H1-rf", {Xnd/9, Ynd/9, ZndReticularDermis + 1}, 
  ReticularDermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];


 Extract boundaries of the mesh;
 {MeshXmin, MeshXmax} = {Min[SMTNodes[[All, 2]]], 
   Max[SMTNodes[[All, 2]]]};
 {MeshYmin, MeshYmax} = {Min[SMTNodes[[All, 3]]], 
   Max[SMTNodes[[All, 3]]]};
 {MeshZmin, MeshZmax} = {Min[SMTNodes[[All, 4]]], 
   Max[SMTNodes[[All, 4]]]};

 Switch[AnalysisSwitch, 1,
  UNIAXIAL COMPRESSION X;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0];
  SMTAddEssentialBoundary["Y" == MeshYmin &, 2 -> 0];
  (*SMTAddEssentialBoundary["Y"\[Equal]MeshYmax&,2\[Rule] 0];*)

  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition X;
  Xdisplacement = -Prestress*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];
  ,
  2,(*This is the step that I only want applied once in one session*)
  EQUIBIAXIAL COMPRESSION;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0];
  SMTAddEssentialBoundary["Y" == MeshYmin &, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition X;
  Xdisplacement = -CompressionStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];

  Displacement condition Y;
  Ydisplacement = -CompressionStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["Y" == MeshYmax &, 2 -> Ydisplacement];
  ,
  3,
  SHEAR ALONG Y;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = ShearStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> 0, 
   2 -> ShearDisplacement];

  ,
  4,
  SHEAR ALONG X;
   Essential boundary conditions;
   SMTAddEssentialBoundary["Y" == MeshYmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = ShearStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["Y" == MeshYmax &, 1 -> ShearDisplacement, 
   2 -> 0];
  ,
  5,
  SHEAR AND COMPRESSION;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = -ShearStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 2 -> ShearDisplacement];

  Displacement condition X;
  Xdisplacement = -CompressionStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];

  ];

 Switch[SolverSwitch,
  In core PARDISO solver;
  1,
  SMTAnalysis["Output" -> AnalysisName <> ".out", 
   "OptimizeDll" -> True];
  SMTSetSolver[5, SolverMatrixType, {{60, 0}}];
  ,
  2,
  Out of core PARDISO solver;
  SMTAnalysis["Output" -> AnalysisName <> ".out", 
   "OptimizeDll" -> True];
  SMTSetSolver[5, SolverMatrixType, {{60, 2}}];
  ];
 EXTRACTION OF NODES AND ASSOCIATED DATA;
 Identify nodes of each skin layer;
 StratumCorneumNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Stratum-Corneum"}]}];
 LivingEpidermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Living-Epidermis"}]}];
 PapillaryDermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Papillary-Dermis"}]}];
 ReticularDermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Reticular-Dermis"}]}];

 Select nodes belonging to SC top lines;
 NodesSelectionYmax = 
  Table[4*i + (Xnd + 1)*4*(Xnd), {i, 1, Xnd + 1}];
 (*NodesSelectionYmid=Table[4*i+(Xnd+1)*4*(Xnd/2),{i,1,Xnd+1}]*) (* \
EVEN number of divisions *)

 NodesSelectionYmid = 
  Table[4*i + (Xnd + 1)*4*(Xnd + 1)/2, {i, 1, Xnd + 1}];
 NodesSelectionYmin = Table[4*i, {i, 1, Xnd + 1}];

 Initialise converged step counter;
 ic = 0;

 Output result file at zero step ;
 SMTDump[AnalysisName <> "_" <> ToString[ic]];


 Leading parameter;
 V0 = Velocity0; (* velocity *)

 Analysis characteristics;
 \[Gamma]0 = 0.01;
 \[CapitalDelta]\[Gamma]Min = 10^-6 V0;
 \[CapitalDelta]\[Gamma]Max = 0.1;
 \[Gamma]Max = 1;

 tolNR = 10^-8;
 tolNR0 = 10^-4; (* relaxed tolerance to go through critical points *)

 maxNR = 25; targetNR = 12;

 \[CapitalDelta]\[Lambda]Min = 10^-6 V0;
 \[CapitalDelta]\[Lambda]LineSearch = 
  100*\[CapitalDelta]\[Lambda]Min;
 stressgraph = {{0, 0}};
 SMTNextStep["t[\[Gamma]]" -> (#/V0 &), 
  "\[Lambda][\[Gamma]]" -> (# &), 
  "\[CapitalDelta]\[Gamma]" -> \[Gamma]0];
 plot\[Sigma]e = {{0, 0}};
 While[
  While[
   step = SMTConvergence[
     Conditional tolerance;
     Conditional tolerance;
     If[SMTRData["Multiplier"] > 0.9975*\[Gamma]Max, tolNR, tolNR0],
     maxNR, {"Adaptive \[Gamma]", 
      targetNR, \[CapitalDelta]\[Gamma]Min, \
\[CapitalDelta]\[Gamma]Max, \[Gamma]Max}]
   , 
   (*If[SMTIData["Iteration"]>4 || SMTRData[
   "MultiplierIncrement"]>\[CapitalDelta]\[Lambda]LineSearch,SMTRData[
   "LineSearchStepLength",1],SMTRData["LineSearchStepLength",0.1]];*)

      SMTNewtonIteration[];
   ];
  If[step[[4]] === "MinBound", 
   SMTStatusReport[
     "\[CapitalDelta]\[Gamma] < \[CapitalDelta]\[Gamma]Min"];];
  step[[3]] 
  , If[step[[1]]
   , SMTStepBack[]
   ];

  Display converged configuration;

  SMTShowMesh["DeformedMesh" -> True, Axes -> True, 
   AxesLabel -> {"X", "Y", "Z"}, "Field" -> "w", "Show" -> "Window"];
  (*AppendTo[plot\[Sigma]e,{SMTPostData["\[Lambda]11",{0,0,MeshZmax}],
  SMTPostData["\[Sigma]11(MPa)",Point[{0,0,MeshZmax}]]}];*)

  Save current converged increment configuration;
  ic = ic + 1;
  SMTDump[AnalysisName <> "_" <> ToString[ic]];

  SMTNextStep["t[\[Gamma]]" -> (#/V0 &), 
   "\[Lambda][\[Gamma]]" -> (# &), 
   "\[CapitalDelta]\[Gamma]" -> step[[2]]]
  ];


 SMTNextStep[100, 0];
 SMTNewtonIteration[];

 SMTDump[AnalysisName <> "_" <> ToString[ic] <> 
   ".After.Dissipation"];

 SMTSimulationReport[];

 ]

This is just the loop that I have at the moment, the actual code is much longer than this. I need step two within AnalysisSwitch to only run once and I'm not really sure how to get this to work. I was thinking of something along the lines of SMSIf[switch == 2 && (AnalysisSwitch = 1), SMSReturn[]]; but I have tried that and its not really working. Any ideas would be much appreciated. I apologize for not posting the whole code but it's quite long.

Answer 1

I am not sure if I have understood the question but let me give a minimal example of how to change the type of boundary conditions during AceFEM analysis. We will squeeze a box in one direction and then stretch it in the other direction.

First we define functions which will make the whole process more readable and organized.

Get["AceFEM`"]

setupPhase[noElements_Integer] := Module[
  {},
  SMTInputData[];
  (* Element libraries are downloaded from AceShare server. *)
  SMTAddDomain["Domain1", "OL:SED3H1DFHYH1NeoHooke", {}];
  SMTAddMesh[
   Raster3D@{
     {{{0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {1, 1, 0}}},
     {{{0, 0, 1}, {1, 0, 1}}, {{0, 1, 1}, {1, 1, 1}}}
     },
   "Domain1", "H1", noElements*{1, 1, 1}
   ];
  SMTAddEssentialBoundary[{{"X" == 0 &, 1 -> 0}, {"Y" == 0 &, 2 -> 0}, {"Z" == 0 &, 3 -> 0}}];
  SMTAnalysis[]
];

(* Handy function for visualization. *)
showMesh[imageSize_] := SMTShowMesh[
  "DeformedMesh" -> True,
  "BoundaryConditions" -> True,
  "Field" -> "u",
  "Contour" -> {False, 0, 0.5, 4},
  Axes -> True,
  AxesStyle -> 12,
  ImageSize -> imageSize,
  PlotRange -> {{-0.1, 2}, {-0.1, 1.1}, {-0.1, 1.1}}
];

(* Increase multiplier from 0 to 1 in specified number of steps.*)
analysisLoop[noSteps_Integer] := Module[
  {Δλ, results},
  Δλ = 1./noSteps;
  results = {};
  Do[
   SMTNextStep["Δλ" -> Δλ];
   While[SMTConvergence[10^-8, 15], SMTNewtonIteration[]];
   (* Collect any results here. *)
   AppendTo[results, showMesh[Medium]],
   {i, 1, noSteps}
   ];
  (* Results are returned from the function. *)
  results
];

Then we call the functions for setup phase and analysis loop. Between two loops we change the boundary conditions and reset multiplier (lambda) to start again from zero. Keep in mind that this is only one of many ways to achieve the same result.

setupPhase[4]

(* Set intensity of essential BC (displacement). *)
SMTAddEssentialBoundary[{"Y" == 1 &, 2 -> -0.25}];
(* Run analysis and collect results. *)
results = analysisLoop[5];
(* Reset intensity of the same essential BC (displacement is zero). *)
SMTAddEssentialBoundary[{"Y" == 1 &, 2 -> 0}, "Set" -> True];

(* Reset BC multiplier value. *)
SMTRData["Multiplier", 0.];
(* Set intensity of essential BC in other direction. *)
SMTAddEssentialBoundary[{"X" == 1 &, 1 -> 0.5}];
results = Join[results, analysisLoop[10]];

ListAnimate[results]