Hello fellow Ace Users.

Currently I'm working on a project to implement Peridynamics. This is a discretization technique in the fashion of a meshless particle method. AceGen/AceFEM provides the feature of arbitrary nodes per element which suits my need perfectly as such a peridynamic particle interacts with arbitrary number of neighbour particles. To use the benefits of this method such as modelling discontinuities I'm aiming to utilize an explicit solution procedure.

I appreciate any thoughts on this! I have some code running in AceGen/AceFEM so far, still struggling on some design decisions which lead to the following specific questions:

    1. Whats exactly prarallelized in AceFEM? My recent experience indicate that SMSStandardModule["Tasks"] is not. Is that correct ? How about SMSStandardModule["Tangent and residual"] (I'm talking about the evaluation of the elements not solving the global equation system in parallel.)?
    1. Is there any known (maybe approximate) limit to the performance regarding arbitrary nodes per element?
    1. Does anyone have experiences with explicit simulations in AceFEM/AceGen?
    1. I expect a lot of data due to the particle discretization. Visualisation in post processing will be a to hard task to do in Mathematica. Does anyone have experiences with exporting the simulation data for use in e.g Paraview? If so, what's most performant way to write these to a file without significantly slowing down the simulation? I'm aware of the SMTPut[] feature by the way, but to my knowledge this binds me to Mathematica again.

As always You have my kudos in advance and I'm excited for your comments and answers !

Thanks for the response so far.

I'm back with a 'minimal' example that shows my main concerns.

The code is provided in SimpePDImplementation on GitHub:

The element contains a vary basic implementation of explicit peridynamics following two steps for each time step:

    1. Compute force density for each node (based on its neighbours)
    1. Integrate in time: acceleration = force density / density (per node)

These two task are implemented twice (using the same code), for once implemented into the SKR subroutine and for once as individual element tasks. !As the code is explicit i do not want or have a system of equation to solve, but i definitely want go over all elements in parallel to gain speedup. The results for both implementations are the same as expected, however the SKR implementation run significantly slower (i guess do to the solution of the linear system which is completely zero in this case). enter image description here

While performing the analysis, I checked my CPU usage.

For the SKR implementation I get: CPU usage for SKR Implementation. While for the Task implementation as reported I have: enter image description here

My conclusion so far is that the parallelization only works on the solution of the linear system and at least does not parallelize the loop over all elements for tasks.

I would be great if one of you guy can confirm, or even better of course tell me what I made wrong so that I know wether AceFEM works for my purpose at all.

Best, S


2 Answers 2


I can mostly help with the implementation aspect of meshfree methods and answer point two. Although Mathematica can plot anything you want, but you might need to program it yourself.

  1. In the new version of AceGen/FEM 6.923, prof Korelc introduced a new SMSIO functions and there is improvement of SMSArray usage and some new global fields. With that the formulation of elements with variable number of nodes is greatly simplified and the efficiency is independent of maximum number of nodes. You can put 1000 nodes if you want, the efficiency will not change more than couple %. To see the Improvement you need to introduce your vector of DOF as an array and not as a list with a fixed length of maximum possible values, additionally there are some fields that will give you actual number of nodes (I haven't tested that part). Here is a simple code extract of important things:

Function additionalNodes needs to pad actual amount of nodes with Null to match the max length for each element:


additionalNodes[nodes_, maxNodes_] := Module[{n},
   If[Length[nodes] > maxNodes,
     SMCError = {"Maximum number of nodes:", maxNodes, "Nodes:", nodes};
     SMCAbort["Given number of nodes exceeded the maximum number.", "",""];
   Join[ConstantArray[Null, maxNodes - Length[nodes]]]

You need to specify maximum no. of nodes, all nodes should be given as dummy-nodes and SMSAdditionalNodes should be assigned the above function:

 "SMSTopology" -> "QX"
 , "SMSNoNodes" -> maxNodes
 , "SMSDOFGlobal" -> 2
 , "SMSNodeID" -> "D -D"
 , "SMSAdditionalNodes" -> Function[additionalNodes[{##},maxNodes]]
 , "SMSMMAInitialisation" -> {{Definition[additionalNodes]}, Null}

Inside the SMSStandardModule we can extract the actual number of nodes/DOF from a new fields (Or if you like from ed$$["Data",_] field)

nNodes ⊢ SMSLastTrueNode[];
nDOF ⊢ SMSLastTrueDOF[];
nDOFNode = SMSNoDimensions;

The vector of all unknowns can then be defined as follows with SMSArray:

pe ⊨ SMSArray[nDOF, 
 Function[{dof}, SMSIO["Nodal DOFs"[(dof - 1)/nDOFNode + 1,SMSMod[dof-1,nDOFNode] + 1]]]

Then we can loop over arbitrary points/edges... we can read the displacement and coordinate of each point with node defined as an SMSInteger or a loop parameter as:

XNode ⊢ Table[SMSReal[nd$$[node, "X", j]], {j, SMSNoDimensions}];
uNode ⊨ Table[SMSPart[pe, (node - 1) nDOFNode + j], {j, nDOFNode}];

After your define your potential, the Residual and tangent can be derived same way as before, but only looped over the actual no of DOF defined by nDOF:

SMSDo[Rgi ⊨ SMSD[W, pe, i, "Constant"->SMSVariables[pseudoWConstants]];
  SMSIO[Rgi, "Add to", "Residual"[i]];
    Kgij ⊨SMSD[Rgi, \[DoubleStruckP]e, j];
    SMSIO[Kgij, "Add to", "Tangent"[i, j]];
  ,{j, i, nDOF}
,{i, 1, nDOF}
  • $\begingroup$ Hey Blaz, thanks for the comments. In my experiences Mathematica just cannot handle really large number of Graphic primitives in 3D very well. The implementation you give here is quite interesting (of course especially the SMCError and SMCAbort functionality ;-)). Nevertheless I couldn't make your approach work for me (see the code I posted on GitHub) because I get an Error due to my need of NodeData fields. Can you specify NodeData in the element you showed snippets above ? $\endgroup$ Apr 30, 2019 at 9:06
  • $\begingroup$ Hi, Your GitHub code just gives me error that your mesh is bad: Elem list is not list of numbers. Error has nothing to do with your AceGen code, since it runs when i take a part of Elem that are numbers. And node history is there. Concerning the efficiency. I switched off the assembly of tangent and residual by using: SMTIData["SkipTangent", 1];SMTIData["SkipResidual", 1];SMTIData["SkipSolver", 1], and then i got much faster time (still slower as task). SkipSolver seems to have no effect in my case. So I guess AceFEM tries to assemble empty Tangent and residual and tries to solve it. $\endgroup$ May 1, 2019 at 17:42
  • $\begingroup$ You should be able to set NoNodeStorage same way as you do now. I think SMSIO doesn't supports it jet since it is not in help like element history data. In our case the efficiency improvement comes from SMSArray, since we can now differentiate potential with respect to variable SMSArray instead of a length list to get residual and tangent. $\endgroup$ May 1, 2019 at 17:58
  • $\begingroup$ Thanks for surveying me code, I really appreciate the effort! However, do you agree that it is weird that the procedure by tasks only uses one core ? Or in other words: Am I doing smth. wrong or is this part just not parallelised, however still faster than using SKR for my purpose ? $\endgroup$ May 1, 2019 at 18:30
  • $\begingroup$ It seems the Tasks are not parallelized, but they are still faster, probably because when you call SMTNewtonIteration[], it doesn't just evaluate the tangent, it might do some other things related to solving FEM and preparing matrices that are preformed in the C-language. $\endgroup$ May 2, 2019 at 19:14

I have no experience with meshless methods, but I will try to answer/comment on your questions.

  1. In AceFEM assembling of global matrices and vectors ("Tangent and residual" subroutine) is parallelized. Solving the linear system is also parallelized (Intel MKL PARDISO). "Tasks" subroutine is not parallelized, but in general such procedures should not be computationally too intensive.

  2. Virtual element method (VEM) has been implemented in AceGen/AceFEM and it seems many (up to 100) nodes per element are possible. See Aldakheel et. al., 2018 for example. See even better explanation in BHudobivnik answer.

  3. According to comment by prof. Korelc, explicit simulations are possible. If this helps, Schmied et. al, 2013 have implemented element assembly subroutines in AceGen and used them in LS-Dyna software with explicit time integration.

  4. This really depends what you would like to visualize. I would perform visualization out of Mathematica only if there was really no other option. For that purpose you can always save a limited subset of results in some efficient common data format. Maybe something like HDF? As always, I would recommend to start with small example and deal with problems of large example when/if it happens.

  • $\begingroup$ Pinti ur awesome! I'll check about my issue with parallelization in Tasks and return with a small example when I'm sure then that it is not a problem of my implementation. Thanks also for the paper recommendations. To your last point though, don't you think the IO from AceFEM to Mathematica and then Mathematica to a file is somewhat like a detour? $\endgroup$ Apr 24, 2019 at 8:15
  • $\begingroup$ You are right, it is a detour, but Mathematica gives you all the flexibility you want. AceFEM exports data in HDF format (.h5) so maybe you can already import that to some other software such as Paraview. $\endgroup$
    – Pinti
    Apr 24, 2019 at 8:22
  • $\begingroup$ Now that you mention it, i have problems reading in the .hdf output from AceFEM in Mathematica (Im aware how it works WITH AceFEM but I would like to do post processing while running a simulation occupying AceFEM). None of the following works for me: Import["..._Output.hdf"] Import["..._Output.hdf", "HDF"] Import["..._Output.hdf", "HDF5"] $\endgroup$ Apr 24, 2019 at 8:42
  • $\begingroup$ Try Import["your_file.h5", "Data"] or replace "Data" with "Elements" to see what is available. It works for me with AceFEM 6.921. $\endgroup$
    – Pinti
    Apr 24, 2019 at 9:02
  • $\begingroup$ For me too, on the .h5 file. However it doesn't work on the .hdf file, which contains my actual post processing data if I'm not mistaken. $\endgroup$ Apr 24, 2019 at 9:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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