Standard AceGen Procedure

Load AceGen package

This loads the AceGen package.

<<AceGen`

Description of AceGen Characteristic Steps

Step 1: Initialization

This initializes the AceGen session. FORTRAN is chosen as the final code language. See also SMSInitialize.

SMSInitialize["test", "Language" -> "Fortran"];

Step 2: Definition of Input and Output Parameters

This starts a new subroutine with the name "Test" and four real type parameters. The input parameters of the subroutine are u, x, and L, and parameter g is an output parameter of the subroutine. The input and output parameters of the subroutine are characterized by the double $ sign at the end of the name.

SMSModule["Test", Real[u$$[3], x$$, L$$, g$$[3]]];

Step 3: Definition of Numeric-Symbolic Interface Variables

The standard Mathematica assignment operator = has been replaced by the special AceGen operator ⊨. Operator ⊨ performs stochastic simultaneous optimization of expressions.

{x, L} ⊨ SMSReal[{x$$, L$$}]
ui ⊨ SMSReal[Array[u$$[#1]  & , 3]]

Step 4: Description of the Problem

Ni ⊨ {x/L,  1  -  x/L,  
        x/L*(1  -  x/L)}

u ⊨ Ni  .  ui
f ⊨ u^2
g ⊨ SMSD[f,  ui]

Step 5: Definition of Symbolic - Numeric Interface Variables

This assigns the results to the output parameters of the subroutine. See also SMSExport.

SMSExport[g,  g$$];

Step 6: Code Generation

During the session AceGen generates pseudo-code which is stored into the AceGen database. At the end of the session AceGen translates the code from pseudo-code to the required script or compiled program language and prints out the code to the output file. See also SMSWrite.

SMSWrite[];

This displays the contents of the generated file.

!!test.f

Generation of MathLink code

Here the MathLink version of the source code is generated. The generated code is automatically enhanced by an additional modules necessary for the proper MathLink connection.

<< AceGen`;
SMSInitialize["test", "Environment" -> "MathLink"];
SMSModule["Test", 
        Real[u$$[3], x$$, L$$, 
          g$$[3]], 
"Input" -> {u$$, x$$, 
            L$$}, "Output" -> g$$];
{x, L} ⊨ {SMSReal[x$$], SMSReal[L$$]};
ui ⊨ SMSReal[Array[u$$[#1]  & , 3]];
Ni ⊨ {(x/L), 
          1 - (x/L), (x/L) (1 - (x/L))};
u ⊨ Ni . ui;
f ⊨ (u^2);
g ⊨ SMSD[f, ui];
SMSExport[g, g$$];
SMSWrite[];

!!test.c

Here the MathLink program Test.exe is build from the generated source code and installed so that functions defined in the source code can be called directly from Mathematica. (see also SMSInstallMathLink)

SMSInstallMathLink[]

Here the generated executable is used to calculate gradient for the numerical test example. (see also Verification of Automatically Generated Code ).

Test[{0. , 1. , 7. }, π // N, 10. ]

Generation of heat conduction finite element

Step 1: Initialization

This loads the AceGen code generator.

<< AceGen`;

This initializes the AceGen session.

\(\(SMSInitialize["heatconduction", "Environment" -> "AceFEM"];\)\)

This initializes three-dimensional, eight node, hexahedron element with one degree of freedom per node.

\(\(SMSTemplate["SMSTopology" -> "H1", "SMSDOFGlobal" -> 1, "SMSSymmetricTangent" -> False];\)\)

Step 2: Element subroutine for the evaluation of tangent matrix and residual

SMSStandardModule["Tangent and residual"];

Step 3: Interface to the input data of the element subroutine

Xi ⊨ 
        Array[ SMSReal[nd$$[#, "X", 1]] &, 8];
Yi ⊨ 
        Array[ SMSReal[nd$$[#, "X", 2]] &, 8];
Zi ⊨ 
        Array[ SMSReal[nd$$[#, "X", 3]] &, 8];
ϕi ⊨ 
        Array[SMSReal[nd$$[#, "at", 1]] &, 8];
SMSGroupDataNames = {"Conductivity parameter k0", "Conductivity \
parameter k1", 
"Conductivity parameter k2", "Heat \
source"};
{k0, k1, k2, Q} ⊨ 
        SMSReal[{es$$["Data", 1], es$$["Data", 2], 
            es$$["Data", 3], es$$["Data", 4]}];

Step 4: Numerical integration loop

\(\(SMSDo[IpIndex, 1, 
        SMSInteger[
          es$$["id", "NoIntPoints"]]];\)
\)
{ξ, η, ζ, wGauss} ⊢ 
        Array[SMSReal[es$$["IntPoints", #1, IpIndex]] &, 
          4];

Step 5: Definition of the trial functions

ϕ = (Subscript[N, i]) . (Subscript[ϕ, i]). 
The (Subscript[ϕ, i]) are unknown parameters of the variational problem.

Ni ⊨ 
        MapThread[ 
          1/8 (1 + ξ #1) (1 + η #2) (1 + ζ #3) &, 
Transpose[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 
                1, -1}, 
{-1, -1, 
                1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}]];
X ⊢ SMSFreeze[Ni . Xi];
Y ⊢ SMSFreeze[Ni . Yi];
Z ⊢ SMSFreeze[Ni . Zi];
Jm ⊨ 
        SMSD[{X, Y, Z}, {ξ, η, ζ}];
SMSDefineDerivative[{ξ, η, ζ}, {X, Y, Z}, 
        SMSInverse[Jm]];
ϕ ⊨ Ni . ϕi;

Step 5: Definition of the governing equations

k ⊨ 
        k0 + k1 ϕ + k2 (ϕ^2);
δϕ ⊨ 
        SMSD[ϕ, ϕi];
λ ⊨ 
        SMSReal[rdata$$["Multiplier"]];
Ψi ⊨ 
        Det[Jm] wGauss (k SMSD[δϕ, {X, Y, Z}] . 
                  SMSD[ϕ, {X, Y, 
                      Z}] - δϕ λ Q);

Element contribution to global residual vector Ψi is exported into the p$$ output parameter of the "Tangent and residual" subroutine.

SMSExport[SMSResidualSign  Ψi, 
        p$$, "AddIn" -> True];

Step 6: Definition of the Jacobian matrix

This evaluates the explicit form of the Jacobian (tangent) matrix and exports result into the s$$ output parameter of the user subroutine.

Kij ⊨ 
        SMSD[Ψi, ϕi];
SMSExport[Kij, s$$, "AddIn" -> True];

SMSEndDo[];

Step 7: Code Generation

SMSWrite[];

!! heatconduction.c

Minimization of free energy

Free energy of the heat conduction problem is formulated as
\(TraditionalForm\`Π = ∫\_Ω\((\(1\/2\) k\ Δϕ . Δϕ - ϕ\ Q)\) dΩ\)
where a ϕ indicates temperature, a k is the conductivity and a Q is the heat generation per unit volume and Ω is the domain of the problem. 
The domain of the example is a cube filled with water:

The problem is formulated using various approaches using FindMinimum function:
A.	Trial polynomial interpolation if temperature with fifth order Lagrange polynomials in x y and z direction
	M.G	Gradient method of optimization + Mathematica directly
	M.N	Newton method of optimization + Mathematica directly
	A.G	Gradient method of optimization + AceGen+MathLink
	A.N	Newton method of optimization + AceGen+MathLink

B.	Finite difference interpolation (large sparse Hessian)
	M.G	Gradient method of optimization + Mathematica directly
	M.N	Newton method of optimization + Mathematica directly
	A.G	Gradient method of optimization + AceGen+MathLink
	A.N	Newton method of optimization + AceGen+MathLink

C.AceFEM	Finite element method

A. Trial Lagrange polynomial interpolation

A trial function for temperature ϕ is constructed as a fifth order Lagrange polynomials in x y and z direction. The chosen trial function is constructed in a way that  satisfies boundary conditions.

<< AceGen`;
Clear[x, y, z, a];
kcond = 0.58; Q = 500;
order = 5;
nterm = (order - 1) (order - 1) (order)

Here the fifth order Lagrange polynomials are constructed in three dimensions.

toc = Table[{x, 0}, {x, -0.5, 0.5, 1/order}]; 
    xp = MapIndexed[ 
        InterpolatingPolynomial[ReplacePart[toc, 1, {#, 2}], x] &, 
        Range[2, order]];
yp = 
        MapIndexed[ 
          InterpolatingPolynomial[ReplacePart[toc, 1, {#, 2}], y] &, 
          Range[2, order]];
toc = 
        Table[{x, 0}, {x, 0. , 1. , 1/order}];
zp = 
        MapIndexed[ 
          InterpolatingPolynomial[ReplacePart[toc, 1, {#, 2}], z] &, 
          Range[2, order + 1]];
ϕi = Array[a, nterm];
poly = 
        Flatten[Outer[Times, xp, yp, zp]] // 
          Chop;
ϕ = poly . ϕi;

poly[[28]]
Plot3D[
        poly[[28]] /. z -> 0.5, {x, -0.5, 0.5}, {y, -0.5, 
          0.5}, PlotRange -> All];

Here the Gauss points and weights are calculated for ngp×ngp×ngp Gauss numerical integration of the free energy over the domain [-.0.5m,0.5m]×[-0.5m,0.5m]×[0,1m].

ngp = 6;
<< NumericalMath`GaussianQuadrature`;
g1 = GaussianQuadratureWeights[ngp, -0.5, 0.5];
g2 = GaussianQuadratureWeights[ngp, -0.5, 0.5];
g3 = GaussianQuadratureWeights[ngp, 0, 1];
gp = {g1[[#1[[1]], 1]], g2[[#1[[2]], 1]], 
              g3[[#1[[3]], 1]], 
              g1[[#1[[1]], 2]]*g2[[#1[[2]], 2]]*
                g3[[#1[[3]], 2]]} & /@ 
          Flatten[Array[{#3, #2, #1} &, {ngp, ngp, ngp}], 2];

The subsection Definitions has to be executed before the current subsection.

start = SessionTime[];
Δϕ = {D[ϕ, x], D[ϕ, y], 
          D[ϕ, z]};
Π = 
        1/2 kcond Δϕ . Δϕ - ϕ Q;
Πi = 
        Total[Map[(#[[4]] Π /. {x -> #[[1]], 
                    y -> #[[2]], z -> #[[3]]}) &, 
            gp]];
derivation = SessionTime[] - start

G. Gradient based optimization

start = SessionTime[]; ii = 0;
sol = 
        FindMinimum[Πi, Array[{a[#] , 0. } &, nterm], 
          Method -> "Gradient", 
          EvaluationMonitor :> (ii++;)];
{ii, ϕ /. sol[[2]] /. {x -> 0, y -> 0, 
          z -> 0.5}}
SessionTime[] - start

N. Newton method based optimization

start = SessionTime[]; ii = 0;
sol = 
        FindMinimum[Πi, Array[{a[#] , 0. } &, nterm], 
          Method -> "Newton", 
          EvaluationMonitor :> (ii++;)];
{ii, ϕ /. sol[[2]] /. {x -> 0, y -> 0, 
          z -> 0.5}}
SessionTime[] - start

The subsection Definitions has to be executed before the current subsection.

\(\(\(start = SessionTime[]; 
    SMSInitialize["Thermal", "Environment" -> "MathLink", "Mode" -> "Prototype", "ADMethod" -> "Forward"]\)\(
\)
    \)\)
Πf[
        i_] := (
ai ⊨ 
          SMSReal[Array[a$$, nterm]]; 
ag ⊨ 
          SMSArray[ai]; 
{xa, ya, za, wa} ⊨ 
          Map[SMSArray, Transpose[gp]]; 
{xi, yi, 
            zi} ⊨ 
          SMSFreeze[{SMSPart[xa, i], SMSPart[ya, i], 
              SMSPart[za, i]}]; 
{xpr, ypr, 
            zpr} ⊨ {xp /. x -> xi, yp /. y -> yi, 
            zp /. z -> zi}; 
poly ⊨ 
          SMSArray[
            Flatten[Outer[Times, xpr, ypr, 
                zpr]]]; 
ϕt ⊨ 
          SMSDot[poly, 
            ag]; 
Δϕ ⊨ 
          SMSD[ϕt, {xi, yi, 
              zi}]; 
wi ⊨ 
          SMSPart[wa, 
            i]; 
wi (1/
                  2 kcond Δϕ . Δϕ - ϕt Q)
)

SMSModule["FThermal", 
        Real[a$$[nterm], f$$], "Input" -> a$$, "Output" -> 
          f$$];
SMSExport[0, f$$];
SMSDo[i, 1, gp // Length];
\(	\(Π ⊨ Πf[i];\)\)
\(	\(SMSExport[\ Π, 
        f$$, "AddIn" -> True];\)\)
SMSEndDo[];

SMSModule["GThermal", 
        Real[a$$[nterm], g$$[nterm]], "Input" -> 
          a$$, "Output" -> g$$];
SMSExport[Table[0, {nterm}], g$$];
SMSDo[i, 1, gp // Length];
\(	\(Π ⊨ Πf[i];\)\)
\(	\(SMSDo[j, 1, nterm];\)\)
\(		\(δΠ ⊨ 
        SMSD[Π, ag, j, "Method" -> "Forward"];\)\)
\(		\(SMSExport[δΠ, 
        g$$[j], "AddIn" -> True];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

derivation = SessionTime[] - start

SMSModule["HThermal", 
        Real[a$$[nterm], h$$[nterm, nterm]], "Input" -> 
          a$$, "Output" -> h$$];
SMSDo[i, 1, nterm];
\(	\(SMSDo[j, 1, nterm];\)\)
\(		\(SMSExport[0, h$$[i, j]];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];
SMSDo[i, 1, gp // Length];
\(	\(Π ⊨ Πf[i];\)\)
\(	\(SMSDo[j, 1, nterm];\)\)
\(		\(δΠ ⊨ 
        SMSD[Π, ag, j, "Method" -> "Forward"];\)\)
\(		\(SMSDo[k, 1, nterm];\)\)
\(			\(hij ⊨ 
        SMSD[δΠ, ag, 
          k, "Method" -> "Forward"];\)\)
\(			\(SMSExport[hij, 
        h$$[j, k], "AddIn" -> True];\)\)
\(		\(SMSEndDo[];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

SMSWrite[];

SMSInstallMathLink["Optimize" -> 
        False]
derivation = SessionTime[] - start

G. Gradient based optimization

start = SessionTime[]; ii = 0;
sol = 
        FindMinimum[FThermal[ϕi], {ϕi, Table[0, {nterm}]}, 
          Method -> "Gradient", Gradient -> GThermal[ϕi], 
          EvaluationMonitor :> (ii++;)];
{ii, ϕ /. 
          MapThread[Rule, List @@ sol[[2, 1]]] /. {x -> 0, 
          y -> 0, z -> 0.5}, SessionTime[] - start}

N. Newton method based optimization

start = SessionTime[]; ii = 0;
\(\(sol = 
        FindMinimum[
          FThermal[ϕi], {ϕi, 
            Table[0, {nterm}]}, 
Method -> {"Newton", \ 
              Hessian -> 
                HThermal[ϕi]}, 
Gradient -> 
            GThermal[ϕi], 
          EvaluationMonitor :> \((\(\(ii++\);\))\)];\)\)
{ii, ϕ /. 
          MapThread[Rule, List @@ sol[[2, 1]]] /. {x -> 0, 
          y -> 0, z -> 0.5}, SessionTime[] - start}

B) Finite difference interpolation

The central difference approximation of derivatives is used for the points inside the cube and backward or forward difference for the points on the boundary.

<< AceGen`;
Clear[a, i, j, k];
nx = ny = nz = 11;
dlx = 1. /(nx - 1);
dly = 1. /(ny - 1);
dlz = 1. /(nz - 1);
bound = {0};
nboun = 1;
kcond = 0.58; Q = 500;

nterm = 0; dofs = {};
index = 
        Table[Which[
i ≤ 2  ||  
                i  ≥   nx + 1   ||  
                j ≤ 2  ||  
                j  ≥   ny + 1   ||  k ≤  2 , 
              b[1]
, 
              k == nz + 2, 
If[
                FreeQ[dofs, a[i, j, k - 1]]
, ++nterm; 
                AppendTo[dofs, a[i, j, k - 1] -> nterm]; 
                nterm
, 
                a[i, j, k - 1] /. 
                  dofs
]
, 
              True, 
If[
                FreeQ[dofs, a[i, j, k]]
, ++nterm; 
                AppendTo[dofs, a[i, j, k] -> nterm]; 
                nterm
, 
                a[i, j, k] /. 
                  dofs
]
], 
{i, 1, nx + 2}, {j, 1, ny + 2}, {k, 1, nz + 2}] /. 
          b[i_] :> nterm + i;
ϕi = Array[a, nterm];
nterm

The subsection Definitions have to be executed before the current subsection.

start = SessionTime[];
Πi = 
        Sum[
dlxt = 
            If[i == 2  ||  i == nx + 1 , dlxt = dlx/2, 
              dlx]; 
dlyt = 
            If[j == 2  ||  j == ny + 1 , dlyt = dly/2, 
              dly]; 
dlzt = 
            If[k == 2  ||  k == nz + 1 , dlzt = dlz/2, 
              dlz]; 
vol = 
            dlxt dlyt dlzt; 
aijk = 
            Map[If[# > nterm, bound[[# - nterm]], 
                  a[#]] &, 
Extract[
                index, {{i, j, k}, {i - 1, j, k}, {i + 1, j, k}, {i, j - 1, 
                    k}, 
{i, j + 1, k}, {i, j, k - 1}, {i, 
                    j, k + 1}}]]; 
grad = {((aijk[[3]] - aijk[[2]])/(2 dlxt)), ((aijk[[5]] - aijk[[4]])/(2 dlyt)), ((aijk[[7]] - aijk[[6]])/(2 dlzt))}; 

vol (1/2 kcond grad . grad - 
                Q aijk[[1]])
, {i, 2, nx + 1}, {j, 
            2, ny + 1}, {k, 2, 
            nz + 1}
];
derivation = SessionTime[] - start

G. Gradient based optimization

start = SessionTime[]; ii = 0;
sol = 
        FindMinimum[Πi, Array[{a[#] , 0. } &, nterm], 
          Method -> "Gradient", 
          EvaluationMonitor :> (ii++;)];
{ii, 
      a[index[[(nx + 3)/2, (ny + 3)/2, (nz + 3)/2]]] /. 
        sol[[2]], SessionTime[] - start}

N. Newton method based optimization

start = SessionTime[]; ii = 0;
sol = 
        FindMinimum[Πi, Array[{a[#] , 0. } &, nterm], 
          Method -> "Newton", 
          EvaluationMonitor :> (ii++;)];
{ii, 
      a[index[[(nx + 3)/2, (ny + 3)/2, (nz + 3)/2]]] /. 
        sol[[2]], SessionTime[] - start}

The subsection Definitions have to be executed before the current subsection.

\(\(\(start = SessionTime[]; 
    SMSInitialize["Thermal", "Environment" -> "MathLink", "Mode" -> "Prototype", "ADMethod" -> "Backward"]\)\(
\)
    \)\)
\(Πf[i_, j_, 
        k_] := \((
indexp ⊨ 
          SMSInteger[
            Map[
index$$[\((#[\([1]\)] - 1)\)*\((nyp + 
                          2)\) \((nzp + 2)\) + \((#[\([2]\)] - 1)\)*\((nzp + 
                          2)\) + #[\([3]\)]] &, 
{{i, j, 
                  k}, {i - 1, j, k}, {i + 1, j, k}, {i, j - 1, 
                  k}, 
{i, j + 1, k}, {i, j, k - 1}, {i, j, 
                  k + 1}}]]; 
aijk ⊨ 
          SMSReal[Map[a$$[#] &, indexp]]; 
{dx, dy, dz, kc, 
            Qt} ⊨ 
          SMSReal[Array[mc$$, 5]]; 
SMSIf[
          i == 2\  || \ i == nxp + 1\ ]; 
	
        dlxt ⫤ 
          dx/2; 
SMSElse[]; 
	
        dlxt ⊣ dx; 
SMSEndIf[
          dlxt]; 
SMSIf[
          j == 2\  || \ j == nyp + 1\ ]; 
	
        dlyt ⫤ 
          dy/2; 
SMSElse[]; 
	
        dlyt ⊣ dy; 
SMSEndIf[
          dlyt]; 
SMSIf[
          k == 2\  || \ k == nzp + 1\ ]; 
	
        dlzt ⫤ 
          dz/2; 
SMSElse[]; 
	
        dlzt ⊣ dz; 
SMSEndIf[
          dlzt]; 
vol ⊨ 
          dlxt\ dlyt\ dlzt; 
grad = {\(aijk[\([3]\)] - aijk[\([2]\)]\)\/\(2\ dlxt\), \(aijk[\([5]\)] - aijk[\([4]\)]\)\/\(2\ dlyt\), \(aijk[\([7]\)] - aijk[\([6]\)]\)\/\(2\ dlzt\)}; 
vol\ \((1/
                  2\ kc\ grad . grad - 
              Qt\ aijk[\([1]\)])\)
)\)\)

SMSModule["FThermal", Integer[ndof$$, nt$$[3], index$$["*"]], 
        Real[a$$["*"], mc$$["*"], 
          f$$], "Input" -> {ndof$$, nt$$, index$$, a$$, 
            mc$$}, "Output" -> f$$];
SMSExport[0, f$$];
{nxp, nyp, nzp} ⊨ 
        SMSInteger[Array[nt$$, 3]];
SMSDo[i, 2, nxp + 1];
\(	\(SMSDo[j, 2, nyp + 1];\)\)
\(		\(SMSDo[k, 2, nzp + 1];\)\)
\(			\(Π ⊨ Πf[i, j, 
          k];\)\)
\(			\(SMSExport[Π, 
        f$$, "AddIn" -> True];\)\)
\(		\(SMSEndDo[];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

SMSModule["GThermal", Integer[ndof$$, nt$$[3], index$$["*"]], 
        Real[a$$["*"], mc$$["*"], 
          g$$[ndof$$]], "Input" -> {ndof$$, nt$$, index$$, a$$, 
            mc$$}, "Output" -> g$$];
ndof ⊨ 
        SMSInteger[ndof$$];
{nxp, nyp, nzp} ⊨ 
        SMSInteger[Array[nt$$, 3]];
SMSDo[i, 1, ndof];
\(	\(SMSExport[0, g$$[i]];\)\)
SMSEndDo[];

    
SMSDo[i, 2, nxp + 1];
\(	\(SMSDo[j, 2, nyp + 1];\)\)
\(		\(SMSDo[k, 2, nzp + 1];\)\)
\(			\(Π ⊨ Πf[i, j, 
          k];\)\)
\(			\(SMSDo[i1, 1, indexp // Length];\)\)
\(				\(dof ⊨ 
        SMSPart[indexp, i1];\)\)
\(				\(SMSIf[dof <= ndof];\)\)
\(					\(gi ⊨ 
        SMSD[Π, aijk, i1];\)\)
\(					\(SMSExport[gi, 
        g$$[dof], "AddIn" -> True];\)\)
\(				\(SMSEndIf[];\)\)
\(			\(SMSEndDo[];\)\)
\(		\(SMSEndDo[];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

derivation = SessionTime[] - start

SMSModule["HThermal", Integer[ndof$$, nt$$[3], index$$["*"]], 
        Real[a$$["*"], mc$$["*"], 
          h$$[ndof$$, ndof$$]], "Input" -> {ndof$$, nt$$, index$$, 
            a$$, mc$$}, "Output" -> 
          h$$];
ndof ⊨ 
        SMSInteger[ndof$$];
{nxp, nyp, nzp} ⊨ 
        SMSInteger[Array[nt$$, 3]];
SMSDo[i, 1, ndof];
\(	\(SMSDo[j, 1, ndof];\)\)
\(		\(SMSExport[0, h$$[i, j]];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

    
SMSDo[i, 2, nxp + 1];
\(	\(SMSDo[j, 2, nyp + 1];\)\)
\(		\(SMSDo[k, 2, nzp + 1];\)\)
\(			\(Π ⊨ Πf[i, j, 
          k];\)\)
\(			\(SMSDo[i1, 1, indexp // Length];\)\)
\(				\(dofi ⊨ 
        SMSPart[indexp, i1];\)\)
\(				\(SMSIf[dofi <= ndof];\)\)
\(					\(gi ⊨ 
        SMSD[Π, aijk, i1];\)\)
\(					\(SMSDo[j1, 1, 
        indexp // Length];\)\)
\(						\(dofj ⊨ 
        SMSPart[indexp, j1];\)\)
\(						\(SMSIf[dofj <= ndof];\)\)
\(							\(hij ⊨ 
        SMSD[gi, aijk, j1];\)\)
\(							\(SMSExport[hij, 
        h$$[dofi, dofj], "AddIn" -> True];\)\)
\(						\(SMSEndIf[];\)\)
\(					\(SMSEndDo[];\)\)
\(				\(SMSEndIf[];\)\)
\(			\(SMSEndDo[];\)\)
\(		\(SMSEndDo[];\)\)
\(	\(SMSEndDo[];\)\)
SMSEndDo[];

SMSWrite[];

SMSInstallMathLink["Optimize" -> 
        True]
derivation = SessionTime[] - start

G. Gradient based optimization

start = SessionTime[]; ii = 0;
indexb = Flatten[index];
\(sol = 
      FindMinimum[
FThermal[nterm, {nx, ny, nz}, indexb, 
          Join[ϕi, bound], {dlx, dly, dlz, kcond, 
            Q}]
, {ϕi, 
          Table[0, {nterm}]}, 
Method -> "Gradient", 
Gradient -> 
          GThermal[nterm, {nx, ny, nz}, indexb, 
            Join[ϕi, bound], {dlx, dly, dlz, kcond, 
              Q}]
, 
        EvaluationMonitor :> \((\(\(ii++\);\))\)]; {ii, 
      a[index[\([\((nx + 3)\)/2, \((ny + 3)\)/2, \((nz + 3)\)/2]\)]] /. 
        MapThread[Rule, List @@ sol[\([2, 1]\)]], 
      SessionTime[] - start}\)

N. Newton method based optimization

start = SessionTime[]; ii = 0;
indexb = 
        Flatten[index /. 
            b[i_] :> nterm + i];
\(sol = 
      FindMinimum[
FThermal[nterm, {nx, ny, nz}, indexb, 
          Join[ϕi, bound], {dlx, dly, dlz, kcond, 
            Q}]
, {ϕi, 
          Table[0, {nterm}]}, 
Method -> {"Newton", \ 
            Hessian -> 
              HThermal[nterm, {nx, ny, nz}, indexb, 
                Join[ϕi, bound], {dlx, dly, dlz, kcond, 
                  Q}]}, 
Gradient -> 
          GThermal[nterm, {nx, ny, nz}, indexb, 
            Join[ϕi, bound], {dlx, dly, dlz, kcond, 
              Q}]
, 
        EvaluationMonitor :> \((\(\(ii++\);\))\)]; {ii, 
      a[index[\([\((nx + 3)\)/2, \((ny + 3)\)/2, \((nz + 3)\)/2]\)]] /. 
        MapThread[Rule, List @@ sol[\([2, 1]\)]], 
      SessionTime[] - start}\)

The tangent matrix is in the case of finite difference approximation extremely sparse.

MatrixPlot[
        HThermal[nterm, {nx, ny, nz}, index, 
          Join[0 ϕi, bound], {dlx, dly, dlz, kcond, Q}]];

C) Finite element method

First the finite element mesh 30×30×30 is used to obtain convergence solution at the central point of the cube. The procedure to generate heat-conduction element that is used in this example is explained in AceGen manual section  Description of FE Characteristic Steps.

<< AceFEM`;
start = SessionTime[];
SMTInputData[];
k = 0.58; Q = 500;
nn = 30;
SMTAddDomain["cube", "heatconduction", {k, 0, 0, 
          Q}];
\(\(SMTAddEssentialBoundary[{"X" == \(-0.5\) || \ "X" == 0.5\  || \ "Y" == \(-0.5\)\  || \ "Y" == 
                0.5\  || \ "Z" == 0.  &, 0}];\)\)
SMTMesh["cube", "H1", {nn, nn, 
          nn}, {
{{{-0.5, -0.5, 0}, {0.5, -0.5, 
                0}}, {{-0.5, 0.5, 0}, {0.5, 0.5, 
                0}}}, {{{-0.5, -0.5, 1}, {0.5, -0.5, 
                1}}, {{-0.5, 0.5, 1}, {0.5, 0.5, 
                1}}}
}];
SMTAnalysis["Solver" -> 5];

SMTNextStep[0, 1];
SMTNewtonIteration[];

SMTPointValues[{0, 0, 0.5}, SMTPost[1]]
SessionTime[] - start

SMTShowMesh["Mesh" -> True, "Elements" -> 
          True, "Field" -> SMTPost[1], "Contour" -> 
          True];

Coupled thermo-mechanical problem

Four noded solid element -solidQ1

<< AceGen`; << AceFEM`;
\(\(SMSInitialize[
        ctfile = "solidQ1", "Environment" -> "AceFEM", "VectorLength" -> 3500];\)\)
\(\(SMSTemplate["SMSTopology" -> "Q1", "SMSNoNodes" -> 4, "SMSDOFGlobal" -> 
          2, "SMSNodeID" -> "D"];\)\)
SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex, 1, 
        SMSInteger[es$$["id", "NoIntPoints"]]];
{ξ, η, ζ, wGauss} ⊢ 
        Array[SMSReal[es$$["IntPoints", #1, IpIndex]] &, 4];
XYZ ⊨ 
        SMSReal[Table[
            nd$$[i, "X", j], {i, SMSNoNodes},  {j, 
              SMSNoDimensions}]];
{Xi, Yi} ⊨ 
        Transpose[XYZ];
Ni ⊨ 
        (1/4) {(1 - ξ) (1 - η), (1 + ξ) (1 - η), (1 + ξ) (1 + η), (1 - ξ) (1 + η)};
X ⊢ SMSFreeze[Ni . Xi];
Y ⊢ SMSFreeze[Ni . Yi];
Z ⊢ SMSFreeze[ζ];
Jm ⊨ 
        SMSD[{X, Y, Z}, {ξ, η, ζ}];
SMSDefineDerivative[{ξ, η, ζ}, {X, Y, Z}, 
        SMSInverse[Jm]];

    
dof ⊨ 
        SMSReal[Table[
            nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
𝕒t ⊨ 
        Flatten[dof];
ui ⊨ 
        Table[dof[[i, 1]], {i, 
            SMSNoNodes}];
vi ⊨ 
        Table[
          dof[[i, 2]], {i, 
            SMSNoNodes}];
u ⊨ Ni . ui; 
    v ⊨ 
      Ni . vi;

SMSGroupDataNames = {"Elastic modulus", "Poisson \
ration"};
{𝔼, ν} ⊨ 
        SMSReal[{es$$["Data", 1], 
            es$$["Data", 2]}];
{λ, μ} ⊨ 
        SMSHookeToLame[𝔼, ν];

    
𝔻 ⊨ 
        SMSD[{u, v, 0}, {X, Y, Z}];
ϵ ⊨  (𝔻 + 
              Transpose[𝔻])/
          2;
Π ⊨ 
        Simplify[λ/2 
               Tr[ϵ]^2 + μ ((ϵ ϵ) /. 
                  List -> Plus)];
Ψi ⊨ 
        wGauss Det[Jm] SMSD[Π, 𝕒t];
SMSExport[SMSResidualSign Ψi, 
        p$$, "AddIn" -> True];
Kij ⊨ 
        SMSD[Ψi, 𝕒t];
SMSExport[Kij, s$$, "AddIn" -> True];
SMSEndDo[]; 

SMSStandardModule["Postprocessing"];
dof ⊨ 
        SMSReal[Table[
            nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
ui ⊨ 
        Table[dof[[i, 1]], {i, 
            SMSNoNodes}];
vi ⊨ 
        Table[dof[[i, 2]], {i, 
            SMSNoNodes}];
SMSNPostNames = {"DeformedMeshX", "DeformedMeshY"};
SMSExport[{ui, vi} // Transpose, npost$$];

    
SMSWrite[];

Four noded solid element with thermal expansion - solidthermalQ1

<< AceGen`; << AceFEM`;
\(\(SMSInitialize[
        ctfile = "solidthermalQ1", "Environment" -> "AceFEM", "VectorLength" -> 3500];\)\)
\(\(SMSTemplate["SMSTopology" -> "Q1", "SMSNoNodes" -> 8, 
"SMSDOFGlobal" -> {2, 2, 2, 2, 1, 1, 
            1, 1}, 
"SMSAdditionalNodes" -> "{#1,#2,#3,#4}&", 
"SMSNodeID" -> {"D", "D", "D", "D", "T", "T", "T", "T"}, 
"SMSSegments" -> {{1, 5}, {2, 6}, {3, 7}, {4, 
              8}, {1, 2, 3, 4}}, "SMSSymmetricTangent" -> 
          False];\)\)
SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex, 1, 
        SMSInteger[es$$["id", "NoIntPoints"]]];
{ξ, η, ζ, wGauss} ⊢ 
        Array[SMSReal[es$$["IntPoints", #1, IpIndex]] &, 4];
XYZ ⊨ 
        SMSReal[Table[
            nd$$[i, "X", j], {i, SMSNoNodes},  {j, 
              SMSNoDimensions}]];
{Xi, Yi} ⊨ 
        Transpose[Take[XYZ, 4]];
Ni ⊨ 
        (1/4) {(1 - ξ) (1 - η), (1 + ξ) (1 - η), (1 + ξ) (1 + η), (1 - ξ) (1 + η)};
X ⊢ SMSFreeze[Ni . Xi];
Y ⊢ SMSFreeze[Ni . Yi];
Z ⊢ SMSFreeze[ζ];
Jm ⊨ 
        SMSD[{X, Y, Z}, {ξ, η, ζ}];
SMSDefineDerivative[{ξ, η, ζ}, {X, Y, Z}, 
        SMSInverse[Jm]];

    
dof ⊨ 
        SMSReal[
          Table[nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
𝕒t ⊨ 
        Flatten[dof];
ui ⊨ 
        Table[dof[[i, 1]], {i, 
            4}];
vi ⊨ 
        Table[dof[[i, 2]], {i, 
            4}];
ϕi ⊨ 
        Table[dof[[i, 1]], {i, 5, 
            8}];
u ⊨ Ni . ui; 
    v ⊨ Ni . vi; ϕ ⊨ 
      Ni . ϕi;

\(\(SMSGroupDataNames = {"Elastic modulus", "Poisson ration", "Thermal expansion coefficient"};\)\)
{𝔼, ν, αt} ⊨ 
        SMSReal[{es$$["Data", 1], es$$["Data", 2], 
            es$$["Data", 3]}];
{λ, μ} ⊨ 
        SMSHookeToLame[𝔼, ν];

    
𝔻 ⊨ 
        SMSD[{u, v, 0}, {X, Y, Z}];
ϵ ⊨  (𝔻 + 
                Transpose[𝔻])/
            2 + αt IdentityMatrix[
              3] ϕ;
Π ⊨ 
        Simplify[λ/2 
               Tr[ϵ]^2 + μ ((ϵ ϵ) /. 
                  List -> Plus)];
Ψi ⊨ 
        wGauss Det[Jm] Join[
            SMSD[Π, Take[𝕒t, 8]], {0, 0, 0, 
              0}];
SMSExport[SMSResidualSign Ψi, 
        p$$, "AddIn" -> True];
Kij ⊨ 
        SMSD[Ψi, 𝕒t];
SMSExport[Kij, s$$, "AddIn" -> True];
SMSEndDo[]; 

SMSStandardModule["Postprocessing"];
dof ⊨ 
        SMSReal[Table[
            nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
ui ⊨ 
        Table[dof[[i, 1]], {i, 
            4}];
vi ⊨ 
        Table[dof[[i, 2]], {i, 
            4}];
ϕi ⊨ 
        Table[dof[[i, 1]], {i, 5, 
            8}];
\(\(SMSNPostNames = {"DeformedMeshX", "DeformedMeshY", "T"};\)\)
SMSExport[
        Join[Transpose[{ui, vi, ϕi}], Map[{0, 0, #} & , ϕi]], 
        npost$$];
SMSWrite[];

6 noded heat conduction element - heatT2

<< AceGen`; << AceFEM`;
\(\(SMSInitialize[
        ctfile = "heatT2", "Environment" -> "AceFEM", "VectorLength" -> 3500];\)\)
\(\(SMSTemplate["SMSTopology" -> "T2", "SMSDOFGlobal" -> 1, "SMSNodeID" -> "T"];\)\)
SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex, 1, 
        SMSInteger[es$$["id", "NoIntPoints"]]];
{ξ, η, ζ, wGauss} ⊢ 
        Array[SMSReal[es$$["IntPoints", #1, IpIndex]] &, 4];
XYZ ⊨ 
        SMSReal[Table[
            nd$$[i, "X", j], {i, SMSNoNodes},  {j, 
              SMSNoDimensions}]];
{Xi, Yi} ⊨ 
        Transpose[XYZ];
ζc = 1 - ξ - η;
Ni ⊨ {(2 ξ - 1) ξ, (2 η - 
                1) η, (2 ζc - 1) ζc,  
          4 ξ η, 4  η ζc, 
          4 ζc ξ};
X ⊢ SMSFreeze[Ni . Xi];
Y ⊢ SMSFreeze[Ni . Yi];
Z ⊢ SMSFreeze[ζ];
Jm ⊨ 
        SMSD[{X, Y, Z}, {ξ, η, ζ}];
SMSDefineDerivative[{ξ, η, ζ}, {X, Y, Z}, 
        SMSInverse[Jm]];

    
dof ⊨ 
        SMSReal[Table[
            nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
𝕒t ⊨ 
        Flatten[dof];
ϕi ⊨ 
        Table[dof[[i, 1]], {i, 
            SMSNoNodes}];
\(\(\(ϕ ⊨ Ni . ϕi;\)\(
\)
    \)\)
SMSGroupDataNames = {"Conductivity parameter", "Heat \
source"};
{k, Q} ⊨ 
        SMSReal[{es$$["Data", 1], 
            es$$["Data", 2]}];
δϕ ⊨ SMSD[ϕ, ϕi];
Ψi ⊨ 
        Det[Jm] wGauss (k SMSD[δϕ, {X, Y}] . 
                  SMSD[ϕ, {X, Y}] - δϕ Q);
SMSExport[SMSResidualSign Ψi, 
        p$$, "AddIn" -> True];
Kij ⊨ SMSD[Ψi, ϕi];
SMSExport[Kij, s$$, "AddIn" -> True];
SMSEndDo[]; 

SMSStandardModule["Postprocessing"];
dof ⊨ 
        SMSReal[Table[
            nd$$[i, "at", j], {i, SMSNoNodes}, {j,  
              SMSDOFGlobal[[i]]}]];
ϕi ⊨ 
        Table[dof[[i, 1]], {i, 
            SMSNoNodes}];
\(\(SMSNPostNames = {"DeformedMeshX", "DeformedMeshY", "T"};\)\)
SMSExport[Map[{0, 0, #} & , ϕi], 
        npost$$];

    
SMSWrite[];

Example: 2D ribbon

Solid

<< AceFEM`;
SMTInputData[];

    
SMTAddDomain["D", "solidQ1", {3000,  .2}];
SMTMesh["D", "Q1", {80, 5}, 
        Array[ {#2, #2 Sin[#2 // N]/20 + 4 #1 } &, {2, 40}]];
SMTAddNaturalBoundary[
        Abs["X" Sin["X" // N]/20 + 8 - "Y"] < 
              0.1  &&  "ID" == "D" &, 0, -5];
SMTAddEssentialBoundary["X" == 
              1  &&  "ID" == "D"  &, 0, 
        0];
SMTAddEssentialBoundary["X" == 
              40 &&  "ID" == "D"  &, 0, 0];
SMTAnalysis[];

SMTNextStep[1, 1];
SMTNewtonIteration[]
SMTNewtonIteration[]

Show[
        SMTShowMesh["Show" -> 
            False, "BoundaryConditions" -> True, "Mesh" -> 
            True, "Domains" -> "D", "NodeMarks" -> 
            True], 
SMTShowMesh["Show" -> 
            False, "BoundaryConditions" -> 
            True, "DeformedMesh" -> True, "Field" -> 
            SMTPost[2], "Legend" -> 
            False], 
PlotRange -> All];

Solid - thermal

<< AceFEM`;
SMTInputData[];

    
\(\(SMTAddDomain["D", "solidthermalQ1", {3000,  .2,  .00001}];\)\)
SMTMesh["D", "Q1", {40, 5}, 
        Array[ {#2, #2 Sin[#2 // N]/20 + 4 #1 } &, {2, 40}]];
SMTAddNaturalBoundary[
        Abs["X" Sin["X" // N]/20 + 8 - "Y"] < 
              0.1  &&  "ID" == "D" &, 0, -5];
SMTAddEssentialBoundary["X" == 
              1  &&  "ID" == "D"  &, 0, 
        0];
SMTAddEssentialBoundary["X" == 
              40 &&  "ID" == "D"  &, 0, 
        0];

SMTAddDomain["H", "heatT2", {0.5, 
          300}];
SMTMesh["H", "T2", {40, 5}, 
        Array[ {#2, #2 Sin[#2 // N]/20 + 4 #1 } &, {2, 40}]];
\(\(SMTAddEssentialBoundary[\(("X" == 1\  || \ 
                Abs["X"\ Sin["X" // N]/20 + 8 - "Y"] < 
                  0.1)\)\  && \ "ID" == "T"\  &, 0];\)
\)
SMTAnalysis[];

SMTNextStep[1, 1];
SMTNewtonIteration[]
SMTNewtonIteration[]

SMTShowMesh["Show" -> True, "BoundaryConditions" -> 
        True, "NodeMarks" -> 
        True, 
"Mesh" -> 
        True, "Domains" -> "D"]

\(SMTShowMesh["Show" -> True, "BoundaryConditions" -> 
        True, "NodeMarks" -> True, "Mesh" -> 
        True, "DeformedMesh" -> False, "Legend" -> 
        False, "Domains" -> "H", "Field" -> "T", "Contour" -> 100]\)

def = 
        Transpose[
          SMTNodeData["NodeID"] /. {"D" -> {0, 
                  0}, "T" -> {5, 5}}];
Show[
SMTShowMesh["Show" -> 
            False, "Mesh" -> 
            False, 
"DeformedMesh" -> 
            True, "Scale" -> 2, "Mesh" -> 
            True, "Domains" -> "D", "NodeMarks" -> 
            False], 
SMTShowMesh["Show" -> 
            False, "BoundaryConditions" -> 
            False, 
"DeformedMesh" -> 
            def, "Mesh" -> 
            True, "Domains" -> "D", "NodeMarks" -> 
            True], 
SMTShowMesh["Show" -> 
            False, "BoundaryConditions" -> False, "Mesh" -> 
            True, "Domains" -> "H", "NodeMarks" -> 
            True, "Field" -> "T", "DeformedMesh" -> 
            def, "Legend" -> False, "Contour" -> 
            100], 
PlotRange -> All];