A simple plot helper for visualizing FElupe results (2D case)

[ZAARAOUI999]

Hello Mr. Andreas, I think that it can be a good idea if FElupe has a plot helper, this can be useful so that any user can visualize some data directly with his python IDE without the need to export his results to another file and use another software to look what he had, at least for some cases. The code below uses the matplotlib library to visualize results.

The code:

import felupe as fem
import numpy as np

import matplotlib.pyplot as plt
from matplotlib import rc, tri, collections
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

#3D Mesh case (not yet implemented)
# m1 = fem.Cube(b=(10,10,10), n=10)
# r1 = fem.RegionHexahedron(m1)

#2D Mesh case
#get some data to plot
m1 = fem.Rectangle(b=(10,10), n=20)
r1 = fem.RegionQuad(m1)
f1 = fem.FieldsMixed(r1, n=1)
mat1 = fem.NeoHooke(mu=1, bulk=1)
s1 = fem.SolidBody(mat1, f1)
bounds, _ = fem.dof.uniaxial(f1,right=10, clamped=True)
move = fem.math.linsteps([0,10],10)
ramp = {bounds['move']:move}
step = fem.Step([s1], ramp, bounds)
job = fem.Job([step])
job.evaluate(verbose=False)

###################### -Plot helper- ############################
def quatplot(u, v, cells, ax=None, **kwargs): #stackoverflow
    """plot meshes"""
    if not ax: ax = plt.gca()
    uv = np.c_[u,v]
    verts= uv[cells]
    pc = collections.PolyCollection(verts, **kwargs)
    ax.add_collection(pc)
    ax.autoscale()

def plot(field, values, axes=['x','y','z'], label:str='', component:int=0,
         offset:float=0., show_mesh:bool=True, deformed:bool=True,
         show_min_max:bool=True):
    """
    A simple plot helper.
    """
    
    region = field.region
    points = region.mesh.copy().points
    cells = region.mesh.cells
    dim = points.shape[-1]
    
    #get displacements
    du_ = field.fields[0].extract(grad=False)
    du = fem.project(du_, region)
    if deformed:
        #add dispalacements to the reference points coordinates
        points[:,:dim]+=du[:,:dim]
        
    if dim==3:
        # split points array into separated axes
        u, v, w = points[:,0], points[:,1], points[:,2]
        #make a 3d figure
        ax = plt.figure().add_subplot(projection='3d')
        #setting plot labels
        ax.set_xlabel(axes[0]+'$\longrightarrow$')
        ax.set_ylabel(axes[1]+'$\longrightarrow$')
        ax.set_zlabel(axes[2]+'$\longrightarrow$')
        #get specified component from the given data
        values = values[:,component]
        triang = tri.Triangulation(u, v)  # Delaunay triangulation of all points.
        #remove unwanted triangles
        mask = tri.TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=.1,rescale=False)
        triang.set_mask(mask)
        #plot interpolated data 
        pcm = ax.tricontourf(triang, values, cmap='turbo', levels=128, zdir='z', offset=offset)
    else:
        # split points array into separated axes
        u, v = points[:,0], points[:,1]
        #make a 2d figure
        ax = plt.figure().add_subplot()
        #setting plot labels
        ax.set_xlabel(axes[0]+'$\longrightarrow$')
        ax.set_ylabel(axes[1]+'$\longrightarrow$')
        #get specified component from the given data
        values = values[:,component]
        triang = tri.Triangulation(u, v)  # Delaunay triangulation of all points.
        #remove unwanted triangles
        mask = tri.TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=.1,rescale=False)
        triang.set_mask(mask)
        #plot interpolated data 
        pcm = ax.tricontourf(triang, values, cmap='turbo', levels=128)
        if show_mesh:
            #plot mesh edges
            quatplot(u, v, cells, ax=ax, color='k', facecolor='none')
    
    #add colorbar
    c_axis = inset_axes(ax, width='5%', height='50%',loc='center left', borderpad=-10)
    c = plt.gcf().colorbar(pcm, cax=c_axis, pad=1)
    c.ax.yaxis.set_ticks_position('right')
    #set title
    if show_min_max:
        label += '\nmin: %.1e\nmax: %.1e\n'%(values.min(), values.max()) 
    c.ax.set_title(x = 0, label=label,verticalalignment='top',
                   horizontalalignment='left', fontproperties={'size':10})

#################################################################

Some examples of results:

[adtzlr]

Hi, thanks for your efforts. What about creating your own additional package on PyPI with stuff like that, e.g. feluplot or felupetools? I fear that I do not have the time right now (writing a PhD thesis, having a full-time job, my family and this hobby project) to go through the code, generalize it here and there, ensure that matplotlib is not a must-dependency and write the docs. You could simply publish your own package with small enhancements around felupe. What about that?

[adtzlr]

But your code is fantastic, the result looks amazing! I'll think on what to do with it.

[adtzlr]

I always had in mind to keep the core of FElupe as small as possible. Well, sort of 😄.

[ZAARAOUI999]

😁 Okey then FEluPlot looks like a good choice. Thanks.

[ZAARAOUI999]

FEplot is available now as a PyPI package, I hope that it will be helpful.

[adtzlr]

Great! I'll add a list with community extensions to the README of felupe.