'eaf' is an initial deformation which provide 'divsig0' According to its expression

eaf[kk_] := 
  Table[If[pcontour[[i]][[1]] + 
      uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] > 
     pillarscenters[[i]][[1]], t/100, 0], {i, 
    Dimensions[alpha][[1]]}];

I would like that such an initial deformation stops once the the sum

pcontour[[i]][[1]] + 
      uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] 

is less than

pillarscenters[[i]][[1]]

By doing so, when plotting the displacement of a point in the mesh, this should increase until a certain time point and than stay constant, whereas when plotting 'eaf' this should increase until the same time point and than drop to zero according to the If function

However, when I do 1 iteration it works for eaf but not for the displacement which continously increase over time and I do not understand why since if I test the condition

pcontour[[i]][[1]] + 
          uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] > 
         pillarscenters[[i]][[1]]

it evaluates to true for the first time steps and false for the rest of the interval, than the displacement should be stop too...

What problem are you solving? What do all these If[] mean?

Check $Version. My version "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)" works good. I calculated 8 iterations. The process branches at t>10.4. I attach the file.

Attachments:

uv8.nb

Remove all unnecessary. The last lines of code should be such

ci = {u[0, x, y] == 0., v[0, x, y] == 0.};
cid = {Derivative[1, 0, 0][u][0, x, y] == 0, 
   Derivative[1, 0, 0][v][0, x, y] == 0};
uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 2; tf = 15;
Do[{uif[i], vif[i]} = 
   NDSolveValue[{Activate[
      divsig - (divsig0 /. uif -> uif[i - 1]) == inertia + fadh], ci, 
     cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, 
    Method -> {"PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> {"MeshOrder" -> 2, 
            "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, 
    MaxStepSize -> 0.5];, {i, 1, kf}]

Run the code and make sure the solution on the first and second iteration is different at t>10.4. Therefore it is necessary to increase kf.

Table[ContourPlot[uif[kf][t, x, y], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", 
  Contours -> 20, PlotLabel -> Row[{"t=", t}]], {t, 2, tf, 2}]
Plot[Evaluate[Table[uif[i][t, -2*10^-6, -2*10^-6], {i, 1, kf}]], {t, 
  0, tf}, AxesLabel -> {"t", "u"}]

It is bizarre, when I open the file from scratch, it is not working anymore... I'll try to fix it, but actually even when it worked, it was very slow...was it the same on your machine?

Thanks rachele

here the whole code

Clear["Global`*"]
Needs["NDSolve`FEM`"]

c1hs[x_, scale_] = 
  If[(x/scale > -1), 1, 0] If[(x/scale < 1), 1, 
     0] (0.5 + x/scale (0.75 - 0.25 (x/scale)^2)) + 
   If[(x/scale >= 1), 1, 0];
Id = IdentityMatrix[2];
Tens[a_, b_] := Outer[Times, a, b];

ix = {1, 0};
iy = {0, 1};
teta = 0;
ndirx = Cos[teta];
ndiry = Sin[teta];
Vloc = Array[Subscript[v, ##] &, {2, 2}];
Vloc[[1 ;; 2, 1]] = {ndirx, ndiry};
Vloc[[1 ;; 2, 2]] = {-ndiry, ndirx};
V = Vloc;
rcell = 5*10^(-6);
cell = ImplicitRegion[((x)^2 + (y)^2 <= rcell^2), {x, y}];
substrate = 
  ImplicitRegion[-40*10^(-6) < x < 40*10^(-6) && -40*10^(-6) < y < 
     40*10^(-6), {x, y}];
cellcenter = {0, 0};

rigid = c1hs[-x + 8*10^(-5), 10^(-14)];
soft = c1hs[x - 8*10^(-5), 10^(-14)];
pillarsdist = 10*10^(-6);
rpillars = 1 10^(-6);
lspillars = 
  N[Table[(x - i)^2 + (y - j)^2 - rpillars^2, {i, -pillarsdist, 
     pillarsdist, pillarsdist}, {j, -pillarsdist, pillarsdist, 
     pillarsdist}]];
pillarscenters = 
  Table[{i, j}, {i, -pillarsdist, pillarsdist, 
    pillarsdist}, {j, -pillarsdist, pillarsdist, pillarsdist}];
pillars2cell = 
  Table[pillarscenters[[i, j]] - cellcenter, {i, 
    Dimensions[pillarscenters][[1]]}, {j, 
    Dimensions[pillarscenters][[2]]}];
Normpillars2cell = 
  Table[Norm[pillars2cell[[i, j]]], {i, 
    Dimensions[pillarscenters][[1]]}, {j, 
    Dimensions[pillarscenters][[2]]}];
pillarsangles = 
  Table[ArcTan[pillars2cell[[i, j]][[1]], 
    pillars2cell[[i, j]][[2]]], {i, 
    Dimensions[pillarscenters][[1]]}, {j, 
    Dimensions[pillarscenters][[2]]}];
pillars = 
  Table[c1hs[-lspillars[[i, j]], 1 10^(-30)], {i, 
    Dimensions[lspillars][[1]]}, {j, Dimensions[lspillars][[2]]}];

alpha = pillarsangles;
beta = Pi/10;
nppodx = Table[
   Cos[alpha[[i, j]]], {i, Dimensions[alpha][[1]]}, {j, 
    Dimensions[alpha][[2]]}];
nppody = Table[
   Sin[alpha[[i, j]]], {i, Dimensions[alpha][[1]]}, {j, 
    Dimensions[alpha][[2]]}];
nppod = Table[{nppodx[[i, j]], nppody[[i, j]]}, {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
xd = Table[
   x*nppodx[[i, j]] + y*nppody[[i, j]], {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
lscone = N[
   Table[Sqrt[(x - xd[[i, j]]*nppodx[[i, j]])^2 + (y - 
          xd[[i, j]]*nppody[[i, j]])^2] - Tan[beta]*xd[[i, j]], {i, 
     Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]];
ppods = Table[
   c1hs[-lscone[[i, j]], 10^(-14)], {i, Dimensions[alpha][[1]]}, {j, 
    Dimensions[alpha][[2]]}];

xcontour = 
  Table[rcell Cos[alpha[[i, j]]], {i, Dimensions[alpha][[1]]}, {j, 
    Dimensions[alpha][[2]]}];
ycontour = 
  Table[rcell Sin[alpha[[i, j]]], {i, Dimensions[alpha][[1]]}, {j, 
    Dimensions[alpha][[2]]}];
pcontour = 
  Table[{xcontour[[i, j]], ycontour[[i, j]]}, {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
phi = ArcTan[x, y];
irx = Cos[phi];
iry = Sin[phi];
eax = Table[
   If[pcontour[[i, j]][[1]] != 0, 
    Norm[(pillarscenters[[i, j]][[1]] - pcontour[[i, j]][[1]])/rcell],
     0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
eay = Table[
   If[pcontour[[i, j]][[2]] != 0, 
    Norm[(pillarscenters[[i, j]][[2]] - pcontour[[i, j]][[2]])/rcell],
     0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
alphaf = 0.04;
alphar = 0.01;
T = 60;
ea = Sin[2*Pi*t/T];
dea = D[ea, t];
eaf = Table[
   If[pcontour[[i, j]][[1]] + 
      uif[t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100 , 0], {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
ear = Table[{eax[[i, j]] ea c1hs[ea, 10^(-6)], 
    eay[[i, j]] ea c1hs[ea, 10^(-6)]}, {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];

mesh = ToElementMesh[cell, MaxCellMeasure -> {"Length" -> 1 10^(-6)}, 
  "BoundaryMeshGenerator" -> "Continuation"]

Ecell = 10^4;
EY1 = Ecell(**rigid+(Ecell/10)*soft*);
EY2 = Ecell/7;
nu12 = 0.3;
nu21 = 0.3;
G12 = EY1/(2*(1 + nu12));
ro = 1000;

lzonef = 2*10^(-6);
lzoner = -2*10^(-6);
zonef = c1hs[(x nppodx[[1, 1]] + y nppody[[1, 1]]) - lzonef, 10^(-12)];
zoner = c1hs[-(x nppodx[[1, 1]] + y nppody[[1, 1]]) + lzoner, 
   10^(-12)];
muv = 10^(16);
mustab = 10^(3);
mupil = 10^(9);
fadhf = muv*zonef*
   c1hs[-dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]};
fadhr = muv*zoner*
   c1hs[dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]};
fstab = mustab {D[u[t, x, y], t], D[v[t, x, y], t]};
fpillars = mupil {D[u[t, x, y], t], D[v[t, x, y], t]}*soft;
fadh = fadhf + fadhr + fstab;

mechprop = {Y -> EY1, \[Nu] -> nu12, rho -> ro};
\[Lambda]dp = Y \[Nu]/(1 + \[Nu])/(1 - 2 \[Nu]);
\[Lambda] = Y \[Nu]/(1 - \[Nu]^2);
\[Mu] = Y/2/(1 + \[Nu]);
\[CurlyEpsilon] = {{ix.Inactive[Grad][u[t, x, y], {x, y}], 
    1/2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + 
       ix.Inactive[Grad][v[t, x, y], {x, y}])}, {1/
     2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + 
       ix.Inactive[Grad][v[t, x, y], {x, y}]), 
    iy.Inactive[Grad][v[t, x, y], {x, y}]}};
\[CurlyEpsilon]0 = -ppods[[1, 1]] (eaf[[1, 
     1]]) ((x nppodx[[1, 1]] + y nppody[[1, 1]])/
      10^(-6))^2 {{nppodx[[1, 1]] nppodx[[1, 1]], 
     nppodx[[1, 1]] nppody[[1, 1]]}, {nppodx[[1, 1]] nppody[[1, 1]], 
     nppody[[1, 1]] nppody[[1, 1]]}};

\[Sigma] = \[Lambda] Tr[\[CurlyEpsilon]] Id + 
   2 \[Mu] (\[CurlyEpsilon]);
\[Sigma]0 = \[Lambda] Tr[\[CurlyEpsilon]0] Id + 
   2 \[Mu] (\[CurlyEpsilon]0);
divsig = {Inactive[Div][(\[Sigma].ix), {x, y}], 
    Inactive[Div][(\[Sigma].iy), {x, y}]} /. mechprop;
divsig0 = {Inactive[Div][(\[Sigma]0.ix), {x, y}], 
    Inactive[Div][(\[Sigma]0.iy), {x, y}]} /. mechprop;
inertia = 
  rho {D[u[t, x, y], {t, 2}], D[v[t, x, y], {t, 2}]} /. mechprop;

ci = {u[0, x, y] == 0., v[0, x, y] == 0.};
cid = {Derivative[1, 0, 0][u][0, x, y] == 0, 
   Derivative[1, 0, 0][v][0, x, y] == 0};

Button["Stop", stop = True]
stop = False;
currentTime = "Initialization";
interval = T;
SetOptions[EvaluationNotebook[], 
  WindowStatusArea -> Dynamic["t = " <> ToString[CForm[currentTime]]]];
{uif, vif} = 
  NDSolveValue[{Activate[divsig - divsig0 == inertia + fadh], ci, cid,
     WhenEvent[stop, end = t; "StopIntegration"]}, {u, 
    v}, {x, y} \[Element] mesh, {t, 0, interval}, Method -> {
     "PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5},
          "IntegrationOrder" -> 4}}}, 
   EvaluationMonitor :> (currentTime = t)(* ,AccuracyGoal\[Rule]5*), 
   MaxStepSize -> 0.5];

Hello

Thank you very much for the suggestion. However, when I run the code, it computes correctly up to tf = 10 but it does not stop and cannot post process the results.

Do you have any idea why?

Thanks rachele