Fixed several free surface evolution bugs

This commit is contained in:
Jose Luis Cercós Pita 2012-08-06 18:24:18 +02:00 committed by Yorik van Havre
parent 76c641bed8
commit 87c42d512c
4 changed files with 35 additions and 16 deletions

View File

@ -54,6 +54,11 @@ class simFSEvolution:
# Integrate variables
for i in range(0,nx):
for j in range(0,ny):
# Get value at pos using characteristics method
gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j])
gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + \
dt*np.linalg.norm(grad[i*ny+j])
# Free surface points position
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*grad[i*ny+j][2]
# Velocity potential
@ -67,6 +72,9 @@ class simFSEvolution:
for i in range(0,nx):
for j in [0,ny-1]:
self.boundaryCondition(i,j, waves, dt, t)
for j in range(0,ny):
for i in [0,nx-1]:
self.boundaryCondition(i,j, waves, dt, t)
def evaluateGradient(self):
""" Evaluate potential gradients over free surface.
@ -89,7 +97,7 @@ class simFSEvolution:
"""
nx = self.fs['Nx']
ny = self.fs['Ny']
grad = np.ndarray(3, dtype=np.float32)
grad = np.zeros(3, dtype=np.float32)
for i in range(0,nx):
for j in range(0,ny):
# Get source position (desingularized)
@ -100,6 +108,8 @@ class simFSEvolution:
# Get distance between points
d = pos-srcPos
grad = grad + d/np.dot(d,d)*src*area
# Discard Z induced effect by desingularization
grad[2] = 0.
return grad
def boundaryCondition(self, i,j, waves, dt, t):
@ -113,6 +123,8 @@ class simFSEvolution:
"""
pos = self.fs['pos'][i,j]
pos[2] = 0.
self.fs['velPot'][i,j] = 0.
self.fs['accPot'][i,j] = 0.
for w in waves['data']:
A = w[0]
T = w[1]
@ -125,11 +137,7 @@ class simFSEvolution:
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
amp = A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + amp
amp = frec*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['vel'][i,j][2] = self.fs['vel'][i,j][2] - amp
amp = frec*frec*A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['acc'][i,j][2] = self.fs['acc'][i,j][2] - amp
amp = grav/frec*A*np.sin(k*l - frec*(t+dt) + phase)
amp = - grav/frec*A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + amp
amp = grav*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + amp

View File

@ -43,7 +43,8 @@ class simInitialization:
# Compute time step
self.dt = 0.1
for w in self.waves['data']:
self.dt = np.min(self.dt, w[1]/200.0)
if(self.dt > w[1]/200.0):
self.dt = w[1]/200.0
def loadData(self, FSmesh, waves):
""" Convert data to numpy format.
@ -99,6 +100,7 @@ class simInitialization:
ny = self.fs['Ny']
for i in range(0,nx):
for j in range(0,ny):
self.fs['pos'][i,j][2] = 0.
for w in self.waves['data']:
A = w[0]
T = w[1]
@ -111,11 +113,7 @@ class simInitialization:
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
amp = A*np.sin(k*l + phase)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + amp
amp = frec*A*np.cos(k*l + phase)
self.fs['vel'][i,j][2] = self.fs['vel'][i,j][2] - amp
amp = frec*frec*A*np.sin(k*l + phase)
self.fs['acc'][i,j][2] = self.fs['acc'][i,j][2] - amp
amp = grav/frec*A*np.sin(k*l + phase)
amp = - grav/frec*A*np.sin(k*l + phase)
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + amp
amp = grav*A*np.cos(k*l + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + amp

View File

@ -48,11 +48,12 @@ class Singleton(type):
class FreeCADShipSimulation(threading.Thread):
__metaclass__ = Singleton
def __init__ (self, device, endTime, output, FSmesh, waves):
def __init__ (self, device, endTime, output, simInstance, FSmesh, waves):
""" Thread constructor.
@param device Device to use.
@param endTime Maximum simulation time.
@param output [Rate,Type] Output rate, Type=0 if FPS, 1 if IPF.
@param simInstance Simulaation instance.
@param FSmesh Free surface mesh faces.
@param waves Waves parameters (A,T,phi,heading)
"""
@ -70,6 +71,7 @@ class FreeCADShipSimulation(threading.Thread):
# Storage data
self.endTime = endTime
self.output = output
self.sim = simInstance
self.FSmesh = FSmesh
self.waves = waves
@ -93,9 +95,11 @@ class FreeCADShipSimulation(threading.Thread):
waves = init.waves
dt = init.dt
t = 0.0
nx = FS['Nx']
ny = FS['Ny']
msg = Translator.translate("\t[Sim]: Iterating...\n")
FreeCAD.Console.PrintMessage(msg)
while self.active:
while self.active and t < self.endTime:
msg = Translator.translate("\t\t[Sim]: Generating linear system matrix...\n")
FreeCAD.Console.PrintMessage(msg)
matGen.execute(FS, A)
@ -106,7 +110,16 @@ class FreeCADShipSimulation(threading.Thread):
FreeCAD.Console.PrintMessage(msg)
fsEvol.execute(FS, waves, dt, t)
t = t + dt
FreeCAD.Console.PrintMessage('t = %g s' % (t))
FreeCAD.Console.PrintMessage('t = %g s\n' % (t))
# Update FreeCAD
"""
pos = self.sim.FS_Position[:]
for i in range(0, nx):
for j in range(0, ny):
pos[i*ny+j].z = float(FS['pos'][i,j][2])
self.sim.FS_Position = pos[:]
FreeCAD.ActiveDocument.recompute()
"""
# Set thread as stopped (and prepare it to restarting)
self.active = False
threading.Event().set()

View File

@ -69,7 +69,7 @@ class TaskPanel:
msg = Translator.translate("Launching simulation...\n")
App.Console.PrintMessage(msg)
# Build simulation thread
simulator = Sim(device, endTime, output, FSMesh, waves)
simulator = Sim(device, endTime, output, self.sim, FSMesh, waves)
simulator.start()
msg = Translator.translate("Done!\n")
App.Console.PrintMessage(msg)