3
$\begingroup$

How can I get value of stress (or some other field) in each integration point of finite element with AceFEM package? I need this information to show mesh with values of stress averaged over each element as shown on image bellow (source of example: Krysl 2015, pg 360).

average stress

AceFEM documentation shows only examples of already recovered stress that is continuous over elements (extrapolated to nodes). How can I create image similar to above on the following simple example of 3D analysis?

<< AceFEM`;

(* A example function. Its argument is the number of elements on one edge. *)
Clear[example]
example[n_Integer] := Module[{},
  SMTInputData[];
  SMTAddDomain["test", "OL:SED3H1DFHYH1NeoHooke", {}];
  SMTAddMesh[
   Hexahedron[{{0, 0, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 0}, {1, 1, 0}, {1, 1, 1}, {1, 0, 1}}],
     "test", "H1", {n, n, n}
   ];
  SMTAddEssentialBoundary[{ "X" == 0 &, 1 -> 0, 2 -> 0, 3 -> 0}];
  SMTAddNaturalBoundary[ Line[{{1, 0, 0}, {1, 0, 1}}], 3 -> Line[{-100}]];
  SMTAnalysis[];
  SMTNextStep["\[Lambda]" -> 1];
  While[SMTConvergence[10^-8, 10], SMTNewtonIteration[];];
  ]

example[4]
SMTShowMesh["BoundaryConditions" -> True, "DeformedMesh" -> True,"Field" -> "Sxx", 
  ImageSize -> 200, Axes -> True,AxesLabel -> {"X", "Y", "Z"}, Ticks -> None, "Legend" -> False
]

continuous stress

$\endgroup$
4
$\begingroup$

There is supposed to be option in the SMSShowMesh, where you can define the "Field" as a discontinuous field. We can obtain discontinuous field by using SMTPost[_, "Smooth" -> False] command, but the mesh output is not what one would expect:

pp=SMTPost["Sxx", "Smooth" -> False];
SMTShowMesh["Field" -> (Table[Mean[#], {i, #}] & /@pp),
ImageSize -> 260, Axes -> True]

enter image description here

However there is a simple but slower workaround, where we have to call SMTShowMesh for each element separately and then we can use Show to merge all elements together:

p1 = SMTShowMesh["BoundaryConditions" -> True, "DeformedMesh" -> True,
   "Field" -> "Sxx", ImageSize -> 260, Axes -> True, 
  AxesLabel -> {"X", "Y", "Z"}, Ticks -> None, "Legend" -> True, 
  "FillElements" -> False];
range = MinMax[Mean /@ pp];
Show[
 p1 /. (GraphicsComplex[__] -> {}),
 MapIndexed[
  SMTShowMesh["DeformedMesh" -> True, "Field" -> Table[#, SMTNoNodes],
    "ZoomElements" -> #2, 
    "Contour" -> {range[[1]], range[[2]], 20}] &, Mean/@ pp
  ], PlotRange -> All
 ]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Interesting solution that keeps the original layout of SMTShowMesh. BTW, you could avoid using some pure functions in your code and make it more elegant, e.g. range = MinMax[Mean /@ pp]. $\endgroup$ – Pinti Apr 21 '17 at 9:07
  • $\begingroup$ I have edited the answer with your suggestions. I didnt knew before that MinMax function even exists. The first way seems to be a bug or unfinished option, because it sets the colors of elements to nodes, so we get the white color because there are more noodes than elements. For unsmooth mesh, you have to have the NoElements*NoNodesPerElement number of nodes, and specify the values at each node of element seperately and also interpolate them seperately. $\endgroup$ – BHudobivnik Apr 21 '17 at 9:32
3
$\begingroup$

After @BHudobivnik showed the undocumented function SMTPost[_,"Smooth"->False] that gives values of field in each integration point, I came up with alternative way to plot the mesh.

(* integration point values *)
igPts = SMTPost["Sxx", "Smooth" -> False];

(* average value per element *)
elPts = Mean /@ igPts;

(* Rescaling to interval {0,2/3} in neccesseary because Hue ranges from Red to Red *)
colors = Hue /@ Rescale[elPts, Reverse@MinMax[elPts], {0, 2/3}];

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

This method of showing the mesh is fast, but it doesn't preserve some nice features of SMTShowMesh, such as showing boundary condition symbols.

Graphics3D[
 Thread[{FaceForm /@ colors, Hexahedron /@ nodesPerElement}],
 ImageSize -> 200
 ]

mesh

$\endgroup$

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.