diff --git a/HySoP/hysop/operator/fct2op.py b/HySoP/hysop/operator/fct2op.py
index fe085a5e5dea5e4950c93b64bf6db7565f58d863..af527f902da1954595814dabcbbae454117c7d7a 100644
--- a/HySoP/hysop/operator/fct2op.py
+++ b/HySoP/hysop/operator/fct2op.py
@@ -39,17 +39,26 @@ class Fct2Op(object):
We used Fct2Op.apply like a def function
"""
- if self.choice == 'gradU' :
+ if self.choice == 'gradU' : # field2 = gradU
#Make Scalar product
temp0 = self.field2[0][:,:,:]* y[0][:,:,:] +self.field2[1][:,:,:]* y[1][:,:,:]+self.field2[2][:,:,:]* y[2][:,:,:]
temp1 = self.field2[3][:,:,:]* y[0][:,:,:] +self.field2[4][:,:,:]* y[1][:,:,:]+self.field2[5][:,:,:]* y[2][:,:,:]
temp2 = self.field2[6][:,:,:]* y[0][:,:,:] +self.field2[7][:,:,:]* y[1][:,:,:]+self.field2[8][:,:,:]* y[2][:,:,:]
result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
return result
- if self.choice == 'divWU' :
+ if self.choice == 'divWU' : # field2 = velocity field
print 'taille field22:', self.field2.shape
result = DifferentialOperator_d(y, self.field2, self.choice, self.topology).apply()
return result
+ if self.choice == 'hessianU' : # field2 = velocity field
+ gradU, maxgersh = DifferentialOperator_d(self.field1, self.field2, choice='gradU', topology=self.topology).apply()
+ GradUline1 = np.concatenate((gradU[0],gradU[1],gradU[2]))
+ GradUline2 = np.concatenate((gradU[3],gradU[4],gradU[5]))
+ GradUline3 = np.concatenate((gradU[6],gradU[7],gradU[8]))
+ result1, maxgersh = DifferentialOperator_d(self.field1, GradUline1, choice='gradU', topology=self.topology).apply()
+ result2, maxgersh = DifferentialOperator_d(self.field1, GradUline2, choice='gradU', topology=self.topology).apply()
+ result3, maxgersh = DifferentialOperator_d(self.field1, GradUline3, choice='gradU', topology=self.topology).apply()
+ return result1, result2, result3
else :
print 'choice', choice , 'is not define'
diff --git a/HySoP/hysop/physics/compute_forces.py b/HySoP/hysop/physics/compute_forces.py
index a0594688b3840f6c3cfce90e3d691cf1d8abdd55..cbf38cc2124f02e27ea52a58f51f4679daa9020b 100644
--- a/HySoP/hysop/physics/compute_forces.py
+++ b/HySoP/hysop/physics/compute_forces.py
@@ -5,6 +5,7 @@ Compute forces
"""
from ..constants import *
import numpy as np
+from ..operator.fct2op import Fct2Op
import mpi4py.MPI as MPI
import time
@@ -38,7 +39,7 @@ class Compute_forces(object):
"""
Compute indicator functions for the control box (including the obstacle)
"""
- self.normal=np.zeros([self.dim*2, self.dim])
+ self.normal=np.zeros([self.dim*2, self.dim], dtype=PARMES_REAL)
self.Up = 0
self.Down = 1
self.East = 2
@@ -52,7 +53,7 @@ class Compute_forces(object):
self.normal[self.North][2] = 1.0
self.normal[self.South][2] = -1.0
- tmp_ind = np.zeros([self.dim, 2])
+ tmp_ind = np.zeros([self.dim, 2], dtype=PARMES_INTEGER)
distMin = self.boxMin - self.coordMin
distMax = self.boxMax - self.coordMax
@@ -74,7 +75,7 @@ class Compute_forces(object):
# otherwise, there is no volume control inside the local domain
- ind = np.zeros(self.dim*2)
+ ind = np.zeros(self.dim*2, dtype=PARMES_INTEGER)
ind[self.Down] = tmp_ind[0][0]
ind[self.Up] = tmp_ind[0][1]
ind[self.West] = tmp_ind[1][0]
@@ -88,7 +89,7 @@ class Compute_forces(object):
ind[self.North] = ind[self.North]-1
# Count the number of points on each face and inside the domain
- nbPoints = np.zeros(2*self.dim)
+ nbPoints = np.zeros(2*self.dim, dtype=PARMES_INTEGER)
nbPoints[self.Down] = (ind[self.East]-ind[self.West]+1)*(ind[self.North]-ind[self.South]+1)
nbPoints[self.Up] = (ind[self.East]-ind[self.West]+1)*(ind[self.North]-ind[self.South]+1)
nbPoints[self.West] = (ind[self.Up]-ind[self.Down]+1)*(ind[self.North]-ind[self.South]+1)
@@ -110,13 +111,13 @@ class Compute_forces(object):
if (distMin[2]>0. and distMax[2]>0. and self.boxMin[2]> self.coordMax[2]):
nbPoints[self.North] = 0
- boxDim = np.zeros(self.dim)
+ boxDim = np.zeros(self.dim, dtype=PARMES_INTEGER)
boxDim[0] = ind[self.Up] - ind[self.Down] + 1
boxDim[1] = ind[self.East] - ind[self.West] + 1
boxDim[2] = ind[self.North] - ind[self.South] + 1
nbPointsBox = boxDim[0] * boxDim[1] * boxDim[2]
- tmp_box = np.zeros([self.dim,nbPointsBox])
+ tmp_box = np.zeros([self.dim,nbPointsBox], dtype=PARMES_INTEGER)
count_box=0
if(boxDim>0):
for k in xrange (ind[self.South],ind[self.North]): ## enlever les boucles !!
@@ -134,12 +135,12 @@ class Compute_forces(object):
self.chi_box = tmp_box[:][0:count_box]
- self.chi_up = np.zeros([self.dim,nbPoints[self.Up]])
- self.chi_down = np.zeros([self.dim,nbPoints[self.Down]])
- self.chi_east = np.zeros([self.dim,nbPoints[self.East]])
- self.chi_west = np.zeros([self.dim,nbPoints[self.West]])
- self.chi_north = np.zeros([self.dim,nbPoints[self.North]])
- self.chi_south = np.zeros([self.dim,nbPoints[self.South]])
+ self.chi_up = np.zeros([self.dim,nbPoints[self.Up]], dtype=PARMES_INTEGER)
+ self.chi_down = np.zeros([self.dim,nbPoints[self.Down]], dtype=PARMES_INTEGER)
+ self.chi_east = np.zeros([self.dim,nbPoints[self.East]], dtype=PARMES_INTEGER)
+ self.chi_west = np.zeros([self.dim,nbPoints[self.West]], dtype=PARMES_INTEGER)
+ self.chi_north = np.zeros([self.dim,nbPoints[self.North]], dtype=PARMES_INTEGER)
+ self.chi_south = np.zeros([self.dim,nbPoints[self.South]], dtype=PARMES_INTEGER)
if(nbPoints[self.South]!=0.): # South boundary is in the domain
count = 0
@@ -197,7 +198,7 @@ class Compute_forces(object):
self.bufferForce = 0.
# Compute coef used to compute the drag coefficient
- self.coef = 2./(uinf**2*np.pi*self.obstacle.radius**2)
+ self.coef = 2./(uinf**2*np.pi*self.obstacle.radius**2) # attention, uinf !!!!!!!!!
def apply(self, dt, velo, vort, Re):
@@ -205,21 +206,28 @@ class Compute_forces(object):
Computation of the drag according to the "impulsion" formula presented by
Noca et. al in Journal of Fluids and Structures, 1999
"""
- localForce = 0.
- force = 0.
+ localForce = np.zeros(self.dim, dtype=PARMES_REAL)
+ force = np.zeros(self.dim, dtype=PARMES_REAL)
+ # computation of jacobian (du) and hessian (d2u) matrix
+ hessian1, hessian2, hessian3 = Fct2Op(vort, velo, choice = 'hessianU', topology=self.topo).apply(None,None)
+ self.jacobian, maxgersh = DifferentialOperator_d(vort, velo, choice='gradU', topology=self.topo).apply()
+ self.a = 2.*hessian1[0]+hessian1[4]+hessian2[1]+hessian1[8]+hessian3[2]
+ self.b = hessian2[0]+hessian1[1]+2.*hessian2[4]+hessian2[8]+hessian3[5]
+ self.c = hessian3[0]+hessian1[6]+hessian3[4]+hessian2[7]+2.*hessian3[8]
+
# Downstream and upstream surfaces
dsurf = step[1] * step[2]
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_down,self.normal[Down][:],X,Re,dsurf)
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_up,self.normal[Up][:],X,Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_down,self.normal[Down][:],Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_up,self.normal[Up][:],Re,dsurf)
# East and West surfaces
dsurf = step[0] * step[2]
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_east,self.normal[East][:],Y,Re,dsurf)
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_west,self.normal[West][:],Y,Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_east,self.normal[East][:],Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_west,self.normal[West][:],Re,dsurf)
# North and south surfaces
dsurf = step[0] * step[1]
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_south,self.normal[South][:],Z,Re,dsurf)
- localForce = integrateOnSurface(localForce,velo,vort,self.chi_north,self.normal[North][:],Z,Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_south,self.normal[South][:],Re,dsurf)
+ localForce = integrateOnSurface(localForce,velo,vort,self.chi_north,self.normal[North][:],Re,dsurf)
# over the volume
dvol = step[0] * step[1] * step[2]
localForce = integrateOnBox(localForce,vort,self.chi_box,dvol,dt)
@@ -236,7 +244,7 @@ class Compute_forces(object):
fact = -dvol/((self.dim-1)*dt)
int1 = 0.
# For all points in the box
- for ind in xrange (int(chi[1].shape)): # ?????
+ for ind in xrange (chi.shape[1]):
i = chi[0][ind]
j = chi[0][ind]
k = chi[0][ind]
@@ -248,7 +256,7 @@ class Compute_forces(object):
force = force + fact * (int1 - self.bufferForce)
self.bufferForce = int1 # Save for next time step ...
- return force #??????????
+ return force
def integrateOnSurface(self,force,velo,vort,chi,NormalVec,direction,Re,dsurf):
"""
@@ -259,10 +267,11 @@ class Compute_forces(object):
fact = 1./(self.dim-1)
# For each point of the current plane
- int1 = 0.
- int2 = 0.
- diff_dir = 0.
- for ind in xrange (int(chi[1].shape)): # ?????
+ int1 = np.zeros(self.dim, dtype=PARMES_REAL)
+ int2 = np.zeros(self.dim, dtype=PARMES_REAL)
+ X_n_DivT = np.zeros(self.dim, dtype=PARMES_REAL)
+ n_T = np.zeros(self.dim, dtype=PARMES_REAL)
+ for ind in xrange (chi.shape[1]):
i = chi[0][ind]
j = chi[0][ind]
k = chi[0][ind]
@@ -281,6 +290,26 @@ class Compute_forces(object):
- fact * n_u * np.cross(coords,vort[:][i,j,k]) \
+ fact * n_w * np.cross(coords,velo[:][i,j,k])
+ # part 2 of the force, the one concerning T = nu(nabla u + nabla uT)
+ if (chi == self.chi_down):
+ X_n_DivT[0] = -coords[1]*self.b-coords[2]*self.c
+ X_n_DivT[1] = -coords[0]*self.b
+ X_n_DivT[2] = -coords[0]*self.c
+ n_T[0] = -1./Re*2*self.jacobian[0]
+ n_T[1] = -1./Re*self.jacobian[3]+self.jacobian[1]
+ n_T[2] = -1./Re*self.jacobian[6]+self.jacobian[2]
+
+ if (chi == self.chi_up):
+ X_n_DivT[0] = coords[1]*self.b+coords[2]*self.c
+ X_n_DivT[1] = coords[0]*self.b
+ X_n_DivT[2] = coords[0]*self.c
+ n_T[0] = 1./Re*2*self.jacobian[0]
+ n_T[1] = 1./Re*self.jacobian[3]+self.jacobian[1]
+ n_T[2] = 1./Re*self.jacobian[6]+self.jacobian[2]
+
+ int2 = int2 + X_n_DivT+n_T
+ # Product with element of surface and sum to the total (input) force
+ force = force + (int1 + int2) * dsurf
def __str__(self):