AceFEM: Efficient calculation of volume of axisymmetric elements

Question

I am looking for the most efficient calculation of volume of axisymmetric elements in finite element analysis with AceFEM package. Ideally I would like to have this implemented in AceGen element code as part of standard subroutine ("Postprocessing" or "Tasks").

This calculation will be performed at in every load/time step, for large number of elements, therefore I want it to be as fast as possible to avoid creating an unnecessary bottleneck. Bellow is a simple example of axisymmetric domain and we know its total volume is Pi (3.14159...).

<< AceFEM`;
(* Simple domain with 4 elements, assuming axial symmetry around Y axis. Element code can be downloaded from AceShare. *)
SMTInputData[];
SMTAddDomain["test","OL:SEAXQ1DFHYQ1NeoHooke", {"E *" -> 1000., "\[Nu] *" -> 0.3}];
SMTAddMesh[
  {{0, 0}, {0.5, 0}, {1, 0}, {0, 0.5}, {0.4, 0.4}, {1, 0.5}, {0, 1}, {0.5, 1}, {1, 1}},
  {"test" -> {{1, 2, 5, 4}, {2, 3, 6, 5}, {4, 5, 8, 7}, {5, 6, 9, 8}}}
  ];
SMTAnalysis[];
SMTShowMesh["Marks"->True,ImageSize->100,AxesLabel->{"X","Y"},Axes->True,Ticks -> None]

This is calculation of node coordinates for each element. Displacements (SMTNodeData["at"]) are currently 0, but they could be arbitrary.

elementCoordinates := With[
  {nodeCoordinates = SMTNodeData["X"] + SMTNodeData["at"]},
  nodeCoordinates[[#]] & /@ SMTElementData["Nodes"]
  ]
elementCoordinates

The simplest high level function for calculating volume of solid of revolution. It is short, but very slow!

volumeFun = 2*Pi*First[RegionCentroid[#]]*Area[#] &;
Total[
  Composition[volumeFun, Polygon] /@ elementCoordinates
  ] // RepeatedTiming
(* {0.014, 3.14159} *)

This is an alternative function with low level arithmetic that can be efficiently compiled. It is much faster and can be even improved by compiling to C. But are there even faster methods?

(* Function splits rectangular polygon to two triangles and calculates their area and centroids with simple formulas and then sums two volumes together.*)
volumeFun2 = Compile[{{pts, _Real, 2}},
   Module[{triangles, area, centroid},

    triangles = {pts[[{1, 4, 3}]], pts[[{1, 3, 2}]]};
    area = (Abs[(#[[1, 1]] - #[[3, 1]]) (#[[2, 2]] - #[[1, 2]]) - (#[[1, 1]] - #[[2, 1]]) (#[[3, 2]] - #[[1, 2]])]/2) & /@ triangles;
    centroid = Apply[Plus, ({#[[All, 1]], #[[All, 2]]}/3) & /@ triangles, {2}];

    Total[2.*Pi* centroid[[All, 1]]*area]
    ]
   ];

Total[
  volumeFun2 /@ elementCoordinates
  ] // RepeatedTiming
(* {0.00017, 3.14159} *)

Question: How can I speed up this procedure by implementing it inside element code in AceGen?

---

EDIT:

Another example with larger number of elements in initial and deformed configuration. I have downloaded element source code ("OL:SEAXQ1DFHYQ1NeoHooke") from AceShare and complied it locally. Then I can also modify source code with new tasks and postprocessing fields.

ClearAll[example];
(* syntax: example[type of configuration, number of elements per edge]*)
example["Initial", n_Integer] := (
   SMTInputData[];
   SMTAddDomain["test", "MySEAXQ1DFHYQ1NeoHooke", {}];
   SMTAddMesh[Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], "test","Q1", {n, n}];
   SMTAddEssentialBoundary[{{"X" == 0 &, 1 -> 0, 2 -> 0}, {"Y" == 0 &, 2 -> 0}, {"X" == 1 &, 1 -> 0.5}}];
   SMTAnalysis[];
   );

example["Deformed", n_Integer] := (
   example["Initial", n];
   SMTNextStep["\[CapitalDelta]\[Lambda]" -> 1];
   While[
    SMTConvergence[10^-8, 15],
    SMTNewtonIteration[];
    ];
   );

Volume of initial configuration, calculated with simplest and slowest function.

example["Initial", 10]
Total[Composition[volumeFun, Polygon] /@ elementCoordinates]
(* 3.14159 *)

Volume of deformed configuration. This example can be used to compare different methods to calculate volume of elements/domain.

example["Deformed", 10]
SMTShowMesh["BoundaryConditions" -> True, "DeformedMesh" -> True, ImageSize -> 100]

Total[Composition[volumeFun, Polygon] /@ elementCoordinates]
(* 4.2285 *)

Answer 1

The best way to get the Volume of the mesh is to use the "Tasks" Standard module. You have some choices what to input/output of the function. There are 6 task types available (see the AceGen help on "User defined Tasks" for more details). The Task type 1 is what you want in your case:

SMSStandardModule["Tasks"];
task \[DoubleRightTee] SMSInteger[Task$$];

You have to choose the names of the task which you will call:

SMSCharSwitch = {"V"};

You have to define the type on the first place (1), the number of real/integer input/output parameters. For volume we need type 1 and 1 real output, so we write: {1,0,0,0,1} to TasksData$$ for the first task (-1!):

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

We can then export volume. You can just loop over all the int. points points and export the volume 2 \[Pi] X Jed wgp at each point to RealOutput$$[1], where Jed=Det[SMSD[\[DoubleStruckCapitalX],\[CapitalXi]]], wgp is the int. point weight. "AddIn"->True is here needed, because we do not want to overwrite the volume each time:

SMSDo[Ig, 1, SMSInteger[es$$["id", "NoIntPoints"]]];
    ElementDefinitions[];
    SMSSwitch[task
    , 1, SMSExport[2 \[Pi] X Jed wgp , RealOutput$$[1], "AddIn" -> True];
    , _, 0;
    ];
SMSEndDo[];

Alternative to task is also to export the raw postprocessing data at each integration point by using the SMTPost command (one needs to export the 2 \[Pi] X Jed wgp in PP):

SMSDo[Ig, 1, SMSInteger[es$$["id", "NoIntPoints"]]];
    ElementDefinitions[];
    GPlotQuantities={...};
    AppendTo[GPlotQuantities, {"V", 2 \[Pi] X Jed wgp}]
    SMSGPostNames = GPlotQuantities[[All, 1]];
    SMSExport[GPlotQuantities[[All, 2]], gpost$$[Ig, #1] &];
SMSEndDo[];

For test I have used higher density mesh to check the efficiency, because the RepeatedTiming does not return representative measure of time on low density meshes, since the summation is then neglected, and for high density, you can use Timing or AbsoluteTiming:

...
SMTAddMesh[Polygon[{{0,0},{1,0},{1,1},{0,1}}],"test","Q1",{40,40}];    (*replaces SMTAddMesh*)
...
SMTAnalysis[];

Then You can use different volume calculations the test (here are your examples):

Total[Composition[volumeFun,Polygon]/@elementCoordinates]//AbsoluteTiming
(*{6.12715,3.14159}*)
Total[volumeFun2/@elementCoordinates]//AbsoluteTiming
(*{0.0106032,3.14159}*)

The Tasks values are caled here:

SMTTask["V","Summation"->True]//AbsoluteTiming 
(*{0.00105146,3.14159}*)
Total@SMTTask["V","Summation"->False]//AbsoluteTiming
(*{0.0040994,3.14159}*)

The unsmoothened values of postrocessing quantities are exported for each element, whoose summ equals to total volume:

Total[Flatten[#]&@SMTPost["V","Smooth"->False]]//AbsoluteTiming
(*{0.00463692,3.14159}*)

Note the SMTTask and SMTPost commands are comparable, if you suppress the summation in SMTTask and use Mathematica´s Total and it takes it takes 4 times longer.

EDIT: If you need to export the volume of the deformed configuration, you can use identical procedure as for undeformed, and just multiply the nodal point voulume with JE i.e. JE 2 \[Pi] X Jed wgp. JF is the determinant of deformation gradient, which tells the relative volume change of element:

\[DoubleStruckCapitalF = IdentityMatrix[3] + \[DoubleStruckCapitalH], 
JF = Det[\[DoubleStruckCapitalF]];

I added export of deformed volume (2) and also export both voulumes at same time (3):

SMSCharSwitch = {"V0","V","Vall"};

SMSIf[task < 0
    , SMSSwitch[task
    , -1, SMSExport[{1, 0, 0, 0, 1}, TasksData$$];
    , -2, SMSExport[{1, 0, 0, 0, 1}, TasksData$$];
    , -3, SMSExport[{1, 0, 0, 0, 2}, TasksData$$];
    , _, 0
   ];
   SMSReturn[];
];
...
fGauss \[DoubleRightTee]  2 \[Pi] X Jed wgp
SMSSwitch[task
    , 1, SMSExport[fGauss , RealOutput$$[1], "AddIn" -> True];
    , 2, SMSExport[JF fGauss , RealOutput$$[1], "AddIn" -> True];
    , 3, SMSExport[{1, JF} fGauss, RealOutput$$[#1] &, "AddIn" -> True];
    , _, 0;
];

The efficiency of each is simmilar (I used the 200x200 mesh):

SMTTask["Vall"] // AbsoluteTiming
(*{0.0192612, {3.14159, 4.19879}}*)
SMTTask["V"] // AbsoluteTiming
(*{0.0185199, 4.19879}*)
SMTTask["V0"] // AbsoluteTiming
(*{0.0175194, 3.14159}*)