3
$\begingroup$

I have been struggling with this problem for a long time. I know that this question is not specifically related to AceGen or AceFEM, but I would like to know if anyone had a similar experience in postprocessing the results.

I have exported the results of my field in undeformed and deformed configurations. First, I will start with undeformed configuration data. The undeformed data can be found here. It is constructed in the form {x,y,z,f}, where f is the field values. The intention here is to plot the data using RegionPlot3D.

First, I import the data and use the interpolation command to construct an interpolation out of it

   dataUnDef = Import["DataUnDef.csv"];
intdataUnDef = Interpolation[dataUnDef];

Then by using the following command I readily plot the function in the specified range.

RegionPlot3D[
   intdataUnDef[x,y,z]>= 0.52, {x,0,25}, {y,0,25},{z,0,25}, Mesh -> None, 
   PlotStyle -> Red,PlotPoints->100,BoundaryStyle ->None, Boxed ->False, 
   Axes -> False]

The plot is somehow the same as I observe if I use SMTShowMesh to plot the field.

The field plotted in undeformed state

Now in the next step, I would like to use the same procedure to plot the values in the deformed configuration. The deformed data can be found here. And it is in the form {x,y,z,f}. Like the undeformed case, I import the data first and construct an interpolation function out of it. However, since the mesh is irregular now, thus, I am obliged to use the Interpolation order of one.

dataDef = Import["DataDef.csv"];
intdataDef=Interpolation[dataDef,InterpolationOrder->1];

Next, I apply the interpolation function to the coordinates of my regular mesh (which corresponds to the coordinates of my undeformed mesh). In this case, I can have the values of my field in the deformed configuration in a regular mesh and then I can use the interpolation function of higher order to construct a nice interpolation function out of it.

First I define the coordinates

Xnew = Table[
    Table[Table[
      x, {x, 0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {i, 
      0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {j, 0, 
     25 + 0.3012048192771104`*3, 0.3012048192771104`}] // Flatten;
Ynew = Table[
    Table[Table[
      i, {x, 0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {i, 
      0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {j, 0, 
     25 + 0.3012048192771104`*3, 0.3012048192771104`}] // Flatten;
Znew = Table[
    Table[Table[
      j, {x, 0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {i, 
      0, 25 + 0.3012048192771104`, 0.3012048192771104`}], {j, 0, 
     25 + 0.3012048192771104`*3, 0.3012048192771104`}] // Flatten;

NewCoord = Transpose@{Xnew, Ynew, Znew};

Note that I have added two additional layers to my grid in the Z direction in order to cover any pile-ups in the deformed configuration of my domain.

Then, I apply the interpolation function to it.

newDataDef = Apply[intdataDef, NewCoord, 1];

And finally, I use the interpolation function (which can have an order of 3 now) and plot the data

        newDataDef2 = 
          Insert[NewCoord // Transpose, newDataDef// Flatten, 4] // 
           Transpose;
    newintDataDef = Interpolation[newDataDef2];
RegionPlot3D[
 newintDataDef[x,y,z] >= 0.52, {x, 0, 25}, {y, 0, 25}, {z,
   0, 25.6024}, Mesh -> None, PlotStyle -> Red, PlotPoints -> 100, 
 BoundaryStyle -> None, Boxed -> False, Axes -> False]

The deformed configuration.

Now, this deformed configuration is not what I was looking for. First, it is showing some spurious values of my field on the top surface of the domain. Secondly, the deformed shape is not that obvious here. I was wondering if anyone has some suggestions on how to improve the plot of the deformed configuration.

Thanks in advance

$\endgroup$
  • $\begingroup$ If I understand correctly, you would like to plot a 3D contour over deformed configuration (based on FEM mesh)? $\endgroup$ – Pinti Feb 5 at 14:37
  • $\begingroup$ @Pinti Exactly! But The above problem exists. $\endgroup$ – Msen Rezaee Feb 5 at 14:38
  • $\begingroup$ Unrelated to the core of the question, but regular coordinate arrays can be elegantly created with CoordinateBoundingBoxArray[{{0, 0, 0}, {25, 25, 25}}, 5 {1, 1, 1}]; instead of fiddling with nested Table. $\endgroup$ – Pinti Feb 6 at 8:45
3
$\begingroup$

Finite element solution is based on interpolating function and it is the best to use it also for visualization, since it is the "true" solution.

Lets start with an example from AceFEM documentation "Mixed 3D Solid FE, Elimination of Local Unknowns", because it represents a nice 3D problem with large deformations.

Get["AceFEM`"]

SMTInputData[];
SMTAddDomain["A", "ExamplesHypersolid3D", {"E *" -> 1000., "\[Nu] *" -> .3}];
SMTAddEssentialBoundary[{ "X" == 0 &, 1 -> 0, 2 -> 0, 3 -> 0}, { "X" == 10 &, 3 -> -1}];
SMTAddMesh[
  Hexahedron[{{0, 0, 0}, {10, 0, 0}, {10, 2, 0}, {0, 2, 0}, {0, 0, 3}, {10, 0, 2}, {10, 2, 2}, {0, 2, 3}}],
   "A", "H1", {15, 6, 6}
];
SMTAnalysis[];

SMTNextStep["λ" -> 1];
While[
  While[
   (step = SMTConvergence[10^-8, 10, {"Adaptive BC", 8, .001, 1, 5}]), 
   SMTNewtonIteration[];
  ];
  SMTStatusReport[];
  If[step[[4]] === "MinBound", SMTStatusReport["Analyze"]; SMTStepBack[]];
  step[[3]] 
  ,
  If[step[[1]], SMTStepBack[]];
  SMTNextStep["Δλ" -> step[[2]] ]
];

Show[
 SMTShowMesh["FillElements" -> False],
 SMTShowMesh["DeformedMesh" -> True, "Mesh" -> False, "Field" -> Map[Norm, SMTNodeData["at"]]]
 ]

AceFEM_mesh

crds = SMTNodeData["X"];
displacements = SMTNodeData["at"];
connectivity = SMTElementData["Nodes"];
field = SMTPostData["Exx"]; (* Choose field (values at nodes) here. *)

Up to here we just generated some data to work with and procedure doesn't need AceFEM` any more. We will use utility functions from NDSolve`FEM` context to create ElementMesh object with information about node coordinates and element connectivity. Take care about node ordering because for some element topology AceFEM representation and ElementMesh are not directly compatible.

mesh = ToElementMesh[
  "Coordinates" -> crds,
  "MeshElements" -> {HexahedronElement[connectivity]}
  ]

Create also deformed representation of ElementMesh by using convenient function ElementMeshDeformation and observe the results.

meshDef = ElementMeshDeformation[mesh, Transpose@displacements]
meshDef["Wireframe"["MeshElementStyle" -> FaceForm@LightBlue]]

Deformed_ElementMesh

Now comes the crucial part. We create InterpolatingFunction objects based on ElementMesh, with special options to "ignore" values outside the solution domain. There will be two functions, one for undeformed configuration and one for deformed, because they are so different.

int = ElementMeshInterpolation[{mesh}, field,
  "ExtrapolationHandler" -> {Function[Indeterminate], "WarningMessage" -> False}
  ]
intDef = ElementMeshInterpolation[{defMesh}, field,
  "ExtrapolationHandler" -> {Function[Indeterminate], "WarningMessage" -> False}
  ]

This InterpolatingFunction object can be easily used in ContourPlot3D and related functions. The following is a function that combines contour plot of chosen FEM solution with mesh representation. Small details are explained with inline comments.

MeshContourPlot3D // ClearAll
MeshContourPlot3D[int_InterpolatingFunction, opts : OptionsPattern[]] := Module[
  {mesh, plotRange, contourPlot, x, y, z},

  (* Extract ElementMesh from solution. *)
  mesh = int["ElementMesh"];
  (* Automatically determine plot range from mesh size and expand it for 5% for nicer plot. *)
  plotRange = MinMax[#, Scaled[0.05]] & /@ mesh["Bounds"];

  contourPlot = ContourPlot3D[
    int[x, y, z],
    Evaluate[Sequence @@ MapThread[Prepend, {plotRange, {x, y, z}}]],
    Evaluate@FilterRules[{opts}, Options@ContourPlot3D],
    PerformanceGoal -> "Speed",
    Mesh -> None,
    (* Default contour value is 0, but can be changed through options. *)
    Contours -> {0}
    ];

  Show[
   mesh["Wireframe"[
     Sequence @@ FilterRules[{opts}, Options@ElementMeshWireframe3D]
     ]],
   contourPlot
   ]
  ]

This is undeformed configuration with default options.

MeshContourPlot3D[int]

Plot1

And deformed configuration with some styling.

MeshContourPlot3D[intDef,
 "MeshElementStyle" -> EdgeForm@Blue,
 Contours -> {0},
 ContourStyle -> Red,
 PerformanceGoal -> "Quality"
]

Plot2

$\endgroup$
  • $\begingroup$ +1 Thanks for your comprehensive explanation. It takes time for me to apply it myself. $\endgroup$ – Msen Rezaee Feb 6 at 8:53

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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