Commit a2d43b17 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

copy/paste error corrected

parent f40defd5
Pipeline #6952 passed with stage
in 1 minute and 21 seconds
......@@ -329,7 +329,7 @@
name="Script"
command="SetScript"
number_of_elements="1"
default_values="import vtk
import numpy as np
from migflow import scontact

p = scontact.ParticleProblem(2)
all_ts = self.GetOutputInformation(0).Get(vtk.vtkStreamingDemandDrivenPipeline.TIME_STEPS())
ts = self.GetOutputInformation(0).Get(vtk.vtkStreamingDemandDrivenPipeline.UPDATE_TIME_STEP())
it = all_ts.index(ts)
f = inputs[0] 
p.read_vtk(directory,it)
xmin = np.amin(f.Points[:,0])
xmax = np.amax(f.Points[:,0])
ymin = np.amin(f.Points[:,1])
ymax = np.amax(f.Points[:,1])
zmin = np.amin(f.Points[:,2])
zmax = np.amax(f.Points[:,2])
step = np.minimum(np.minimum(np.abs(xmax-xmin)/2,step),np.abs(ymax-ymin)/2)
if zmax-zmin == 0:
 z = [0]
else:
 step = np.minimum(step,np.abs(zmax-zmin)/2)
 z = np.arange(zmin,zmax+(zmax-zmin)*1e-8,step)
x = np.arange(xmin,xmax+(xmax-xmin)*1e-8,step)
y = np.arange(ymin,ymax+(ymax-ymin)*1e-8,step)
nx = x.size
ny = y.size
x,y,z = np.meshgrid(x,y,z)
points = np.column_stack([x.flat,y.flat,z.flat])
stress = p.computeStressTensor(points,radius)
tensor = np.zeros((stress.shape[0],6))
tensor[:,0] = stress[:,0]
tensor[:,1] = stress[:,3]
tensor[:,3] = stress[:,1]
a = stress[:,0]
d = stress[:,3]
b = stress[:,1]
de = ((a-d)**2+4*b**2)**0.5
l1 = (a+d+de)/2 
l2 = (a+d-de)/2
output.SetPoints(points)
ncells = (nx-1)*(ny-1)
types = np.full([ncells],9,dtype=np.uint8)
location = np.arange(0,5*ncells,5,np.uint32)
cells = np.zeros((nx-1,ny-1,5),dtype=np.uint32)
ix = np.arange(0,nx-1,dtype=int)
iy = np.arange(0,ny-1,dtype=int)
cells[:,:,0] = 4
cells[:,:,1] = ix[:,None]+iy[None,:]*nx
cells[:,:,2] = (ix[:,None]+1)+iy[None,:]*nx
cells[:,:,3] = (ix[:,None]+1)+(iy[None,:]+1)*nx
cells[:,:,4] = ix[:,None]+(iy[None,:]+1)*nx
output.SetCells(types,location,cells.reshape([-1]))
output.PointData.append(tensor,"tensor")
lmax = np.maximum(np.abs(l1),np.abs(l2))
lmin = np.minimum(np.abs(l1),np.abs(l2))
output.PointData.append(lmax,"lmax")
output.PointData.append(lmax,"lmin")
output.PointData.append(l1,"l1")
output.PointData.append(l2,"l2")
"
default_values="import vtk
import numpy as np
from migflow import scontact

p = scontact.ParticleProblem(2)
all_ts = self.GetOutputInformation(0).Get(vtk.vtkStreamingDemandDrivenPipeline.TIME_STEPS())
ts = self.GetOutputInformation(0).Get(vtk.vtkStreamingDemandDrivenPipeline.UPDATE_TIME_STEP())
it = all_ts.index(ts)
f = inputs[0] 
p.read_vtk(directory,it)
xmin = np.amin(f.Points[:,0])
xmax = np.amax(f.Points[:,0])
ymin = np.amin(f.Points[:,1])
ymax = np.amax(f.Points[:,1])
zmin = np.amin(f.Points[:,2])
zmax = np.amax(f.Points[:,2])
step = np.minimum(np.minimum(np.abs(xmax-xmin)/2,step),np.abs(ymax-ymin)/2)
if zmax-zmin == 0:
 z = [0]
else:
 step = np.minimum(step,np.abs(zmax-zmin)/2)
 z = np.arange(zmin,zmax+(zmax-zmin)*1e-8,step)
x = np.arange(xmin,xmax+(xmax-xmin)*1e-8,step)
y = np.arange(ymin,ymax+(ymax-ymin)*1e-8,step)
nx = x.size
ny = y.size
x,y,z = np.meshgrid(x,y,z)
points = np.column_stack([x.flat,y.flat,z.flat])
stress = p.computeStressTensor(points,radius)
tensor = np.zeros((stress.shape[0],6))
tensor[:,0] = stress[:,0]
tensor[:,1] = stress[:,3]
tensor[:,3] = stress[:,1]
a = stress[:,0]
d = stress[:,3]
b = stress[:,1]
de = ((a-d)**2+4*b**2)**0.5
l1 = (a+d+de)/2 
l2 = (a+d-de)/2
output.SetPoints(points)
ncells = (nx-1)*(ny-1)
types = np.full([ncells],9,dtype=np.uint8)
location = np.arange(0,5*ncells,5,np.uint32)
cells = np.zeros((nx-1,ny-1,5),dtype=np.uint32)
ix = np.arange(0,nx-1,dtype=int)
iy = np.arange(0,ny-1,dtype=int)
cells[:,:,0] = 4
cells[:,:,1] = ix[:,None]+iy[None,:]*nx
cells[:,:,2] = (ix[:,None]+1)+iy[None,:]*nx
cells[:,:,3] = (ix[:,None]+1)+(iy[None,:]+1)*nx
cells[:,:,4] = ix[:,None]+(iy[None,:]+1)*nx
output.SetCells(types,location,cells.reshape([-1]))
output.PointData.append(tensor,"tensor")
lmin = np.copy(l2)
lmax = np.copy(l1)
index = np.where(np.abs(l1)<np.abs(l2))
lmin[index] = l1[index]
lmax[index] = l2[index]
d = lmax-lmin
p = lmax+lmin

d_o_p = d/p
d_o_p[np.isnan(d_o_p)] = 0
output.PointData.append(np.abs(lmax),"lmax")
output.PointData.append(np.abs(lmin),"lmin")
output.PointData.append(l1,"l1")
output.PointData.append(l2,"l2")
output.PointData.append(d,"deviatoric")
output.PointData.append(p,"pressure")
output.PointData.append(p_o_d,"p_o_d")
"
panel_visibility="advanced">
<Hints>
<Widget type="multi_line"/>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment