strange artifacts seen in 3D view of script-animated high-poly mesh

Question

I'm animating a high-poly mesh using an algorithm which is supposed to be challenging to visualize. See sect. 8.1 of this paper. I'm using a key_block method to animate the distortion of a high-poly icosphere, and the image shown below looks weird. Is this an artifact of the 3D view rendering, or is there something just bad about the mesh (which could be quite true)?

Since the application will be technical, and getting surfaces in the right place is important (there will be other things in close proximity), just going to fewer polygons or similar things, are not the kind of answer I can use.

The image was made using n_subdivisions = 7, but I have lowered it to 6 in the script to start with. It is a little resource-heavy and takes about a minute.

Yes the script is a little long, but it's the easiest way to allow others to generate the problem.

def get_ico_data(n_subdiv=3, size=1.0, location=(0,0,1.0)):
    """get icosphere data from ops primitive because 
    It's easier than figuring it all out myself"""
    make_ico = bpy.ops.mesh.primitive_ico_sphere_add
    make_ico(subdivisions=n_subdiv, size=size, location=location)
    ico = bpy.context.active_object
    ico_data = ico.data.copy()
    bpy.ops.object.delete()
    return ico_data

def make_new_meshobject_from_data(data, name=None, scene=None):
    if scene == None:
        scene = bpy.context.scene
    if name == None:
        name = "bob"
    me = bpy.data.meshes.new(name)
    ob = bpy.data.objects.new(name, me)
    ob.data = data.copy()
    scene.objects.link(ob)
    ob.select = True
    return ob

def velo1(X, t, T):  # written for convenience, not speed
    """from here: http://geuz.org/gmsh/doc/preprints/gmsh_visu_preprint.pdf"""
    pi, twopi = np.pi, 2.0*np.pi
    g = np.cos(pi * t / T)    
    x, y, z = X.T
    vx = 2.*np.sin(   pi*x)**2 * np.sin(twopi*y)    * np.sin(twopi*z)    * g
    vy =   -np.sin(twopi*x)    * np.sin(   pi*y)**2 * np.sin(twopi*z)    * g
    vz =   -np.sin(twopi*x)    * np.sin(twopi*y)    * np.sin(   pi*z)**2 * g    
    V = np.vstack((vx, vy, vz)).T
    return V

def RK4a(x, t, n, h, F, T):
    for i in range(n):  # written for readability, not speed
        ho2, ho6 = h/2.0, h/6.0
        k1 = F(x[i]         ,  t      , T)
        k2 = F(x[i] + k1*ho2,  t + ho2, T)
        k3 = F(x[i] + k2*ho2,  t + ho2, T)
        k4 = F(x[i] + k3*h  ,  t + h  , T)
        x[i+1] = x[i] + ho6 * (k1 + 2.*(k2 + k3) + k4)
        t += h

import bpy
import numpy as np

pi, twopi = np.pi, 2.0*np.pi
radius, center = 0.15, [0.35, 0.35, 0.35]

ico_data = get_ico_data(n_subdiv=6, size=1.0, location=(0,0,0))
print(len(ico_data.vertices))
for v in ico_data.vertices:
    for i in range(3):
        v.co[i] = radius*v.co[i] + center[i]

X0 = np.array([v.co for v in ico_data.vertices])
print("X0.shape: ", X0.shape)
print("X0: ", X0[:,2].max(), X0[:,2].min())

T = 3.0  # how far to go

n, dt = 601, 0.01
t = 0.

# data using RK4
s1, s2 = X0.shape
X_RK = np.zeros((n+2, s1, s2))
X_RK[0] = X0
RK4a(X_RK, t, n, dt, velo1, T)
ddata_RK = X_RK[:n+1]

X = X0.copy()
data = []
data.append(X.copy())
for i in range(n):
    X += dt * velo1(X, t, T)
    data.append(X.copy())
    t += dt
ddata = np.array(data)

np.save("make_my_dayta7",    ddata   )
np.save("make_my_dayta_RK7", ddata_RK)

ico = make_new_meshobject_from_data(ico_data)
ve = ico.data.vertices
ico.scale = [6.0]*3
ico.location = (0,0,0)
print("ve: ", max([v.co[2] for v in ve]), min([v.co[2] for v in ve]))

ico_RK = make_new_meshobject_from_data(ico_data)
ico_RK.scale = [6.0]*3
ico_RK.location = (6,0,0)

dayta    = np.load("make_my_dayta7.npy"   )
dayta_RK = np.load("make_my_dayta_RK7.npy")

n_frames = dayta.shape[0]
print("hey dayta.shape: ", dayta.shape)

bpy.context.scene.frame_end = n_frames

# From here: https://blender.stackexchange.com/a/36915/5334
## ------- part 2  ---set up keyblocks / shape_keys

for i_frame in range(n_frames):
    block    = ico.shape_key_add(name=str(i_frame), from_mix=False)  # returns a key_blocks member
    block_RK = ico_RK.shape_key_add(name=str(i_frame), from_mix=False)  # returns a key_blocks member

    block.value, block_RK.value = 1.0, 1.0
    block.mute,  block_RK.mute  = True, True 

    for (vert, co) in zip(block.data, dayta[i_frame]):
        vert.co = co

    for (vert, co) in zip(block_RK.data, dayta_RK[i_frame]):
        vert.co = co

    # keyframe off on frame zero
    block.mute, block_RK.mute = True, True
    block.keyframe_insert(data_path='mute', frame=0, index=-1)
    block_RK.keyframe_insert(data_path='mute', frame=0, index=-1)

    block.mute, block_RK.mute = False, False
    block.keyframe_insert(data_path='mute', frame=i_frame + 1, index=-1)
    block_RK.keyframe_insert(data_path='mute', frame=i_frame + 1, index=-1)

    block.mute, block_RK.mute = True, True
    block.keyframe_insert(data_path='mute', frame=i_frame + 2, index=-1)
    block_RK.keyframe_insert(data_path='mute', frame=i_frame + 2, index=-1)

Answer 1

OK, after running the script and looking the mesh, I can say that this is a clear case of overlapping faces, created due to heavy stretching of triangles on top of each other, especially in the center area:

You can probably avoid this if you somehow generate more geometry in these heavily stretched areas (perhaps calculate each face's area in each iteration and subdivide any face that's above a certain area?). This will make sure you don't get huge highly stretched tris overlapping each other.

It will be trickier to reduce density in areas that don't need it, though, since in your animation some stretched areas will later have hardly any stretch at all.