diff --git a/CodesEnVrac/CodeGH/src-THI/stretch.f b/CodesEnVrac/CodeGH/src-THI/stretch.f
index 58f817fa6ef9d6832f6171c77044ba55a90314fe..7173c5cda497f93b0098b5bf268236a620d00cea 100644
--- a/CodesEnVrac/CodeGH/src-THI/stretch.f
+++ b/CodesEnVrac/CodeGH/src-THI/stretch.f
@@ -31,32 +31,66 @@ cc
ibb=mod(i-3+nx1,nx1)+1
c
- aux1=omg1(itt,j,k)*vxg(itt,j,k)+8.*(omg1(ib,j,k)*vxg(ib,j,k)-
- 1 omg1(it,j,k)*vxg(it,j,k))-omg1(ibb,j,k)*vxg(ibb,j,k)
- aux2=omg2(i,jtt,k)*vxg(i,jtt,k)+8.*(omg2(i,jb,k)*vxg(i,jb,k)-
- 1 omg2(i,jt,k)*vxg(i,jt,k))-omg2(i,jbb,k)*vxg(i,jbb,k)
- aux3=omg3(i,j,ktt)*vxg(i,j,ktt)+8.*(omg3(i,j,kb)*vxg(i,j,kb)-
- 1 omg3(i,j,kt)*vxg(i,j,kt))-omg3(i,j,kbb)*vxg(i,j,kbb)
-
- strg1(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
-c
- aux1=omg1(itt,j,k)*vyg(itt,j,k)+8.*(omg1(ib,j,k)*vyg(ib,j,k)-
- 1 omg1(it,j,k)*vyg(it,j,k))-omg1(ibb,j,k)*vyg(ibb,j,k)
- aux2=omg2(i,jtt,k)*vyg(i,jtt,k)+8.*(omg2(i,jb,k)*vyg(i,jb,k)-
- 1 omg2(i,jt,k)*vyg(i,jt,k))-omg2(i,jbb,k)*vyg(i,jbb,k)
- aux3=omg3(i,j,ktt)*vyg(i,j,ktt)+8.*(omg3(i,j,kb)*vyg(i,j,kb)-
- 1 omg3(i,j,kt)*vyg(i,j,kt))-omg3(i,j,kbb)*vyg(i,j,kbb)
-
- strg2(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
-!
- aux1=omg1(itt,j,k)*vzg(itt,j,k)+8.*(omg1(ib,j,k)*vzg(ib,j,k)-
- 1 omg1(it,j,k)*vzg(it,j,k))-omg1(ibb,j,k)*vzg(ibb,j,k)
- aux2=omg2(i,jtt,k)*vzg(i,jtt,k)+8.*(omg2(i,jb,k)*vzg(i,jb,k)-
- 1 omg2(i,jt,k)*vzg(i,jt,k))-omg2(i,jbb,k)*vzg(i,jbb,k)
- aux3=omg3(i,j,ktt)*vzg(i,j,ktt)+8.*(omg3(i,j,kb)*vzg(i,j,kb)-
- 1 omg3(i,j,kt)*vzg(i,j,kt))-omg3(i,j,kbb)*vzg(i,j,kbb)
-
- strg3(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
+ aux1=
+omg1(itt,j,k)*vxg(itt,j,k)
++8.*(omg1(ib,j,k)*vxg(ib,j,k)
+-omg1(it,j,k)*vxg(it,j,k))
+-omg1(ibb,j,k)*vxg(ibb,j,k)
+
+
+ aux2=
+omg2(i,jtt,k)*vxg(i,jtt,k)
++8.*(omg2(i,jb,k)*vxg(i,jb,k)
+-omg2(i,jt,k)*vxg(i,jt,k))
+-omg2(i,jbb,k)*vxg(i,jbb,k)
+
+ aux3=
+omg3(i,j,ktt)*vxg(i,j,ktt)
++8.*(omg3(i,j,kb)*vxg(i,j,kb)
+-omg3(i,j,kt)*vxg(i,j,kt))
+-omg3(i,j,kbb)*vxg(i,j,kbb)
+
+strg1(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
+
+ aux1=
+omg1(itt,j,k)*vyg(itt,j,k)
++8.*(omg1(ib,j,k)*vyg(ib,j,k)
+-omg1(it,j,k)*vyg(it,j,k))
+-omg1(ibb,j,k)*vyg(ibb,j,k)
+
+ aux2=
+omg2(i,jtt,k)*vyg(i,jtt,k)
++8.*(omg2(i,jb,k)*vyg(i,jb,k)
+-omg2(i,jt,k)*vyg(i,jt,k))
+-omg2(i,jbb,k)*vyg(i,jbb,k)
+
+ aux3=
+omg3(i,j,ktt)*vyg(i,j,ktt)
++8.*(omg3(i,j,kb)*vyg(i,j,kb)
+-omg3(i,j,kt)*vyg(i,j,kt))
+-omg3(i,j,kbb)*vyg(i,j,kbb)
+
+ strg2(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
+
+ aux1=
+omg1(itt,j,k)*vzg(itt,j,k)
++8.*(omg1(ib,j,k)*vzg(ib,j,k)
+-omg1(it,j,k)*vzg(it,j,k))
+-omg1(ibb,j,k)*vzg(ibb,j,k)
+
+ aux2=
+omg2(i,jtt,k)*vzg(i,jtt,k)
++8.*(omg2(i,jb,k)*vzg(i,jb,k)
+-omg2(i,jt,k)*vzg(i,jt,k))
+-omg2(i,jbb,k)*vzg(i,jbb,k)
+
+ aux3=
+omg3(i,j,ktt)*vzg(i,j,ktt)
++8.*(omg3(i,j,kb)*vzg(i,j,kb)
+-omg3(i,j,kt)*vzg(i,j,kt))
+-omg3(i,j,kbb)*vzg(i,j,kbb)
+
+ strg3(i,j,k)=-aux1*dxinv-aux2*dyinv-aux3*dzinv
enddo
enddo
diff --git a/Examples/mainED.py b/Examples/mainED.py
new file mode 100755
index 0000000000000000000000000000000000000000..21e3fcdf56856fa50f8cc34c81b54643905095d2
--- /dev/null
+++ b/Examples/mainED.py
@@ -0,0 +1,79 @@
+# -*- coding: utf-8 -*-
+import time
+import parmepy as pp
+from math import *
+
+
+def vitesse(x, y, z):
+ vx = 1. + x
+ vy = - x * y
+ vz = x * y * z + 10.
+ return vx, vy, vz
+
+def vorticite(x, y, z):
+ wx = x * y
+ wy = y * z
+ wz = - y
+ return wx, wy, wz
+
+def scalaire(x, y, z):
+ if x < 0.5 and y < 0.5 and z < 0.5:
+ return 1.
+ else:
+ return 0.
+
+
+def run():
+ # Parameters
+ nb = 32
+ timeStep = 0.02
+ finalTime = 1.
+ outputFilePrefix = './res/Stretching_'
+ outputModulo = 1
+
+ t0 = time.time()
+
+ ## Domain
+ box = pp.Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
+
+ ## Fields
+ #scal = pp.ContinuousField(domain=box, name='Scalar')
+ #velo = pp.ContinuousField(domain=box, name='Velocity', vector=True)
+ #vorti = pp.ContinuousField(domain=box, name='Vorticity', vector=True)
+ #scal = pp.AnalyticalField(domain=box, name='Scalar')
+ velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
+ vorti = pp.AnalyticalField(domain=box, formula=vorticite, name='Vorticity', vector=True)
+
+ ## Operators
+ #advec = pp.Transport(velo, scal)
+ stretch = pp.Stretching(velo,vorti)
+
+ ## Solver creation (discretisation of objects is done in solver initialisation)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
+
+ ##Problem
+ #pb = pp.Problem(topo3D, [advec,stretch])
+ pb = pp.Problem(topo3D, [stretch])
+
+ ## Setting solver to Problem
+ ## Problem => ParticularSover = basic.init()
+ pb.setSolver(finalTime, timeStep, solver_type='basic',
+ io=pp.Printer(fields=[vorti, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
+ ## Problem => ParticularSover = basic.initialize()
+ pb.initSolver()
+
+ t1 = time.time()
+
+ ## Solve problem
+ pb.solve()
+
+ tf = time.time()
+
+ print "\n"
+ print "Total time : ", tf - t0, "sec (CPU)"
+ print "Init time : ", t1 - t0, "sec (CPU)"
+ print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+ run()
diff --git a/Examples/main_Shear_2D.py b/Examples/main_Shear_2D.py
index 1b0a7e968ac18e534d36d82f05c4296bd0ad4004..dd1ac45cfb4a563cd5015adfe04a043e5dcec46e 100644
--- a/Examples/main_Shear_2D.py
+++ b/Examples/main_Shear_2D.py
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
import time
-import new_ParMePy as pp
+import parmepy as pp
import math
diff --git a/HySoP/hysop/__init__.py b/HySoP/hysop/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..0c1788e728d7c5938974a7efcbda0eb6fbc83769
--- /dev/null
+++ b/HySoP/hysop/__init__.py
@@ -0,0 +1,65 @@
+"""
+@package parmepy
+
+Python package dedicated to flow simulation using particular methods on hybrid architectures (MPI-GPU)
+
+"""
+__version__ = 1.00
+__all__ = ['Box', 'CartesianTopology', 'ScalarField']
+
+import os
+import site
+import mpi4py.MPI as MPI
+import pyopencl as cl
+
+rank_world = MPI.COMM_WORLD.Get_rank()
+if(rank_world == 0):
+ print "Starting Parmes version " + str(__version__) + ".\n"
+
+import constants
+
+import domain.box
+## Box-type physical domain
+Box = domain.box.Box
+
+## Cartesian grid
+
+## MPI topologies and associated meshes
+import domain.topology
+CartesianTopology = domain.topology.CartesianTopology
+
+## Fields
+import fields.discrete
+import fields.continuous
+import fields.analytical
+ScalarField = fields.discrete.ScalarField
+ContinuousField = fields.continuous.ContinuousField
+AnalyticalField = fields.analytical.AnalyticalField
+
+## Operators
+import operator.transport
+import operator.velocity
+import operator.stretching
+Transport = operator.transport.Transport
+Velocity = operator.velocity.Velocity
+Stretching = operator.stretching.Stretching
+
+## Problem
+import problem.problem
+Problem = problem.problem.Problem
+
+## Time integrator
+#import particular_solvers.integrator.integrator
+import particular_solvers.integrator.euler
+#ODESolver = particular_solvers.integrator.integrator.ODESolver
+Euler = particular_solvers.integrator.euler.Euler
+
+## Solver
+import particular_solvers.basic
+import particular_solvers.gpu
+ParticleSolver = particular_solvers.basic.ParticleSolver
+GPUParticleSolver = particular_solvers.gpu.GPUParticleSolver
+
+## Tools
+import tools.printer
+Printer = tools.printer.Printer
diff --git a/HySoP/hysop/operator/differentialOperator.py b/HySoP/hysop/operator/differentialOperator.py
new file mode 100755
index 0000000000000000000000000000000000000000..db3e5f715965f0a0ac55982f70b9c1cc693bad20
--- /dev/null
+++ b/HySoP/hysop/operator/differentialOperator.py
@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator.Stretching
+
+Penalization operator representation
+"""
+from .continuous import ContinuousOperator
+from .differentialOperator_d import DifferentialOperator_d
+
+
+class DifferentialOperator(ContinuousOperator):
+ """
+ Differential operator
+ """
+
+ def __init__(self, field1, field2, choice=None): # changer commentaires
+ """
+ Constructor.
+ Create a Differential operator from given velocity and vorticity variables.
+
+ @param velocity ContinuousVectorField : velocity variable.
+ @param vorticity ContinuousVectorField : vorticity variable.
+ """
+ ContinuousOperator.__init__(self)
+ self.field1 = field1
+ self.field2 = field2
+ self.choice = choice
+ self.addVariable([field1, field2])
+
+ def discretize(self, field1=None, field2=None, result=None): # changer commentaires
+ """
+ Stretching operator discretization method.
+ Create a discrete Stretching operator from given specifications.
+
+ @param idVelocityD : Index of velocity discretisation to use.
+ @param idVorticityD : Index of vorticity discretisation to use.
+ @param res_vorticity : result stretch (vorticity).
+ @param method : the method to use.
+ """
+ if self.discreteOperator is None:
+ self.discreteOperator = DifferentialOperator_d(self, field1, field2, result, self.choice)
+
+ def __str__(self):
+ """ToString method"""
+ s = " Differencial operator (ContinuousOperator)"
+ if self.discreteOperator is not None:
+ s += "with the following discretization:\n"
+ s += str(self.discreteOperator)
+ else:
+ s += ". Not discretised"
+ return s + "\n"
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : DifferentialOperator"
+ print DifferentialOperator.__doc__
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
new file mode 100755
index 0000000000000000000000000000000000000000..453faaa88fefb26d15bd09efc87422c6f5717594
--- /dev/null
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -0,0 +1,284 @@
+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete stretching representation
+"""
+from ..constants import *
+from .discrete import DiscreteOperator
+import numpy as np
+import time
+
+
+class DifferentialOperator_d(DiscreteOperator):
+ """
+ Operator representation.
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, opDif, field1, field2, result=None, choice=None):
+ """
+ Constructor.
+ Create a Stretching operator on a given continuous domain.
+ Work on a given field2 and field1 to compute the stretching term.
+
+ @param stretch : Stretching operator.
+ """
+ DiscreteOperator.__init__(self)
+ self.field1 = field1
+ self.field2 = field2
+ self.result = result
+ self.choice = choice
+ self.compute_time = 0.
+
+ def apply(self):
+ """
+ Apply operator.
+ """
+
+ if self.field2.topology != self.field1.topology:
+ print "Error, not the same topology for field1 and field2"
+ else:
+ dim = self.field1.topology.dim
+ coord_min = self.field1.topology.mesh.start
+ mesh_size = self.field1.topology.mesh.size
+ resolution = self.field1.topology.mesh.resolution
+ aux = [np.zeros((dim), dtype=PARMES_REAL, order=ORDER) for d in xrange(dim)]
+ #self.result = [np.zeros((resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(dim)]
+
+ if self.choice == 'divergenceProduct':
+
+# for k in xrange (resolution[2]):
+# kp =(k + 1 + resolution[2]) % resolution[2]
+# kpp=(k + 2 + resolution[2]) % resolution[2]
+# km =(k - 1 + resolution[2]) % resolution[2]
+# kmm=(k - 2 + resolution[2]) % resolution[2]
+#
+# for j in xrange (resolution[1]):
+# jp =(j + 1 + resolution[1]) % resolution[1]
+# jpp=(j + 2 + resolution[1]) % resolution[1]
+# jm =(j - 1 + resolution[1]) % resolution[1]
+# jmm=(j - 2 + resolution[1]) % resolution[1]
+#
+# for i in xrange (resolution[0]):
+# ip =(i + 1 + resolution[0]) % resolution[0]
+# ipp=(i + 2 + resolution[0]) % resolution[0]
+# im =(i - 1 + resolution[0]) % resolution[0]
+# imm=(i - 2 + resolution[0]) % resolution[0]
+#
+## X components of aux and result
+# aux[0][0] = 1.0 * self.field1[0][imm,j,k] * self.field2[0][imm,j,k] -\
+# 8.0 * self.field1[0][im ,j,k] * self.field2[0][im ,j,k] +\
+# 8.0 * self.field1[0][ip ,j,k] * self.field2[0][ip ,j,k] -\
+# 1.0 * self.field1[0][ipp,j,k] * self.field2[0][ipp,j,k]
+#
+# aux[0][1] = 1.0 * self.field1[1][i,jmm,k] * self.field2[0][i,jmm,k] -\
+# 8.0 * self.field1[1][i,jm,k] * self.field2[0][i,jm,k] +\
+# 8.0 * self.field1[1][i,jp,k] * self.field2[0][i,jp,k] -\
+# 1.0 * self.field1[1][i,jpp,k] * self.field2[0][i,jpp,k]
+#
+# aux[0][2] = 1.0 * self.field1[2][i,j,kmm] * self.field2[0][i,j,kmm] -\
+# 8.0 * self.field1[2][i,j,km] * self.field2[0][i,j,km] +\
+# 8.0 * self.field1[2][i,j,kp] * self.field2[0][i,j,kp] -\
+# 1.0 * self.field1[2][i,j,kpp] * self.field2[0][i,j,kpp]
+#
+# self.result[0][i,j,k]= aux[0][0] / (12. * mesh_size[0]) +\
+# aux[0][1] / (12. * mesh_size[1]) +\
+# aux[0][2] / (12. * mesh_size[2])
+#
+## Y components of aux and result
+# aux[1][0] = 1.0 * self.field1[0][imm,j,k] * self.field2[1][imm,j,k] -\
+# 8.0 * self.field1[0][im,j,k] * self.field2[1][im,j,k] +\
+# 8.0 * self.field1[0][ip,j,k] * self.field2[1][ip,j,k] -\
+# 1.0 * self.field1[0][ipp,j,k] * self.field2[1][ipp,j,k]
+#
+# aux[1][1] = 1.0 * self.field1[1][i,jmm,k] * self.field2[1][i,jmm,k] -\
+# 8.0 * self.field1[1][i,jm,k] * self.field2[1][i,jm,k] +\
+# 8.0 * self.field1[1][i,jp,k] * self.field2[1][i,jp,k] -\
+# 1.0 * self.field1[1][i,jpp,k] * self.field2[1][i,jpp,k]
+#
+# aux[1][2] = 1.0 * self.field1[2][i,j,kmm] * self.field2[1][i,j,kmm] -\
+# 8.0 * self.field1[2][i,j,km] * self.field2[1][i,j,km] +\
+# 8.0 * self.field1[2][i,j,kp] * self.field2[1][i,j,kp] -\
+# 1.0 * self.field1[2][i,j,kpp] * self.field2[1][i,j,kpp]
+#
+# self.result[1][i,j,k]= aux[1][0] / (12. * mesh_size[0]) +\
+# aux[1][1] / (12. * mesh_size[1]) +\
+# aux[1][2] / (12. * mesh_size[2])
+#
+## Z components of aux and result
+# aux[2][0] = 1.0 * self.field1[0][imm,j,k] * self.field2[2][imm,j,k] -\
+# 8.0 * self.field1[0][im,j,k] * self.field2[2][im,j,k] +\
+# 8.0 * self.field1[0][ip,j,k] * self.field2[2][ip,j,k] -\
+# 1.0 * self.field1[0][ipp,j,k] * self.field2[2][ipp,j,k]
+#
+# aux[2][1] = 1.0 * self.field1[1][i,jmm,k] * self.field2[2][i,jmm,k] -\
+# 8.0 * self.field1[1][i,jm,k] * self.field2[2][i,jm,k] +\
+# 8.0 * self.field1[1][i,jp,k] * self.field2[2][i,jp,k] -\
+# 1.0 * self.field1[1][i,jpp,k] * self.field2[2][i,jpp,k]
+#
+# aux[2][1] = 1.0 * self.field1[2][i,j,kmm] * self.field2[2][i,j,kmm] -\
+# 8.0 * self.field1[2][i,j,km] * self.field2[2][i,j,km] +\
+# 8.0 * self.field1[2][i,j,kp] * self.field2[2][i,j,kp] -\
+# 1.0 * self.field1[2][i,j,kpp] * self.field2[2][i,j,kpp]
+#
+# self.result[2][i,j,k]= aux[2][0] / (12. * mesh_size[0]) +\
+# aux[2][1] / (12. * mesh_size[1]) +\
+# aux[2][2] / (12. * mesh_size[2])
+
+
+# for k in xrange (resolution[2]):
+# kk =(k + resolution[2]) % resolution[2]
+# k2m=(k - 2 + resolution[2]) % resolution[2]
+# kp =(k + 1 + resolution[2]) % resolution[2]
+# k3m=(k - 3 + resolution[2]) % resolution[2]
+#
+# for j in xrange (resolution[1]):
+# jj =(j + resolution[1]) % resolution[1]
+# j2m=(j - 2 + resolution[1]) % resolution[1]
+# jp =(j + 1 + resolution[1]) % resolution[1]
+# j3m=(j - 3 + resolution[1]) % resolution[1]
+#
+# for i in xrange (resolution[0]):
+# ii =(i + resolution[0]) % resolution[0]
+# i2m=(i - 2 + resolution[0]) % resolution[0]
+# ip =(i + 1 + resolution[0]) % resolution[0]
+# i3m=(i - 3 + resolution[0]) % resolution[0]
+#
+## X components of aux and result
+# aux[0][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[0][ip,j,k] +\
+# 8.0 * self.field1[0][i2m ,j,k] * self.field2[0][i2m ,j,k] -\
+# 1.0 * self.field1[0][ii ,j,k] * self.field2[0][ii ,j,k] -\
+# 1.0 * self.field1[0][i3m,j,k] * self.field2[0][i3m,j,k]
+#
+# aux[0][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[0][i,jp,k] +\
+# 8.0 * self.field1[1][i,j2m,k] * self.field2[0][i,j2m,k] -\
+# 1.0 * self.field1[1][i,jj,k] * self.field2[0][i,jj,k] -\
+# 1.0 * self.field1[1][i,j3m,k] * self.field2[0][i,j3m,k]
+#
+# aux[0][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[0][i,j,kp] +\
+# 8.0 * self.field1[2][i,j,k2m] * self.field2[0][i,j,k2m] -\
+# 1.0 * self.field1[2][i,j,kk] * self.field2[0][i,j,kk] -\
+# 1.0 * self.field1[2][i,j,k3m] * self.field2[0][i,j,k3m]
+#
+# self.result[0][i,j,k]= aux[0][0] / (12. * mesh_size[0]) +\
+# aux[0][1] / (12. * mesh_size[1]) +\
+# aux[0][2] / (12. * mesh_size[2])
+#
+## Y components of aux and result
+# aux[1][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[1][ip,j,k] +\
+# 8.0 * self.field1[0][i2m ,j,k] * self.field2[1][i2m ,j,k] -\
+# 1.0 * self.field1[0][ii ,j,k] * self.field2[1][ii ,j,k] -\
+# 1.0 * self.field1[0][i3m,j,k] * self.field2[1][i3m,j,k]
+#
+# aux[1][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[1][i,jp,k] +\
+# 8.0 * self.field1[1][i,j2m,k] * self.field2[1][i,j2m,k] -\
+# 1.0 * self.field1[1][i,jj,k] * self.field2[1][i,jj,k] -\
+# 1.0 * self.field1[1][i,j3m,k] * self.field2[1][i,j3m,k]
+#
+# aux[1][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[1][i,j,kp] +\
+# 8.0 * self.field1[2][i,j,k2m] * self.field2[1][i,j,k2m] -\
+# 1.0 * self.field1[2][i,j,kk] * self.field2[1][i,j,kk] -\
+# 1.0 * self.field1[2][i,j,k3m] * self.field2[1][i,j,k3m]
+#
+# self.result[1][i,j,k]= aux[1][0] / (12. * mesh_size[0]) +\
+# aux[1][1] / (12. * mesh_size[1]) +\
+# aux[1][2] / (12. * mesh_size[2])
+#
+## Z components of aux and result
+# aux[2][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[2][ip,j,k] +\
+# 8.0 * self.field1[0][i2m ,j,k] * self.field2[2][i2m ,j,k] -\
+# 1.0 * self.field1[0][ii ,j,k] * self.field2[2][ii ,j,k] -\
+# 1.0 * self.field1[0][i3m,j,k] * self.field2[2][i3m,j,k]
+#
+# aux[2][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[2][i,jp,k] +\
+# 8.0 * self.field1[1][i,j2m,k] * self.field2[2][i,j2m,k] -\
+# 1.0 * self.field1[1][i,jj,k] * self.field2[2][i,jj,k] -\
+# 1.0 * self.field1[1][i,j3m,k] * self.field2[2][i,j3m,k]
+#
+# aux[2][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[2][i,j,kp] +\
+# 8.0 * self.field1[2][i,j,k2m] * self.field2[2][i,j,k2m] -\
+# 1.0 * self.field1[2][i,j,kk] * self.field2[2][i,j,kk] -\
+# 1.0 * self.field1[2][i,j,k3m] * self.field2[2][i,j,k3m]
+#
+# self.result[2][i,j,k]= aux[2][0] / (12. * mesh_size[0]) +\
+# aux[2][1] / (12. * mesh_size[1]) +\
+# aux[2][2] / (12. * mesh_size[2])
+#
+
+ for k in xrange (1,resolution[2]-1):
+ kp =(k + 1 + resolution[2]) % resolution[2]
+ km=(k - 1 + resolution[2]) % resolution[2]
+
+ for j in xrange (1,resolution[1]-1):
+ jm=(j - 1 + resolution[1]) % resolution[1]
+ jp =(j + 1 + resolution[1]) % resolution[1]
+
+ for i in xrange (1,resolution[0]-1):
+ im=(i - 1 + resolution[0]) % resolution[0]
+ ip =(i + 1 + resolution[0]) % resolution[0]
+
+# X components of aux and result
+ aux[0][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[0][ip,j,k] -\
+ 1.0 * self.field1[0][im ,j,k] * self.field2[0][im ,j,k]
+
+
+ aux[0][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[0][i,jp,k] -\
+ 1.0 * self.field1[1][i,jm,k] * self.field2[0][i,jm,k]
+
+ aux[0][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[0][i,j,kp] -\
+ 1.0 * self.field1[2][i,j,km] * self.field2[0][i,j,km]
+
+ self.result[0][i,j,k]= aux[0][0] / (2. * mesh_size[0]) +\
+ aux[0][1] / (2. * mesh_size[1]) +\
+ aux[0][2] / (2. * mesh_size[2])
+
+# Y components of aux and result
+ aux[1][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[1][ip,j,k] -\
+ 1.0 * self.field1[0][im,j,k] * self.field2[1][im,j,k]
+
+ aux[1][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[1][i,jp,k] -\
+ 1.0 * self.field1[1][i,jm,k] * self.field2[1][i,jm,k]
+
+ aux[1][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[1][i,j,kp] -\
+ 1.0 * self.field1[2][i,j,km] * self.field2[1][i,j,km]
+
+ self.result[1][i,j,k]= aux[1][0] / (2. * mesh_size[0]) +\
+ aux[1][1] / (2. * mesh_size[1]) +\
+ aux[1][2] / (2. * mesh_size[2])
+
+# Z components of aux and result
+ aux[2][0] = 1.0 * self.field1[0][ip,j,k] * self.field2[2][ip,j,k] -\
+ 1.0 * self.field1[0][im,j,k] * self.field2[2][im,j,k]
+
+ aux[2][1] = 1.0 * self.field1[1][i,jp,k] * self.field2[2][i,jp,k] -\
+ 1.0 * self.field1[1][i,jm,k] * self.field2[2][i,jm,k]
+
+ aux[2][2] = 1.0 * self.field1[2][i,j,kp] * self.field2[2][i,j,kp] -\
+ 1.0 * self.field1[2][i,j,km] * self.field2[2][i,j,km]
+
+ self.result[2][i,j,k]= aux[2][0] / (2. * mesh_size[0]) +\
+ aux[2][1] / (2. * mesh_size[1]) +\
+ aux[2][2] / (2. * mesh_size[2])
+
+
+ elif self.choice == 'curl':
+ pass
+ else:
+ raise ValueError("Unknown differiential operator: " + str(self.choice))
+
+ def printComputeTime(self):
+ self.timings_info[0] = "\"Differential Operator total\""
+ self.timings_info[1] = str(self.total_time)
+ self.timings_info[1] += str(self.compute_time[0]) + " " + str(self.compute_time[1]) + " " + str(self.compute_time[2]) + " "
+ print "Differential Operator total time : ", self.total_time
+ print "\t Differential Operator :", self.compute_time
+
+
+ def __str__(self):
+ s = "DifferentialOperator_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : DifferentialOperator_d"
+ print DifferentialOperator_d.__doc__
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
new file mode 100755
index 0000000000000000000000000000000000000000..21640cd26ed8b32aad4020e55123bd1664b58765
--- /dev/null
+++ b/HySoP/hysop/operator/stretching.py
@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator.Stretching
+
+Penalization operator representation
+"""
+from .continuous import ContinuousOperator
+from .stretching_d import Stretching_d
+
+
+class Stretching(ContinuousOperator):
+ """
+ Stretching operator representation
+ """
+
+ def __init__(self, velocity, vorticity):
+ """
+ Constructor.
+ Create a Stretching operator from given velocity and vorticity variables.
+
+ @param velocity ContinuousVectorField : velocity variable.
+ @param vorticity ContinuousVectorField : vorticity variable.
+ """
+ ContinuousOperator.__init__(self)
+ ## velocity variable (vector)
+ self.velocity = velocity
+ ## vorticity variable (vector)
+ self.vorticity = vorticity
+ self.addVariable([velocity, vorticity])
+
+ def discretize(self, idVelocityD=0, idVorticityD=0, res_vorticity=None, method=None): ## rajouter un attribut méthode (avec RK4, RK2 etc pour integration en temps) ?
+ """
+ Stretching operator discretization method.
+ Create a discrete Stretching operator from given specifications.
+
+ @param idVelocityD : Index of velocity discretisation to use.
+ @param idVorticityD : Index of vorticity discretisation to use.
+ @param res_vorticity : result stretch (vorticity).
+ @param method : the method to use.
+ """
+ if self.discreteOperator is None:
+ self.discreteOperator = Stretching_d(self, idVelocityD, idVorticityD, res_vorticity, method)
+
+ def __str__(self):
+ """ToString method"""
+ s = " Stretching operator (ContinuousOperator)"
+ if self.discreteOperator is not None:
+ s += "with the following discretization:\n"
+ s += str(self.discreteOperator)
+ else:
+ s += ". Not discretised"
+ return s + "\n"
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Stretching"
+ print Stretching.__doc__
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
new file mode 100755
index 0000000000000000000000000000000000000000..351e8dec8d95d74ddad8ed5e4ca40f95f1909faa
--- /dev/null
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete stretching representation
+"""
+from ..constants import *
+from .discrete import DiscreteOperator
+from .differentialOperator import DifferentialOperator
+import numpy as np
+import time
+
+
+class Stretching_d(DiscreteOperator):
+ """
+ Stretching operator representation.
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, stretch, idVelocityD=0, idVorticityD=0, res_vorticity=None, method=None):
+ """
+ Constructor.
+ Create a Stretching operator on a given continuous domain.
+ Work on a given velocity and vorticity to compute the stretching term.
+
+ @param stretch : Stretching operator.
+ """
+ DiscreteOperator.__init__(self)
+ ## Velocity.
+ self.velocity = stretch.velocity.discreteField[idVelocityD]
+ ## Vorticity.
+ self.vorticity = stretch.vorticity.discreteField[idVorticityD]
+ ## Result vorticity on the grid
+ self.res_vorticity = res_vorticity.discreteField[idVorticityD]
+ ## input fields
+ self.input = [self.velocity, self.vorticity]
+ ## output fields
+ self.output = [self.res_vorticity]
+ self.compute_time = 0.
+ self.method = method
+ ## self.name = "stretching"
+
+ def apply(self, t, dt):
+ """
+ Apply Stretching operator.
+
+ @param t : current time.
+ @param dt : time step.
+
+ Stretching algorithm:
+ @li 1. discretization of the divergence operator (div wu) : \n
+ @li 2. Time integration. Performs a RK4 resolution of dw/dt = div(wu).
+ -define the function of ODEsolver as the result of the first step
+ -integrate with the chosen method (RK4. RK2. euler)
+ """
+ #c_time = 0.
+ if self.numMethod is not None:
+ ## Space discretisation of div (wu)
+ self.stretchTerm = self.DifferencialOperator(stretch.vorticity, stretch.velocity, choice='divergenceProduct')
+ self.stretchTerm.discretize(self.vorticity, self.velocity, self.stretchTerm) # attention au self.stretchTerm ??
+ self.stretchTerm.discreteOperator.apply()
+
+ # Time integration
+ self.method(stretch)
+ for i in range(self.vorticity.topology.dim):
+ self.res_vorticity = self.method.integrate(self.vorticity, self.stretchTerm.discreteOperator.result, t, dt, i) ## method RK4
+ self.vorticity=np.copy(self.res_vorticity)
+
+ #c_time = (evt.profile.end - evt.profile.start) * 1e-9
+ #self.compute_time += c_time
+ #self.total_time += c_time
+ #return c_time
+
+ def printComputeTime(self):
+ self.timings_info[0] = "\"Stretching total\""
+ self.timings_info[1] = str(self.total_time)
+ self.timings_info[1] += str(self.compute_time[0]) + " " + str(self.compute_time[1]) + " " + str(self.compute_time[2]) + " "
+ print "Stretching total time : ", self.total_time
+ print "\t Stretching :", self.compute_time
+
+ def __str__(self):
+ s = "Stretching_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Stretching_d"
+ print Stretching_d.__doc__
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index f1da47c3173650da477ed913bde5dc2a6280c9df..ccf55f6852365909e5d4f09209648859c9b54328 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -4,8 +4,13 @@
Particular solver description.
"""
from solver import Solver
+#from integrator.runge_kutta4 import RK4
+#from integrator.integrator import ODESolver
+from integrator.euler import Euler
+#import integrator.euler.Euler as Euler
from ..fields.continuous import ContinuousField
from ..operator.transport import Transport
+from ..operator.stretching import Stretching
from ..operator.remeshing import Remeshing
from ..operator.splitting import Splitting
from ..tools.timer import Timer
@@ -44,7 +49,8 @@ class ParticleSolver(Solver):
## Advection operator of the problem. Advection need special treatment in particle methods.
self.advection = None
## ODE Solver method.
- self.ODESolver = ODESolver
+ #self.ODESolver = ODESolver
+ self.ODESolver = Euler
## Interpolation method.
self.InterpolationMethod = InterpolationMethod
## Remeshing method.
@@ -64,8 +70,12 @@ class ParticleSolver(Solver):
for op in self.problem.operators:
if isinstance(op, Transport):
self.advection = op
- if self.advection is None:
- raise ValueError("Particular Solver : Cannot create from a problem with no Transport operator")
+ if isinstance(op, Stretching):
+ self.stretch = op
+ #if isinstance(op, Velocity)
+ #self.velocity = op
+ #if self.advection is None:
+ # raise ValueError("Particular Solver : Cannot create from a problem with no Transport operator")
def initialize(self):
"""
@@ -73,8 +83,9 @@ class ParticleSolver(Solver):
Initialize a particle method solver regarding the problem.
"""
self.p_position = ContinuousField(domain=self.problem.domains[0], name='ParticlePosition')
+ self.p_vorticity = ContinuousField(domain=self.problem.domains[0], name='VortiStretch')
self.p_scalar = ContinuousField(domain=self.problem.domains[0], name='ParticleScalar')
- self.problem.addVariable([self.p_position, self.p_scalar])
+ self.problem.addVariable([self.p_position, self.p_scalar, self.p_vorticity])
## Discretise domains
for d in self.problem.domains:
d.discretize(self.problem.topology.resolution)
@@ -89,13 +100,19 @@ class ParticleSolver(Solver):
for op in self.problem.operators:
if op is self.advection:
op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
+ elif op is self.stretch:
+ op.discretize(res_vorticity=self.p_vorticity, method=self.ODESolver)
+ ## Velocity is only program on gpu for instance
+ #elif op is self.velocity:
+ # op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
else:
op.discretize()
- ## Create remeshing operator
- self.remeshing = Remeshing(partPositions=self.advection.discreteOperator.res_position,
- partScalar=self.advection.discreteOperator.res_scalar,
- resscal=self.advection.discreteOperator.scalar,
- method=self.RemeshingMethod)
+ if self.advection is not None:
+ ## Create remeshing operator
+ self.remeshing = Remeshing(partPositions=self.advection.discreteOperator.res_position,
+ partScalar=self.advection.discreteOperator.res_scalar,
+ resscal=self.advection.discreteOperator.scalar,
+ method=self.RemeshingMethod)
for op in self.problem.operators:
if op.needSplitting:
if op is self.advection:
diff --git a/HySoP/hysop/particular_solvers/integrator/__init__.py b/HySoP/hysop/particular_solvers/integrator/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..48c83c107ba4b8011ac4d47b2f8275f204530750 100644
--- a/HySoP/hysop/particular_solvers/integrator/__init__.py
+++ b/HySoP/hysop/particular_solvers/integrator/__init__.py
@@ -0,0 +1,6 @@
+"""
+@package parmepy.particular_solvers.integrator
+
+Everything concerning time integration.
+
+"""
diff --git a/HySoP/hysop/particular_solvers/integrator/euler.py b/HySoP/hysop/particular_solvers/integrator/euler.py
new file mode 100755
index 0000000000000000000000000000000000000000..5b94477b45bb3eaf7d3cefc11c9dae3ad417612a
--- /dev/null
+++ b/HySoP/hysop/particular_solvers/integrator/euler.py
@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+"""
+@package Utils
+Forward Euler method.
+"""
+from integrator import ODESolver
+
+
+#providedClass = "Euler"
+
+
+class Euler(ODESolver):
+ """
+ ODESolver implementation for solving an equation system with forward Euler method.
+ y'(t) = f(t,y)
+
+ y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
+
+ """
+ def __init__(self, f=None, dim=1):
+ """
+ Constructor.
+
+ @param f function f(t,y) : Right hand side of the equation to solve.
+ @param dim : dimensions
+ @param conditions : function to apply boundary conditions.
+ """
+ ODESolver.__init__(self, f)
+
+ def integrate(self, y, f, t, dt, direction):
+ """
+ Integration step for forward Euler method.
+ y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
+
+ @param y : position at time t.
+ @param t : current time.
+ @param dt : time step.
+ @return y at t+dt.
+ """
+ return [y[direction][...] + dt * self.f[direction][...] ]
+
+ def __str__(self):
+ """ToString method"""
+ return "Euler Method"
+
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : " + providedClass
+ print eval(providedClass).__doc__
diff --git a/HySoP/hysop/particular_solvers/integrator/integrator.py b/HySoP/hysop/particular_solvers/integrator/integrator.py
index 5974312a3fc6b29b8f96567dbfbf6e3fbff9d56f..f915d1a96532817cf05f7dc35d8974cba34b6147 100644
--- a/HySoP/hysop/particular_solvers/integrator/integrator.py
+++ b/HySoP/hysop/particular_solvers/integrator/integrator.py
@@ -3,7 +3,7 @@
@package Utils
ODESolver interface.
"""
-from ..Param import *
+from ...constants import *
class ODESolver:
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
new file mode 100755
index 0000000000000000000000000000000000000000..fe8f05794d54266987fa9485635af90eb6ac9856
--- /dev/null
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+"""
+@package Utils
+RK4 method interface.
+"""
+from ..Param import *
+from ODESolver import ODESolver
+
+
+class RK4(ODESolver):
+ """
+ ODESolver implementation for solving an equation system with RK2 method.
+ y'(t) = f(t,y)
+
+ y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
+
+ """
+ def __init__(self, f=None, conditions=lambda x: x, dim=3):
+ """
+ Constructor.
+
+ @param f function f(t,y) : Right hand side of the equation to solve.
+ @param dim : dimensions
+ @param conditions : function to apply boundary conditions.
+ """
+ ODESolver.__init__(self, f)
+ # Boundary conditions function.
+ self.boundaryConditions = conditions
+ self.gpu_kernel = """
+ """
+ self.nb_flop = 5
+
+ def integrate(self, f, yp , t, dt):
+ """
+ Integration step for RK4 Method
+ yp= y + dt*y'[y+dt/2*y'(y)].
+
+ @param y : position at time t.
+ @param fy : y'(y) value at time t.
+ @param yp : result position at time t + dt
+ @param t : current time.
+ @param dt : time step.
+ """
+ K1 = yp + dt * f(t, yp)
+ K2 = yp + dt * f(t + 0.5 * dt, (yp + K1) * 0.5)
+ K3 = yp + dt * f(t + 0.5 * dt, (yp + K2) * 0.5)
+ return (yp + (1/6.0) * dt * ( f(t, yp) +
+ 2. * f(t + 0.5 * dt, (yp + K1) * 0.5) +
+ 2. * f(t + 0.5 * dt, (yp + K2) * 0.5) +
+ f(t + dt, K3 )) )
+
+ def __str__(self):
+ """ToString method"""
+ return "RK4 Method"
+
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : RK2"
+ print RK4.__doc__
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 38e3520a1303a8df93fabedf8ea7df98e7caccf2..2f9c65cf30f39928bd0904d09c0cd2cd5c3812df 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -21,6 +21,7 @@ packages = ['parmepy',
'parmepy.operator',
'parmepy.problem',
'parmepy.particular_solvers',
+ 'parmepy.particular_solvers.integrator',
'parmepy.tools',
'parmepy.test',
'parmepy.test.test_domain',
diff --git a/HySoP/tests/test_Stretching/test_DivProduct.py b/HySoP/tests/test_Stretching/test_DivProduct.py
new file mode 100644
index 0000000000000000000000000000000000000000..a0a114f63a4c6061842574a045cbd2fc532d9891
--- /dev/null
+++ b/HySoP/tests/test_Stretching/test_DivProduct.py
@@ -0,0 +1,131 @@
+# -*- coding: utf-8 -*-
+import unittest
+
+from parmepy.problem.problem import *
+from parmepy.operator.transport import *
+from parmepy.operator.continuous import *
+from parmepy.operator.differentialOperator import *
+from parmepy.operator.stretching import *
+from parmepy.fields.discrete import *
+from parmepy.fields.continuous import *
+from parmepy.fields.analytical import *
+from parmepy.domain.topology import *
+from parmepy.domain.box import *
+from parmepy.constants import *
+from parmepy.particular_solvers.basic import *
+from parmepy.particular_solvers.integrator.euler import *
+from parmepy.particular_solvers.solver import *
+import numpy.testing as npt
+import math
+
+
+
+
+class test_DivProduct(unittest.TestCase):
+ """
+ DiscreteVariable test class
+ """
+ def setUp(self):
+ self.e = 0.2 # Accepted error between result and analytical result
+ self.dim = 3
+ self.boxLength = [1., 1., 1.]
+ self.boxMin = [0., 0., 0.]
+ self.nbElem = [32, 32, 32]
+ self.timeStep = 0.02
+ self.box = Box(dimension=self.dim,
+ length=self.boxLength,
+ origin=self.boxMin)
+
+ def testOperatorDiv(self):
+ # Continuous fields and operator declaration
+ self.velo = AnalyticalField(domain=self.box, formula=self.vitesse, name='Velocity', vector=True)
+ self.vorti = AnalyticalField(domain=self.box, formula=self.vorticite, name='Vorticity', vector=True)
+ self.div = DifferentialOperator(self.vorti, self.velo, choice='divergenceProduct')
+ # Topology definition / Fields and operator discretization
+ self.result = [np.ones((self.nbElem), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+ self.topo3D = CartesianTopology(domain=self.box, resolution=self.nbElem, dim=self.dim)
+ self.vorti.discretize(self.topo3D)
+ self.velo.discretize(self.topo3D)
+ self.vorti.initialize()
+ self.velo.initialize()
+ self.div.discretize(self.vorti.discreteField[0], self.velo.discreteField[0], result=self.result)# a changer
+ self.div.discreteOperator.apply()
+ self.resol = self.topo3D.mesh.resolution
+ self.FinalTime = 0.2
+ self.anal=self.analyticalDivProduct(self.topo3D.mesh.coords[0][:,0,0], \
+ self.topo3D.mesh.coords[1][0,:,0], \
+ self.topo3D.mesh.coords[2][0,0,:])
+ npt.assert_almost_equal(self.anal,self.result,decimal=8)
+# # Comparison with analytical solution
+# aux = np.zeros((4), dtype=PARMES_REAL, order=ORDER)
+# maxVal = 0.
+# for k in xrange (1,self.resol[2]-1):
+# for j in xrange (1,self.resol[1]-1):
+# for i in xrange (1,self.resol[0]-1):
+# aux[0] = abs(self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+# self.topo3D.mesh.coords[1][0,j,0], \
+# self.topo3D.mesh.coords[2][0,0,k])[0] - self.result[0][i,j,k])
+# aux[1] = abs(self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+# self.topo3D.mesh.coords[1][0,j,0], \
+# self.topo3D.mesh.coords[2][0,0,k])[1] - self.result[1][i,j,k])
+# aux[2] = abs(self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+# self.topo3D.mesh.coords[1][0,j,0], \
+# self.topo3D.mesh.coords[2][0,0,k])[2] - self.result[2][i,j,k])
+# aux[3]= aux[0]+aux[1]+aux[2]
+# #print "aux : ", aux[3], "i,j,k", i,j,k, "a0=",aux[0], "a1=",aux[1],"a2=",aux[2],"x=",self.topo3D.mesh.coords[0][i,0,0],"y=", self.topo3D.mesh.coords[1][0,j,0],"z=", self.topo3D.mesh.coords[2][0,0,k],"\n"
+## print "res0 : ", self.result[0][i,j,k], "res1 : ", self.result[1][i,j,k], "res2 : ", self.result[2][i,j,k],"\n"
+## print "Anal0 : ", self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+## self.topo3D.mesh.coords[1][0,j,0], \
+## self.topo3D.mesh.coords[2][0,0,k])[0], "Anal1 : ", self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+## self.topo3D.mesh.coords[1][0,j,0], \
+## self.topo3D.mesh.coords[2][0,0,k])[1], "Anal2 : ", self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+## self.topo3D.mesh.coords[1][0,j,0], \
+## self.topo3D.mesh.coords[2][0,0,k])[2],"\n"
+## print "a0=",aux[0], "a1=",aux[1],"a2=",aux[2],"\n"
+# if (aux[3] > maxVal) :
+# maxVal = aux[3]
+# self.assertLess(maxVal, self.e)
+## if abs(self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+## self.topo3D.mesh.coords[1][0,j,0], \
+## self.topo3D.mesh.coords[2][0,0,k]) - self.res[i,j,k,:]) > maxVal:
+## aux[:] = abs(self.analyticalDivProduct(self.topo3D.mesh.coords[0][i,0,0], \
+## self.topo3D.mesh.coords[1][0,j,0], \
+## self.topo3D.mesh.coords[2][0,0,k]) - self.res[i,j,k,:])
+##
+## maxVal = aux[0] + aux[1] + aux[2]
+## self.assertLess(maxVal, self.e)
+ print "Max error : ", maxVal
+
+ def vitesse(self, x, y, z):
+ vx = 1. + x
+ vy = - x * y
+ vz = x * y * z + 10.
+ return vx, vy, vz
+
+ def vorticite(self, x, y, z):
+ wx = x * y
+ wy = y * z
+ wz = - y
+ return wx, wy, wz
+
+ def analyticalDivProduct(self, x, y, z):
+ sx = y*(1.+2.*x) + z *(1.+x)
+ sy = -2.*x*y*(z+y)
+ sz = x*y*(2.*y*z + 2.*z*z - y) + 10.*(y+z)
+ return sx, sy, sz
+
+ def runTest(self):
+ self.setUp()
+ self.testOperatorDiv()
+
+
+def suite():
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(test_DivProduct))
+ return suite
+
+if __name__ == "__main__":
+ unittest.TextTestRunner(verbosity=2).run(suite())
+ #t = test_DivProduct()
+ #t.setUp()
+ #t.testOperatorDiv()