<jupyter_start><jupyter_text>Meshed up swarmsUsing the shadow particle fanciness to achieve great things<jupyter_code>import underworld as uw
from underworld import function as fn

import glucifer

uw.matplotlib_inline()

# import matplotlib.pyplot as pyplot
# import matplotlib.pylab as pylab
import numpy as np
import math
import os



workdir = os.path.abspath(".")
outputPath = os.path.join(workdir,"MeshSwarms/")

if uw.rank() == 0:
    if not os.path.exists(outputPath):
        os.makedirs(outputPath)
    
uw.barrier()<jupyter_output><empty_output><jupyter_text>Create mesh and finite element variables------Note: the use of a pressure-sensitive rheology suggests that it is important to use a Q2/dQ1 element<jupyter_code>minX  =  0.0;   minY = 0.0
maxX  =  1.0;   maxY = 1.0


resX = 32
resY = 32

elementType="Q2/dPc1"  # This is enough for a test but not to use the code in anger

mesh = uw.mesh.FeMesh_Cartesian( elementType = (elementType), 
                                 elementRes  = ( resX, resY), 
                                 minCoord    = ( minX, minY), 
                                 maxCoord    = ( maxX, maxY),
                                 periodic    = [False, False]  ) 


velocityField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=mesh.dim )
pressureField    = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )

velocityField.data[:]   = [0.,0.]
pressureField.data[:]   = 0.

# Define strain rate tensor

strainRateFn = fn.tensor.symmetric( velocityField.fn_gradient )
strainRate_2ndInvariantFn = fn.tensor.second_invariant(strainRateFn)<jupyter_output><empty_output><jupyter_text>Boundary conditionsPure shear with moving  side walls — all boundaries are zero traction with outflow top and bottomto accommodate changing volume<jupyter_code>iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]

leftI  = mesh.specialSets["MinI_VertexSet"]
rightI = mesh.specialSets["MaxI_VertexSet"]
bothIs = leftI + rightI

baseWall   = mesh.specialSets["MinJ_VertexSet"]
topWall    = mesh.specialSets["MaxJ_VertexSet"]

allWalls = iWalls + jWalls

velocityBCs = uw.conditions.DirichletCondition( variable        = velocityField, 
                                                indexSetsPerDof = (iWalls+topWall, jWalls) )

# for index in topWall:
#     velocityField.data[index,0] = 1.0

velocityField.data[topWall.data,0] = np.sin(np.pi * mesh.data[topWall.data,0])
# Already set to zero everywhere else 
stokesPIC = uw.systems.Stokes( velocityField  = velocityField, 
                               pressureField  = pressureField,
                               conditions     = [velocityBCs,],
                               fn_viscosity   = 1.0, 
                               fn_bodyforce   = (0.0,0.0) )

# And a suitable solver package is already attached to it

solver = uw.systems.Solver( stokesPIC )

solver.solve()
import unsupported.triangulation.meshswarm2D 
reload(unsupported.triangulation.meshswarm2D)

mswarm = unsupported.triangulation.meshswarm2D.meshSwarm2D(mesh, velocityField)
mswarmPoints = mesh.data.copy()
mswarm.lagrSwarm.add_particles_with_coordinates(mswarmPoints)

mswarm.update_triangulation()
trackVar = unsupported.triangulation.meshswarm2D.meshSwarmVariable2D(mswarm, dataType='float', count=1)
trackVar.lagrangianData.data[:,0] = (1.0-trackVar.lagrangianData.swarm.particleCoordinates.data[:,1])

## Create tracking variables

# import numpy.ma as ma
# If we want to add all the non-boundary nodes to the swarm:
# masked = ma.array(trackVariable.data)
# masked.mask = ma.make_mask_none(masked.shape)
# masked.mask[mesh.specialSets["AllWalls_VertexSet"].data] = True
# non_bound_pts = masked.compressed()

# import matplotlib.pyplot as plt


# figure = plt.figure(figsize=(7,6))
# ax = figure.add_subplot(111)

# s = mswarm
# points = mswarm.lagrSwarm.particleCoordinates.data
# m = ax.tripcolor(points[:,0], points[:,1], trackVar.lagrangianData.data.reshape(-1), alpha=0.6)

# ax.triplot(points[:,0], points[:,1], color="Black", 
#            alpha=0.4, zorder=100)

# plt.colorbar(mappable=m)


glfig = glucifer.Figure(quality=3)

glfig.append( glucifer.objects.Mesh(mesh=mesh, opacity=0.5))
glfig.append( glucifer.objects.Points(swarm=mswarm.lagrSwarm, 
                                      fn_colour=trackVar.lagrangianData,
                                      pointSize=5.0))

#time = 0.0
tri = mswarm.triangulation

for timestep in range(0,1000):
    dt = mswarm.get_max_dt()
    mswarm.integrate(dt)
            
    if timestep % 100 == 0:
        glfig.save_image(filename="SwarmImage"+str(timestep)+".png")
        
    pts = mswarm.prune_swarm_min_length(0.005)
    with mswarm.deform_swarm():
        mswarm.lagrSwarm.particleCoordinates.data[pts,0] = -1.0e99
    print "Deleted ", len(pts), "points -> ", mswarm.triangulation.npoints
       
    pts = mswarm.refine_swarm_max_length(0.025)
    newVals = trackVar.evaluate(pts)
    newpts = mswarm.lagrSwarm.add_particles_with_coordinates(pts)
    mypts  = np.where(newpts != -1)[0] 
    trackVar.lagrangianData.data[newpts[mypts],0] = newVals[mypts]

    print "Added ", len(pts), "points -> ", mswarm.triangulation.npoints
    mswarm.update_triangulation()
      
    pts = mswarm.refine_swarm_too_narrow(5.0)
    newVals = trackVar.evaluate(pts)
    newpts = mswarm.lagrSwarm.add_particles_with_coordinates(pts)
    mypts  = np.where(newpts != -1)[0] 
    trackVar.lagrangianData.data[newpts[mypts],0] = newVals[mypts]

    print "Added ", len(pts), "points -> ", mswarm.triangulation.npoints
    mswarm.update_triangulation()

    # "Boundary" conditions
    
    top_pts = np.where(mswarm.lagrSwarm.particleCoordinates.data[:,1] == 1.0)
    trackVar.lagrangianData.data[top_pts] = 0.0
    
    base_pts = np.where(mswarm.lagrSwarm.particleCoordinates.data[:,1] == 0.0)
    trackVar.lagrangianData.data[base_pts] = 1.0
   
print timestep
mswarm.update_triangulation()

figure = plt.figure(figsize=(7,6))
ax = figure.add_subplot(111)

s = mswarm
points = mswarm.lagrSwarm.particleCoordinates.data
m = ax.tripcolor(points[:,0], points[:,1], trackVar.lagrangianData.data.reshape(-1), alpha=0.6)

ax.triplot(points[:,0], points[:,1], color="Black", 
           alpha=0.4, zorder=100)

plt.colorbar(mappable=m)
trackVar.lagrangianData.data.min()
# figure.savefig("testmesh.png", dpi=600)







Ffn = fn.math.sin( np.pi * fn.coord()[0]) + fn.math.sin( np.pi * fn.coord()[1])
Fxfn  = np.pi * fn.math.cos( np.pi * fn.coord()[0]) * fn.math.sin( np.pi * fn.coord()[1])
Fyfn  = np.pi * fn.math.sin( np.pi * fn.coord()[0]) * fn.math.cos( np.pi * fn.coord()[1])
Fxxfn = -np.pi**2 * fn.math.sin( np.pi * fn.coord()[0]) * fn.math.sin( np.pi * fn.coord()[1])
Fyyfn = -np.pi**2 * fn.math.sin( np.pi * fn.coord()[0]) * fn.math.sin( np.pi * fn.coord()[1])
F = np.sin(np.pi * mswarm.triangulation.points[:,0]) * \
    np.sin(np.pi * mswarm.triangulation.points[:,1])
    
FxA   = np.pi * np.cos(np.pi * mswarm.triangulation.points[:,0]) * \
                np.sin(np.pi * mswarm.triangulation.points[:,1])

Fx2A  = -np.pi**2 * np.sin(np.pi * mswarm.triangulation.points[:,0]) * \
                    np.sin(np.pi * mswarm.triangulation.points[:,1])


FyA   = -np.pi * np.sin(np.pi * mswarm.triangulation.points[:,0]) * \
                 np.cos(np.pi * mswarm.triangulation.points[:,1])

Fy2A  = -np.pi**2 * np.sin(np.pi * mswarm.triangulation.points[:,0]) * \
                   np.sin(np.pi * mswarm.triangulation.points[:,1])

    
    
DX, DY, DXX, DXY, DYY   = mswarm.triangulation_compute_gradient(F, second=True)
DXX2, DYX2 = mswarm.triangulation_compute_gradient(DX)
DXY2, DYY2 = mswarm.triangulation_compute_gradient(DY)
FxA.max(), DX.max(), (DXX2-DXX).max()

figure = plt.figure(figsize=(10,10))
ax = figure.add_subplot(111)

s = mswarm

points = mswarm.interpolator.points

m = ax.tripcolor(points[:,0], points[:,1], DY-FyA, alpha=0.6)

ax.triplot(points[:,0], points[:,1], color="Black", 
           alpha=0.4, zorder=100)

plt.colorbar(mappable=m)

# figure.savefig("meshGrad.png", dpi=300)
del2fn = gradXVar.fn_gradient[0] + gradYVar.fn_gradient[1]



gptVariable.data[:,0] = trackVar.evaluate(gaussSwarm.particleCoordinates.data)

fig1 = glucifer.Figure( figsize=(600,600) )
fig1.append( glucifer.objects.Mesh(mesh=imesh, opacity=0.25))
fig1.append( glucifer.objects.VectorArrows(mesh=mesh, fn=velocityField, scaling=0.1 ))
fig1.append( glucifer.objects.Surface(mesh=imesh, fn=Fxxfn ))

# fig1.append( glucifer.objects.Points(swarm=gaussSwarm, fn_colour=gptVariable,
#                                      colours="Blue Red", pointSize=0.0))

fig1.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14120608765304954, 'fraction_numerical': 0.023229202425906856, 'mean_word_length': 3.4342973110096398, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '436786', 'license': 'MIT', 'n_tokens_mistral': 3344, 'n_tokens_neox': 3191, 'n_words': 525, 'path': 'notebooks/__work_in_progress/models_inprogress/MeshSwarm.ipynb', 'repo_name': 'lmoresi/uw-docker-landing-page', 'size': '8739'}