Regenerated branch
This commit is contained in:
Jose Luis Cercós Pita
@@ -1,119 +1,327 @@
|
||||
#***************************************************************************
|
||||
#* *
|
||||
#* Copyright (c) 2011, 2012 *
|
||||
#* Jose Luis Cercos Pita <jlcercos@gmail.com> *
|
||||
#* *
|
||||
#* *
|
||||
#* Copyright (c) 2011, 2012 *
|
||||
#* Jose Luis Cercos Pita <jlcercos@gmail.com> *
|
||||
#* *
|
||||
#* This program is free software; you can redistribute it and/or modify *
|
||||
#* it under the terms of the GNU Lesser General Public License (LGPL) *
|
||||
#* as published by the Free Software Foundation; either version 2 of *
|
||||
#* the License, or (at your option) any later version. *
|
||||
#* for detail see the LICENCE text file. *
|
||||
#* *
|
||||
#* This program is distributed in the hope that it will be useful, *
|
||||
#* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
#* GNU Library General Public License for more details. *
|
||||
#* *
|
||||
#* You should have received a copy of the GNU Library General Public *
|
||||
#* it under the terms of the GNU Lesser General Public License (LGPL) *
|
||||
#* as published by the Free Software Foundation; either version 2 of *
|
||||
#* the License, or (at your option) any later version. *
|
||||
#* for detail see the LICENCE text file. *
|
||||
#* *
|
||||
#* This program is distributed in the hope that it will be useful, *
|
||||
#* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
#* GNU Library General Public License for more details. *
|
||||
#* *
|
||||
#* You should have received a copy of the GNU Library General Public *
|
||||
#* License along with this program; if not, write to the Free Software *
|
||||
#* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *
|
||||
#* USA *
|
||||
#* *
|
||||
#* USA *
|
||||
#* *
|
||||
#***************************************************************************
|
||||
|
||||
# numpy
|
||||
import numpy as np
|
||||
import FreeCAD
|
||||
|
||||
grav=9.81
|
||||
|
||||
class simInitialization:
|
||||
def __init__(self, FSmesh, waves, context=None, queue=None):
|
||||
""" Constructor.
|
||||
@param FSmesh Initial free surface mesh.
|
||||
@param waves Considered simulation waves (A,T,phi,heading).
|
||||
@param context OpenCL context where apply. Only for compatibility,
|
||||
must be None.
|
||||
@param queue OpenCL command queue. Only for compatibility,
|
||||
must be None.
|
||||
"""
|
||||
self.context = context
|
||||
self.queue = queue
|
||||
self.loadData(FSmesh, waves)
|
||||
self.execute()
|
||||
# Compute time step
|
||||
self.dt = 0.1
|
||||
for w in self.waves['data']:
|
||||
if(self.dt > w[1]/200.0):
|
||||
self.dt = w[1]/200.0
|
||||
def __init__(self, h, FSMesh, SeaMesh, waves, context=None, queue=None):
|
||||
""" Constructor.
|
||||
@param h Water height.
|
||||
@param FSMesh Initial free surface mesh.
|
||||
@param waves Considered simulation waves (A,T,phi,heading).
|
||||
@param context OpenCL context where apply. Only for compatibility,
|
||||
must be None.
|
||||
@param queue OpenCL command queue. Only for compatibility,
|
||||
must be None.
|
||||
"""
|
||||
self.context = context
|
||||
self.queue = queue
|
||||
self.loadData(h, FSMesh, SeaMesh, waves)
|
||||
self.execute()
|
||||
# Compute time step
|
||||
self.dt = 0.1
|
||||
for w in self.waves['data']:
|
||||
if(self.dt > w[1]/200.0):
|
||||
self.dt = w[1]/200.0
|
||||
|
||||
def loadData(self, FSmesh, waves):
|
||||
""" Convert data to numpy format.
|
||||
@param FSmesh Initial free surface mesh.
|
||||
@param waves Considered simulation waves (A,T,phi,heading).
|
||||
"""
|
||||
nx = len(FSmesh)
|
||||
ny = len(FSmesh[0])
|
||||
nW = len(waves)
|
||||
# Mesh data
|
||||
p = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
n = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
a = np.ndarray((nx,ny), dtype=np.float32)
|
||||
phi = np.ndarray((nx,ny), dtype=np.float32)
|
||||
Phi = np.ndarray((nx,ny), dtype=np.float32)
|
||||
s = np.ndarray((nx,ny), dtype=np.float32)
|
||||
ss = np.ndarray((nx,ny), dtype=np.float32)
|
||||
for i in range(0, nx):
|
||||
for j in range(0, ny):
|
||||
pos = FSmesh[i][j].pos
|
||||
normal = FSmesh[i][j].normal
|
||||
area = FSmesh[i][j].area
|
||||
p[i,j,0] = pos.x
|
||||
p[i,j,1] = pos.y
|
||||
p[i,j,2] = pos.z
|
||||
n[i,j,0] = normal.x
|
||||
n[i,j,1] = normal.y
|
||||
n[i,j,2] = normal.z
|
||||
a[i,j] = area
|
||||
phi[i,j] = 0.
|
||||
Phi[i,j] = 0.
|
||||
s[i,j] = 0.
|
||||
ss[i,j] = 0.
|
||||
self.fs = {'Nx':nx, 'Ny':ny, 'pos':p, 'normal':n, 'area':a, \
|
||||
'velPot':phi, 'accPot':Phi, 'velSrc':s, 'accSrc':ss}
|
||||
# Waves data
|
||||
w = np.ndarray((nW, 4), dtype=np.float32)
|
||||
for i in range(0,nW):
|
||||
w[i,0] = waves[i][0]
|
||||
w[i,1] = waves[i][1]
|
||||
w[i,2] = waves[i][2]
|
||||
w[i,3] = waves[i][3]
|
||||
self.waves = {'N':nW, 'data':w}
|
||||
# Linear system matrix
|
||||
nF = nx*ny
|
||||
nB = 0 # No body for the moment
|
||||
N = nx*ny + nB
|
||||
self.A = np.ndarray((N, N), dtype=np.float32)
|
||||
def loadData(self, h, FSMesh, SeaMesh, waves):
|
||||
""" Convert data to numpy format.
|
||||
@param FSMesh Initial free surface mesh.
|
||||
@param waves Considered simulation waves (A,T,phi,heading).
|
||||
"""
|
||||
# Data will classified in four groups:
|
||||
# Free surface:
|
||||
# Is a key part of the simulation, so is
|
||||
# separated from the rest of water involved
|
||||
# elements.
|
||||
# Sea:
|
||||
# BEM method required a closed domain, so
|
||||
# water floor and sides must be append, but
|
||||
# are not a key objective of the simulation.
|
||||
# Body:
|
||||
# Is the main objective of the simulation.
|
||||
# Waves:
|
||||
# Data that is append as boundary condition.
|
||||
# BEM:
|
||||
# Used to solve the BEM problem and evolution.
|
||||
|
||||
# --------------------------------------------
|
||||
# Free surface data
|
||||
# N, Nx, Ny = Number of points in each
|
||||
# direction
|
||||
# pos = Positions
|
||||
# vel = Velocities
|
||||
# n = Normals
|
||||
# area = Areas
|
||||
# --------------------------------------------
|
||||
nx = len(FSMesh)
|
||||
ny = len(FSMesh[0])
|
||||
p = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
V = np.zeros((nx,ny, 3), dtype=np.float32)
|
||||
n = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
a = np.ndarray((nx,ny), dtype=np.float32)
|
||||
x1 = np.zeros((nx,ny), dtype=np.float32)
|
||||
x2 = np.zeros((nx,ny), dtype=np.float32)
|
||||
x3 = np.zeros((nx,ny), dtype=np.float32)
|
||||
dx1 = np.zeros((nx,ny), dtype=np.float32)
|
||||
dx2 = np.zeros((nx,ny), dtype=np.float32)
|
||||
dx3 = np.zeros((nx,ny), dtype=np.float32)
|
||||
for i in range(0, nx):
|
||||
for j in range(0, ny):
|
||||
pos = FSMesh[i][j].pos
|
||||
normal = FSMesh[i][j].normal
|
||||
area = FSMesh[i][j].area
|
||||
p[i,j,0] = pos.x
|
||||
p[i,j,1] = pos.y
|
||||
p[i,j,2] = pos.z
|
||||
n[i,j,0] = normal.x
|
||||
n[i,j,1] = normal.y
|
||||
n[i,j,2] = normal.z
|
||||
a[i,j] = area
|
||||
self.fs = {'h': h, 'N':nx*ny, 'Nx':nx, 'Ny':ny, \
|
||||
'pos':p, 'vel':V, 'normal':n, 'area':a, \
|
||||
'x1':x1, 'x2':x2, 'x3':x3,\
|
||||
'dx1':dx1, 'dx2':dx2, 'dx3':dx3}
|
||||
# --------------------------------------------
|
||||
# Sea data (dictionary with components
|
||||
# ['front','back','left','right','bottom'])
|
||||
# N, Nx, Ny = Number of points in each
|
||||
# direction
|
||||
# pos = Positions
|
||||
# vel = Velocities
|
||||
# n = Normals
|
||||
# area = Areas
|
||||
# --------------------------------------------
|
||||
self.sea = {'ids':['front','back','left','right','bottom']}
|
||||
N = 0
|
||||
for index in self.sea['ids']:
|
||||
mesh = SeaMesh[index]
|
||||
nx = len(mesh)
|
||||
ny = len(mesh[0])
|
||||
p = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
V = np.zeros((nx,ny, 3), dtype=np.float32)
|
||||
n = np.ndarray((nx,ny, 3), dtype=np.float32)
|
||||
a = np.ndarray((nx,ny), dtype=np.float32)
|
||||
for i in range(0, nx):
|
||||
for j in range(0, ny):
|
||||
pos = mesh[i][j].pos
|
||||
normal = mesh[i][j].normal
|
||||
area = mesh[i][j].area
|
||||
p[i,j,0] = pos.x
|
||||
p[i,j,1] = pos.y
|
||||
p[i,j,2] = pos.z
|
||||
n[i,j,0] = normal.x
|
||||
n[i,j,1] = normal.y
|
||||
n[i,j,2] = normal.z
|
||||
a[i,j] = area
|
||||
d = {'N':nx*ny, 'Nx':nx, 'Ny':ny, 'pos':p, 'vel':V, 'normal':n, 'area':a}
|
||||
self.sea[index] = d
|
||||
N = N + nx*ny
|
||||
self.sea['N'] = N
|
||||
self.sea['h'] = h
|
||||
# --------------------------------------------
|
||||
# Body data
|
||||
# N, Nx, Ny = Number of points in each
|
||||
# direction
|
||||
# pos = Positions
|
||||
# vel = Velocities
|
||||
# n = Normals
|
||||
# area = Areas
|
||||
# --------------------------------------------
|
||||
self.b = {'N':0, 'pos':None, 'vel':None, 'normal':None, 'area':None}
|
||||
# --------------------------------------------
|
||||
# Waves data
|
||||
# N = Number of waves
|
||||
# data = Waves data
|
||||
# --------------------------------------------
|
||||
nW = len(waves)
|
||||
w = np.ndarray((nW, 4), dtype=np.float32)
|
||||
for i in range(0,nW):
|
||||
w[i,0] = waves[i][0]
|
||||
w[i,1] = waves[i][1]
|
||||
w[i,2] = waves[i][2]
|
||||
w[i,3] = waves[i][3]
|
||||
self.waves = {'h':h, 'N':nW, 'data':w}
|
||||
# --------------------------------------------
|
||||
# BEM data
|
||||
# N = nFS + nSea + nB
|
||||
# A,B,dB = Linear system matrix and vectors
|
||||
# p1,... = Velocity potentials (phi) for
|
||||
# each RK4 step. In reallity are
|
||||
# the independent term of the
|
||||
# BEM linear system, so is the
|
||||
# potential for the free surface,
|
||||
# and the gradient projected over
|
||||
# the normal along all other terms.
|
||||
# dp1,... = Acceleration potentials
|
||||
# (dphi/dt) for each RK4 step.
|
||||
# In reallity are the
|
||||
# independent term of the BEM
|
||||
# linear system, so is the
|
||||
# potential for the free surface,
|
||||
# and the gradient projected over
|
||||
# the normal along all other terms.
|
||||
# Ap,Adp = BEM solution vectors, that
|
||||
# contains the potential gradients
|
||||
# on free surface, and the potential
|
||||
# along all toher surfaces.
|
||||
# --------------------------------------------
|
||||
nFS = self.fs['N']
|
||||
nSea = self.sea['N']
|
||||
nB = self.b['N']
|
||||
N = nFS + nSea + nB
|
||||
A = np.zeros((N, N), dtype=np.float32)
|
||||
B = np.zeros((N), dtype=np.float32)
|
||||
dB = np.zeros((N), dtype=np.float32)
|
||||
p1 = np.zeros((N), dtype=np.float32)
|
||||
p2 = np.zeros((N), dtype=np.float32)
|
||||
p3 = np.zeros((N), dtype=np.float32)
|
||||
p4 = np.zeros((N), dtype=np.float32)
|
||||
Ap = np.zeros((N), dtype=np.float32)
|
||||
dp1 = np.zeros((N), dtype=np.float32)
|
||||
dp2 = np.zeros((N), dtype=np.float32)
|
||||
dp3 = np.zeros((N), dtype=np.float32)
|
||||
dp4 = np.zeros((N), dtype=np.float32)
|
||||
Adp = np.zeros((N), dtype=np.float32)
|
||||
self.bem = {'N':N, 'A':A, 'B':B, 'dB':dB, \
|
||||
'p1':p1, 'p2':p2, 'p3':p3, 'p4':p4, 'Ap':Ap, \
|
||||
'dp1':dp1, 'dp2':dp2, 'dp3':dp3, 'dp4':dp4, 'Adp':Adp }
|
||||
|
||||
def execute(self):
|
||||
""" Compute initial conditions. """
|
||||
# --------------------------------------------
|
||||
# Free surface beach nodes.
|
||||
# Beach nodes are the nodes of the free
|
||||
# surface where the waves are imposed. All
|
||||
# the other nodes are computed allowing non
|
||||
# linear waves due to the ship interaction.
|
||||
# The beach will have enough dimension to
|
||||
# control at least half wave length
|
||||
# --------------------------------------------
|
||||
# Get maximum wave length
|
||||
wl = 0.0
|
||||
for w in self.waves['data']:
|
||||
T = w[1]
|
||||
wl = max(wl, 0.5 * grav / np.pi * T*T)
|
||||
# Get nodes dimensions
|
||||
nx = self.fs['Nx']
|
||||
ny = self.fs['Ny']
|
||||
lx = self.fs['pos'][nx-1,0][0] - self.fs['pos'][0,0][0]
|
||||
ly = self.fs['pos'][0,ny-1][1] - self.fs['pos'][0,0][1]
|
||||
dx = lx / nx
|
||||
dy = ly / ny
|
||||
# Get number of nodes involved
|
||||
wnx = max(1, int(round(0.5*wl / dx)))
|
||||
wny = max(1, int(round(0.5*wl / dy)))
|
||||
wnx = min(wnx, nx)
|
||||
wny = min(wny, ny)
|
||||
self.fs['Beachx'] = wnx
|
||||
self.fs['Beachy'] = wny
|
||||
# --------------------------------------------
|
||||
# Free surface initial condition.
|
||||
# Since RK4 scheme starts on the end of
|
||||
# previous step, we only write on last
|
||||
# stage value (p4 and dp4)
|
||||
# --------------------------------------------
|
||||
nx = self.fs['Nx']
|
||||
ny = self.fs['Ny']
|
||||
h = self.fs['h']
|
||||
for i in range(0,nx):
|
||||
for j in range(0,ny):
|
||||
# Since initial values of the potencial, and this acceleration,
|
||||
# depends on z, we need to compute first the positions.
|
||||
self.fs['pos'][i,j][2] = 0.
|
||||
for w in self.waves['data']:
|
||||
A = w[0]
|
||||
T = w[1]
|
||||
phase = w[2]
|
||||
heading = np.pi*w[3]/180.0
|
||||
wl = 0.5 * grav / np.pi * T*T
|
||||
k = 2.0*np.pi/wl
|
||||
frec = 2.0*np.pi/T
|
||||
pos = self.fs['pos'][i,j]
|
||||
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
|
||||
# hfact = np.sinh(k*(pos[2]+h)) / np.cosh(k*h)
|
||||
hfact = 1.0
|
||||
amp = A*np.sin(k*l + phase)*hfact
|
||||
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + amp
|
||||
amp = - A*frec*np.cos(k*l + phase)*hfact
|
||||
self.fs['vel'][i,j][2] = self.fs['vel'][i,j][2] + amp
|
||||
# And now we can compute potentials.
|
||||
for w in self.waves['data']:
|
||||
A = w[0]
|
||||
T = w[1]
|
||||
phase = w[2]
|
||||
heading = np.pi*w[3]/180.0
|
||||
wl = 0.5 * grav / np.pi * T*T
|
||||
k = 2.0*np.pi/wl
|
||||
frec = 2.0*np.pi/T
|
||||
pos = self.fs['pos'][i,j]
|
||||
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
|
||||
hfact = np.cosh(k*(pos[2]+h)) / np.cosh(k*h)
|
||||
amp = - grav/frec*A*np.cos(k*l + phase)*hfact
|
||||
self.bem['p4'][i*ny+j] = self.bem['p4'][i*ny+j] + amp
|
||||
amp = - grav*A*np.sin(k*l + phase)*hfact
|
||||
self.bem['dp4'][i*ny+j] = self.bem['dp4'][i*ny+j] + amp
|
||||
# --------------------------------------------
|
||||
# Sea initial condition on sides.
|
||||
# 1. Since RK4 scheme starts on the end of
|
||||
# previous step, we only write on last
|
||||
# stage value (p4 and dp4)
|
||||
# 2. In the sea boundaries we are
|
||||
# interested on the gradient of the
|
||||
# potentials projected over the normal,
|
||||
# so we really store this value.
|
||||
# 3. In the floor this value is ever null.
|
||||
# --------------------------------------------
|
||||
ids = ['front','back','left','right','bottom']
|
||||
i0 = self.fs['N']
|
||||
for index in ids:
|
||||
sea = self.sea[index]
|
||||
nx = sea['Nx']
|
||||
ny = sea['Ny']
|
||||
for i in range(0,nx):
|
||||
for j in range(0,ny):
|
||||
for w in self.waves['data']:
|
||||
A = w[0]
|
||||
T = w[1]
|
||||
phase = w[2]
|
||||
heading = np.pi*w[3]/180.0
|
||||
wl = 0.5 * grav / np.pi * T*T
|
||||
k = 2.0*np.pi/wl
|
||||
frec = 2.0*np.pi/T
|
||||
pos = sea['pos'][i,j]
|
||||
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
|
||||
normal = sea['normal'][i,j]
|
||||
hfact = np.cosh(k*(pos[2]+h)) / np.cosh(k*h)
|
||||
factor = np.dot(normal,np.array([np.cos(heading), np.sin(heading), 0.]))
|
||||
amp = frec*A*np.sin(k*l + phase)*hfact
|
||||
self.bem['p4'][i0 + i*ny+j] = self.bem['p4'][i*ny+j] + factor*amp
|
||||
amp = - grav*A*k*np.cos(k*l + phase)*hfact
|
||||
self.bem['dp4'][i0 + i*ny+j] = self.bem['dp4'][i*ny+j] + factor*amp
|
||||
i0 = i0 + sea['N']
|
||||
|
||||
def execute(self):
|
||||
""" Compute initial conditions. """
|
||||
nx = self.fs['Nx']
|
||||
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]
|
||||
phase = w[2]
|
||||
heading = np.pi*w[3]/180.0
|
||||
wl = 0.5 * grav / np.pi * T*T
|
||||
k = 2.0*np.pi/wl
|
||||
frec = 2.0*np.pi/T
|
||||
pos = self.fs['pos'][i,j]
|
||||
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 = - 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
|
||||
|
||||
Reference in New Issue
Block a user