I resend you the code as I updated it from your last suggestion
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]]\[NotEqual]0,Norm[(pillarscenters[[\
i,j]][[1]]-pcontour[[i,j]][[1]])/pcontour[[i,j]][[1]]],0],{i,\
Dimensions[alpha][[1]]},{j,Dimensions[alpha][[2]]}];
eay=Table[If[pcontour[[i,j]][[2]]\[NotEqual]0,Norm[(pillarscenters[[i,\
j]][[2]]-pcontour[[i,j]][[2]])/pcontour[[i,j]][[2]]],0],{i,Dimensions[\
alpha][[1]]},{j,Dimensions[alpha][[2]]}];*)
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[kf][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 ndirx+y ndiry)-lzonef,10^(-12)];
zoner=c1hs[-(x ndirx+y ndiry)+lzoner,10^(-12)];*)
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={Y1\[Rule]EY1,Y2\[Rule]EY2,\[Nu]12\[Rule]nu12,\[Nu]21\
\[Rule]nu21,mu12\[Rule]G12,rho\[Rule]ro};*)
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]]}};
(*\[CurlyEpsilon]0=-ppods[[1,1]]ea((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]]}};*)
(*S0=({{1/Y1,-\[Nu]12/EY1},
{-\[Nu]21/Y2,1/Y2}});
S1=({{\[Infinity],1/(2mu12)},
{1/(2mu12),\[Infinity]}});
IS0=Inverse[S0];
IS1=1/S1;
\[Sigma]=Sum[IS0[[a,b]] \
\[CurlyEpsilon][[b,b]]Tens[V[[;;,a]],V[[;;,a]]],{a,1,2},{b,1,2}]+
Sum[IS1[[a,b]] \
\[CurlyEpsilon][[a,b]]Tens[V[[;;,a]],V[[;;,b]]],{a,1,2},{b,1,2}];
\[Sigma][[2,1]]=\[Sigma][[1,2]];
\[Sigma]0=Sum[IS0[[a,b]] \
\[CurlyEpsilon]0[[b,b]]Tens[V[[;;,a]],V[[;;,a]]],{a,1,2},{b,1,2}]+
Sum[IS1[[a,b]] \
\[CurlyEpsilon]0[[a,b]]Tens[V[[;;,a]],V[[;;,b]]],{a,1,2},{b,1,2}];
\[Sigma]0[[2,1]]=\[Sigma]0[[1,2]];*)
\[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};
uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 1;
tf = 15;
currentTime = "Initialization";
SetOptions[EvaluationNotebook[],
WindowStatusArea -> Dynamic["t = " <> ToString[CForm[currentTime]]]];
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" -> 1,
"MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}},
EvaluationMonitor :> (currentTime = t), MaxStepSize -> 0.5];, {i,
1, kf}]
