6
$\begingroup$

I'm back with an unusual issue. I`d like to generate a finite element for AceFEM, that calls an external (AceGen generated) subroutine, for the material. Upfront, as the code is a little broad, I provide a Link with the Mathematica notebook on the bottom.

First I'm using AceGen to generate this subroutine:

<< AceGen`;
SMSInitialize["mate_stvenant", "Language" -> "C"];
SMSModule["Material", 
  Real[d$$[2], F$$[9], hin$$[1], hnin$$[1], Pvec$$[9], dPdF$$[9, 9], 
   hout$$[1], info$$[1]]];

{Eg, νg} ⊢ Table[SMSReal[d$$[i]], {i, 2}];
{λg, μg} ⊨ SMSHookeToLame[Eg, νg];

Fvec ⊢ Table[SMSReal[F$$[i]], {i, 9}];
SMSFreeze[\[DoubleStruckCapitalF], Partition[Fvec, 3]];

\[DoubleStruckCapitalE] ⊨ 1/2 (\[DoubleStruckCapitalF]\[Transpose].\[DoubleStruckCapitalF] - IdentityMatrix[3]);
Π ⊨ λg/
    2 Tr[\[DoubleStruckCapitalE]]^2 + μg Tr[\[DoubleStruckCapitalE].\[DoubleStruckCapitalE]];
\[DoubleStruckCapitalP] ⊨ SMSD[Π, \[DoubleStruckCapitalF]];

Pvec ⊨ Flatten[\[DoubleStruckCapitalP]];
dPdF ⊨ SMSD[Pvec, Fvec, "Method" -> "Forward"];

(*export quantities*)
SMSExport[Pvec, Pvec$$];
SMSExport[dPdF, dPdF$$];
SMSExport[{0}, hout$$];
SMSExport[{1}, info$$];
SMSWrite[];

Then, I use AceGen to generate the FE-code, by using SMSCall to make use of the just generated Material-Subroutine, which is added into the element source code by the "Splice" option of SMSWrite. The AceGen input reads as follows:

<< AceGen`;
SMSInitialize["mate_elmt", "Environment" -> "AceFEM", "Mode" -> "Prototype"];
SMSTemplate[
  "SMSTopology" -> "O2"
  , "SMSDOFGlobal" -> {3, 3, 3, 3, 3, 3, 3, 3, 3, 3}
  , "SMSNoElementData" -> 1
  , "SMSDefaultIntegrationCode" -> {43}
  , "SMSNoTimeStorage" -> es$$["id", "NoIntPoints"]
  , "SMSDomainDataNames" -> {"E", "\[Nu]"}
  , "SMSDefaultData" -> {21000, 0.3}
  , "SMSSymmetricTangent" -> False
  ];

SMSStandardModule["Tangent and residual"];

matd\[DoubleStruckCapitalB] ⊢ Table[SMSReal[es$$["Data", i]], {i, 2}];

XI ⊢ Table[SMSReal[nd$$[i, "X", j]], {i, SMSNoNodes}, {j, SMSNoDimensions}];
UI ⊢ Table[SMSReal[nd$$[i, "at", j]], {i, SMSNoNodes}, {j, SMSDOFGlobal[[i]]}];
DOFVector ⊨ Flatten[UI];

LocalTolerance ⊢ SMSReal[rdata$$["SubIterationTolerance"]];
GlobalIterationNumber ⊢ SMSInteger[idata$$["Iteration"]];
PrintOn ⊢ SMSInteger[ed$$["Data", 1]];


SMSDo[Ig, 1, SMSInteger[es$$["id", "NoIntPoints"]]];
Ξ = {ξ, η, ζ} ⊢ Table[SMSReal[es$$["IntPoints", i, Ig]], {i, 3}];
ω ⊢ SMSReal[es$$["IntPoints", 4, Ig]];
SHP ⊨ {ξ*(-1 + 2*ξ), η*(-1 + 2*η), ζ*(-1 + 2*ζ), (-1 + ζ + η + ξ)*(-1 + 2*ζ + 2*η + 2*ξ), 4*η*ξ, 4*ζ*η, 4*ζ*ξ, -4*ξ*(-1 + ζ + η + ξ), -4* η*(-1 + ζ + η + ξ), -4*ζ*(-1 + ζ + η + ξ)};

Igd = SMSInteger[(Ig - 1) LengthOfLocalField];
\[DoubleStruckH]gnIO ⊢ Table[SMSReal[ed$$["hp", Igd + i]], {i, 1}];
\[DoubleStruckH]gIO ⊢ Table[SMSReal[ed$$["ht", Igd + i]], {i, 1}];

SMSFreeze[X, SHP.XI];
U ⊨ SHP.UI;
Jm ⊨ SMSD[X, Ξ];
detJ ⊨ SMSDet[Jm];
\[DoubleStruckCapitalH] ⊨ SMSD[U, X, "Dependency" -> {Ξ, X, SMSInverse[Jm]}];
SMSFreeze[\[DoubleStruckCapitalF], IdentityMatrix[3] + \[DoubleStruckCapitalH]];
Fvec ⊨ Flatten[\[DoubleStruckCapitalF]];

(*Call Material Subroutine*)
MaterialSubordinate = SMSCall[
   "Material",
   matd\[DoubleStruckCapitalB],
   Fvec,
   \[DoubleStruckH]gIO,
   \[DoubleStruckH]gnIO,
   Real[Pvec$$[9]],
   Real[dPdF$$[9, 9]],
   Real[hout$$[1]],
   Real[io$$[1]]
   ];

h ⊢ SMSReal[Table[hout$$[i], {i, 1}], "Subordinate" -> MaterialSubordinate];
SMSExport[h, Table[ed$$["ht", Igd + i], {i, 1}]];

matio ⊢ SMSReal[Table[io$$[i], {i, 1}], "Subordinate" -> MaterialSubordinate];

dPdF ⊢ SMSReal[Table[dPdF$$[i, j], {i, 9}, {j, 9}], "Subordinate" -> MaterialSubordinate];
Pvec ⊢ SMSReal[Table[Pvec$$[i], {i, 9}], "Subordinate" -> MaterialSubordinate, "Dependency" -> {Fvec, dPdF}];
\[DoubleStruckCapitalP] ⊨ Partition[Pvec, 3];

W ⊨ Tr[\[DoubleStruckCapitalP].\[DoubleStruckCapitalF]\[Transpose]];
SMSDo[m, 1, Length[DOFVector]];
δΠ ⊨ SMSD[W, DOFVector, m, "Constant" -> {\[DoubleStruckCapitalP]}];
SMSExport[detJ ω δΠ, p$$[m], "AddIn" -> True];
   SMSDo[n, 1, Length[DOFVector]];
   ΔδΠ = 
  SMSD[δΠ, DOFVector, n];
   SMSExport[detJ ω ΔδΠ, 
  s$$[m, n], "AddIn" -> True];
   SMSEndDo[];
SMSEndDo[];

SMSEndDo[];
SMSWrite["Splice" -> {"mate_stvenant.c"}];

So far, this procedure works, but not reliable on my machine. The test example in AceFEM, a simple displacement driven tension test, is given by:

<< AceFEM`;
SMTInputData["Threads" -> 1];
SMTAddDomain["TestCube", "mate_elmt", {"E" -> 210000, "ν" -> 0.3}];
SMTAddMesh[
  Hexahedron[{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 
     1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}}], "TestCube", 
  "O2", {1, 1, 2} 2];
SMTAddEssentialBoundary[{"X" == 0 &, 1 -> 0}];(*Back*)
SMTAddEssentialBoundary[{"Z" == 0 &, 3 -> 0}];(*Bottom*)
SMTAddEssentialBoundary[{"Y" == 0 &, 2 -> 0}];(*Side*)
SMTAddEssentialBoundary[{"X" == 1 &, 1 -> 1}];(*Front*)
SMTAnalysis["Output" -> NotebookDirectory[] <> "out.dat"]; εσ = {};
Monitor[Do[
   SMTNextStep["λ" -> LoadingCurve[t], "t" -> t];
   While[SMTConvergence[], SMTNewtonIteration[]];
   ℬ = 
    Show[SMTShowMesh["DeformedMesh" -> True], 
     SMTShowMesh["DeformedMesh" -> False, "FillElements" -> False], 
     Lighting -> "Neutral"];
   , {t, 0, tmax, .5}], ℬ];
ℬ

For multiple executions of this problem, I get three different, randomly appearing responses:

  1. Correct working.
  2. Divergence, at the first and rarely on other steps.
  3. A complete abortion of AceFem with reference to a memory issue.

To me this looks like a memory related e.g. allocation problem. This is even more reasonable as I was not able to replicate this behaviour on "bigger" machine that I have excess to. Of course I checked my code, but it seems like all the IO fields, passed from the element to the subroutine and vice versa, are given in the correct dimensions. Following the links below, you'll find the element source code file generated on my machine, as well as the AceFem error messages with respect to the last mentioned case.

Element Source Code

Failure Log

Full Mathematica Notebook

As you will notice from the log, my local machine is a Mac. If required I can report more specific version information etc..

I'm very excited about your thought on this issue, Thanks in advance!

$\endgroup$
  • $\begingroup$ I tested you code on Windows and it doesn't work either. Is there a particular reason why are you using this subroutine approach, since your material law is very simple? $\endgroup$ – Pinti May 9 '18 at 8:57
  • $\begingroup$ It is indeed very simple, to make a good example here. I found this approach very convenient when working with more complex once, as the AceGen input in such cases is also very large. Doesn't it work on your machine at all or do you also see this varying behaviour ? $\endgroup$ – Sascha Maassen May 9 '18 at 9:25
  • $\begingroup$ Aha, so it is just a MWE, I thought so. On my machine I always get an error, but sometimes the message is different. It is weird! $\endgroup$ – Pinti May 9 '18 at 10:20
5
$\begingroup$

Bug fixed in AceFEM version 6.901


I have run your code, and it was running fine. I was using AceGen/FEM version 6.816 on Linux. I have then installed the newest version 6.824 and regenerated the code. There I was able to reproduce the error. Then I checked the differences in both "mate_elmt.c" Among minor things, there is additional definition inside code

    "es->id.WorkingVectorSize=5373;"

"WorkingVectorSize" is a new option of SMSWrite. By removing this string and saving the "c" code, your example was running fine after. This seems to be a version problem so I will report this to the Author.

For now you can try doing this manually by using following command that removes it automatically:

    Export[#,
        StringReplace[Import[#, "Text"], 
        "es->id.WorkingVectorSize=" ~~ DigitCharacter .. ~~ ";" -> ""], 
    "Text"]&@"mate_elmt.c";

If it doesn’t work in your case then you can install older version until issue is resolved.

$\endgroup$
  • $\begingroup$ +1, Great troubleshooting! $\endgroup$ – Pinti May 9 '18 at 14:08
  • $\begingroup$ Blaz, this works very well so far, also for more complex Materials. Also credit for Pintis quick response. Thank you guys very much! $\endgroup$ – Sascha Maassen May 10 '18 at 11:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.