# Solve the following delay differential equations?

Posted 6 months ago
899 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
Sort By:
Posted 6 months ago
 There are no definitions for divsig - divsig0 == inertia + fadh, ci, cid
Posted 6 months ago
 here the whole code Clear["Global*"] Needs["NDSolveFEM"] 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]; 
Posted 6 months ago
 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}] 
Posted 6 months ago
 HelloThank 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 6 months 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
Posted 6 months ago
 Check $Version and correct ArcTan::indet: Indeterminate expression ArcTan[0,0] encountered. Answer Posted 6 months 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 Answer Posted 6 months ago  I resend you the code as I updated it from your last suggestion Clear["Global*"] Needs["NDSolveFEM"] 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}]  Answer Posted 6 months ago  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"}]  Answer Posted 6 months ago  HelloI 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.  Answer Posted 6 months ago  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 6 months 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...
Posted 6 months ago
 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}]  Attachments:
Posted 6 months 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
Posted 6 months ago
 What problem are you solving? What do all these If[] mean?
Posted 6 months 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 functionHowever, 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...
Community posts can be styled and formatted using the Markdown syntax.