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];
