1
$\begingroup$

Is there any way in AceFem to use newton iterative scheme without using multiplier λ or time step t? For example, if one is trying to calculate solution to some nonlinear equation no multiplier is needed. In the problem that I am facing, I am trying to use AceFem to calculate geometry of some surface, assembled from surface patches. Every surface patch presents one finite element which has three degrees of freedom in every node. This degrees of freedom are twist vectors and with their variation, shapes of patches can be changed. To determine values of twist vectors I defined energy potential in finite element which has to be minimized. This potential is assembled from mean curvature in combination with Gaussian curvature (2H^2-2K). Node locations and normal vectors in nodes are always fixed therefore no boundary conditions are defined in analysis (also twist vectors cannot be defined as boundary conditions since they are all unknowns in every node).

Here is my code in AceFem:

<< AceFEM`;

(* Geometry - Nodes *)
XIO = {{-10, -10, 10 Sqrt[2]}, {10, -10, 10 Sqrt[2]}, {10, 10, 10 Sqrt[2]}, {-10, 10, 10 Sqrt[2]}}

(* Newton iteration parameters *)
I0max = 15;
tol = 10^-8;

SMTInputData["Threads" -> 1];

(* Finite element *)
elem = "DKQ-5";
a1 = {"tol"};
a2 = {tol};
SMTAddDomain["A", elem, Array[a1[[#]] -> a2[[#]] &, Length[a1]]];

(* Mesh import *)
SMTAddNode[XIO];
SMTAddElement["A", {1, 2, 3, 4}];

SMTAnalysis[];

(* Normal vectors *)
Normals = 
  Array[SMTNodeData[#,"X"]/(Sqrt[SMTNodeData[#, "X"].SMTNodeData[#, "X"]]) &,SMTIData["NoNodes"]];
Array[SMTNodeData[#, "Data", Normals[[#]]] &, SMTIData["NoNodes"]];
SMTUpdatePostData[];

(* Analysis *) 
SMTNextStep[];

While[SMTConvergence[tol, I0max],
  SMTNewtonIteration[];
  SMTStatusReport[];
  TwistNode1 = SMTNodeData[1, "at"];
  Print[TwistNode1];
  Energy = SMTTask["Energy"];
  Print["Energy = ", Energy];
  ];

The problem is that analysis is converging and no optimal solution is found.

Here is also element code:

<< AceGen`;
elem = "DKQ-5";

SMSInitialize[elem, "Language" -> "C", "Mode" -> "Optimal","VectorLength" ->3000];
SMSTemplate["CDriver","SMSTopology" -> "S1","SMSDOFGlobal" -> 3,"SMSNoNodeData" -> 3,"SMSDefaultIntegrationCode" -> {19 + 5, 19 + 5},"SMSCharSwitch" -> {"Energy", "Surface area"},"SMSSymmetricTangent" -> False,"SMSDomainDataNames" -> {"tol"},"SMSPostIterationCall" -> True,"SMSDefaultData" -> {10^-8}];

InputData[] := Block[{},
   XIO⊨Table[SMSReal[nd$$[i, "X", j]], {i, SMSNoNodes},{j, 3}];
   T0x⊨Array[SMSReal[nd$$[#, "at", 1]] &,SMSNoNodes];
   T0y⊨Array[SMSReal[nd$$[#, "at", 2]] &,SMSNoNodes];
   T0z⊨Array[SMSReal[nd$$[#, "at", 3]] &,SMSNoNodes];

   T0⊨{{T0x[[1]],T0y[[1]],T0z[[1]]},{T0x[[2]],T0y[[2]],T0z[[2]]}, {T0x[[3]],T0y[[3]],T0z[[3]]},{T0x[[4]],T0y[[4]],T0z[[4]]}};

   A3i⊨Array[SMSReal[nd$$[#1, "Data", #2]] &, {SMSNoNodes, 3}];

   {E0}⊨SMSReal[Table[es$$["Data", i], {i, Length[SMSDomainDataNames]}]];

   Do[A3i[[i]]⊨A3i[[i]]/SMSSqrt[A3i[[i]].A3i[[i]]], {i, 1, SMSNoNodes}];

   t0⊨-1;
   t1⊨1;

   A1v1⊨(Cross[Cross[A3i[[1]], (XIO[[2]] - XIO[[1]])], A3i[[1]]]);
   A2v1⊨(Cross[Cross[A3i[[1]], (XIO[[4]] - XIO[[1]])], A3i[[1]]]);

   A1v2⊨(Cross[Cross[A3i[[2]], (XIO[[2]] - XIO[[1]])], A3i[[2]]]);
   A2v2⊨(Cross[Cross[A3i[[2]], (XIO[[3]] - XIO[[2]])], A3i[[2]]]);

   A1v3⊨(Cross[Cross[A3i[[3]], (XIO[[3]] - XIO[[4]])], A3i[[3]]]);
   A2v3⊨(Cross[Cross[A3i[[3]], (XIO[[3]] - XIO[[2]])], A3i[[3]]]);

   A1v4⊨(Cross[Cross[A3i[[4]],(XIO[[3]]-XIO[[4]])],A3i[[4]]]);
   A2v4⊨(Cross[Cross[A3i[[4]], (XIO[[4]] - XIO[[1]])], A3i[[4]]]);

   AR12A1v1⊨(6 (-A1v1.XIO[[1]] + A1v1.XIO[[2]]) A1v2.A1v2 + 3 A1v1.A1v2 (A1v2.XIO[[1]] - A1v2.XIO[[2]]))/((t0 - t1) ((A1v1.A1v2)^2 - 4 A1v1.A1v1 A1v2.A1v2));

   AR12A1v2⊨(3 A1v1.A1v2 (A1v1.XIO[[1]] - A1v1.XIO[[2]]) + 6 A1v1.A1v1 (-A1v2.XIO[[1]] + A1v2.XIO[[2]]))/((t0 - t1) ((A1v1.A1v2)^2 - 4 A1v1.A1v1 A1v2.A1v2));

   AR23A2v2⊨(6 (-A2v2.XIO[[2]] + A2v2.XIO[[3]]) A2v3.A2v3 + 3 A2v2.A2v3 (A2v3.XIO[[2]] - A2v3.XIO[[3]]))/((t0 - t1) ((A2v2.A2v3)^2 - 4 A2v2.A2v2 A2v3.A2v3));

   AR23A2v3⊨(3 A2v2.A2v3 (A2v2.XIO[[2]] - A2v2.XIO[[3]]) + 6 A2v2.A2v2 (-A2v3.XIO[[2]] + A2v3.XIO[[3]]))/((t0 - t1) ((A2v2.A2v3)^2 - 4 A2v2.A2v2 A2v3.A2v3));

   AR43A1v4⊨(6 (-A1v4.XIO[[4]] + A1v4.XIO[[3]]) A1v3.A1v3 + 3 A1v4.A1v3 (A1v3.XIO[[4]] - A1v3.XIO[[3]]))/((t0 - t1) ((A1v4.A1v3)^2 - 4 A1v4.A1v4 A1v3.A1v3)); 

   AR43A1v3⊨(3 A1v4.A1v3 (A1v4.XIO[[4]] - A1v4.XIO[[3]]) + 6 A1v4.A1v4 (-A1v3.XIO[[4]] + A1v3.XIO[[3]]))/((t0 - t1) ((A1v4.A1v3)^2 - 4 A1v4.A1v4 A1v3.A1v3));

   AR14A2v1⊨(6 (-A2v1.XIO[[1]] + A2v1.XIO[[4]]) A2v4.A2v4 + 3 A2v1.A2v4 (A2v4.XIO[[1]] - A2v4.XIO[[4]]))/((t0 - t1) ((A2v1.A2v4)^2 - 4 A2v1.A2v1 A2v4.A2v4)); 

   AR14A2v4⊨(3 A2v1.A2v4 (A2v1.XIO[[1]] - A2v1.XIO[[4]]) + 6 A2v1.A2v1 (-A2v4.XIO[[1]] + A2v4.XIO[[4]]))/((t0 - t1) ((A2v1.A2v4)^2 - 4 A2v1.A2v1 A2v4.A2v4));


   A10⊨{A1v1*AR12A1v1, A1v2*AR12A1v2, A1v3*AR43A1v3, A1v4*AR43A1v4};
   A20⊨{A2v1*AR14A2v1, A2v2*AR23A2v2, A2v3*AR23A2v3, A2v4*AR14A2v4};

   ];

Discretization[] := Block[{},
   N1⊨1/16*{(u - 1)^2*(v - 1)^2*(u + 2)*(v + 2), (-u - 1)^2*(v - 1)^2*(-u + 2)*(v + 2), (-u - 1)^2*(-v - 1)^2*(-u + 2)*(-v + 2), (u - 1)^2*(-v - 1)^2*(u + 2)*(-v + 2)};
   N2⊨1/16*{(u - 1)^2*(v - 1)^2*(u + 1)*(v + 2), -(-u - 1)^2*(v - 1)^2*(-u + 1)*(v + 2), -(-u - 1)^2*(-v - 1)^2*(-u + 1)*(-v + 2), (u - 1)^2*(-v - 1)^2*(u + 1)*(-v + 2)};
   N3⊨1/16*{(u - 1)^2*(v - 1)^2*(u + 2)*(v + 1), (-u - 1)^2*(v - 1)^2*(-u + 2)*(v + 1), -(-u - 1)^2*(-v - 1)^2*(-u + 2)*(-v + 1), -(u - 1)^2*(-v - 1)^2*(u + 2)*(-v + 1)};
   N4⊨1/16*{(-1 + v)^2*(1 + v)*(-1 + u)^2*(1 + u), (-1 + v)*(1 + v)^2*(-1 + u)^2*(1 + u), (-1 + v)^2*(1 + v)*(-1 + u)*(1 + u)^2, (-1 + v)*(1 + v)^2*(-1 + u)*(1 + u)^2};

   X⊨N1.XIO + N2.A10 + N3.A20 + N4.T0;

   A1⊨SMSD[X, u];
   A2⊨SMSD[X, v];

   A11⊨A1.A1;
   A22⊨A2.A2;
   A12⊨A1.A2;

   A⊨A11*A22 - A12^2;

   A3⊨1/SMSSqrt[A]*Cross[A1, A2];

   Anad11⊨A22/A;
   Anad22⊨A11/A;
   Anad12⊨-A12/A;

   B11⊨A3.SMSD[A1, u];
   B22⊨A3.SMSD[A2, v];
   B12⊨(A3.SMSD[A2, u] + A3.SMSD[A1, v])/2;

   B1nad1⊨Anad11*B11 + Anad12*B12;
   B2nad2⊨Anad12*B12 + Anad22*B22;
   B2nad1⊨Anad11*B12 + Anad12*B22;
   B1nad2⊨Anad12*B12 + Anad22*B12;

   H⊨1/2*(B1nad1 + B2nad2);

   K⊨B1nad1*B2nad2 - B2nad1*B1nad2;

   ];

constant = {};
SMSStandardModule["Tangent and residual"];

{nA, nB}⊢SMSInteger[{es$$["id", "NoIntPointsA"], es$$["id","NoIntPointsB"]}];

InputData[];

SMSDo[iu, 1, nA];
    u⊢SMSReal[es$$["IntPoints", 1, iu]];
    SMSDo[iv, 1, nB];
       v⊢SMSReal[es$$["IntPoints", 2, (iv - 1)*nA + 1]];
       wGauss⊢SMSReal[es$$["IntPoints", 4, iu + (iv - 1) nA]];

       Discretization[];

       W⊨4*H^2 - 2*K;

       p0={T0x, T0y, T0z} // Transpose // Flatten;
       SMSDo[i, 1, Length[p0]];
          Rg⊨wGauss*SMSD[SMSSqrt[A]*W, p0, i,"Constant" -> constant];
          SMSExport[Rg, p$$[i], "AddIn" -> True];

          SMSDo[j, 1, Length[p0]];
             Kg⊨SMSD[Rg, p0, j];
             SMSExport[Kg, s$$[i, j], "AddIn" -> True];
          SMSEndDo[];
       SMSEndDo[];

    SMSEndDo[];

SMSEndDo[];

constant = {};
SMSStandardModule["Tasks"];
task⊨SMSInteger[Task$$];

{nA, nB}⊢SMSInteger[{es$$["id", "NoIntPointsA"], es$$["id","NoIntPointsB"]}];

InputData[];

SMSIf[task < 0
  , SMSSwitch[task
  , -1,
   SMSExport[{1, 0, 0, 0, 1}, TasksData$$];
  , -2,
   SMSExport[{1, 0, 0, 0, 1}, TasksData$$];
  ];

  ,
  SMSDo[iu, 1, nA];
     u⊢SMSReal[es$$["IntPoints", 1, iu]];
     SMSDo[iv, 1, nB];
        v⊢SMSReal[es$$["IntPoints", 2, (iv - 1)*nA + 1]];

        wGauss⊢SMSReal[es$$["IntPoints", 4, iu + (iv - 1) nA]];

        Discretization[];

        W⊨4*H^2 - 2*K;

        SMSSwitch[task
        , 1,
        SMSExport[wGauss*SMSSqrt[A]*W, RealOutput$$[1],"AddIn" -> True]
        , 2,
        SMSExport[wGauss*SMSSqrt[A], RealOutput$$[1], "AddIn" -> True]
        ];

    SMSEndDo[];

  SMSEndDo[];

  ];

SMSWrite[];

Every help is welcome. Thank you. Tomo

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Pinti Sep 29 '19 at 17:32
  • $\begingroup$ It is much easier to answer your question if you provide a minimal self-contained code example of your problem. In general you can use SMTNewtonIteration without changing the BC multiplier (lambda). You could, for example, in every iteration modify SMTDomainData or SMTElementData and that way increment the solution. $\endgroup$ – Pinti Sep 29 '19 at 17:39
  • $\begingroup$ Thank you for your recommendations regarding the usage of this forum. I will check how to navigate trough it. In my problem I do not have to modify SMTDomainData or SMTElementData. I only have to calculate degrees of freedom at the nodes with Newton method where no external forces are applied on the geometry. $\endgroup$ – TVeldin Sep 30 '19 at 14:34
  • $\begingroup$ Could you please make your example self-contained and include all code one would need to rerun it and investigate it? Currently code for element routine ("DKQ-5") is missing. If element routine is very long and complicated, consider making it shorter, while still demonstrating its behavior. It is very difficult, for anyone, to figure out what is the problem, without a working example. $\endgroup$ – Pinti Oct 1 '19 at 13:42
  • $\begingroup$ Thank you for adding element routine code. Unfortunately code cannot be run/compiled, because special AceGen syntax is missing (e.g. \[DoubleRightTee]). Please edit the code so that anyone can copy it and run on their own computer. You can also use these tools for nice formatting of special characters. $\endgroup$ – Pinti Oct 3 '19 at 11:32

Browse other questions tagged or ask your own question.