Search code examples
pythonmatplotlibcontourf

matplotlib tricontourf ploblem when I give more data point

I have got a problem when i try to plot the stress.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import matplotlib.cm as cm

def plot(x_plot, y_plot, a_plot):
    x = np.array(x_plot)
    y = np.array(y_plot)
    a = np.array(a_plot)

    triang = mtri.Triangulation(x, y)
    refiner = mtri.UniformTriRefiner(triang)
    tri_refi, z_test_refi = refiner.refine_field(a, subdiv=4)

    plt.figure(figsize=(18, 9))
    plt.gca().set_aspect('equal')
    #     levels = np.arange(23.4, 23.7, 0.025)
    levels = np.linspace(a.min(), a.max(), num=1000)
    cmap = cm.get_cmap(name='jet')
    plt.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
    plt.scatter(x, y, c=a, cmap=cmap)
    plt.colorbar()

    plt.title('stress plot')
    plt.show()

First I have plot by using only 8 points:

x = [2.3384750000000003, 3.671702, 0.3356813, 3.325298666666667, 2.660479, 1.3271675666666667, 1.6680919666666665, 0.6659845666666667]
y = [0.614176, 0.5590579999999999, 0.663329, 0.24002166666666666, 0.26821433333333333, 0.31229233333333334, 0.6367503333333334, 0.3250663333333333]
a = [2.572, 0.8214, 5.689, -0.8214, -2.572, -4.292, 4.292, -5.689]
plot(x, y, a)

version a

then i try to give an information of the bound of rectangle:

x = [2.3384750000000003, 1.983549, 3.018193, 2.013683, 3.671702, 3.978008, 4.018905, 0.3356813, 0.0, 0.0, 1.0070439, 3.325298666666667, 2.979695, 2.660479, 1.3271675666666667, 0.9909098, 1.6680919666666665, 0.6659845666666667]
y = [0.614176, -0.038322, 0.922264, 0.958586, 0.5590579999999999, -0.1229, 0.87781, 0.663329, 1.0, 0.0, 0.989987, 0.24002166666666666, -0.079299, 0.26821433333333333, 0.31229233333333334, -0.014787999999999999, 0.6367503333333334, 0.3250663333333333]
a = [2.572, 2.572, 2.572, 2.572, 0.8214, 0.8214, 0.8214, 5.689, 5.689, 5.689, 5.689, -0.8214, -0.8214, -2.572, -4.292, -4.292, 4.292, -5.689]
plot(x, y, a)

version b

I don't know how to fix it and why this happen. The figure that i want is:

photoshoped version

I have do the scatter plot of each point in second figure and there are correct but why the color is not contour.

Thank you very much.

Solution

  • The field returned by the UniformTriRefiner does not interpolate well in the case of the additional points. Instead it introduces new minima and maxima with values up to 20 times larger than the original points.

    The plot below shows what is happening.

    import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import matplotlib.cm as cm

def plot(x_plot, y_plot, a_plot, ax=None):
    if ax == None: ax = plt.gca()
    x = np.array(x_plot)
    y = np.array(y_plot)
    a = np.array(a_plot)

    triang = mtri.Triangulation(x, y)
    refiner = mtri.UniformTriRefiner(triang)
    tri_refi, z_test_refi = refiner.refine_field(a, subdiv=2)

    levels = np.linspace(z_test_refi.min(), z_test_refi.max(), num=100)
    cmap = cm.get_cmap(name='jet')

    tric = ax.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
    ax.scatter(x, y, c=a, cmap=cmap, vmin= z_test_refi.min(),vmax= z_test_refi.max())
    fig.colorbar(tric, ax=ax)

    ax.set_title('stress plot')

fig, (ax, ax2) = plt.subplots(nrows=2, sharey=True,sharex=True, subplot_kw={"aspect":"equal"} )
    
x = [2.3384750000000003, 3.671702, 0.3356813, 3.325298666666667, 2.660479, 1.3271675666666667, 1.6680919666666665, 0.6659845666666667]
y = [0.614176, 0.5590579999999999, 0.663329, 0.24002166666666666, 0.26821433333333333, 0.31229233333333334, 0.6367503333333334, 0.3250663333333333]
a = [2.572, 0.8214, 5.689, -0.8214, -2.572, -4.292, 4.292, -5.689]
plot(x, y, a, ax)


x = [2.3384750000000003, 1.983549, 3.018193, 2.013683, 3.671702, 3.978008, 4.018905, 0.3356813, 0.0, 0.0, 1.0070439, 3.325298666666667, 2.979695, 2.660479, 1.3271675666666667, 0.9909098, 1.6680919666666665, 0.6659845666666667]
y = [0.614176, -0.038322, 0.922264, 0.958586, 0.5590579999999999, -0.1229, 0.87781, 0.663329, 1.0, 0.0, 0.989987, 0.24002166666666666, -0.079299, 0.26821433333333333, 0.31229233333333334, -0.014787999999999999, 0.6367503333333334, 0.3250663333333333]
a = [2.572, 2.572, 2.572, 2.572, 0.8214, 0.8214, 0.8214, 5.689, 5.689, 5.689, 5.689, -0.8214, -0.8214, -2.572, -4.292, -4.292, 4.292, -5.689]
plot(x, y, a, ax2)

plt.show()

    enter image description here

    As can be seen, the values of the "interpolated" field overshoot the original values by a large amount.
    The reason is that by default, UniformTriRefiner.refine_field uses a cubic interpolation (a CubicTriInterpolator). The documentation states

    The interpolation is based on a Clough-Tocher subdivision scheme of the triangulation mesh (to make it clearer, each triangle of the grid will be divided in 3 child-triangles, and on each child triangle the interpolated function is a cubic polynomial of the 2 coordinates). This technique originates from FEM (Finite Element Method) analysis; the element used is a reduced Hsieh-Clough-Tocher (HCT) element. Its shape functions are described in 1. The assembled function is guaranteed to be C1-smooth, i.e. it is continuous and its first derivatives are also continuous (this is easy to show inside the triangles but is also true when crossing the edges).

    In the default case (kind ='min_E'), the interpolant minimizes a curvature energy on the functional space generated by the HCT element shape functions - with imposed values but arbitrary derivatives at each node.

    While this is truely very technical, I emphazised some important bits, namely that the interpolation is smooth and continuous with a defined derivative. To guarantee this behaviour, overshoots are inevitable when the data is very sparse but with large amplitude fluctuations.

    Here, the data is simply not suitable for cubic interpolation. One would either try to aquire denser data, or use a linear interpolation instead.