Density plot

[manarshhjot]

Hello,

I want to have a plot of the material density of the elements. Can someone please guide me how to do it.

[mikerife]

Hi @manarshhjot I think that using PyVista here would be the most straight-forward way to get this type of plot. We can use mesh.grid to create mesh that PyVista can plot, and have an array of density for each element that can be used as cell data to be plotted. See this:

https://docs.pyvista.org/user-guide/what-is-a-mesh.html

A normally 'hard' part would be to gather the actual density values for each element. I noticed from your last discussion that your model may already be set up to gather this easily as it had logic to define a Material ID for each element, each with a unique density. Are you still using that model with respect to this question?

Mike

[akaszynski]

If possible, if you could upload a minimum working example we might be able to generate a simple script to pull out material values. I could probably make one, just busy with development...

[mikerife]

Hi @manarshhjot & @akaszynski Here is an example with two meshed blocks having densities equal to 2 and 4. The Jupyter Lab notebook is attached. But here is the important code...I'm guessing that someone more experienced with Python than I am (I'm pretty new) can probably make this more efficient for a large number of elements and/or material models.

mats = mapdl.mesh.material_type
mats_unique = np.unique(mats)

dens = []

for material in np.unique(mats):
    dens.append(mapdl.get_value('dens',material,'temp',0))

res = dict(zip(mats_unique, dens))

def replace(list1, dictionary):
    return [item if not item in dictionary else dictionary[item] for item in list1]

mats_dens = replace(mats, res)

mesh = mapdl.mesh.grid

mesh.cell_data['density'] = mats_dens
mesh.plot(scalars='density', show_edges=True)

Mike
PlotDensities.zip

p.s. I have updated this from my original posting. The original used looping and assumed the elements used a list of material IDs with no gaps in the material ID numbers. This new method will account for any numbering gaps in the materials. Still not sure if it is optimal.

[akaszynski]

@germa89 threw down the gauntlet saying that he didn't think we could do better without resorting to C. Read it and weep @germa89.

First, let's make the mesh bigger by changing the element size in the above example to mapdl.esize(0.01). With that you get 3 million elements. That's better.

Now, here's the timings. The top is from @mikerife's implementation using a dictionary as a lookup table. Bottom is mine using an array as a lookup table. Effectively I am resorting to C, numpy doesn't incur a GIL lock for each individual lookup so there's a much smaller penalty in the array lookup.

Still, even with 3M elements, the total time is just ~300 milliseconds. Not as good as ~7 milliseconds, but it's not bad and the majority of the time on a 3M element mesh will be in the solution, not post-processing.

>>> mapdl.mesh
ANSYS Mesh
  Number of Nodes:              3080702
  Number of Elements:           3000000
  Number of Element Types:      1
  Number of Node Components:    0
  Number of Element Components: 0

>>> timeit pretty_fast(mats)
333 ms ± 892 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

>>> timeit wicked_fast(mats)
6.59 ms ± 50.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

mats = mapdl.mesh.material_type

def pretty_fast(mats):
    mats = mapdl.mesh.material_type
    mats_unique = np.unique(mats)

    dens = []
    for material in np.unique(mats):
        dens.append(mapdl.get_value('dens',material,'temp',0))

    res = dict(zip(mats_unique, dens))


    def replace(lst, dictionary):
        return [item if not item in dictionary else dictionary[item] for item in lst]

    mats_dens = replace(mats, res)
    return mats_dens


def wicked_fast(mats):
    import re

    # get all material numbers from the output from mplist
    output = mapdl.mplist()
    mats_unique = np.array(re.findall('MATERIAL NUMBER\s+([[0-9]+)', output), dtype=int)
    dens = [mapdl.get_value('dens', mat, 'temp', 0) for mat in mats_unique]

    # create a trivial lookup table
    # this will be faster than a dictionary lookup and take up (a tiny bit)
    # more memory. Typical memory vs speed tradeoff, but the table is so
    # small. FYI, tables just like these are used within MAPDL.
    dens_lookup = np.empty(mats_unique.max() + 1)
    dens_lookup[mats_unique] = dens

    mats_dens = dens_lookup[mats]
    return mats_dens

[germa89]

Nice benchmark @akaszynski !!

[mikerife]

@akaszynski Two orders of magnitude improvement! Wow, that's great.

[mikerife]

Hi All
The answer above with the example Jupyter Lab zipped file will work for densities defined by the MP command. The TB command can also be used especially for creating density (single value), and as a function of a field variable like time or x/y/z coordinates. I don't see a way of using *GET to get the field variable values, so we would need to know those...which if defined in a PyMAPDL script should not be too difficult! Mike

[germa89]

Are those material properties stored in the RST file? I guess they are.

We could resource to parse some command output (TBLIST??)