diff --git a/Examples/NSDebug_faux2D.py b/Examples/NSDebug_faux2D.py
new file mode 100755
index 0000000000000000000000000000000000000000..9479608ce5ca919c87ada9eb42322636493d0b27
--- /dev/null
+++ b/Examples/NSDebug_faux2D.py
@@ -0,0 +1,439 @@
+#!/usr/bin/python
+
+"""
+Navier Stokes 3D : flow past bluff bodies (using penalization).
+
+All parameters are set and defined in python module dataNS_bb.
+
+"""
+
+import parmepy as pp
+from parmepy.f2py import fftw2py
+import numpy as np
+from parmepy.problem.simulation import Simulation
+from parmepy.domain.obstacle.controlBox import ControlBox
+from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
+ Remesh, Support, Splitting, dtCrit, SpaceDiscretisation, GhostUpdate
+from parmepy.numerics.integrators.runge_kutta2 import RK2 as RK2
+from parmepy.numerics.integrators.runge_kutta3 import RK3 as RK3
+from parmepy.numerics.integrators.runge_kutta4 import RK4 as RK4
+from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
+from parmepy.numerics.interpolation import Linear
+from parmepy.numerics.remeshing import L6_4 as rmsh
+import parmepy.tools.numpywrappers as npw
+from parmepy.constants import HDF5
+import parmepy.tools.io_utils as io
+
+
+## ----------- A 3d problem -----------
+print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
+
+## pi constant
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+# ====== Flow constants =======
+uinf = 1.0
+VISCOSITY = 1. / 300.
+
+# ========= Geometry of the domain =========
+dim = 3
+boxlength = npw.realarray([6.0, 6.0, 0.08])
+boxorigin = npw.realarray([-2.0, -3.0, -0.04])
+#boxlength = npw.realarray([5.12, 5.12, 5.12])
+#boxorigin = npw.realarray([-2.56, -2.56, -2.56])
+# The domain
+box = pp.Box(dim, length=boxlength, origin=boxorigin)
+
+# Sphere inside the domain
+RADIUS = 0.5
+obst_pos = [0., 0., 0.]
+from parmepy.domain.obstacle.sphere import Sphere, HemiSphere
+from parmepy.domain.obstacle.cylinder import Cylinder, HemiCylinder
+#sphere = Sphere(box, position=obst_pos, radius=RADIUS)
+sphere = Cylinder(box, position=obst_pos, radius=RADIUS)
+
+# ========= Discretisation parameters =========
+nbElem = [301, 301, 5]
+#nbElem = [65, 65, 65]
+
+# ========= Vector Fields and variable parameters =========
+
+# Function to compute initial velocity
+def computeVel(res, x, y, z, t):
+ res[0][...] = uinf
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+# Function to compute initial vorticity
+def computeVort(res, x, y, z, t):
+ res[0][...] = 0.
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+from parmepy import VariableParameter, Field
+velo = Field(domain=box, formula=computeVel,
+ name='Velocity', isVector=True)
+vorti = Field(domain=box, formula=computeVort,
+ name='Vorticity', isVector=True)
+dt = VariableParameter(data=0.0125, name='dt')
+
+# ========= Operators that describe the problem =========
+op = {}
+
+# Adaptative timestep method : dt = min(values(dtCrit))
+# where dtCrit is a list of criterions on which the computation
+# of the adaptative time step is based
+# ex : dtCrit = ['gradU', 'cfl', 'stretch'], means :
+# dt = min (dtAdv, dtCfl, dtStretch), where dtAdv is equal to LCFL / |gradU|
+# For dtAdv, the possible choices are the following:
+# 'vort' (infinite norm of vorticity) : dtAdv = LCFL / |vort|
+# 'gradU' (infinite norm of velocity gradient), dtAdv = LCFL / |gradU|
+# 'deform' (infinite norm of deformation tensor),
+# dtAdv = LCFL / (0.5(gradU + gradU^T))
+# The operator is defined on a Cartesian topology with ghost points, that will
+# be dedicated to operators using finite differences.
+from parmepy.operator.adapt_timestep import AdaptTimeStep
+from parmepy.mpi.topology import Cartesian
+NBGHOSTS = 2
+ghosts = np.ones((box.dimension)) * NBGHOSTS
+topostr = Cartesian(box, box.dimension, nbElem, ghosts=ghosts)
+
+op['dtAdapt'] = AdaptTimeStep(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ dt_adapt=dt,
+ method={TimeIntegrator: RK3,
+ SpaceDiscretisation: FD_C_4,
+ dtCrit: ['vort', 'stretch']},
+ topo=topostr,
+ io_params={},
+ lcfl=0.125,
+ cfl=0.5)
+
+from parmepy.operator.advection import Advection
+op['advection'] = Advection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ method={Scales: 'p_M4',
+ Splitting: 'classic'}) # Scales advection
+
+from parmepy.operator.stretching import Stretching
+op['stretching'] = Stretching(velo, vorti,
+ resolutions={velo: nbElem, vorti: nbElem},
+ topo=topostr)
+
+from parmepy.operator.diffusion import Diffusion
+op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
+
+from parmepy.operator.poisson import Poisson
+op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem})
+
+from parmepy.operator.differential import Curl
+op['curl'] = Curl(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
+ method={SpaceDiscretisation: FD_C_4})
+
+# Discretisation of computational operators and link
+# to the associated topologies.
+for ope in op.values():
+ ope.discretize()
+topofft = op['poisson'].discreteFields[op['poisson'].vorticity].topology
+topocurl = op['curl'].discreteFields[op['curl'].invar].topology
+topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
+topostr = op['stretching'].discreteFields[op['stretching'].velocity].topology
+
+# Penalization of the velocity on a sphere inside the domain.
+# We enforce penalization on the same topology as for fft operators.
+from parmepy.operator.penalization import Penalization
+op['penalization'] = Penalization(velo, [sphere], coeff=[1e8],
+ topo=topofft, resolutions={velo: nbElem})
+op['penalization'].discretize()
+
+# Kill the vorticity at the inlet.
+# We enforce penalization on the same topology as for fft operators.
+ref_step_fft = topofft.mesh.space_step
+ibpos = npw.zeros(dim)
+iblength = npw.zeros(dim)
+ibpos[...] = boxorigin[...]
+iblength[...] = boxlength[...]
+iblength[0] = 10 * ref_step_fft[0]
+inlet_band = ControlBox(box, ibpos, iblength)
+op['kill_vort'] = Penalization(vorti, [inlet_band], coeff=[1e10],
+ topo=topofft, resolutions={vorti: nbElem})
+op['kill_vort'].discretize()
+
+# Correction of the velocity, used as post-process for Poisson solver,
+# based on a fixed flowrate through the entrance of the domain.
+# Time-dependant required-flowrate (Variable Parameter)
+def computeFlowrate(simu):
+ # === Time-dependant flow rate ===
+ t = simu.tk
+ Tstart = 3.0
+ flowrate = np.zeros(3)
+ flowrate[0] = uinf * box.length[1] * box.length[2]
+ if t >= Tstart and t <= Tstart + 1.0:
+ flowrate[1] = sin(pi * (t - Tstart)) * \
+ box.length[1] * box.length[2]
+ # === Constant flow rate ===
+ # flowrate = np.zeros(3)
+ # flowrate[0] = uinf * box.length[1] * box.length[2]
+ return flowrate
+req_flowrate = VariableParameter(formula=computeFlowrate)
+
+frate = npw.zeros(3)
+frate[0] = uinf * box.length[1] * box.length[2]
+req_flowrate = VariableParameter(data=frate, name='flowRate')
+
+from parmepy.operator.velocity_correction import VelocityCorrection
+op['correction'] = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ req_flowrate=req_flowrate, topo=topofft)
+op['correction'].discretize()
+
+# ========= Bridges between the different topologies =========
+from parmepy.operator.redistribute import Redistribute
+distr = {}
+distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
+distr['fft2advec'] = Redistribute([velo, vorti],
+ op['poisson'], op['advection'])
+distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
+distr['adv2str'] = Redistribute([velo, vorti],
+ op['advection'], op['stretching'])
+distr['str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
+distr['curl2adv'] = Redistribute([velo, vorti], op['curl'], op['advection'])
+distr['curl2str'] = Redistribute([velo, vorti], op['curl'], op['stretching'])
+distr['fft2str'] = Redistribute([velo, vorti], op['poisson'], op['stretching'])
+distr['str2curl'] = Redistribute([velo], op['stretching'], op['curl'])
+for ope in distr.values():
+ ope.discretize()
+
+# ========= Monitoring operators =========
+
+from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
+ Reprojection_criterion
+monitors = {}
+monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
+ frequency=100, formattype=HDF5, prefix='fields',
+ xmfalways=True)
+
+monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
+ formattype=HDF5, prefix='curl_io')
+
+monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
+ formattype=HDF5, prefix='advec_io',
+ xmfalways=True)
+
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
+ io_params={},
+ viscosity=VISCOSITY, isNormalized=False)
+
+reprojection_constant = 0.04
+reprojection_rate = 1
+monitors["reprojection"] = Reprojection_criterion(vorti, reprojection_constant,
+ reprojection_rate, topofft,
+ io_params={})
+# Add reprojection for poisson operator
+#op["poisson"].activateProjection(monitors["reprojection"])
+
+# Compute and save drag and lift on a
+# 3D Control box
+ref_step = topostr.mesh.space_step
+cbpos = npw.zeros(dim)
+cblength = npw.zeros(dim)
+cbpos[...] = boxorigin[...]
+cbpos[0] += 10 * ref_step[0]
+cblength[...] = boxlength[...]
+cblength[0] -= 20 * ref_step[0]
+cb = ControlBox(box, cbpos, cblength)
+coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
+
+monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+ topostr, cb, obstacles=[sphere], io_params={})
+
+step_dir = ref_step[0]
+#thickSliceXY = ControlBox(box, origin=[-2.56, -2.56, -2.0 * step_dir],
+# lengths=[5.12, 5.12, 4.0 * step_dir])
+#monitors['printerSliceXY'] = Printer(variables=[velo, vorti], topo=topofft,
+# frequency=1, formattype=HDF5, prefix='sliceXY',
+# subset=thickSliceXY, xmfalways=True)
+
+#thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
+# lengths=[10.24, 4.0 * step_dir, 5.12])
+#monitors['printerSliceXZ'] = Printer(variables=[velo, vorti], topo=topofft,
+# frequency=100, formattype=HDF5, prefix='sliceXZ',
+# subset=thickSliceXZ, xmfalways=True)
+
+#subBox = ControlBox(box, origin=[-0.7, -1.25, -1.25], lengths=[7.0, 2.5, 2.5])
+#monitors['printerSubBox'] = Printer(variables=[velo, vorti], topo=topofft,
+# frequency=100, formattype=HDF5, prefix='subBox',
+# subset=subBox, xmfalways=True)
+
+# ========= Setup for all declared operators =========
+for ope in op.values():
+ ope.setUp()
+for ope in distr.values():
+ ope.setUp()
+# Set up for monitors
+for monit in monitors.values():
+ monit.setUp()
+
+allops = dict(op.items() + distr.items() + monitors.items())
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=100., timeStep=0.0125, iterMax=1000000)
+#simu = Simulation(tinit=0.0, tend=75., timeStep=dt['dt'], iterMax=6)
+
+ind = sphere.discretize(topofft)
+vdfft = velo.discreteFields[topofft].data
+wdfft = vorti.discreteFields[topofft]
+
+## ========= Fields initialization =========
+# Initialisation, mode 1:
+# - compute velocity on topofft
+# - penalize
+# - compute vorticity with curl
+def initFields_mode1():
+ # The velocity is initialized on the fftw topology
+ # penalized and distributed on curl topo
+ velo.initialize(topo=topofft)
+ op['penalization'].apply(simu)
+ distr['fft2curl'].apply(simu)
+ distr['fft2curl'].wait()
+ # Vorticity is initialized after penalization of velocity
+ op['curl'].apply(simu)
+ distr['curl2adv'].apply(simu)
+ distr['curl2adv'].wait()
+
+## Call initialization
+initFields_mode1()
+
+##### Check initialize process results #####
+norm_vel_fft = velo.norm(topofft)
+norm_vel_curl = velo.norm(topocurl)
+norm_vort_fft = vorti.norm(topofft)
+norm_vort_curl = vorti.norm(topocurl)
+
+#assert np.allclose(norm_vort_curl, norm_vort_fft)
+assert np.allclose(norm_vel_curl, norm_vel_fft)
+
+# Print initial state
+for mon in monitors.values():
+ mon.apply(simu)
+
+## Note Franck : tests init ok, same values on both topologies, for v and w.
+# === After init ===
+#
+# v init on topocurl (which may be equal to topofft)
+#
+# === Definition of the 'one-step' sequence of operators ===
+#
+# 1 - w = curl(v), topocurl
+# 2 - w = advection(v,w), topo-advec
+# 3 - w = stretching(v,w), topostr
+# 4 - w = diffusion(w), topofft
+# 5 - v = poisson(w), topofft
+# 6 - v = correction(v, w), topofft
+# 7 - v = penalization(v), topofft
+#
+# Required bridges :
+# 1 --> 2 : (v,w) on topo-advec
+# 2 --> 3 : (v,w) on topostr
+# 3 --> 4 : w on topofft
+# 7 --> 1 : v on topocurl
+
+#fullseq = ['advection', # in: v,w out: w, on topoadvec
+# 'adv2str', 'stretching', # in: v,w out: w on topostr
+# # in: w, out: v, w on topofft
+# 'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
+# 'penalization',
+# 'fft2curl', 'curl', 'forces',
+# 'curl2adv', 'energy', 'profiles', 'printerAdvec',
+# ]
+
+fullseq = ['stretching',
+ 'str2diff', 'diffusion', 'poisson', 'correction',
+ 'penalization',
+ 'fft2curl', 'curl',
+ 'curl2fft', 'poisson', 'correction',
+ 'fft2adv',
+ 'advection',
+ #'printerSliceXY', #'printerSliceXZ', 'printerSubBox',
+ #'energy',
+ 'adv2str',
+ 'forces',
+ #'dtAdapt'
+ ]
+
+cb2 = op['correction'].cb
+#cb.discretize(topo=topofft)
+def pseudo_norm(topo, var):
+ sloc = npw.zeros(3)
+ gloc = npw.zeros(3)
+ sl = topo.mesh.iCompute
+ for i in xrange(3):
+ sloc[i] = np.sum((var[i][sl]))
+ # sfft[i] = la.norm((vfft[i][slfft]))
+ topo.comm.Allreduce(sloc, gloc)
+ return gloc
+
+def run(sequence):
+ ## for name in sequence:
+ ## assert allops.keys().count(name) > 0 or distr.keys().count(name) > 0\
+ ## or monitors.keys().count(name) > 0, 'unknow key:' + name
+ ## allops[name].apply(simu)
+ op['advection'].apply(simu)
+# distr['adv2str'].apply(simu)
+# distr['adv2str'].wait()
+# op['stretching'].apply(simu)
+# distr['str2diff'].apply(simu)
+# distr['str2diff'].wait()
+ op['diffusion'].apply(simu)
+ op['kill_vort'].apply(simu)
+ op['poisson'].apply(simu)
+ op['correction'].apply(simu)
+ op['penalization'].apply(simu)
+ distr['fft2curl'].apply(simu)
+ distr['fft2curl'].wait()
+ op['curl'].apply(simu)
+ distr['curl2fft'].apply(simu)
+ distr['curl2fft'].wait()
+ op['poisson'].apply(simu)
+ op['correction'].apply(simu)
+ monitors['printerFFT'].apply(simu)
+ monitors['energy'].apply(simu)
+ distr['fft2advec'].apply(simu)
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ monitors['forces'].apply(simu)
+ distr['fft2advec'].wait()
+
+seq = fullseq
+
+simu.initialize()
+
+while not simu.isOver:
+ if topocurl.rank == 0:
+ simu.printState()
+ run(seq)
+ simu.advance()
+
+
+## print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
+
+# === Finalize for all declared operators ===
+
+# Note FP : bug in fft finalize. To be checked ...
+# --> clean_fftw must be called only once
+#for ope in op.values():
+# ope.finalize()
+## Clean memory buffers
+
+fftw2py.clean_fftw_solver(box.dimension)
+for ope in distr.values():
+ ope.finalize()
+for monit in monitors.values():
+ monit.finalize()