Using Ansys pyMAPDL,

pyDPF-Post and more

Part 1: Running and post-processing an ANSYS
model outside of ANSYS

Alex Grishin, PhD

We Make Innovation Work
[Link]
PADT, Inc
Preliminaries: What is pyANSYS?
• In this article, we’ll explore a relatively new ANSYS API for python formerly called
pyANSYS (see the link here: [Link] )
• It’s main use is to give external programs written in Python access to ANSYS model
information
• There are hooks to both launch ANSYS (at least MAPDL –through pyMAPDL) and to
access model information from previous runs
• The main framework for accessing model information from previous runs is called
pyDPF (Data Processing Framework). This is a general purpose (physics-agnostic) tool
for accessing the data. A more specific module for structural and thermal models is
called pyDPF-Post.
• In this article, we’ll focus on pyMAPDL and the pyDPF-Post

We Make Innovation Work

[Link]
 Preliminaries: What are the benefits of such an interface?

• Most serious ANSYS users would benefit from learning at least the basics of how the DPF-Post and pyMAPDL
work for a couple reasons:

1. To extend their current simulation capabilities. Although we won’t dive deeply into it, in this article
we’ll explore how users may export a model’s system equations and post-process a model outside of
Ansys These two tasks may seem unrelated, but a future article will tie them together (and show users
how to solve the system equations outside of ANSYS!)
2. To extend their current post-processing capabilities. This is actually the main ‘killer-app’ of the DPF-
Core product. Although ANSYS users have access to the ACT customization process automation tools,
that framework is limited to what’s currently available in the .NET version that ships with ANSYS. To
take a simple example: Suppose users want to perform a custom spectrum analysis (perhaps an
industry-specific mode-combination technique) using the eigenvectors from an ANSYS modal solution.
In such a case, the ACT framework is not ideal as it has no native capability for manipulating large
matrix systems*

• Read more about the DPF-Core and view a quick example here:
[Link]

*at least none that are exposed to users

We Make Innovation Work
[Link]
 Preliminaries: Installing. From Scratch
• pyDPF-Core, pyDPF-Post and pyMADPL have a complex web of dependencies.
• The best and easiest way to ensure that everything goes smoothly is to allow these installations to install
everything they need ‘from scratch’ –that is to say, outside and independently of any package manager
• This is true even for users who already have a powerful package manager like Anaconda (if you have a
standard, default installation). There should be no reason you can’t have both an external python
installation AND Anaconda. And we’re suggesting this is going to be the cleanest solution!*
• When you visit the two github download links (here for DPF-Core, and here for pyMADPL), you’ll see the
following version support:

• So, we suggest installing Python 3.8 if you dont’ already have it. And again: We mean a
‘standalone’ installation. Users may download an installer for Windows here *veteran python users should also
know that the following useful
modules ship WITH the pyANSYS
products:
• numpy
• scipy
• pyVista
We Make Innovation Work • matplotlib
[Link]
• and more...
 Preliminaries: Installing. From Scratch
• After downloading the installer, right-click on it and “Run as administrator” (you really should have
administrative privileges to install software)

• Check “Add Python 3.8 to PATH” –but only if you don’t already
have a path variable set for another installation (or you’re ok
with overriding it. Or you simply specify the Scripts folder each
time you want to launch an editor).
• An alternative (if you insist on having multiple versions of
python installed) is to set custom environment variables –one
for each installation. Then execute the environment variable
corresponding to the installation you want to access at the
command prompt
• Click “Install Now”
We Make Innovation Work
[Link]
 Preliminaries: Installing. From Scratch
• Now, launch a command prompt as adminstrator (right-click on Command Prompt from the main menu and
select Run as administrator)
• If everything went smoothly so far, you should be able to type “python –version” at the prompt and see the
following:

• Install pymadl by typing the following at the prompt (followed by <enter>)

• Install pypdf-post by typing the following at the prompt

We Make Innovation Work

[Link]
 Preliminaries: Installing. From Scratch

• Try a simple test of the installation at the command prompt. Open a python shell by typing “python”
• Then, at the prompt, type “from ansys-mapdl-core import launch_mapdl”
• If successful, you wont’ see any errors or warnings (like below)
• However, some users may get an error related to the version of Google’s protobuf installed
• If you fall into this category, uninstall protobuf, and then install protobuf version 3.20* (first type “pip
uninstall protobuf <enter>” followed by “pip install protobuf == 3.20 <enter>”

• From here, you can

install an editor of
your choice (or use
the command shell)
• We’ll use the
standard python
shell

*WARNING: This is true as of this writing

on 7/13/2022, but will probably not be
true a year from now)
We Make Innovation Work
[Link]
The Example Problem: 2D NAFEMS benchmark problem
• We chose the following problem for demonstration
• This is a NAFEMS “Challenge Problem” from 2015 (too late to submit our solution!). The
link may be found here: [Link]
• We thought it might help users understand the problem more intimately by looking at
the solution in different frameworks (that post-process element results slightly
differently) and understanding how the post-processing is done
• In a future article, even more insight will be gained by actually looking at the underlying
system equations

We Make Innovation Work

[Link]
 The Example Problem: In Workbench...
• The challenge problem asks us to build both a single linear element model, as well as a single quadratic
element model. The “challenge” is to explain the differences and make statements about which is
“correct”. These models are represented in our Workbench project in systems A and B respectively
• In addtion, we’ve added another model for completeness in system C: A 64-element quadratic model.
We should be able to spot a trend as to what’s going on as we go from left to right..

We Make Innovation Work

[Link]
 The Example Problem: What we’ll be doing...

• As we mentioned at the outset, we chose a simple model as a benchmark to demonstrate pyMAPDL

and pydpf-post’s capabilities
• The setup is this: We built the model in Workbench, and so we’ll be accessing the project’s solution
folder to open the model in pyMAPDL and export the matrices.
• We’ll also run the model in pyMAPDL and view the results using pydpf-post
• The python files will reside in a python project folder
dp0\SYS\MECH
[Link],
[Link], Single Linear
[Link]
WB Project [Link]
Rectangle
[Link] Model
python project folder [Link]

.py files,
dp0\SYS-1\MECH
Single [Link],
[Link],
Quadratic [Link]
Rectangle [Link]
[Link]
Model [Link]

We Make Innovation Work

[Link]
 The Example Problem: Single Linear Rectangle Model
• In a challenges like this, the authors intentionally omit details. For example, how should we constrain
the model? Well, that’s part of the challenge
• Without exploiting the symmetry (or anti-symmetry in this case), the mininal constraint condition
suggests that we can only apply constraints where the applied loads are zero (locations 1 and 2 below)
• And since we don’t know how to apply this sort of combined non-uniform traction in Workbench, we
script it in*

• We’ll fix
uy here
1

2
• We’ll fix
ux and uy
here

*The details of this are beyond the scope of this article.

Contact the author if you have questions...
We Make Innovation Work
[Link]
 The Example Problem: Single Linear Rectangle Model
file: [Link]
from [Link] import launch_mapdl
• We’ve already run this model in Workbench, but from [Link] import post
before we show the result, we’ll do the following: import numpy as np
import sys,os
• Re-read the model in the solution folder and run it
using pyMAPDL def runandsolve(mapdl,path):
#mapdl = launch_mapdl(run_location=rundir,override=True)
• Export the matrices using MADPL’s using APDL [Link]()
Math functions specifically designed for this [Link]()
[Link](path)
purpose (see Eric Miller’s introductory article on [Link]('/solu')
APDL Math here: [Link]('yes')
[Link]()
[Link]
math-access-to-the-ansys-solver-matrices-with- def matrix_export(mapdl,kname='kmat',fname='fmat'):
[Link]()
apdl/ (again, we won’t be using the exported [Link]('*smat,vks,d,import,full,[Link],stiff')
matrices in this article –that will have to wait for a [Link]('*dmat,f,d,import,full,[Link],rhs')
[Link]('*smat,usr2solv,d,import,full,[Link],usr2solv')
future article) [Link]('*dim,snums,,usr2solv_rowdim')
• View the results in pypdf-post [Link]('*vfill,snums(1),ramp,1,1')
• to make things a little more compact, we’ve [Link]('*vec,vsnums,d,import,apdl,snums')
bundled most of the necessary pyMADPL code [Link]('*mult,usr2solv,tran,vsnums,,vmap')
runstr = '*export,vks,mmf,' + kname + '.mmf'
into a python file called [Link] shown at right [Link](runstr)
runstr = '*export,f,mmf,' + fname + '.mmf'
[Link](runstr)
We Make Innovation Work
[Link]('*export,vmap,csv,[Link]')
[Link]
 The Example Problem: Single Linear Rectangle Model
• Ok. Now, all we have to do is open a python shell and execute the following lines of
python (see next slide)

file: [Link]
import os
import sys
• Change this path to wpath = r"C:\Users\[Link]\exporting_matrices" • Specify the python working
your location before [Link](wpath) directory
execution! [Link](wpath)
• import [Link]
from kmat import *

solverfilesloc = wpath + "\\plane_stress_files\\dp0\\SYS\\MECH"

dspath = solverfilesloc + "\\[Link]"
resultpath = solverfilesloc + "\\[Link]"
• Export a stiffness matrix called k1 mapdl = launch_mapdl(run_location=solverfilesloc,override=True)
runandsolve(mapdl,dspath)
• This function is defined in [Link]
and a load vector called f1
matrix_export(mapdl,kname='k1',fname='f1') • This function is defined in [Link]

solution = post.load_solution(resultpath)
disp = [Link]()
stress = [Link]() • Plot ux
disp.x.plot_contour(cpos='xy',cmap='rainbow',n_colors=9) • Plot sx
[Link].plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
[Link]() • close the MAPDL connection
We Make Innovation Work
[Link]
 The Example Problem: Single Linear Rectangle Model
• From the Start Menu, open a python shell and simply cut-and-paste the
lines from the previous slide into the shell
• You’ll see the displacement plot first pop up in a ‘blocking’ window
• You’ll have to explicitly close this window before you can see the stress
result, which will pop up the same way
• Just make sure to close both windows to execute all the code
• Note that this behavior is a feature of the python shell we’re using to execute this
code. Different editors will have behave differently when confronted with an
interactive blocking process σx
ux

We Make Innovation Work

[Link]
The Example Problem: Single Linear Rectangle Model

• Just verify that these plots are the same as what’s shown in Workbench
• This verifies the NAFEMS solution for the linear rectangle case: A null result. Dont’
worry, there’s nothing wrong here. The authors of the challenge problem
specifically designed it to produce a null result over a linear rectangle!
ux σx

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

file: [Link]
• Following the NAFEMS challenge
import os
problem directions, next verify the
import sys
single quadratic element by wpath = r"C:\Users\[Link]\exporting_matrices"
cutting-and-pasting the code at [Link](wpath)
right into a new python shell [Link](wpath)
• Notice, this code is nearly identical from kmat import *
to the code we previously
executed for the linear rectangle, solverfilesloc = wpath + "\\plane_stress_files\\dp0\\SYS-1\\MECH"
but this time we’re pointing to a dspath = solverfilesloc + "\\[Link]"
resultpath = solverfilesloc + "\\[Link]"
different folder for the [Link] and
mapdl = launch_mapdl(run_location=solverfilesloc,override=True)
results file (see slide 10) runandsolve(mapdl,dspath)
matrix_export(mapdl,kname='k2',fname='f2')

solution = post.load_solution(resultpath)
disp = [Link]()
stress = [Link]()
disp.x.plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
[Link].plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
We Make Innovation Work
[Link]()
[Link]
 The Example Problem: Single Quadratic Rectangle Model
• And again: compare

ux (Workbench) σx (Workbench)

ux (dpf-post) σx (dpf-post)

We Make Innovation Work

[Link]
The Example Problem: Single Quadratic Rectangle Model
• Everything looks identical except for the displacement plot. Why?
• Well, first verify that the nodal values ARE the same
• And this line tells us what nodes those
• This line in DPF-Post returns a numpy array of ux values... values correspond to...

ux
ux
• They ARE the same
values
• It’s just that the DPF-
post plot seems to
just be interpolating
the 4 corner node
values (or something.
We’ll get to the
We Make Innovation Work bottom of this!)
[Link]
The Example Problem: Single Quadratic Rectangle Model
• It’s worth spending a little time contemplating what’s going on here
• To help with this, we’ll first introduce a VERY useful tool which came to us silently,
without charge and without any mention when we downloaded the pyANSYS
modules: pyVista
• All the mesh-based contour plots we see in pyMAPDL, pyDPF-Post, and DPF-Core
are provided using this tool ‘under the hood’
• To prove this, type the following in the same open python session (we want the
‘mapdl’ and ‘solution’ objects to be instantiated before doing this)
import pyvista as pv <enter>
[Link] <enter>
[Link] <enter>

• You should see something like this:

• So, both the [Link] and

the [Link] use a pyvista
core [Link])

We Make Innovation Work

[Link]
The Example Problem: Single Quadratic Rectangle Model
• Now, copy the grid to a new unstructured grid called ‘copygrid’ by typing:
copygrid = [Link]() <enter>
• Then type copygrid <enter>
• This is what you should see:

• This makes sense, right? We should be looking at a single quadrilateral rectangle, right?
• Note that if you copied the [Link] instead, you would see 5 elements. This is
because that mesh differs in that it also contains the surface effect elements (and these are
just line elements coincident with the quadrilateral edges) which are used to apply the
surface tractions...

We Make Innovation Work

[Link]
The Example Problem: Single Quadratic Rectangle Model
• Ok. Now, we’re going to use pyVista to create a useful and pretty plot of this mesh
and it’s nodes. Just cut-and-paste the following into your open python session:
nodes = [Link]([Link]) • create a pv PolyData object to render the nodes
plabels = [str(i+1) for i in range(myquad.n_points)] • create node labels
p = [Link]() • create a Plotter instance
p.add_mesh(copygrid,color='brown') • add the mesh to the plotter
p.add_point_labels(nodes,plabels,font_size=20, • add the node (PolyData) object
point_color='gold',
point_size=40,render_points_as_spheres=True,
always_visible=True) • set the camera for 2D viewing
p.camera_position='xy' • show the mesh
[Link]()

• Here’s what you

should see!

We Make Innovation Work

[Link]
The Example Problem: Single Quadratic Rectangle Model
• Returning to our problem. From slide 18, we learned the following about the nodal
displacement x-component values:

• According to
Workbench, the solution
should be as shown as
shown at the right:

• This solution can be

verified in pyDPF-Post
by creating a numpy
array to store the data
like so:
ux = [Link](8) • Initialize a numpy array to store the ux values
ids = [Link](disp.x.get_scoping_at_field())
• store the nodal ids (like in slide 18)
ux[ids-1] = disp.x.get_data_at_field()
• get the solution but store it ine correct
• pyvista arrays are zero-based, while ANSYS numbering is always order
ones-based (a FORTRAN legacy) so we have to subtract 1
We Make Innovation Work
[Link]
 The Example Problem: Single Quadratic Rectangle Model

• We’ll use our copied pyVista mesh to plot these values!

• This is remarkably easy to do with pyVista. All we have to do is create a ‘point_data’ array with
these values. And pyVista supplies us with a shortcut for doing so. Such arrays are managed by
pyVista within python dictionaries.

copygrid['ux'] = ux

• This line adds the result data to the grid for plotting.
Note the dictionary-like syntax
• Plot the displacement like so:

p = [Link]()
p.add_mesh(copygrid,scalars='ux',
cmap='rainbow',n_colors=9)
p.camera_position='xy'
• This is what you should see!
[Link]()
• This is identical to the pyDPF-Post
plot
We Make Innovation Work
[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Ok. We’ve verified the pyDPF solution. But we still neeed to explain why it differs from the WB
solution
• Here’s where things get really interesting (for readers interested in the guts of FEA)
• We menionted (on slide 18) that pyDPF-Post seems to be only interpolating the corner nodes,
and so it does –even for the plot we just created OUTSIDE of both ANSYS and pyDPF-Post where
we were quite careful to map all the nodal results to their correct place.
• We can verify that pyDPF-Post is doing that by querying the mesh’s connectivity:

• Remember that in pyVista (and python generally), arrays are indexed by zero. So when we query this
mesh’s node connectivity, we expect to see (according to the ANSYS convention and mesh we see in
Workbench):
• [0,1,2,3,4,6,7,5]
• What’s happening is that pyDPF-Post is copying the mid-side nodes to lie on top of the corner nodes!

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• We can use pyVista to generate a contour plot like the one shown in Workbench
• Hopefully, doing so will reveal much about how modern FEA software works AND the NAFEMS
challenge problem
• First, let’s create a better displacement plot. We’ll build on what we’ve learned so far and create a
quadratic quadrilateral mesh ‘from scratch’. This time, we’ll properly define the midside nodes
• pyVista is a high-level python wrapper around a larger and older framework called VTK
• When we define a mesh (‘grid’ in VTK parlance) in pyVista, we’re using the names and conventions
of VTK. So, for example, to define a quadratic quadrilateral, we first have to look that up in VTK
(see the link here).
• The type we’re looking for is called VTK_QUADRATIC_QUAD (id = 23)
• To create this grid type in pyVista with the nodes we already have, simply cut-and-paste the
following into the python session:

import vtk
cells = [Link]([8,0,1,2,3,4,6,7,5])
cell_type = [Link]([vtk.VTK_QUADRATIC_QUAD])
myquad = [Link](cells,cell_type,[Link])

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Now add the nodal displacement data as before:

myquad['ux'] = ux

• And plot the displacement:

p = [Link]()
p.add_mesh(myquad,scalars='ux',cmap='
rainbow',n_colors=9)
• Compare!
p.camera_position='xy'
[Link]()

• This is what we get!

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Well, that’s MUCH better. At least values along edges look right.
• It looks like the Workbench solution still has an advantage, though.
• Notice that we can query the center of the mesh, even though we DON’T HAVE A NODE THERE
• Remember, a finite element solution is smooth throughout it’s domain (sometimes called ‘the
support’). Though we don’t have nodes within that domain, we should still be able to query values
ANYWHERE (via the shape functions). This fact is leveraged when creating finite element contour
plots (even if users can’t query any location they want. In a future article, we’ll show how DPF-Post
and pyVista will allow us to overcome obstacles like that)

• To capture this additional information, we’ll use a hgher order

grid called VTK_BIQUADRATIC_QUAD (id=28)
• add a point at the center
import vtk
eightpts = [Link]([Link])
ninepts = [Link]((eightpts,[Link]([0.5,0.5,0.])))
cells = [Link]([9,0,1,2,3,4,6,7,5,8])
cell_type = [Link]([vtk.VTK_BIQUADRATIC_QUAD])
mybiquad = [Link](cells,cell_type,ninepts)

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Now, append the center ‘ux’ value from Workbench to our ux solution and add it the mesh:

ux = [Link](ux,0.5)
mybiquad['ux'] = ux
• And plot as before...
p = [Link]()
p.add_mesh(mybiquad,scalars='ux',cmap
='rainbow',n_colors=9)
p.camera_position='xy'
[Link]()

• And compare...
• This is as far as we can go with a
single cell in pyVista. We can verify
that every node has the correct value
• Differences in the two plots are due
to different plotting algorithms. We
will return to this in a future article

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Let’s return again to the challenge problem. Recall that we get a null solution (both displacements
and stresses) with the linear rectangle. We’re NOT getting a null solution now
• As we mentioned, this problem was designed to produce a null solution for a linear rectangle, but
the exact solution is only zero over a region of the rectangle as reflected in the (exact) component
stress solutions shown below (we’ve already verified σx)

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• Let’s follow the author and reproduce some of the “low fidelity” stress results he reports.
• In particular, we want to show how CS1 and CS2 below can come about with the solution post-
processing techniques we’ve just seen

We Make Innovation Work

[Link]
 The Example Problem: Single Quadratic Rectangle Model

• In the same python session we’ve been using (we need the ‘stress’ object active), execute the following:
stress.von_mises.plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
• This should result in the following:

• This is identical to
the plot the
NAFEMS author
calls CS1 (and it’s
also what you get in
MAPDL)

• Compare to the Workbench von Mises stress

plot of the quadratic rectangle (above). This
is identical to what the NAFEMS author calls
CS2
We Make Innovation Work
[Link]
The Example Problem: Single Quadratic Rectangle Model

• Based on the work we’ve already done, we can easily check to see if the quadratic or
bi-quadratic pyVista cell can produce solutions similar to the Workbench result.
• We’ll re-use the pyVista grids we’ve already created
• Just as we obtained the ux result, stored it in an array with correct node ordering and
plotted it over the pyVista cell vertices, we can do the same with the von Mises stress
• Try the following:
seqv = [Link](8)
ids = [Link](stress.von_mises.get_scoping_at_field())
seqv[ids-1] = stress.von_mises.get_data_at_field()
myquad['seqv'] = seqv
mybiquad['seqv'] =[Link](seqv,1.0) • add the WB stress at center
• To plot the eight-node • To plot the nine-node
quadrilateral von Mises bi-quadrilateral von
stress: Mises stress:
p = [Link]() p = [Link]()
p.add_mesh(myquad,scalars='seqv', p.add_mesh(mybiquad,scalars='seqv',
cmap='rainbow',n_colors=9) cmap='rainbow',n_colors=9)
p.camera_position='xy' p.camera_position='xy'
[Link]() [Link]()
We Make Innovation Work
[Link]
The Example Problem: Single Quadratic Rectangle Model
• And this is what results

• quadratic grid • bi-quadratic

grid

• Compare to the Workbench (CS2) stress result

• Again, we can be confident that all node values match
• Only the contour algorithm differs
• A more thorough discussion will have to wait for another article.
• Suffice it so say that what we’ve been exploring are different ways
to post-process a result. pyVista offers one set of assumptions,
while WB is using a different graphical tesselation
We Make Innovation Work
[Link]
The Example Problem: Conclusions
• so, what can we say about this model?
• In slide 29, we showed the exact solution, but let’s
forget we saw that.
• We can’t believe the single linear element because
loads are only applied at corner nodes in such an
element. And the applied tractions all sum to zero at
corner nodes (although the stress estimate of the
center of that element yields the correct result due
to the careful planning of problem’s author!)
• But let’s look again at what the author calls CS1. We
KNOW why that pattern emerges in pyDPF-Post. But
it happes in MAPDL too. This is because the MAPDL
convention is to extrapolate stress values from
integration points to the corner nodes only (this fact
may explain why pyDPF-Post simply copies corner
nodes when considering midside node values)
• MAPDL stores no stress values for midside nodes,
while in pyDPF-Post, midside node values are
averaged with corner node values at those corners
We Make Innovation Work
[Link]
The Example Problem: Conclusions

• We defintely believe the CS2 models as a trend (at least at the edges if not the center).
It’s not really clear if the Workbench plotting algorithm adds any more information than
the one we used with pyVista (botom left. In other words: is the Workbench solution
‘closer’ to the exact result?)
σeqv σ eqv

pyVista CS2: Workbench

We Make Innovation Work

[Link]
 The Example Problem: Conclusions
• We conclude by comparing the algorithms on an 8 x 8 mesh of quadratic elements (the third
Workbench model in the project)
• First, let’s look at the pyDPF-Post plot. Run the following...
import os
import sys file: [Link]
wpath = r"C:\Users\[Link]\exporting_matrices"
[Link](wpath)
[Link](wpath)
from kmat import *

solverfilesloc = wpath + "\\plane_stress_files\\dp0\\SYS-2\\MECH"

dspath = solverfilesloc + "\\[Link]" • This is what you should see...
resultpath = solverfilesloc + "\\[Link]" σeqv
mapdl = launch_mapdl(run_location=solverfilesloc,override=True)
runandsolve(mapdl,dspath)
matrix_export(mapdl,kname='k3',fname='f3’)el = el

solution = post.load_solution(resultpath)
disp = [Link]()
stress = [Link]()
#commenting the displacement plots because we don't need them for demo...
#disp.x.plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
stress.von_mises.plot_contour(cpos='xy',cmap='rainbow',n_colors=9)
[Link]()

We Make Innovation Work

[Link]
 The Example Problem: Conclusions
• Now, we’ll create our quadratic pyVista solution. This time, we’ll use a slightly different syntax (it’s easier
just to read in an array of cell connectivities than to create the conectivity format we used for a single cell)
• First, we’ll create a cell connectivity array like so:

cells = [Link]((0,8),int) #initialize

for i in range([Link].n_cells):
el = [Link].element_by_id(i+1)
carray = [Link]([Link]).astype(int)
cells = [Link](cells,[carray],axis=0)

• And here’s the new syntax:

newquad = [Link]({vtk.VTK_QUADRATIC_QUAD:cells},[Link])

• Now, get the ids and seqv data as before...

seqv = [Link](newquad.n_points)
ids = [Link](stress.von_mises.get_scoping_at_field())
seqv[ids-1] = stress.von_mises.get_data_at_field()
newquad['seqv'] = seqv
We Make Innovation Work
[Link]
 The Example Problem: Conclusions
• Finally, plot as before (but this time, set the ‘show_edges’ option to ‘yes’)
p = [Link]()
p.add_mesh(newquad,scalars='seqv',
cmap='rainbow',n_colors=9, • show edges to revealt the
show_edges='yes') VTK_QUADRILATERAL_QUAD cell edges
p.camera_position='xy'
[Link]()

We Make Innovation Work

[Link]
 The Example Problem: Conclusions

• Finaly, comparing all three algorithms on an 8 x 8 mesh of quadratic elements reveals that the
Workbench algorithm (CS2) and our pyVista solution (middle figure) are both smoother (converge
faster) than the CS1 solution (DPF-Post and MAPDL)
σeqv σeqv σeqv

CS1: pyDPF-Post on 8 x 8 quadratic grid pyVista on 8 x 8 quadratic grid CS2: Workbench

(highlighting the unique edge with 8 x 8 quadratic
pattern of this pyvista cell type) mesh

We Make Innovation Work

[Link]
The Example Problem: Conclusions

• One last thing: though the NAFEMS author neglects displacements, we plot the
distorted displacement of the refined model just to show readers how this contrived
model behaves (compare to single quadratic element on slide 28)

ux uy

We Make Innovation Work

[Link]
 Final Notes
• The author of the NAFEMS challenge problem points out that different commercial FE packages may employ
different strategies for graphically rendering nodal solutions
• The differences betweeen them may be especially accute for coarse, high-order element types (quadratic and
above)
• The author makes a compelling case that there’s nothing to be gained by what he calls ‘low fidelity’
approaches (see slides 30 and 31 –although that wasn’t the intent)
• While good engineering practice would guide most analysts to refine a mesh before making conclusions or
identifying trends, a little care in rendering results with more accuracy could potentially reduce the number of
such refinements an analyst has to make
• The central challenge posed by the author was to provide an estimate for the von Mises stress at the center fo
the plate for both the single linear and quadratic rectangles and choose which one is ‘correct’. This was a
misdirection from the issue of mesh refinement and the relation of each model to the exact solution. The
single linear rectangle ‘accidentally’ yields the exact solution (from the perspective of analyst refining the mesh
–not the problem’s author). This is misleading because that solution is only exact in that the loads cancel to
within numerical accuracy
• Issues of stress contour rendering (which consumed this article) were only raised by the author in his solution
• We’ll return to this problem and explore the issue of contour rendering in a little more detail in a future article

We Make Innovation Work

[Link]