From 8014f1a9621ebc31d621dacbed19cc4117566f15 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Mon, 12 Mar 2012 17:04:35 +0000
Subject: [PATCH] Fix sequential code (back to grid.points array)

---
 parmepy/examples/mainJM.py                    |  4 +--
 parmepy/new_ParMePy/Operator/RemeshingDOp.py  |  8 +++--
 .../new_ParMePy/Utils/GPUParticularSolver.py  |  2 --
 parmepy/new_ParMePy/Utils/ParticularSolver.py | 34 ++++++++++++++-----
 4 files changed, 34 insertions(+), 14 deletions(-)

diff --git a/parmepy/examples/mainJM.py b/parmepy/examples/mainJM.py
index c405cf17d..018181072 100644
--- a/parmepy/examples/mainJM.py
+++ b/parmepy/examples/mainJM.py
@@ -17,7 +17,7 @@ def scalaire(x):
 def run():
     # Parameter
     dim = 3
-    nb = 128
+    nb = 32
     boxLength = [1., 1., 1.]
     boxMin = [0., 0., 0.]
     nbElem = [nb, nb, nb]
@@ -48,7 +48,7 @@ def run():
     ## cpu
     #solver = pp.ParticularSolver(ODESolver, InterpolationMethod, RemeshingMethod)
     ## gpu
-    solver = pp.GPUParticularSolver(ODESolver, InterpolationMethod, RemeshingMethod)
+    solver = pp.GPUParticularSolver(ODESolver, InterpolationMethod, RemeshingMethod, device_type='gpu')
     cTspPb.setSolver(solver)
     #cTspPb.setPrinter(pp.Printer(frequency=outputModulo, outputPrefix=outputFilePrefix, problem=cTspPb))
 
diff --git a/parmepy/new_ParMePy/Operator/RemeshingDOp.py b/parmepy/new_ParMePy/Operator/RemeshingDOp.py
index 79712bfb7..9efb9a8b3 100644
--- a/parmepy/new_ParMePy/Operator/RemeshingDOp.py
+++ b/parmepy/new_ParMePy/Operator/RemeshingDOp.py
@@ -17,7 +17,7 @@ class RemeshingDOp(DiscreteOperator):
     DiscreteOperator.DiscreteOperator specialization.
     """
 
-    def __init__(self, partPositions, partScalar, partBloc, partTag, resscal, method):
+    def __init__(self, partPositions, partScalar, resscal, method, partBloc=None, partTag=None):
         """
         Constructor.
         Create a Remeshing operator on a given discrete domain.
@@ -40,7 +40,11 @@ class RemeshingDOp(DiscreteOperator):
         self.ptag = partTag
         ## Result scalar
         self.resultScalar = resscal
-        self.addVariable([self.ppos, self.pscal, self.pbloc, self.ptag, self.resultScalar])
+        self.addVariable([self.ppos, self.pscal, self.resultScalar])
+        if not self.pbloc is None:
+            self.addVariable([self.pbloc])
+        if not self.ptag is None:
+            self.addVariable([self.ptag])
         self.numMethod = method
         self.gpu_kernel = None
         self.compute_time = 0.
diff --git a/parmepy/new_ParMePy/Utils/GPUParticularSolver.py b/parmepy/new_ParMePy/Utils/GPUParticularSolver.py
index 0f4d5b537..2afe12432 100644
--- a/parmepy/new_ParMePy/Utils/GPUParticularSolver.py
+++ b/parmepy/new_ParMePy/Utils/GPUParticularSolver.py
@@ -98,8 +98,6 @@ class GPUParticularSolver(ParticularSolver):
         ## Remeshing : particles positions, particles scalar, particles tag, particles bloc -> grid scalar
         remesh = RemeshingDOp(partPositions=self.ppos,
                               partScalar=self.pscal,
-                              partBloc=self.pscal,
-                              partTag=self.pscal,
                               resscal=discreteProblem.advecOp.scalar,
                               method=self.RemeshingMethod)
         discreteProblem.addOperator([remesh])
diff --git a/parmepy/new_ParMePy/Utils/ParticularSolver.py b/parmepy/new_ParMePy/Utils/ParticularSolver.py
index f7a880e05..77778aa34 100644
--- a/parmepy/new_ParMePy/Utils/ParticularSolver.py
+++ b/parmepy/new_ParMePy/Utils/ParticularSolver.py
@@ -48,13 +48,33 @@ class ParticularSolver(Solver):
         @param discreteProblem : Problem to initialize.
         """
         dim = discreteProblem.domains[0].dimension
-        self.parts = ParticleField(discreteProblem.domains[0])
-        self.ppos = DiscreteVariable(domain=self.parts, dimension=discreteProblem.advecOp.velocity.dimension)
+        self.grid = discreteProblem.domains[0]
+        self.gvelo = discreteProblem.advecOp.velocity
+        self.gscal = discreteProblem.advecOp.scalar
+        ### Create grid.points TODO: A enlever
+        self.grid.points = np.empty_like(self.gvelo.values)
+        if dim == 1:
+            for x in xrange(self.grid.elementNumber[0]):
+                self.grid.points[x] = self.grid.axes[0][x]
+        if dim == 2:
+            for x in xrange(self.grid.elementNumber[0]):
+                for y in xrange(self.grid.elementNumber[1]):
+                    self.grid.points[x, y][0] = self.grid.axes[0][x]
+                    self.grid.points[x, y][1] = self.grid.axes[1][y]
+        if dim == 3:
+            for x in xrange(self.grid.elementNumber[0]):
+                for y in xrange(self.grid.elementNumber[1]):
+                    for z in xrange(self.grid.elementNumber[2]):
+                        self.grid.points[x, y, z][0] = self.grid.axes[0][x]
+                        self.grid.points[x, y, z][1] = self.grid.axes[1][y]
+                        self.grid.points[x, y, z][2] = self.grid.axes[2][z]
+        self.parts = ParticleField(self.grid)
+        self.ppos = DiscreteVariable(domain=self.parts, dimension=self.gvelo.dimension)
         self.parts.setPositions(self.ppos)
-        self.pvelo = DiscreteVariable(domain=self.parts, dimension=discreteProblem.advecOp.velocity.dimension)
-        self.pscal = DiscreteVariable(domain=self.parts, dimension=discreteProblem.advecOp.scalar.dimension, initialValue=lambda x: 0., scalar=True)
-        self.pbloc = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
-        self.ptag = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
+        self.pvelo = DiscreteVariable(domain=self.parts, dimension=self.gvelo.dimension)
+        self.pscal = DiscreteVariable(domain=self.parts, dimension=self.gscal.dimension, scalar=True)
+        #self.pbloc = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
+        #self.ptag = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
         # Interpolation : particles positions, grid velocity -> particles velocity
         interpol = InterpolationDOp(position=self.ppos,
                                     velocityField=discreteProblem.advecOp.velocity,
@@ -70,8 +90,6 @@ class ParticularSolver(Solver):
         # Remeshing : particles positions, particles scalar, particles tag, particles bloc -> grid scalar
         remesh = RemeshingDOp(partPositions=self.ppos,
                               partScalar=self.pscal,
-                              partBloc=self.pbloc,
-                              partTag=self.ptag,
                               resscal=discreteProblem.advecOp.scalar,
                               method=self.RemeshingMethod)
         # Advection : grid position, grid velocity, grid scalar -> particles position, particles scalar
-- 
GitLab