Commit 7686246c authored by Nathan Coppin's avatar Nathan Coppin
Browse files

updating

parents 5adbf949 4caffdf3
......@@ -63,7 +63,7 @@
name="Script"
command="SetScript"
number_of_elements="1"
default_values="import numpy as np

in_particles,in_bnd,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
vals_t = []
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
 disks = in_bnd.Cells[1:2*n_disks:2]
 vals.append(in_contacts.FieldData["particle_disk"])
 vals_t.append(in_contacts.FieldData["particle_disk_t"])
 contacts = in_contacts.FieldData["particle_disk_idx"]
 points_d = np.ndarray((contacts.shape[0],2,3))
 points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
 points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
 points.append(points_d)

#particle-segments contacts

if n_segments :
 segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
 vals.append(in_contacts.FieldData["particle_segment"])
 vals_t.append(in_contacts.FieldData["particle_segment_t"])
 contacts = in_contacts.FieldData["particle_segment_idx"]
 points_s = np.ndarray((contacts.shape[0],2,3))
 points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
 s = in_bnd.Points[segments[contacts[:,0],:]]
 t = s[:,1,:]-s[:,0,:]
 t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
 d = points_s[:,0,:]-s[:,0,:]
 l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
 points_s[:,1,:] = s[:,0,:]+l[:,None]*t
 points.append(points_s)

#particle-triangle contacts
if n_triangles :
 triangles = in_bnd.Cells[2*n_disks+3*n_segments,2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
 vals.append(in_contacts.FieldData["particle_triangle"])
 vals_t.append(in_contacts.FieldData["particle_triangle_t"])
 vals.append(in_contacts.FieldData["particle_triangle_t"])
 contacts = in_contacts.FieldData["particle_triangle_idx"]
 points_t = np.ndarray((contacts.shape[0],2,3))
 points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
 points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
 points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
 +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
 +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
 j = (i+1)%nf
 pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
default_values="import numpy as np

in_particles,in_bnd,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
 vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
 disks = in_bnd.Cells[1:2*n_disks:2]
 vals.append(in_contacts.FieldData["particle_disk"])
 vals_t.append(in_contacts.FieldData["particle_disk_t"])
 if vals_s is not None:
 vals_s.append(in_contacts.FieldData["particle_disk_s"])
 contacts = in_contacts.FieldData["particle_disk_idx"]
 points_d = np.ndarray((contacts.shape[0],2,3))
 points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
 points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
 points.append(points_d)

#particle-segments contacts

if n_segments :
 segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
 vals.append(in_contacts.FieldData["particle_segment"])
 vals_t.append(in_contacts.FieldData["particle_segment_t"])
 if vals_s is not None:
 vals_s.append(in_contacts.FieldData["particle_segment_s"])
 contacts = in_contacts.FieldData["particle_segment_idx"]
 points_s = np.ndarray((contacts.shape[0],2,3))
 points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
 s = in_bnd.Points[segments[contacts[:,0],:]]
 t = s[:,1,:]-s[:,0,:]
 t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
 d = points_s[:,0,:]-s[:,0,:]
 l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
 points_s[:,1,:] = s[:,0,:]+l[:,None]*t
 points.append(points_s)

#particle-triangle contacts
if n_triangles :
 triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
 vals.append(in_contacts.FieldData["particle_triangle"])
 vals_t.append(in_contacts.FieldData["particle_triangle_t"])
 if vals_s is not None:
 vals_s.append(in_contacts.FieldData["particle_triangle_s"])
 contacts = in_contacts.FieldData["particle_triangle_idx"]
 points_t = np.ndarray((contacts.shape[0],2,3))
 points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
 points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
 points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
 vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
 +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
 +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
 output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
 j = (i+1)%nf
 pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
panel_visibility="advanced">
<Hints>
<Widget type="multi_line"/>
......
......@@ -43,11 +43,14 @@ def RequestData():
vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []
#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
......@@ -55,6 +58,8 @@ def RequestData():
disks = in_bnd.Cells[1:2*n_disks:2]
vals.append(in_contacts.FieldData["particle_disk"])
vals_t.append(in_contacts.FieldData["particle_disk_t"])
if vals_s is not None:
vals_s.append(in_contacts.FieldData["particle_disk_s"])
contacts = in_contacts.FieldData["particle_disk_idx"]
points_d = np.ndarray((contacts.shape[0],2,3))
points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
......@@ -67,6 +72,8 @@ def RequestData():
segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
vals.append(in_contacts.FieldData["particle_segment"])
vals_t.append(in_contacts.FieldData["particle_segment_t"])
if vals_s is not None:
vals_s.append(in_contacts.FieldData["particle_segment_s"])
contacts = in_contacts.FieldData["particle_segment_idx"]
points_s = np.ndarray((contacts.shape[0],2,3))
points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
......@@ -80,14 +87,21 @@ def RequestData():
#particle-triangle contacts
if n_triangles :
triangles = in_bnd.Cells[2*n_disks+3*n_segments,2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
vals.append(in_contacts.FieldData["particle_triangle"])
vals_t.append(in_contacts.FieldData["particle_triangle_t"])
vals.append(in_contacts.FieldData["particle_triangle_t"])
if vals_s is not None:
vals_s.append(in_contacts.FieldData["particle_triangle_s"])
contacts = in_contacts.FieldData["particle_triangle_idx"]
points_t = np.ndarray((contacts.shape[0],2,3))
points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
d0 = in_bnd.Points[triangles[contacts[:,0],:]][:,1,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
d1 = in_bnd.Points[triangles[contacts[:,0],:]][:,2,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
N = np.cross(d0,d1)
N /= np.linalg.norm(N,axis=1)[:,newaxis]
dd = in_bnd.Points[triangles[contacts[:,0],:]][:,0,:] - in_particles.Points[contacts[:,1],:]
dist = einsum('ij,ij->i', N, dd)
points_t[:,1,:] = in_particles.Points[contacts[:,1],:] + N*dist[:,newaxis]
points.append(points_t)
#merge everything
......@@ -95,6 +109,8 @@ def RequestData():
points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
......@@ -110,6 +126,8 @@ def RequestData():
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
......
......@@ -208,8 +208,8 @@ class ParticleProblem :
"""Returns the contact forces on each grain."""
return (self._get_matrix("ContactForces",self._dim))
def get_boundary_forces(self,tag="default") :
"""Returns the net normal and tangential forces acting on a boundary because of the contacts.
def get_boundary_impulses(self,tag="default") :
"""Returns the net normal and tangential impulses acting on a boundary because of the contacts.
Keyword arguments:
tag -- Name of the boundary
"""
......@@ -360,7 +360,6 @@ class ParticleProblem :
nameb = name.encode()
fdata.append((nameb,v[:,[0]]))
fdata.append((nameb+b"_t",v[:,[1]]))
fdata.append((nameb+b"_dir_n",basis[:,0,:self._dim-1]))
if self._dim == 3:
fdata.append((nameb+b"_s",v[:,[2]]))
fdata.append((nameb+b"_idx",o))
......@@ -413,25 +412,10 @@ class ParticleProblem :
v = []
v.append(np.vstack([fdata[k] for k in ks]))
v.append(np.vstack([fdata[k+"_t"] for k in ks]))
basis = []
cheat = np.vstack([fdata[k+"_dir_n"] for k in ks])[:,:2]
if self._dim == 2:
cheat = np.append(cheat,((1 - cheat[:,0]**2)**0.5).reshape((cheat.shape[0],1)),axis=1)
cheat[np.isnan(cheat[:,1]),1] = 0
basis.append(cheat)
basis.append(cheat[:,[1,0]]*np.array([-1, 1]))
else:
cheat = np.append(cheat,((1-cheat[:,0]**2-cheat[:,1]**2)**0.5).reshape((cheat.shape[0],1)),axis=1)
cheat[np.isnan(cheat[:,2]),2] = 0
basis.append(cheat)
dirt = cheat[:,[2,0,1]]*np.array([-1,1,-1])-cheat*np.einsum("ij,ij->i",cheat,cheat[:,[2,0,1]]*np.array([-1, 1, -1]))[:,np.newaxis]
nnorm = np.linalg.norm(dirt,axis=1)
nnorm[nnorm==0] = 1
basis.append(dirt/nnorm[:,np.newaxis])
basis.append(np.cross(cheat,dirt/nnorm[:,np.newaxis]))
if self._dim == 3:
v.append(np.vstack([fdata[k+"_s"] for k in ks]))
basis = np.column_stack(basis)
v = np.column_stack(v)*contact_factor
basis = np.zeros((v.shape[0],self._dim,self._dim))
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(basis),_np2c(t,np.int32))
......
......@@ -85,6 +85,8 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
vp0 = np.copy(particles.velocity())
pp0 = np.copy(particles.position())
cp0 = np.copy(particles.contact_forces())
if particles._rotation_enabled :
op0 = np.copy(particles.omega())
# Predictor
#
......@@ -113,6 +115,8 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.velocity()[:,:] = vp0
if particles._rotation_enabled :
particles.omega()[:,:] = op0
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces())
......
......@@ -1121,13 +1121,6 @@ int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion, double *externalForces)
{
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv *= dt;
p->contacts[j].ct *= dt;
#if DIMENSION==3
p->contacts[j].cs *= dt;
#endif
}
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
......@@ -1149,13 +1142,6 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
}
}
free(oldVelocity);
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
return 0;
}
for (size_t i = 0; i < vectorSize(p->position); ++i){
......@@ -1183,13 +1169,6 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
p->contactForces[j*DIMENSION+i] = (p->velocity[j*DIMENSION+i]-p->contactForces[j*DIMENSION+i])*p->particles[j].m/dt;
}
}
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
free(oldVelocity);
return 1;
}
......
......@@ -31,8 +31,13 @@ import os
import time
import shutil
import random
import sys
listeVit = [14,36,64,95,130,169,211,255,303,354,408,467]
inter = int(sys.argv[1])
index = int(inter/5)
numDep = inter%5
# Physical parameters
vit = -0.05 #stream velocity
vit = -0.001*listeVit[index] #stream velocity
muwall = 0.3 #friction coefficient between the walls and the particles
mupart = 0.15 #friction coefficient between particles
rho = 1000 #fluid density
......@@ -41,9 +46,9 @@ g = np.array([0,-9.81,0]) #gravity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
tEnd = 30 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
outf = 2500 #iterations between data frames
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
......@@ -54,12 +59,12 @@ filename = outputdir + "/Drag.csv"
#
#Injected particles
p = scontact.ParticleProblem(3,True,True)
p.read_vtk("depot/output",10)
p.read_vtk("depot/"+ str(numDep),2)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
p.velocity()[:,1] = vit
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
p2.read_vtk("depot/output",10)
p.read_vtk("depot/"+ str(numDep),2)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
......@@ -69,17 +74,17 @@ p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0,0])
fluid.set_wall_boundary("Right", velocity = [0,0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
fluid.set_wall_boundary("Front", velocity = [0,0,0])
fluid.set_wall_boundary("Rear", velocity = [0,0,0])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
fluid.export_vtk(outputdir, 0.,0)
#fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
#fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,0,0])
#fluid.set_wall_boundary("Right", velocity = [0,0,0])
#fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
#fluid.set_wall_boundary("Front", velocity = [0,0,0])
#fluid.set_wall_boundary("Rear", velocity = [0,0,0])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir, 0.,0)
#
# COMPUTATION LOOP
#
......@@ -93,9 +98,10 @@ while t < tEnd :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.52))
F = np.zeros(3)
time_integration.particle_changed(fluid,p2)
time_integration.iterate(fluid,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
with open(filename,"a") as file1:
F /= dt
F += fluid.boundary_force()[0,:]
file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
t += dt
......@@ -103,6 +109,6 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p2.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
#fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
......@@ -124,7 +124,7 @@ class Weight(unittest.TestCase) :
ioutput = int(ii/outf) + 1
#p.write_vtk(outputdir, ioutput, t)
print("%i : %.2g/%.2g" % (ii, t, tEnd))
error = np.sqrt(sum(sum(forces) - bnd_forces)**2)
error = np.sqrt(sum(sum(forces) - bnd_forces/dt)**2)
error /= np.sqrt(sum(sum(forces)**2))
print(error*100,"%")
self.assertLess(error,.5/100., "error is too large in weight")
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