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:
- Correct working.
- Divergence, at the first and rarely on other steps.
- 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.
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!
