Message Boards Message Boards

GROUPS:

Solve the following delay differential equations?

Posted 1 month ago
298 Views
|
16 Replies
|
0 Total Likes
|

Dear all,

I am solving a PDE such as :

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

Now, the variable 'divsig0' depends on a previous defined variable, 'eaf', which in turn depends itself on the solution of the PDE through 'uif' like this

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

It seems I should used a Delay Differential Equation according to the documentation, but I have trouble to specify the initial hystory function...could anyone help me?

I can send the whole code if necessary...

Best RA

16 Replies

There are no definitions for divsig - divsig0 == inertia + fadh, ci, cid

Posted 1 month ago

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

Here you can use the iteration method. I checked that the solution quickly converged. An example with three iterations

uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 3; tf = 10;

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

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

fig1

Posted 1 month ago

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

Posted 1 month ago

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

Check $Version and correct ArcTan::indet: Indeterminate expression ArcTan[0,0] encountered.

Posted 1 month ago

This normally shouldn't cause any problem for the NDSolve command, it is just a warning and it is additionally not used in the NDSolve

Posted 1 month ago

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

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"}]

fig2

Posted 1 month ago

Hello

I tried again, but still not working even by simplifying as you said the code. I get this error

NDSolveValue::femcnsd: The PDE coefficient 7.14286*10^15 <<3>>+3.84615*10^15 <<1>>^2 <<1>> (<<1>>)-1.55408*10^16 (-0.707107 x-0.707107 y) If[-(1/(200000 Sqrt[2]))+uif[0][3][t,-(1/(200000 Sqrt[2])),-(1/(200000 Sqrt[2]))]>-(1/100000),t/100,0] (-1. If[100000000000000 (Times[<<2>>]+Times[<<2>>])>=1,1,0]-1. (0.5 +1.*10^14 (0.32492 Plus[<<2>>]-1. Power[<<2>>]) (0.75 -2.5*10^27 Power[<<2>>])) If[100000000000000 (Times[<<2>>]+Times[<<2>>])>-1,1,0] If[100000000000000 (Times[<<2>>]+Times[<<2>>])<1,1,0]) does not evaluate to a numeric scalar at the coordinate {-2.34045*10^-6,-4.24011*10^-6}; it evaluated to 0. +7.23139*10^10 If[-3.53553*10^-6+uif[0.][3.][0.,-3.53553*10^-6,-3.53553*10^-6]>-0.00001,0. 0.01,0.] instead.

Set::shape: Lists {uif[i],vif[i]} and NDSolveValue[{<<1>>},{u,v},{x,y}\[Element]ElementMesh[{{5.*10^-6,1.046*10^-21},{4.95417*10^-6,6.75457*10^-7},{4.80783*10^-6,1.37285*10^-6},<<45>>,{-4.32484*10^-6,1.63214*10^-7},{-1.63214*10^-7,-4.32484*10^-6},<<1155>>},{<<1>>},{<<1>>},{<<1>>}],<<1>>,Method->{PDEDiscretization->{MethodOfLines,SpatialDiscretization->{FiniteElement,Rule[<<2>>],Rule[<<2>>]}}},MaxStepSize->0.5] are not the same shape.

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:
Posted 1 month ago

Actually it works like this, but the problem is that the line

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

is wrong because one needs to write

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

since with the iteration method, uif depends now on kf too. Am I right? Then withis change, it does not work anymore...

No, you're not right, we made this replacement in (divsig0 /. uif -> uif[i - 1]) in line Do[…]. If you want to make this replacement earlier, then make it so that there is no error due to the coincidence of the names of the iterator

eaf[kk_] := 
  Table[If[pcontour[[i, j]][[1]] + 
      uif[kk][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[
   kk_] := -ppods[[1, 
     1]] (eaf[kk][[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[kk_] := \[Lambda] Tr[\[CurlyEpsilon]0[kk]] Id + 
   2 \[Mu] (\[CurlyEpsilon]0[kk]);
divsig = {Inactive[Div][(\[Sigma].ix), {x, y}], 
    Inactive[Div][(\[Sigma].iy), {x, y}]} /. mechprop;
divsig0[kk_] := {Inactive[Div][(\[Sigma]0[kk].ix), {x, y}], 
    Inactive[Div][(\[Sigma]0[kk].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 = 4; tf = 15;

Do[{uif[kk], vif[kk]} = 
   NDSolveValue[{Activate[divsig - divsig0[kk - 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];, {kk, 1, kf}]
Table[ContourPlot[uif[i][tf, x, y], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, Contours -> 20, 
  ColorFunction -> "TemperatureMap", PlotLabel -> Row[{"i=", i}]], {i,
   1, kf}]

fig3

Attachments:
Posted 1 month ago

Actually it runs like this, but I am not sure that this war 'uif' is correctly taken into accoutn then. Normally, eaf[kk] should be equal to 0, and therefore the displacement, once

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

is not satisfied anymore. But when you plot

Plot[{uif[kf][t, pcontour[[1, 1]][[1]], pcontour[[1, 1]][[2]]], 
  vif[kf][t, pcontour[[1, 1]][[1]], pcontour[[1, 1]][[2]]]}, {t, 0, 
  tf}, PlotRange -> All]

this does not occur. Then I am wondering if the 'If' in the definition of 'eaf[kk] is correctly taken into account?

best rachele

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

Posted 1 month ago

'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...

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract