unused plugin files removed

8 files changed, 0 insertions, 2297 deletions

diff --git a/source/blender/python/manta_pp/source/advection.cpp b/source/blender/python/manta_pp/source/advection.cpp
deleted file mode 100644
index 54e8820d307..00000000000
--- a/source/blender/python/manta_pp/source/advection.cpp
+++ /dev/null
@@ -1,324 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Plugins for pressure correction:
- * - solve_pressure
- *
- ******************************************************************************/
-
-#include "../util/vectorbase.h"
-#include "grid.h"
-#include "kernel.h"
-
-using namespace std;
-
-namespace Manta { 
-
-static inline bool isNotFluid(FlagGrid& flags, int i, int j, int k)
-{
-	if ( flags.isFluid(i,j,k)   ) return false;
-	if ( flags.isFluid(i-1,j,k) ) return false; 
-	if ( flags.isFluid(i,j-1,k) ) return false; 
-	if ( flags.is3D() ) {
-		if ( flags.isFluid(i,j,k-1) ) return false;
-	}
-	return true;
-}
-
-//! Semi-Lagrange interpolation kernel
-KERNEL(bnd=1) template<class T> 
-void SemiLagrange (FlagGrid& flags, MACGrid& vel, Grid<T>& dst, Grid<T>& src, Real dt, bool isLevelset) 
-{
-	if (flags.isObstacle(i,j,k)) {
-		dst(i,j,k) = 0;
-		return;
-	}
-	if (!isLevelset && isNotFluid(flags,i,j,k) ) {
-		dst(i,j,k) = src(i,j,k);
-		return;
-	}
-	
-	// SL traceback
-	Vec3 pos = Vec3(i+0.5f,j+0.5f,k+0.5f) - vel.getCentered(i,j,k) * dt;
-	dst(i,j,k) = src.getInterpolated(pos);
-}
-
-static inline bool isNotFluidMAC(FlagGrid& flags, int i, int j, int k)
-{
-	if ( flags.isFluid(i,j,k)   ) return false;
-	return true;
-}
-
-//! Semi-Lagrange interpolation kernel for MAC grids
-KERNEL(bnd=1)
-void SemiLagrangeMAC(FlagGrid& flags, MACGrid& vel, MACGrid& dst, MACGrid& src, Real dt) 
-{
-	if (flags.isObstacle(i,j,k)) {
-		dst(i,j,k) = 0;
-		return;
-	}
-	if ( isNotFluidMAC(flags,i,j,k) ) {
-		dst(i,j,k) = src(i,j,k);
-		return;
-	}
-	
-	// get currect velocity at MAC position
-	// no need to shift xpos etc. as lookup field is also shifted
-	Vec3 xpos = Vec3(i+0.5f,j+0.5f,k+0.5f) - vel.getAtMACX(i,j,k) * dt;
-	Real vx = src.getInterpolatedComponent<0>(xpos);
-	Vec3 ypos = Vec3(i+0.5f,j+0.5f,k+0.5f) - vel.getAtMACY(i,j,k) * dt;
-	Real vy = src.getInterpolatedComponent<1>(ypos);
-	Vec3 zpos = Vec3(i+0.5f,j+0.5f,k+0.5f) - vel.getAtMACZ(i,j,k) * dt;
-	Real vz = src.getInterpolatedComponent<2>(zpos);
-	
-	dst(i,j,k) = Vec3(vx,vy,vz);
-}
-
-//! Kernel: Correct based on forward and backward SL steps (for both centered & mac grids)
-KERNEL(idx) template<class T> 
-void MacCormackCorrect(FlagGrid& flags, Grid<T>& dst, Grid<T>& old, Grid<T>& fwd,  Grid<T>& bwd, 
-					   Real strength, bool isLevelSet, bool isMAC=false )
-{
-	// note, replacement for isNotFluidMAC and isNotFluid
-	bool skip = false;
-
-	if (!flags.isFluid(idx)) skip = true;
-	if(!isMAC) {
-	if( (idx>=flags.getStrideX()) && 	(!flags.isFluid(idx-flags.getStrideX()) )) skip = true; 
-	if( (idx>=flags.getStrideY()) && 	(!flags.isFluid(idx-flags.getStrideY()) )) skip = true; 
-	if ( flags.is3D() ) {
-		if( (idx>=flags.getStrideZ()) &&(!flags.isFluid(idx-flags.getStrideZ()) )) skip = true;
-	} }
-	if ( skip ) {
-		dst[idx] = isLevelSet ? fwd[idx] : (T)0.0;
-		return;
-	}
-	
-	// note, strenth of correction can be modified here
-	dst[idx] = fwd[idx] + strength * 0.5 * (old[idx] - bwd[idx]);
-}
-
-// Helper to collect min/max in a template
-template<class T> inline void getMinMax(T& minv, T& maxv, const T& val) {
-	if (val < minv) minv = val;
-	if (val > maxv) maxv = val;
-}
-template<> inline void getMinMax<Vec3>(Vec3& minv, Vec3& maxv, const Vec3& val) {
-	getMinMax(minv.x, maxv.x, val.x);
-	getMinMax(minv.y, maxv.y, val.y);
-	getMinMax(minv.z, maxv.z, val.z);
-}
-
-	
-//! Helper function for clamping non-mac grids
-template<class T>
-inline T doClampComponent(const Vec3i& upperClamp, Grid<T>& orig, T dst, const Vec3i& posFwd) {
-	// clamp forward lookup to grid
-	const int i0 = clamp(posFwd.x, 0, upperClamp.x-1);
-	const int j0 = clamp(posFwd.y, 0, upperClamp.y-1);
-	const int k0 = clamp(posFwd.z, 0, (orig.is3D() ? (upperClamp.z-1) : 1) );
-	const int i1 = i0+1, j1 = j0+1, k1= (orig.is3D() ? (k0+1) : k0);
-	
-	if (!orig.isInBounds(Vec3i(i0,j0,k0),1)) {
-	return dst;
-	}
-
-		// find min/max around fwd pos
-		T minv = orig(i0,j0,k0), maxv = minv;
-		getMinMax(minv, maxv, orig(i1,j0,k0));
-		getMinMax(minv, maxv, orig(i0,j1,k0));
-		getMinMax(minv, maxv, orig(i1,j1,k0));
-		getMinMax(minv, maxv, orig(i0,j0,k1));
-		getMinMax(minv, maxv, orig(i1,j0,k1));
-		getMinMax(minv, maxv, orig(i0,j1,k1));
-		getMinMax(minv, maxv, orig(i1,j1,k1)); 
-		
-		// write clamped value
-		return clamp(dst, minv, maxv);
-}
-	
-//! Helper function for clamping MAC grids
-template<int c> 
-inline Real doClampComponentMAC(const Vec3i& upperClamp, MACGrid& orig, Real dst, const Vec3i& posFwd) {
-	// clamp forward lookup to grid
-	const int i0 = clamp(posFwd.x, 0, upperClamp.x-1);
-	const int j0 = clamp(posFwd.y, 0, upperClamp.y-1);
-	const int k0 = clamp(posFwd.z, 0, (orig.is3D() ? (upperClamp.z-1) : 1) );
-	const int i1 = i0+1, j1 = j0+1, k1= (orig.is3D() ? (k0+1) : k0);
-	if (!orig.isInBounds(Vec3i(i0,j0,k0),1)) 
-		return dst;
-	
-	// find min/max around fwd pos
-	Real minv = orig(i0,j0,k0)[c], maxv = minv;
-	getMinMax(minv, maxv, orig(i1,j0,k0)[c]);
-	getMinMax(minv, maxv, orig(i0,j1,k0)[c]);
-	getMinMax(minv, maxv, orig(i1,j1,k0)[c]);
-	getMinMax(minv, maxv, orig(i0,j0,k1)[c]);
-	getMinMax(minv, maxv, orig(i1,j0,k1)[c]);
-	getMinMax(minv, maxv, orig(i0,j1,k1)[c]);
-	getMinMax(minv, maxv, orig(i1,j1,k1)[c]);
-	
-	return clamp(dst, minv, maxv);    
-}
-
-//! Kernel: Clamp obtained value to min/max in source area, and reset values that point out of grid or into boundaries
-//          (note - MAC grids are handled below)
-KERNEL(bnd=1) template<class T>
-void MacCormackClamp(FlagGrid& flags, MACGrid& vel, Grid<T>& dst, Grid<T>& orig, Grid<T>& fwd, Real dt)
-{
-	if (flags.isObstacle(i,j,k))
-		return;
-	if ( isNotFluid(flags,i,j,k) ) {
-		dst(i,j,k) = fwd(i,j,k);
-		return;
-	}
-
-	T     dval       = dst(i,j,k);
-	Vec3i upperClamp = flags.getSize() - 1;
-	
-	// lookup forward/backward
-	Vec3i posFwd = toVec3i( Vec3(i,j,k) - vel.getCentered(i,j,k) * dt );
-	Vec3i posBwd = toVec3i( Vec3(i,j,k) + vel.getCentered(i,j,k) * dt );
-	
-	dval = doClampComponent<T>(upperClamp, orig, dval, posFwd );
-	
-	// test if lookups point out of grid or into obstacle
-	if (posFwd.x < 0 || posFwd.y < 0 || posFwd.z < 0 ||
-		posBwd.x < 0 || posBwd.y < 0 || posBwd.z < 0 ||
-		posFwd.x > upperClamp.x || posFwd.y > upperClamp.y || ((posFwd.z > upperClamp.z)&&flags.is3D()) ||
-		posBwd.x > upperClamp.x || posBwd.y > upperClamp.y || ((posBwd.z > upperClamp.z)&&flags.is3D()) ||
-		flags.isObstacle(posFwd) || flags.isObstacle(posBwd) ) 
-	{        
-		dval = fwd(i,j,k);
-	}
-	dst(i,j,k) = dval;
-}
-
-//! Kernel: same as MacCormackClamp above, but specialized version for MAC grids
-KERNEL(bnd=1) 
-void MacCormackClampMAC (FlagGrid& flags, MACGrid& vel, MACGrid& dst, MACGrid& orig, MACGrid& fwd, Real dt)
-{
-	if (flags.isObstacle(i,j,k))
-		return;
-	if ( isNotFluidMAC(flags,i,j,k) ) {
-		dst(i,j,k) = fwd(i,j,k);
-		return;
-	}
-	
-	Vec3  pos(i,j,k);
-	Vec3  dval       = dst(i,j,k);
-	Vec3i upperClamp = flags.getSize() - 1;
-	
-	// get total fwd lookup
-	Vec3i posFwd = toVec3i( Vec3(i,j,k) - vel.getCentered(i,j,k) * dt );
-	Vec3i posBwd = toVec3i( Vec3(i,j,k) + vel.getCentered(i,j,k) * dt );
-	
-	// clamp individual components
-	dval.x = doClampComponentMAC<0>(upperClamp, orig, dval.x, toVec3i( pos - vel.getAtMACX(i,j,k) * dt) );
-	dval.y = doClampComponentMAC<1>(upperClamp, orig, dval.y, toVec3i( pos - vel.getAtMACY(i,j,k) * dt) );
-	dval.z = doClampComponentMAC<2>(upperClamp, orig, dval.z, toVec3i( pos - vel.getAtMACZ(i,j,k) * dt) );
-	
-	// test if lookups point out of grid or into obstacle
-	if (posFwd.x < 0 || posFwd.y < 0 || posFwd.z < 0 ||
-		posBwd.x < 0 || posBwd.y < 0 || posBwd.z < 0 ||
-		posFwd.x > upperClamp.x || posFwd.y > upperClamp.y || ((posFwd.z > upperClamp.z)&&flags.is3D()) ||
-		posBwd.x > upperClamp.x || posBwd.y > upperClamp.y || ((posBwd.z > upperClamp.z)&&flags.is3D()) 
-		//|| flags.isObstacle(posFwd) || flags.isObstacle(posBwd)  // note - this unfortunately introduces asymmetry... TODO update
-		) 
-	{        
-		dval = fwd(i,j,k);
-	}
- 
-	// writeback
-	dst(i,j,k) = dval;
-}
-
-//! template function for performing SL advection
-template<class GridType> 
-void fnAdvectSemiLagrange(FluidSolver* parent, FlagGrid& flags, MACGrid& vel, GridType& orig, int order, Real strength) {
-	typedef typename GridType::BASETYPE T;
-	
-	Real dt = parent->getDt();
-	bool levelset = orig.getType() & GridBase::TypeLevelset;
-	
-	// forward step
-	GridType fwd(parent);
-	SemiLagrange<T> (flags, vel, fwd, orig, dt, levelset);
-	
-	if (order == 1) {
-		orig.swap(fwd);
-	}
-	else if (order == 2) { // MacCormack
-		GridType bwd(parent);
-		GridType newGrid(parent);
-	
-		// bwd <- backwards step
-		SemiLagrange<T> (flags, vel, bwd, fwd, -dt, levelset);
-		
-		// newGrid <- compute correction
-		MacCormackCorrect<T> (flags, newGrid, orig, fwd, bwd, strength, levelset);
-		
-		// clamp values
-		MacCormackClamp<T> (flags, vel, newGrid, orig, fwd, dt);
-		
-		orig.swap(newGrid);
-	}
-}
-
-//! template function for performing SL advection: specialized version for MAC grids
-template<> 
-void fnAdvectSemiLagrange<MACGrid>(FluidSolver* parent, FlagGrid& flags, MACGrid& vel, MACGrid& orig, int order, Real strength) {
-	Real dt = parent->getDt();
-	
-	// forward step
-	MACGrid fwd(parent);    
-	SemiLagrangeMAC (flags, vel, fwd, orig, dt);
-	
-	if (order == 1) {
-		orig.swap(fwd);
-	}
-	else if (order == 2) { // MacCormack
-		MACGrid bwd(parent);
-		MACGrid newGrid(parent);
-		
-		// bwd <- backwards step
-		SemiLagrangeMAC (flags, vel, bwd, fwd, -dt);
-		
-		// newGrid <- compute correction
-		MacCormackCorrect<Vec3> (flags, newGrid, orig, fwd, bwd, strength, false, true);
-		
-		// clamp values
-		MacCormackClampMAC (flags, vel, newGrid, orig, fwd, dt);
-		
-		orig.swap(newGrid);
-	}
-}
-
-//! Perform semi-lagrangian advection of target Real- or Vec3 grid
-PYTHON void advectSemiLagrange (FlagGrid* flags, MACGrid* vel, GridBase* grid, 
-						   int order = 1, Real strength = 1.0)
-{    
-	assertMsg(order==1 || order==2, "AdvectSemiLagrange: Only order 1 (regular SL) and 2 (MacCormack) supported");
-	
-	// determine type of grid    
-	if (grid->getType() & GridBase::TypeReal) {
-		fnAdvectSemiLagrange< Grid<Real> >(flags->getParent(), *flags, *vel, *((Grid<Real>*) grid), order, strength);
-	}
-	else if (grid->getType() & GridBase::TypeMAC) {    
-		fnAdvectSemiLagrange< MACGrid >(flags->getParent(), *flags, *vel, *((MACGrid*) grid), order, strength);
-	}
-	else if (grid->getType() & GridBase::TypeVec3) {    
-		fnAdvectSemiLagrange< Grid<Vec3> >(flags->getParent(), *flags, *vel, *((Grid<Vec3>*) grid), order, strength);
-	}
-	else
-		errMsg("AdvectSemiLagrange: Grid Type is not supported (only Real, Vec3, MAC, Levelset)");    
-}
-
-} // end namespace DDF 
-
diff --git a/source/blender/python/manta_pp/source/extforces.cpp b/source/blender/python/manta_pp/source/extforces.cpp
deleted file mode 100644
index e761e0a6ccf..00000000000
--- a/source/blender/python/manta_pp/source/extforces.cpp
+++ /dev/null
@@ -1,144 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Set boundary conditions, gravity
- *
- ******************************************************************************/
-
-#include "vectorbase.h"
-#include "grid.h"
-#include "commonkernels.h"
-
-using namespace std;
-
-namespace Manta { 
-
-//! add Forces between fl/fl and fl/em cells
-KERNEL(bnd=1) void KnAddForceField(FlagGrid& flags, MACGrid& vel, Grid<Vec3>& force) {
-	bool curFluid = flags.isFluid(i,j,k);
-	bool curEmpty = flags.isEmpty(i,j,k);
-	if (!curFluid && !curEmpty) return;
-	
-	if (flags.isFluid(i-1,j,k) || (curFluid && flags.isEmpty(i-1,j,k))) 
-		vel(i,j,k).x += 0.5*(force(i-1,j,k).x + force(i,j,k).x);
-	if (flags.isFluid(i,j-1,k) || (curFluid && flags.isEmpty(i,j-1,k))) 
-		vel(i,j,k).y += 0.5*(force(i,j-1,k).y + force(i,j,k).y);
-	if (vel.is3D() && (flags.isFluid(i,j,k-1) || (curFluid && flags.isEmpty(i,j,k-1))))
-		vel(i,j,k).z += 0.5*(force(i,j,k-1).z + force(i,j,k).z);
-}
-
-//! add Forces between fl/fl and fl/em cells
-KERNEL(bnd=1) void KnAddForce(FlagGrid& flags, MACGrid& vel, Vec3 force) {
-	bool curFluid = flags.isFluid(i,j,k);
-	bool curEmpty = flags.isEmpty(i,j,k);
-	if (!curFluid && !curEmpty) return;
-	
-	if (flags.isFluid(i-1,j,k) || (curFluid && flags.isEmpty(i-1,j,k))) 
-		vel(i,j,k).x += force.x;
-	if (flags.isFluid(i,j-1,k) || (curFluid && flags.isEmpty(i,j-1,k))) 
-		vel(i,j,k).y += force.y;
-	if (vel.is3D() && (flags.isFluid(i,j,k-1) || (curFluid && flags.isEmpty(i,j,k-1))))
-		vel(i,j,k).z += force.z;
-}
-
-//! add gravity forces to all fluid cells
-PYTHON void addGravity(FlagGrid& flags, MACGrid& vel, Vec3 gravity) {    
-	Vec3 f = gravity * flags.getParent()->getDt() / flags.getDx();
-	KnAddForce(flags, vel, f);
-}
-
-//! add Buoyancy force based on smoke density
-KERNEL(bnd=1) void KnAddBuoyancy(FlagGrid& flags, Grid<Real>& density, MACGrid& vel, Vec3 strength) {    
-	if (!flags.isFluid(i,j,k)) return;
-	if (flags.isFluid(i-1,j,k))
-		vel(i,j,k).x += (0.5 * strength.x) * (density(i,j,k)+density(i-1,j,k));
-	if (flags.isFluid(i,j-1,k))
-		vel(i,j,k).y += (0.5 * strength.y) * (density(i,j,k)+density(i,j-1,k));
-	if (vel.is3D() && flags.isFluid(i,j,k-1))
-		vel(i,j,k).z += (0.5 * strength.z) * (density(i,j,k)+density(i,j,k-1));    
-}
-
-//! add Buoyancy force based on smoke density
-PYTHON void addBuoyancy(FlagGrid& flags, Grid<Real>& density, MACGrid& vel, Vec3 gravity) {
-	Vec3 f = - gravity * flags.getParent()->getDt() / flags.getParent()->getDx();
-	KnAddBuoyancy(flags,density, vel, f);
-}
-
-		
-//! set no-stick wall boundary condition between ob/fl and ob/ob cells
-KERNEL void KnSetWallBcs(FlagGrid& flags, MACGrid& vel) {
-	bool curFluid = flags.isFluid(i,j,k);
-	bool curObstacle = flags.isObstacle(i,j,k);
-	if (!curFluid && !curObstacle) return;
-	
-	// we use i>0 instead of bnd=1 to check outer wall
-	if (i>0 && (flags.isObstacle(i-1,j,k) || (curObstacle && flags.isFluid(i-1,j,k))))
-		vel(i,j,k).x = 0;
-	if (j>0 && (flags.isObstacle(i,j-1,k) || (curObstacle && flags.isFluid(i,j-1,k))))
-		vel(i,j,k).y = 0;
-	if (vel.is2D() || (k>0 && (flags.isObstacle(i,j,k-1) || (curObstacle && flags.isFluid(i,j,k-1)))))
-		vel(i,j,k).z = 0;
-		
-	if (curFluid) {
-		if ((i>0 && flags.isStick(i-1,j,k)) || (i<flags.getSizeX()-1 && flags.isStick(i+1,j,k)))
-			vel(i,j,k).y = vel(i,j,k).z = 0;
-		if ((j>0 && flags.isStick(i,j-1,k)) || (j<flags.getSizeY()-1 && flags.isStick(i,j+1,k)))
-			vel(i,j,k).x = vel(i,j,k).z = 0;
-		if (vel.is3D() && ((k>0 && flags.isStick(i,j,k-1)) || (k<flags.getSizeZ()-1 && flags.isStick(i,j,k+1))))
-			vel(i,j,k).x = vel(i,j,k).y = 0;
-	}
-}
-
-//! set no-stick boundary condition on walls
-PYTHON void setWallBcs(FlagGrid& flags, MACGrid& vel) {
-	KnSetWallBcs(flags, vel);
-} 
-
-//! Kernel: gradient norm operator
-KERNEL(bnd=1) void KnConfForce(Grid<Vec3>& force, const Grid<Real>& grid, const Grid<Vec3>& curl, Real str) {
-	Vec3 grad = 0.5 * Vec3(        grid(i+1,j,k)-grid(i-1,j,k), 
-								   grid(i,j+1,k)-grid(i,j-1,k), 0.);
-	if(grid.is3D()) grad[2]= 0.5*( grid(i,j,k+1)-grid(i,j,k-1) );
-	normalize(grad);
-	force(i,j,k) = str * cross(grad, curl(i,j,k));
-}
-
-PYTHON void vorticityConfinement(MACGrid& vel, FlagGrid& flags, Real strength) {
-	Grid<Vec3> velCenter(flags.getParent()), curl(flags.getParent()), force(flags.getParent());
-	Grid<Real> norm(flags.getParent());
-	
-	GetCentered(velCenter, vel);
-	CurlOp(velCenter, curl);
-	GridNorm(norm, curl);
-	KnConfForce(force, norm, curl, strength);
-	KnAddForceField(flags, vel, force);
-}
-
-//! enforce a constant inflow/outflow at the grid boundaries
-KERNEL void KnSetInflow(MACGrid& vel, int dim, int p0, const Vec3& val) {
-	Vec3i p(i,j,k);
-	if (p[dim] == p0 || p[dim] == p0+1)
-		vel(i,j,k) = val;
-}
-
-//! enforce a constant inflow/outflow at the grid boundaries
-PYTHON void setInflowBcs(MACGrid& vel, string dir, Vec3 value) {
-	for(size_t i=0; i<dir.size(); i++) {
-		if (dir[i] >= 'x' && dir[i] <= 'z') { 
-			int dim = dir[i]-'x';
-			KnSetInflow(vel,dim,0,value);
-		} else if (dir[i] >= 'X' && dir[i] <= 'Z') {
-			int dim = dir[i]-'X';
-			KnSetInflow(vel,dim,vel.getSize()[dim]-1,value);
-		} else 
-			errMsg("invalid character in direction string. Only [xyzXYZ] allowed.");
-	}
-}
-
-} // namespace
diff --git a/source/blender/python/manta_pp/source/initplugins.cpp b/source/blender/python/manta_pp/source/initplugins.cpp
deleted file mode 100644
index 69579a1327a..00000000000
--- a/source/blender/python/manta_pp/source/initplugins.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Tools to setup fields and inflows
- *
- ******************************************************************************/
-
-#include "vectorbase.h"
-#include "shapes.h"
-#include "commonkernels.h"
-#include "particle.h"
-#include "noisefield.h"
-#include "mesh.h"
-
-using namespace std;
-
-namespace Manta {
-	
-//! Apply noise to grid
-KERNEL 
-void KnApplyNoise(FlagGrid& flags, Grid<Real>& density, WaveletNoiseField& noise, Grid<Real>& sdf, Real scale, Real sigma) 
-{
-	if (!flags.isFluid(i,j,k) || sdf(i,j,k) > sigma) return;
-	Real factor = clamp(1.0-0.5/sigma * (sdf(i,j,k)+sigma), 0.0, 1.0);
-	
-	Real target = noise.evaluate(Vec3(i,j,k)) * scale * factor;
-	if (density(i,j,k) < target)
-		density(i,j,k) = target;
-}
-
-//! Init noise-modulated density inside shape
-PYTHON void densityInflow(FlagGrid& flags, Grid<Real>& density, WaveletNoiseField& noise, Shape* shape, Real scale=1.0, Real sigma=0)
-{
-	Grid<Real> sdf = shape->computeLevelset();
-	KnApplyNoise(flags, density, noise, sdf, scale, sigma);
-}
-
-	
-//! Init noise-modulated density inside mesh
-PYTHON void densityInflowMesh(FlagGrid& flags, Grid<Real>& density, WaveletNoiseField& noise, Mesh* mesh, Real scale=1.0, Real sigma=0)
-{
-	FluidSolver dummy(density.getSize());
-	LevelsetGrid sdf(&dummy, false);
-	mesh->meshSDF(*mesh, sdf, 1.,3);
-	KnApplyNoise(flags, density, noise, sdf, scale, sigma);
-}
-//! sample noise field and set pdata with its values (for convenience, scale the noise values)
-KERNEL(pts) template<class T>
-void knSetPdataNoise(BasicParticleSystem& parts, ParticleDataImpl<T>& pdata, WaveletNoiseField& noise, Real scale) {
-	pdata[idx] = noise.evaluate( parts.getPos(idx) ) * scale;
-}
-KERNEL(pts) template<class T>
-void knSetPdataNoiseVec(BasicParticleSystem& parts, ParticleDataImpl<T>& pdata, WaveletNoiseField& noise, Real scale) {
-	pdata[idx] = noise.evaluateVec( parts.getPos(idx) ) * scale;
-}
-PYTHON void setNoisePdata    (BasicParticleSystem& parts, ParticleDataImpl<Real>& pd, WaveletNoiseField& noise, Real scale=1.) { knSetPdataNoise<Real>(parts, pd,noise,scale); }
-PYTHON void setNoisePdataVec3(BasicParticleSystem& parts, ParticleDataImpl<Vec3>& pd, WaveletNoiseField& noise, Real scale=1.) { knSetPdataNoiseVec<Vec3>(parts, pd,noise,scale); }
-PYTHON void setNoisePdataInt (BasicParticleSystem& parts, ParticleDataImpl<int >& pd, WaveletNoiseField& noise, Real scale=1.) { knSetPdataNoise<int> (parts, pd,noise,scale); }
-
-//! SDF gradient from obstacle flags
-PYTHON Grid<Vec3> obstacleGradient(FlagGrid& flags) {
-	LevelsetGrid levelset(flags.getParent(),false);
-	Grid<Vec3> gradient(flags.getParent());
-	
-	// rebuild obstacle levelset
-	FOR_IDX(levelset) {
-		levelset[idx] = flags.isObstacle(idx) ? -0.5 : 0.5;
-	}
-	levelset.reinitMarching(flags, 6.0, 0, true, false, FlagGrid::TypeReserved);
-	
-	// build levelset gradient
-	GradientOp(gradient, levelset);
-	
-	FOR_IDX(levelset) {
-		Vec3 grad = gradient[idx];
-		Real s = normalize(grad);
-		if (s <= 0.1 || levelset[idx] >= 0) 
-			grad=Vec3(0.);        
-		gradient[idx] = grad * levelset[idx];
-	}
-	
-	return gradient;
-}
-
-PYTHON LevelsetGrid obstacleLevelset(FlagGrid& flags) {
-   LevelsetGrid levelset(flags.getParent(),false);
-	Grid<Vec3> gradient(flags.getParent());
-
-	// rebuild obstacle levelset
-	FOR_IDX(levelset) {
-		levelset[idx] = flags.isObstacle(idx) ? -0.5 : 0.5;
-	}
-	levelset.reinitMarching(flags, 6.0, 0, true, false, FlagGrid::TypeReserved);
-
-	return levelset;
-}    
-
-
-} // namespace
diff --git a/source/blender/python/manta_pp/source/kepsilon.cpp b/source/blender/python/manta_pp/source/kepsilon.cpp
deleted file mode 100644
index 955d05d53d3..00000000000
--- a/source/blender/python/manta_pp/source/kepsilon.cpp
+++ /dev/null
@@ -1,183 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Turbulence modeling plugins
- *
- ******************************************************************************/
- 
-#include "grid.h"
-#include "commonkernels.h"
-#include "vortexsheet.h"
-#include "conjugategrad.h"
-
-using namespace std;
-
-namespace Manta {
-
-// k-epsilon model constants
-const Real keCmu = 0.09;
-const Real keC1 = 1.44;
-const Real keC2 = 1.92;
-const Real keS1 = 1.0;
-const Real keS2 = 1.3;
-
-// k-epsilon limiters
-const Real keU0 = 1.0;
-const Real keImin = 2e-3;
-const Real keImax = 1.0;
-const Real keNuMin = 1e-3;
-const Real keNuMax = 5.0;
-
-//! clamp k and epsilon to limits    
-KERNEL(idx) 
-void KnTurbulenceClamp(Grid<Real>& kgrid, Grid<Real>& egrid, Real minK, Real maxK, Real minNu, Real maxNu) {
-	Real eps = egrid[idx];
-	Real ke = clamp(kgrid[idx],minK,maxK);
-	Real nu = keCmu*square(ke)/eps;
-	if (nu > maxNu) 
-		eps = keCmu*square(ke)/maxNu;
-	if (nu < minNu) 
-		eps = keCmu*square(ke)/minNu;
-
-	kgrid[idx] = ke;
-	egrid[idx] = eps;
-}
-
-//! Compute k-epsilon production term P = 2*nu_T*sum_ij(Sij^2) and the turbulent viscosity nu_T=C_mu*k^2/eps
-KERNEL(bnd=1) 
-void KnComputeProduction(const MACGrid& vel, const Grid<Vec3>& velCenter, const Grid<Real>& ke, const Grid<Real>& eps, 
-							   Grid<Real>& prod, Grid<Real>& nuT, Grid<Real>* strain, Real pscale = 1.0f) 
-{
-	Real curEps = eps(i,j,k);
-	if (curEps > 0) {
-		// turbulent viscosity: nu_T = C_mu * k^2/eps
-		Real curNu = keCmu * square(ke(i,j,k)) / curEps;
-		
-		// compute Sij = 1/2 * (dU_i/dx_j + dU_j/dx_i)
-		Vec3 diag = Vec3(vel(i+1,j,k).x, vel(i,j+1,k).y, vel(i,j,k+1).z) - vel(i,j,k);
-		Vec3 ux = 0.5*(velCenter(i+1,j,k)-velCenter(i-1,j,k));
-		Vec3 uy = 0.5*(velCenter(i,j+1,k)-velCenter(i,j-1,k));
-		Vec3 uz = 0.5*(velCenter(i,j,k+1)-velCenter(i,j,k-1));
-		Real S12 = 0.5*(ux.y+uy.x);
-		Real S13 = 0.5*(ux.z+uz.x);
-		Real S23 = 0.5*(uy.z+uz.y);
-		Real S2 = square(diag.x) + square(diag.y) + square(diag.z) +
-				  2.0*square(S12) + 2.0*square(S13) + 2.0*square(S23);
-		
-		// P = 2*nu_T*sum_ij(Sij^2)
-		prod(i,j,k) = 2.0 * curNu * S2 * pscale;
-		nuT(i,j,k) = curNu;
-		if (strain) (*strain)(i,j,k) = sqrt(S2);
-	} 
-	else {
-		prod(i,j,k) = 0;
-		nuT(i,j,k) = 0;
-		if (strain) (*strain)(i,j,k) = 0;
-	}
-}
-	
-//! Compute k-epsilon production term P = 2*nu_T*sum_ij(Sij^2) and the turbulent viscosity nu_T=C_mu*k^2/eps
-PYTHON void KEpsilonComputeProduction(MACGrid& vel, Grid<Real>& k, Grid<Real>& eps, Grid<Real>& prod, Grid<Real>& nuT, Grid<Real>* strain=0, Real pscale = 1.0f) 
-{
-	// get centered velocity grid
-	Grid<Vec3> vcenter(k.getParent());
-	GetCentered(vcenter, vel);
-	FillInBoundary(vcenter,1);
-	
-	// compute limits
-	const Real minK = 1.5*square(keU0)*square(keImin);
-	const Real maxK = 1.5*square(keU0)*square(keImax);    
-	KnTurbulenceClamp(k, eps, minK, maxK, keNuMin, keNuMax);
-	
-	KnComputeProduction(vel, vcenter, k, eps, prod, nuT, strain, pscale);    
-}
-
-//! Integrate source terms of k-epsilon equation
-KERNEL(idx) 
-void KnAddTurbulenceSource(Grid<Real>& kgrid, Grid<Real>& egrid, const Grid<Real>& pgrid, Real dt) {
-	Real eps = egrid[idx], prod = pgrid[idx], ke = kgrid[idx];
-	if (ke <= 0) ke = 1e-3; // pre-clamp to avoid nan
-	
-	Real newK = ke + dt*(prod - eps);
-	Real newEps = eps + dt*(prod * keC1 - eps * keC2) * (eps / ke);
-	if (newEps <= 0) newEps = 1e-4; // pre-clamp to avoid nan
-
-	kgrid[idx] = newK;
-	egrid[idx] = newEps;
-}
-
-
-//! Integrate source terms of k-epsilon equation
-PYTHON void KEpsilonSources(Grid<Real>& k, Grid<Real>& eps, Grid<Real>& prod) {
-	Real dt = k.getParent()->getDt();
-		
-	KnAddTurbulenceSource(k, eps, prod, dt);
-	
-	// compute limits
-	const Real minK = 1.5*square(keU0)*square(keImin);
-	const Real maxK = 1.5*square(keU0)*square(keImax);
-	KnTurbulenceClamp(k, eps, minK, maxK, keNuMin, keNuMax);    
-}
-
-//! Initialize the domain or boundary conditions
-PYTHON void KEpsilonBcs(FlagGrid& flags, Grid<Real>& k, Grid<Real>& eps, Real intensity, Real nu, bool fillArea) {
-	// compute limits
-	const Real vk = 1.5*square(keU0)*square(intensity);
-	const Real ve = keCmu*square(vk) / nu;
-	
-	FOR_IDX(k) {
-		if (fillArea || flags.isObstacle(idx)) {
-			k[idx] = vk;
-			eps[idx] = ve;
-		}
-	}
-}
-
-//! Gradient diffusion smoothing. Not unconditionally stable -- should probably do substepping etc.
-void ApplyGradDiff(const Grid<Real>& grid, Grid<Real>& res, const Grid<Real>& nu, Real dt, Real sigma) {
-	// should do this (but requires better boundary handling)
-	/*MACGrid grad(grid.getParent());
-	GradientOpMAC(grad, grid);
-	grad *= nu;
-	DivergenceOpMAC(res, grad);
-	res *= dt/sigma;  */
-	
-	LaplaceOp(res, grid);
-	res *= nu;
-	res *= dt/sigma;
-}
-
-//! Compute k-epsilon turbulent viscosity
-PYTHON void KEpsilonGradientDiffusion(Grid<Real>& k, Grid<Real>& eps, Grid<Real>& nuT, Real sigmaU=4.0, MACGrid* vel=0) {
-	Real dt = k.getParent()->getDt();
-	Grid<Real> res(k.getParent());
-	
-	// gradient diffusion of k
-	ApplyGradDiff(k, res, nuT, dt, keS1);
-	k += res;
-
-	// gradient diffusion of epsilon
-	ApplyGradDiff(eps, res, nuT, dt, keS2);
-	eps += res;
-	
-	// gradient diffusion of velocity
-	if (vel) {
-		Grid<Real> vc(k.getParent());
-		for (int c=0; c<3; c++) {
-			GetComponent(*vel, vc, c);
-			ApplyGradDiff(vc, res, nuT, dt, sigmaU);
-			vc += res;
-			SetComponent(*vel, vc, c);    
-		}
-	}
-}
-
-
-
-} // namespace
\ No newline at end of file
diff --git a/source/blender/python/manta_pp/source/meshplugins.cpp b/source/blender/python/manta_pp/source/meshplugins.cpp
deleted file mode 100644
index 8037447d299..00000000000
--- a/source/blender/python/manta_pp/source/meshplugins.cpp
+++ /dev/null
@@ -1,623 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Smoothing etc. for meshes
- *
- ******************************************************************************/
-
-/******************************************************************************/
-// Copyright note:
-//
-// These functions (C) Chris Wojtan
-// Long-term goal is to unify with his split&merge codebase
-//
-/******************************************************************************/
-
-#include <queue>
-#include <algorithm>
-#include "mesh.h"
-#include "kernel.h"
-#include "edgecollapse.h"
-#include <mesh.h>
-#include <stack>
-
-using namespace std;
-
-namespace Manta { 
-
-//! Mesh smoothing 
-/*! see Desbrun 99 "Implicit fairing of of irregular meshes using diffusion and curvature flow"*/
-PYTHON void smoothMesh(Mesh& mesh, Real strength, int steps = 1, Real minLength=1e-5) {
-	const Real dt = mesh.getParent()->getDt();
-	const Real str = min(dt * strength, (Real)1);
-	mesh.rebuildQuickCheck(); 
-	
-	// calculate original mesh volume
-	Vec3 origCM;
-	Real origVolume = mesh.computeCenterOfMass(origCM);
-	
-	// temp vertices
-	const int numCorners = mesh.numTris() * 3;
-	const int numNodes= mesh.numNodes();            
-	vector<Vec3> temp(numNodes);
-	vector<bool> visited(numNodes);
-	
-	for (int s = 0; s<steps; s++) {
-		// reset markers
-		for(size_t i=0; i<visited.size(); i++) visited[i] = false;
-		
-		for (int c = 0; c < numCorners; c++) {
-			const int node = mesh.corners(c).node;
-			if (visited[node]) continue;
-			
-			const Vec3 pos = mesh.nodes(node).pos;
-			Vec3 dx(_0);
-			Real totalLen = 0;
-				
-			// rotate around vertex
-			set<int>& ring = mesh.get1Ring(node).nodes;
-			for(set<int>::iterator it=ring.begin(); it!=ring.end(); it++) {
-				Vec3 edge = mesh.nodes(*it).pos - pos;
-				Real len = norm(edge);
-
-				if (len > minLength) {
-					dx += edge * (_1/len);
-					totalLen += len;
-				} else {
-					totalLen = _0;
-					break;
-				}                
-			}
-			visited[node] = true;
-			temp[node] = pos;            
-			if (totalLen != 0)
-				temp[node] += dx * (str / totalLen);
-		}
-		
-		// copy back
-		for (int n=0; n<numNodes; n++)
-			if (!mesh.isNodeFixed(n))
-				mesh.nodes(n).pos = temp[n];
-	}
-	
-	// calculate new mesh volume
-	Vec3 newCM;
-	Real newVolume = mesh.computeCenterOfMass(newCM);
-	
-	// preserve volume : scale relative to CM
-	Real beta;
-#if defined(WIN32) || defined(_WIN32)
-	beta = pow( (Real)abs(origVolume/newVolume), (Real)(1./3.) );
-#else
-	beta = cbrt( origVolume/newVolume );
-#	endif
-
-	for (int n=0; n<numNodes; n++)
-		if (!mesh.isNodeFixed(n))
-			mesh.nodes(n).pos = origCM + (mesh.nodes(n).pos - newCM) * beta;
-}
-
-//! Subdivide and edgecollapse to guarantee mesh with edgelengths between
-//! min/maxLength and an angle below minAngle
-PYTHON void subdivideMesh(Mesh& mesh, Real minAngle, Real minLength, Real maxLength, bool cutTubes = false) {
-	// gather some statistics
-	int edgeSubdivs = 0, edgeCollsAngle = 0, edgeCollsLen = 0, edgeKill = 0;
-	mesh.rebuildQuickCheck(); 
-
-	vector<int> deletedNodes;
-	map<int,bool> taintedTris;
-	priority_queue<pair<Real,int> > pq;
-	
-	//////////////////////////////////////////
-	// EDGE COLLAPSE                        //
-	//    - particles marked for deletation //
-	//////////////////////////////////////////
-	
-	for (int t=0; t<mesh.numTris(); t++) {
-		if(taintedTris.find(t)!=taintedTris.end())
-			continue;
-		
-		// check if at least 2 nodes are marked for delete
-		bool k[3];
-		int numKill = 0;
-		for (int i=0; i<3; i++) {
-			k[i] = mesh.nodes(mesh.tris(t).c[i]).flags & Mesh::NfKillme;
-			if (k[i]) numKill++;
-		}        
-		if (numKill<2) continue;
-		
-		if (k[0] && k[1])
-			CollapseEdge(mesh, t, 2, mesh.getEdge(t,0), mesh.getNode(t,0), deletedNodes, taintedTris, edgeKill, cutTubes);
-		else if (k[1] && k[2])
-			CollapseEdge(mesh, t, 0, mesh.getEdge(t,1), mesh.getNode(t,1), deletedNodes, taintedTris, edgeKill, cutTubes);
-		else if (k[2] && k[0])
-			CollapseEdge(mesh, t, 1, mesh.getEdge(t,2), mesh.getNode(t,2), deletedNodes, taintedTris, edgeKill, cutTubes);
-	}
-			
-	//////////////////////////////////////////
-	// EDGE COLLAPSING                      //
-	//      - based on small triangle angle //
-	//////////////////////////////////////////
-	
-	if (minAngle > 0) {
-		for(int t=0; t<mesh.numTris(); t++) {
-			// we only want to run through the edge list ONCE.
-			// we achieve this in a method very similar to the above subdivision method.
-
-			// if this triangle has already been deleted, ignore it
-			if(taintedTris.find(t)!=taintedTris.end())
-				continue;
-			
-			// first we find the angles of this triangle
-			Vec3 e0 = mesh.getEdge(t,0), e1 = mesh.getEdge(t,1), e2 = mesh.getEdge(t,2);
-			Vec3 ne0 = e0;
-			Vec3 ne1 = e1;
-			Vec3 ne2 = e2;
-			normalize(ne0);
-			normalize(ne1);
-			normalize(ne2);
-
-			//Real thisArea = sqrMag(cross(-e2,e0));
-			// small angle approximation says sin(x) = arcsin(x) = x,
-			// arccos(x) = pi/2 - arcsin(x),
-			// cos(x) = dot(A,B),
-			// so angle is approximately 1 - dot(A,B).
-			Real angle[3];
-			angle[0] = 1.0-dot(ne0,-ne2);
-			angle[1] = 1.0-dot(ne1,-ne0);
-			angle[2] = 1.0-dot(ne2,-ne1);
-			Real worstAngle = angle[0];
-			int which = 0;
-			if(angle[1]<worstAngle) {
-				worstAngle = angle[1];
-				which = 1;
-			}
-			if(angle[2]<worstAngle) {
-				worstAngle = angle[2];
-				which = 2;
-			}
-
-			// then we see if the angle is too small
-			if(worstAngle<minAngle) {
-				Vec3 edgevect;
-				Vec3 endpoint;
-				switch(which) {
-				case 0:
-					endpoint = mesh.getNode(t,1);
-					edgevect = e1;
-					break;
-				case 1:
-					endpoint = mesh.getNode(t,2);
-					edgevect = e2;
-					break;
-				case 2:
-					endpoint = mesh.getNode(t,0);
-					edgevect = e0;
-					break;
-				default:
-					break;
-				}
-
-				CollapseEdge(mesh, t,which,edgevect,endpoint,deletedNodes,taintedTris, edgeCollsAngle, cutTubes);
-			}
-		}
-	}
-	
-	//////////////////////
-	// EDGE SUBDIVISION //
-	//////////////////////
-
-	Real maxLength2 = maxLength*maxLength;
-	for (int t=0; t<mesh.numTris(); t++) {
-		// first we find the maximum length edge in this triangle
-		Vec3 e0 = mesh.getEdge(t,0), e1 = mesh.getEdge(t,1), e2 = mesh.getEdge(t,2);
-		Real d0 = normSquare(e0);
-		Real d1 = normSquare(e1);
-		Real d2 = normSquare(e2);
-
-		Real longest = max(d0,max(d1,d2));        
-		if(longest > maxLength2) {
-			pq.push(pair<Real,int>(longest,t));
-		}
-	}
-	if (maxLength > 0) {        
-		
-		while(!pq.empty() && pq.top().first>maxLength2) {
-			// we only want to run through the edge list ONCE
-			// and we want to subdivide the original edges before we subdivide any newer, shorter edges,
-			// so whenever we subdivide, we add the 2 new triangles on the end of the SurfaceTri vector
-			// and mark the original subdivided triangles for deletion.
-			//  when we are done subdividing, we delete the obsolete triangles
-		
-			int triA = pq.top().second;
-			pq.pop();
-
-			if(taintedTris.find(triA)!=taintedTris.end())
-				continue;
-
-			// first we find the maximum length edge in this triangle
-			Vec3 e0 = mesh.getEdge(triA,0), e1 = mesh.getEdge(triA,1), e2 = mesh.getEdge(triA,2);
-			Real d0 = normSquare(e0);
-			Real d1 = normSquare(e1);
-			Real d2 = normSquare(e2);
-
-			Vec3 edgevect;
-			Vec3 endpoint;
-			int which;
-			if(d0>d1) {
-				if(d0>d2) {
-					edgevect = e0;
-					endpoint = mesh.getNode(triA, 0);;
-					which = 2;  // 2 opposite of edge 0-1
-				} else {
-					edgevect = e2;
-					endpoint = mesh.getNode(triA, 2);
-					which = 1;  // 1 opposite of edge 2-0
-				}
-			} else {
-				if(d1>d2) {
-					edgevect = e1;
-					endpoint = mesh.getNode(triA, 1);
-					which = 0;  // 0 opposite of edge 1-2
-				} else {
-					edgevect = e2;
-					endpoint = mesh.getNode(triA, 2);
-					which = 1;  // 1 opposite of edge 2-0
-				}
-			}
-			// This edge is too long, so we split it in the middle
-
-			//         *
-			//        / \.
-			//       /C0 \.
-			//      /     \.
-			//     /       \.
-			//    /    B    \.
-			//   /           \.
-			//  /C1        C2 \.
-			// *---------------*
-			//  \C2        C1 /
-			//   \           /
-			//    \    A    /
-			//     \       /
-			//      \     /
-			//       \C0 /
-			//        \ /
-			//         *
-			//
-			//      BECOMES
-			//
-			//         *
-			//        /|\.
-			//       / | \.
-			//      /C0|C0\.
-			//     /   |   \.
-			//    / B1 | B2 \.
-			//   /     |     \.
-			//  /C1  C2|C1 C2 \.
-			// *-------*-------*
-			//  \C2  C1|C2  C1/
-			//   \     |     /
-			//    \ A2 | A1 /
-			//     \   |   /
-			//      \C0|C0/
-			//       \ | /
-			//        \|/
-			//         *
-	
-			int triB = -1; bool haveB = false;
-			Corner ca_old[3],cb_old[3];
-			ca_old[0] = mesh.corners(triA, which);
-			ca_old[1] = mesh.corners(ca_old[0].next);
-			ca_old[2] = mesh.corners(ca_old[0].prev);
-			if (ca_old[0].opposite>=0) {
-				cb_old[0] = mesh.corners(ca_old[0].opposite);
-				cb_old[1] = mesh.corners(cb_old[0].next);
-				cb_old[2] = mesh.corners(cb_old[0].prev);
-				triB = cb_old[0].tri;
-				haveB = true;
-			} 
-			//else throw Error("nonmanifold");
-			
-			// subdivide in the middle of the edge and create new triangles
-			Node newNode;
-			newNode.flags = 0;
-			
-			newNode.pos = endpoint + 0.5*edgevect; // fallback: linear average            
-			// default: use butterfly
-			if (haveB)
-				newNode.pos = ModifiedButterflySubdivision(mesh, ca_old[0], cb_old[0], newNode.pos);
-
-			// find indices of two points of 'which'-edge
-			// merge flags
-			int P0 = ca_old[1].node;
-			int P1 = ca_old[2].node;
-			newNode.flags = mesh.nodes(P0).flags | mesh.nodes(P1).flags;
-			
-			Real len0 = norm(mesh.nodes(P0).pos - newNode.pos);
-			Real len1 = norm(mesh.nodes(P1).pos - newNode.pos);
-			
-			// remove P0/P1 1-ring connection
-			mesh.get1Ring(P0).nodes.erase(P1);
-			mesh.get1Ring(P1).nodes.erase(P0);
-			mesh.get1Ring(P0).tris.erase(triA);
-			mesh.get1Ring(P1).tris.erase(triA);
-			mesh.get1Ring(ca_old[0].node).tris.erase(triA);
-			if (haveB) {
-				mesh.get1Ring(P0).tris.erase(triB);
-				mesh.get1Ring(P1).tris.erase(triB);
-				mesh.get1Ring(cb_old[0].node).tris.erase(triB);
-			}
-			
-			// init channel properties for new node
-			for(int i=0; i<mesh.numNodeChannels(); i++) {
-				mesh.nodeChannel(i)->addInterpol(P0, P1, len0/(len0+len1));
-			}
-			
-			// write to array
-			mesh.addTri(Triangle(ca_old[0].node, ca_old[1].node, mesh.numNodes()));
-			mesh.addTri(Triangle(ca_old[0].node, mesh.numNodes(), ca_old[2].node));
-			if (haveB) {
-				mesh.addTri(Triangle(cb_old[0].node, cb_old[1].node, mesh.numNodes()));
-				mesh.addTri(Triangle(cb_old[0].node, mesh.numNodes(), cb_old[2].node));
-			}
-			mesh.addNode(newNode);
-			
-			const int nt = haveB ? 4 : 2;
-			int triA1 = mesh.numTris()-nt;
-			int triA2 = mesh.numTris()-nt+1;            
-			int triB1=0, triB2=0;
-			if (haveB) {
-				triB1 = mesh.numTris()-nt+2;
-				triB2 = mesh.numTris()-nt+3;
-			}
-			mesh.tris(triA1).flags = mesh.tris(triA).flags;
-			mesh.tris(triA2).flags = mesh.tris(triA).flags;
-			mesh.tris(triB1).flags = mesh.tris(triB).flags;
-			mesh.tris(triB2).flags = mesh.tris(triB).flags;
-			
-			// connect new triangles to outside triangles,
-			// and connect outside triangles to these new ones
-			for (int c=0; c<3; c++) mesh.addCorner(Corner(triA1,mesh.tris(triA1).c[c]));
-			for (int c=0; c<3; c++) mesh.addCorner(Corner(triA2,mesh.tris(triA2).c[c]));
-			if (haveB) {
-				for (int c=0; c<3; c++) mesh.addCorner(Corner(triB1,mesh.tris(triB1).c[c]));
-				for (int c=0; c<3; c++) mesh.addCorner(Corner(triB2,mesh.tris(triB2).c[c]));            
-			}
-			
-			int baseIdx = 3*(mesh.numTris()-nt);
-			Corner* cBase = &mesh.corners(baseIdx);
-			
-			// set next/prev
-			for (int t=0; t<nt; t++)
-				for (int c=0; c<3; c++) {
-					cBase[t*3+c].next = baseIdx+t*3+((c+1)%3);
-					cBase[t*3+c].prev = baseIdx+t*3+((c+2)%3);                        
-				}
-				
-			// set opposites
-			// A1
-			cBase[0].opposite = haveB ? (baseIdx+9) : -1;
-			cBase[1].opposite = baseIdx+5;
-			cBase[2].opposite = -1;
-			if (ca_old[2].opposite>=0) {
-				cBase[2].opposite = ca_old[2].opposite;
-				mesh.corners(cBase[2].opposite).opposite = baseIdx+2;
-			}
-			// A2
-			cBase[3].opposite = haveB ? (baseIdx+6) : -1;
-			cBase[4].opposite = -1;
-			if (ca_old[1].opposite>=0) {
-				cBase[4].opposite = ca_old[1].opposite;
-				mesh.corners(cBase[4].opposite).opposite = baseIdx+4;
-			}
-			cBase[5].opposite = baseIdx+1;
-			if (haveB) {
-				// B1
-				cBase[6].opposite = baseIdx+3;
-				cBase[7].opposite = baseIdx+11;
-				cBase[8].opposite = -1;
-				if (cb_old[2].opposite>=0) {
-					cBase[8].opposite = cb_old[2].opposite;
-					mesh.corners(cBase[8].opposite).opposite = baseIdx+8;
-				}
-				// B2
-				cBase[9].opposite = baseIdx+0;
-				cBase[10].opposite = -1;
-				if (cb_old[1].opposite>=0) {
-					cBase[10].opposite = cb_old[1].opposite;
-					mesh.corners(cBase[10].opposite).opposite = baseIdx+10;
-				}
-				cBase[11].opposite = baseIdx+7;
-			}
-			
-			////////////////////
-			// mark the two original triangles for deletion
-			taintedTris[triA] = true;
-			mesh.removeTriFromLookup(triA);
-			if (haveB) {
-				taintedTris[triB] = true;
-				mesh.removeTriFromLookup(triB);
-			}
-			
-			Real areaA1 = mesh.getFaceArea(triA1), areaA2 = mesh.getFaceArea(triA2);
-			Real areaB1=0, areaB2=0;
-			if (haveB) {
-				areaB1 = mesh.getFaceArea(triB1);
-				areaB2 = mesh.getFaceArea(triB2);                
-			}
-			
-			// add channel props for new triangles
-			for(int i=0; i<mesh.numTriChannels(); i++) {                
-				mesh.triChannel(i)->addSplit(triA, areaA1/(areaA1+areaA2));
-				mesh.triChannel(i)->addSplit(triA, areaA2/(areaA1+areaA2));
-				if (haveB) {
-					mesh.triChannel(i)->addSplit(triB, areaB1/(areaB1+areaB2));
-					mesh.triChannel(i)->addSplit(triB, areaB2/(areaB1+areaB2));
-				}
-			}
-			
-			// add the four new triangles to the prority queue
-			for(int i=mesh.numTris()-nt; i<mesh.numTris(); i++) {
-				// find the maximum length edge in this triangle
-				Vec3 ne0 = mesh.getEdge(i, 0), ne1 = mesh.getEdge(i, 1), ne2 = mesh.getEdge(i, 2);                
-				Real nd0 = normSquare(ne0);
-				Real nd1 = normSquare(ne1);
-				Real nd2 = normSquare(ne2);
-				Real longest = max(nd0,max(nd1,nd2));
-				//longest = (int)(longest * 1e2) / 1e2; // HACK: truncate
-				pq.push(pair<Real,int>(longest,i));                
-			}
-			edgeSubdivs++;            
-		}
-	}
-	
-	//////////////////////////////////////////
-	// EDGE COLLAPSING                      //
-	//      - based on short edge length    //
-	//////////////////////////////////////////
-	if (minLength > 0) {
-		const Real minLength2 = minLength*minLength;
-		for(int t=0; t<mesh.numTris(); t++) {
-			// we only want to run through the edge list ONCE.
-			// we achieve this in a method very similar to the above subdivision method.
-
-			// NOTE:
-			// priority queue does not work so great in the edge collapse case,
-			// because collapsing one triangle affects the edge lengths 
-			// of many neighbor triangles,
-			// and we do not update their maximum edge length in the queue.
-
-			// if this triangle has already been deleted, ignore it
-			//if(taintedTris[t])
-			//  continue;
-
-			if(taintedTris.find(t)!=taintedTris.end())
-				continue;
-			
-			// first we find the minimum length edge in this triangle
-			Vec3 e0 = mesh.getEdge(t,0), e1 = mesh.getEdge(t,1), e2 = mesh.getEdge(t,2);
-			Real d0 = normSquare(e0);
-			Real d1 = normSquare(e1);
-			Real d2 = normSquare(e2);
-
-			Vec3 edgevect;
-			Vec3 endpoint;
-			Real dist2;
-			int which;
-			if(d0<d1) {
-				if(d0<d2) {
-					dist2 = d0;
-					edgevect = e0;
-					endpoint = mesh.getNode(t,0);
-					which = 2;  // 2 opposite of edge 0-1
-				} else {
-					dist2 = d2;
-					edgevect = e2;
-					endpoint =  mesh.getNode(t,2);
-					which = 1;  // 1 opposite of edge 2-0
-				}
-			} else {
-				if(d1<d2) {
-					dist2 = d1;
-					edgevect = e1;
-					endpoint =  mesh.getNode(t,1);
-					which = 0;  // 0 opposite of edge 1-2
-				} else {
-					dist2 = d2;
-					edgevect = e2;
-					endpoint =  mesh.getNode(t,2);
-					which = 1;  // 1 opposite of edge 2-0
-				}
-			}
-			// then we see if the min length edge is too short
-			if(dist2<minLength2) {
-				CollapseEdge(mesh, t,which,edgevect,endpoint, deletedNodes,taintedTris, edgeCollsLen, cutTubes);
-			}
-		}
-	}
-	// cleanup nodes and triangles marked for deletion
-	
-	//  we run backwards through the deleted array, 
-	//  replacing triangles with ones from the back
-	//          (this avoids the potential problem of overwriting a triangle
-	//              with a to-be-deleted triangle)    
-	std::map<int,bool>::reverse_iterator tti = taintedTris.rbegin();
-	for(;tti!=taintedTris.rend(); tti++)
-		mesh.removeTri(tti->first);
-	
-	mesh.removeNodes(deletedNodes);
-	cout << "Surface subdivision finished with " << mesh.numNodes() << " surface nodes and " << mesh.numTris();
-	cout << " surface triangles, edgeSubdivs:" << edgeSubdivs << ", edgeCollapses: " << edgeCollsLen;
-	cout << " + " << edgeCollsAngle << " + " << edgeKill << endl;
-	//mesh.sanityCheck();
-	
-}
-	
-PYTHON void killSmallComponents(Mesh& mesh, int elements = 10) {
-	const int num = mesh.numTris();
-	vector<int> comp(num);
-	vector<int> numEl;
-	vector<int> deletedNodes;
-	vector<bool> isNodeDel(mesh.numNodes());
-	map<int,bool> taintedTris;
-	// enumerate components
-	int cur=0;
-	for (int i=0; i<num; i++) {
-		if (comp[i]==0) {
-			cur++;
-			comp[i] = cur;
-			
-			stack<int> stack;
-			stack.push(i);
-			int cnt = 1;
-			while(!stack.empty()) {
-				int tri = stack.top();
-				stack.pop();
-				for (int c=0; c<3; c++) {
-					int op = mesh.corners(tri,c).opposite;
-					if (op < 0) continue;
-					int ntri = mesh.corners(op).tri;
-					if (comp[ntri]==0) {
-						comp[ntri] = cur;
-						stack.push(ntri);
-						cnt++;
-					}
-				}
-			}
-			numEl.push_back(cnt);
-		}
-	}
-	// kill small components
-	for (int j=0; j<num; j++) {
-		if (numEl[comp[j]-1] < elements) {
-			taintedTris[j] = true;
-			for (int c=0; c<3; c++) {
-				int n=mesh.tris(j).c[c];
-				if (!isNodeDel[n]) {
-					isNodeDel[n] = true;
-					deletedNodes.push_back(n);
-				}
-			}
-		}
-	}
-	
-	std::map<int,bool>::reverse_iterator tti = taintedTris.rbegin();
-	for(;tti!=taintedTris.rend(); tti++)
-		mesh.removeTri(tti->first);
-	
-	mesh.removeNodes(deletedNodes);
-	
-	if (!taintedTris.empty())
-		cout << "Killed small components : " << deletedNodes.size() << " nodes, " << taintedTris.size() << " tris deleted." << endl;
-}
-   
-	
-} //namespace
-
diff --git a/source/blender/python/manta_pp/source/pressure.cpp b/source/blender/python/manta_pp/source/pressure.cpp
deleted file mode 100644
index 3c826b1d519..00000000000
--- a/source/blender/python/manta_pp/source/pressure.cpp
+++ /dev/null
@@ -1,310 +0,0 @@
-/******************************************************************************

- *

- * MantaFlow fluid solver framework

- * Copyright 2011 Tobias Pfaff, Nils Thuerey 

- *

- * This program is free software, distributed under the terms of the

- * GNU General Public License (GPL) 

- * http://www.gnu.org/licenses

- *

- * Plugins for pressure correction: solve_pressure, and ghost fluid helpers

- *

- ******************************************************************************/

-#include "vectorbase.h"

-#include "kernel.h"

-#include "conjugategrad.h"

-

-using namespace std;

-namespace Manta {

-

-//! Kernel: Construct the right-hand side of the poisson equation

-KERNEL(bnd=1, reduce=+) returns(int cnt=0) returns(double sum=0)

-void MakeRhs (FlagGrid& flags, Grid<Real>& rhs, MACGrid& vel, 

-			  Grid<Real>* perCellCorr) 

-{

-	if (!flags.isFluid(i,j,k)) {

-		rhs(i,j,k) = 0;

-		return;

-	}

-	   

-	// compute divergence 

-	// no flag checks: assumes vel at obstacle interfaces is set to zero

-	Real set =          vel(i,j,k).x - vel(i+1,j,k).x + 

-						vel(i,j,k).y - vel(i,j+1,k).y; 

-	if(vel.is3D()) set+=vel(i,j,k).z - vel(i,j,k+1).z;

-	

-	// per cell divergence correction

-	if(perCellCorr) 

-		set += perCellCorr->get(i,j,k);

-	

-	// obtain sum, cell count

-	sum += set;

-	cnt++;

-	

-	rhs(i,j,k) = set;

-}

-

-

-//! Kernel: Apply velocity update from poisson equation

-KERNEL(bnd=1) 

-void CorrectVelocity(FlagGrid& flags, MACGrid& vel, Grid<Real>& pressure) 

-{

-	int idx = flags.index(i,j,k);

-	if (flags.isFluid(idx))

-	{

-		if (flags.isFluid(i-1,j,k)) vel[idx].x -= (pressure[idx] - pressure(i-1,j,k));

-		if (flags.isFluid(i,j-1,k)) vel[idx].y -= (pressure[idx] - pressure(i,j-1,k));

-		if (flags.is3D() && flags.isFluid(i,j,k-1)) vel[idx].z -= (pressure[idx] - pressure(i,j,k-1));

- 

-		if (flags.isEmpty(i-1,j,k)) vel[idx].x -= pressure[idx];

-		if (flags.isEmpty(i,j-1,k)) vel[idx].y -= pressure[idx];

-		if (flags.is3D() && flags.isEmpty(i,j,k-1)) vel[idx].z -= pressure[idx];

-	}

-	else if (flags.isEmpty(idx))

-	{

-		if (flags.isFluid(i-1,j,k)) vel[idx].x += pressure(i-1,j,k);

-		else                        vel[idx].x  = 0.f;

-		if (flags.isFluid(i,j-1,k)) vel[idx].y += pressure(i,j-1,k);

-		else                        vel[idx].y  = 0.f;

-		if (flags.is3D() ) {

-		if (flags.isFluid(i,j,k-1)) vel[idx].z += pressure(i,j,k-1);

-		else                        vel[idx].z  = 0.f;

-		}

-	}

-}

-

-//! Kernel: Set matrix stencils and velocities to enable open boundaries

-KERNEL void SetOpenBound(Grid<Real>& A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak, MACGrid& vel,

-						 Vector3D<bool> lowerBound, Vector3D<bool> upperBound)

-{    

-	// set velocity boundary conditions

-	if (lowerBound.x && i == 0) vel(0,j,k) = vel(1,j,k);

-	if (lowerBound.y && j == 0) vel(i,0,k) = vel(i,1,k);

-	if (upperBound.x && i == maxX-1) vel(maxX-1,j,k) = vel(maxX-2,j,k);

-	if (upperBound.y && j == maxY-1) vel(i,maxY-1,k) = vel(i,maxY-2,k);

-	if(vel.is3D()) {

-		if (lowerBound.z && k == 0)      vel(i,j,0)      = vel(i,j,1);

-		if (upperBound.z && k == maxZ-1) vel(i,j,maxZ-1) = vel(i,j,maxZ-2); 

-	}

-	

-	// set matrix stencils at boundary

-	if ((lowerBound.x && i<=1) || (upperBound.x && i>=maxX-2) ||

-		(lowerBound.y && j<=1) || (upperBound.y && j>=maxY-2) ||

-		(lowerBound.z && k<=1) || (upperBound.z && k>=maxZ-2)) {

-		A0(i,j,k) = vel.is3D() ? 6.0 : 4.0;

-		Ai(i,j,k) = -1.0;

-		Aj(i,j,k) = -1.0;

-		if (vel.is3D()) Ak(i,j,k) = -1.0;

-	}

-}

-

-//! Kernel: Set matrix rhs for outflow

-KERNEL void SetOutflow (Grid<Real>& rhs, Vector3D<bool> lowerBound, Vector3D<bool> upperBound, int height)

-{

-	if ((lowerBound.x && i < height) || (upperBound.x && i >= maxX-1-height) ||

-		(lowerBound.y && j < height) || (upperBound.y && j >= maxY-1-height) ||

-		(lowerBound.z && k < height) || (upperBound.z && k >= maxZ-1-height))

-		rhs(i,j,k) = 0;

-}

-

-

-// *****************************************************************************

-// Ghost fluid helpers

-

-// calculate fraction filled with liquid (note, assumes inside value is < outside!)

-inline static Real thetaHelper(Real inside, Real outside)

-{

-	Real denom = inside-outside;

-	if (denom > -1e-04) return 0.5; // should always be neg, and large enough...

-	return std::max(Real(0), std::min(Real(1), inside/denom));

-}

-

-// calculate ghost fluid factor, cell at idx should be a fluid cell

-inline static Real ghostFluidHelper(int idx, int offset, const Grid<Real> &phi, Real gfClamp)

-{

-	Real alpha = thetaHelper(phi[idx], phi[idx+offset]);

-	if (alpha < gfClamp) return alpha = gfClamp;

-	return (1-(1/alpha)); 

-}

-

-//! Kernel: Adapt A0 for ghost fluid

-KERNEL(bnd=1) 

-void ApplyGhostFluidDiagonal(Grid<Real> &A0, const FlagGrid &flags, const Grid<Real> &phi, Real gfClamp)

-{

-	const int X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ();

-	int idx = flags.index(i,j,k);

-	if (!flags.isFluid(idx)) return;

-

-	if (flags.isEmpty(i-1,j,k)) A0[idx] -= ghostFluidHelper(idx, -X, phi, gfClamp);

-	if (flags.isEmpty(i+1,j,k)) A0[idx] -= ghostFluidHelper(idx, +X, phi, gfClamp);

-	if (flags.isEmpty(i,j-1,k)) A0[idx] -= ghostFluidHelper(idx, -Y, phi, gfClamp);

-	if (flags.isEmpty(i,j+1,k)) A0[idx] -= ghostFluidHelper(idx, +Y, phi, gfClamp);

-	if (flags.is3D()) {

-		if (flags.isEmpty(i,j,k-1)) A0[idx] -= ghostFluidHelper(idx, -Z, phi, gfClamp);

-		if (flags.isEmpty(i,j,k+1)) A0[idx] -= ghostFluidHelper(idx, +Z, phi, gfClamp);

-	}

-}

-

-//! Kernel: Apply velocity update: ghost fluid contribution

-KERNEL(bnd=1)

-void CorrectVelocityGhostFluid(MACGrid &vel, const FlagGrid &flags, const Grid<Real> &pressure, const Grid<Real> &phi, Real gfClamp)

-{

-	const int X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ();

-	const int idx = flags.index(i,j,k);

-	if (flags.isFluid(idx))

-	{

-		if (flags.isEmpty(i-1,j,k)) vel[idx][0] += pressure[idx] * ghostFluidHelper(idx, -X, phi, gfClamp);

-		if (flags.isEmpty(i,j-1,k)) vel[idx][1] += pressure[idx] * ghostFluidHelper(idx, -Y, phi, gfClamp);

-		if (flags.is3D() && flags.isEmpty(i,j,k-1)) vel[idx][2] += pressure[idx] * ghostFluidHelper(idx, -Z, phi, gfClamp);

-	}

-	else if (flags.isEmpty(idx))

-	{

-		if (flags.isFluid(i-1,j,k)) vel[idx][0] -= pressure(i-1,j,k) * ghostFluidHelper(idx-X, +X, phi, gfClamp);

-		else                        vel[idx].x  = 0.f;

-		if (flags.isFluid(i,j-1,k)) vel[idx][1] -= pressure(i,j-1,k) * ghostFluidHelper(idx-Y, +Y, phi, gfClamp);

-		else                        vel[idx].y  = 0.f;

-		if (flags.is3D() ) {

-		if (flags.isFluid(i,j,k-1)) vel[idx][2] -= pressure(i,j,k-1) * ghostFluidHelper(idx-Z, +Z, phi, gfClamp);

-		else                        vel[idx].z  = 0.f;

-		}

-	}

-}

-

-

-// improve behavior of clamping for large time steps:

-

-inline static Real ghostFluidWasClamped(int idx, int offset, const Grid<Real> &phi, Real gfClamp)

-{

-	Real alpha = thetaHelper(phi[idx], phi[idx+offset]);

-	if (alpha < gfClamp) return true;

-	return false;

-}

-

-KERNEL(bnd=1)

-void ReplaceClampedGhostFluidVels(MACGrid &vel, FlagGrid &flags, 

-		const Grid<Real> &pressure, const Grid<Real> &phi, Real gfClamp )

-{

-	const int X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ();

-	const int idx = flags.index(i,j,k);

-	if (flags.isFluid(idx))

-	{

-		if( (flags.isEmpty(i-1,j,k)) && (ghostFluidWasClamped(idx, -X, phi, gfClamp)) )

-			vel[idx-X][0] = vel[idx][0];

-		if( (flags.isEmpty(i,j-1,k)) && (ghostFluidWasClamped(idx, -Y, phi, gfClamp)) )

-			vel[idx-Y][1] = vel[idx][1];

-		if( flags.is3D() && 

-		   (flags.isEmpty(i,j,k-1)) && (ghostFluidWasClamped(idx, -Z, phi, gfClamp)) )

-			vel[idx-Z][2] = vel[idx][2];

-	}

-	else if (flags.isEmpty(idx))

-	{

-		if( (i>-1) && (flags.isFluid(i-1,j,k)) && ( ghostFluidWasClamped(idx-X, +X, phi, gfClamp) ) )

-			vel[idx][0] = vel[idx-X][0];

-		if( (j>-1) && (flags.isFluid(i,j-1,k)) && ( ghostFluidWasClamped(idx-Y, +Y, phi, gfClamp) ) )

-			vel[idx][1] = vel[idx-Y][1];

-		if( flags.is3D() &&

-		  ( (k>-1) && (flags.isFluid(i,j,k-1)) && ( ghostFluidWasClamped(idx-Z, +Z, phi, gfClamp) ) ))

-			vel[idx][2] = vel[idx-Z][2];

-	}

-}

-

-

-// *****************************************************************************

-// Main pressure solve

-

-inline void convertDescToVec(const string& desc, Vector3D<bool>& lo, Vector3D<bool>& up) {

-	for(size_t i=0; i<desc.size(); i++) {

-		if (desc[i] == 'x') lo.x = true;

-		else if (desc[i] == 'y') lo.y = true;

-		else if (desc[i] == 'z') lo.z = true;

-		else if (desc[i] == 'X') up.x = true;

-		else if (desc[i] == 'Y') up.y = true;

-		else if (desc[i] == 'Z') up.z = true;

-		else errMsg("invalid character in boundary description string. Only [xyzXYZ] allowed.");

-	}

-}

-

-//! Perform pressure projection of the velocity grid

-PYTHON void solvePressure(MACGrid& vel, Grid<Real>& pressure, FlagGrid& flags,

-					 Grid<Real>* phi = 0, 

-					 Grid<Real>* perCellCorr = 0, 

-					 Real gfClamp = 1e-04, 

-					 Real cgMaxIterFac = 1.5,

-					 Real cgAccuracy = 1e-3,

-					 string openBound = "",

-					 string outflow = "",

-					 int outflowHeight = 1,

-					 bool precondition = true,

-					 bool enforceCompatibility = false,

-					 bool useResNorm = true )

-{

-	// parse strings

-	Vector3D<bool> loOpenBound, upOpenBound, loOutflow, upOutflow;

-	convertDescToVec(openBound, loOpenBound, upOpenBound);

-	convertDescToVec(outflow, loOutflow, upOutflow);

-	if (vel.is2D() && (loOpenBound.z || upOpenBound.z))

-		errMsg("open boundaries for z specified for 2D grid");

-	

-	// reserve temp grids

-	FluidSolver* parent = flags.getParent();

-	Grid<Real> rhs(parent);

-	Grid<Real> residual(parent);

-	Grid<Real> search(parent);

-	Grid<Real> A0(parent);

-	Grid<Real> Ai(parent);

-	Grid<Real> Aj(parent);

-	Grid<Real> Ak(parent);

-	Grid<Real> tmp(parent);

-	Grid<Real> pca0(parent);

-	Grid<Real> pca1(parent);

-	Grid<Real> pca2(parent);

-	Grid<Real> pca3(parent);

-		

-	// setup matrix and boundaries

-	MakeLaplaceMatrix (flags, A0, Ai, Aj, Ak);

-	SetOpenBound (A0, Ai, Aj, Ak, vel, loOpenBound, upOpenBound);

-	

-	if (phi) {

-		ApplyGhostFluidDiagonal(A0, flags, *phi, gfClamp);

-	}

-	

-	// compute divergence and init right hand side

-	MakeRhs kernMakeRhs (flags, rhs, vel, perCellCorr);

-	

-	if (!outflow.empty())

-		SetOutflow (rhs, loOutflow, upOutflow, outflowHeight);

-	

-	if (enforceCompatibility)

-		rhs += (Real)(-kernMakeRhs.sum / (Real)kernMakeRhs.cnt);

-	

-	// CG setup

-	// note: the last factor increases the max iterations for 2d, which right now can't use a preconditioner 

-	const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);

-	GridCgInterface *gcg;

-	if (vel.is3D())

-		gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );

-	else

-		gcg = new GridCg<ApplyMatrix2D>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );

-	

-	gcg->setAccuracy( cgAccuracy ); 

-	gcg->setUseResNorm( useResNorm );

-

-	// optional preconditioning

-	gcg->setPreconditioner( precondition ? GridCgInterface::PC_mICP : GridCgInterface::PC_None, &pca0, &pca1, &pca2, &pca3);

-

-	for (int iter=0; iter<maxIter; iter++) {

-		if (!gcg->iterate()) iter=maxIter;

-	} 

-	debMsg("FluidSolver::solvePressure iterations:"<<gcg->getIterations()<<", res:"<<gcg->getSigma(), 1);

-	delete gcg;

-	

-	CorrectVelocity(flags, vel, pressure);

-	if (phi) {

-		CorrectVelocityGhostFluid (vel, flags, pressure, *phi, gfClamp);

-		// improve behavior of clamping for large time steps:

-		ReplaceClampedGhostFluidVels (vel, flags, pressure, *phi, gfClamp);

-	}

-}

-

-} // end namespace

-

diff --git a/source/blender/python/manta_pp/source/vortexplugins.cpp b/source/blender/python/manta_pp/source/vortexplugins.cpp
deleted file mode 100644
index fed847ca2bd..00000000000
--- a/source/blender/python/manta_pp/source/vortexplugins.cpp
+++ /dev/null
@@ -1,312 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Plugins for using vortex sheet meshes 
- *
- ******************************************************************************/
- 
-#include <iostream>
-#include "vortexsheet.h"
-#include "vortexpart.h"
-#include "shapes.h"
-#include "commonkernels.h"
-#include "conjugategrad.h"
-#include "randomstream.h"
-#include "levelset.h"
-
-using namespace std;
-
-namespace Manta {
-	
-//! Mark area of mesh inside shape as fixed nodes. 
-//! Remove all other fixed nodes if 'exclusive' is set
-PYTHON void markAsFixed(Mesh& mesh, Shape* shape, bool exclusive=true) 
-{
-	for (int i=0; i<mesh.numNodes(); i++) {
-		if (shape->isInside(mesh.nodes(i).pos))
-			mesh.nodes(i).flags |= Mesh::NfFixed;
-		else if (exclusive)
-			mesh.nodes(i).flags &= ~Mesh::NfFixed;
-	}
-}
-
-//! Adapt texture coordinates of mesh inside shape
-//! to obtain an effective inflow effect
-PYTHON void texcoordInflow(VortexSheetMesh& mesh, Shape* shape, MACGrid& vel) 
-{
-	static Vec3 t0 = Vec3::Zero;
-	
-	// get mean velocity
-	int cnt=0;
-	Vec3 meanV(_0);
-	FOR_IJK(vel) {
-		if (shape->isInsideGrid(i,j,k)) {
-			cnt++;
-			meanV += vel.getCentered(i,j,k);
-		}
-	}
-	meanV /= (Real) cnt;
-	t0 -= mesh.getParent()->getDt() * meanV;
-	mesh.setReferenceTexOffset(t0);
-
-	// apply mean velocity
-	for (int i=0; i<mesh.numNodes(); i++) {
-		if (shape->isInside(mesh.nodes(i).pos)) {
-			Vec3 tc = mesh.nodes(i).pos + t0;
-			mesh.tex1(i) = tc;
-			mesh.tex2(i) = tc;
-		}
-	}
-};
-
-//! Init smoke density values of the mesh surface inside source shape
-PYTHON void meshSmokeInflow(VortexSheetMesh& mesh, Shape* shape, Real amount) 
-{
-	for (int t=0; t<mesh.numTris(); t++) {
-		if (shape->isInside(mesh.getFaceCenter(t)))
-			mesh.sheet(t).smokeAmount = amount;
-	}    
-}
-
-KERNEL(idx) 
-void KnAcceleration(MACGrid& a, const MACGrid& v1, const MACGrid& v0, const Real idt) { 
-	a[idx] = (v1[idx]-v0[idx])*idt; 
-}
-
-//! Add vorticity to vortex sheets based on buoyancy
-PYTHON void vorticitySource(VortexSheetMesh& mesh, Vec3 gravity, 
-							MACGrid* vel=NULL, MACGrid* velOld=NULL,
-							Real scale = 0.1, Real maxAmount = 0, Real mult = 1.0)
-{
-	Real dt = mesh.getParent()->getDt();
-	Real dx = mesh.getParent()->getDx();
-	MACGrid acceleration(mesh.getParent());
-	if (vel)
-		KnAcceleration(acceleration, *vel, *velOld, 1.0/dt);
-	const Real A= -1.0;
-	Real maxV = 0, meanV = 0;
-	
-	for (int t=0; t<mesh.numTris(); t++) {
-		Vec3 fn = mesh.getFaceNormal(t);
-		Vec3 source;
-		if (vel) {
-			Vec3 a = acceleration.getInterpolated(mesh.getFaceCenter(t));        
-			source = A*cross(fn, a-gravity) * scale;
-		} else {
-			source = A*cross(fn, -gravity) * scale;
-		}
-		
-		if (mesh.isTriangleFixed(t)) source = 0;
-	
-		mesh.sheet(t).vorticity *= mult;
-		mesh.sheet(t).vorticity += dt * source / dx;
-		// upper limit
-		Real v = norm(mesh.sheet(t).vorticity);
-		if (maxAmount>0 && v > maxAmount)
-			mesh.sheet(t).vorticity *= maxAmount/v;
-		
-		//stats
-		if (v > maxV) maxV = v;
-		meanV += v;
-	}
-	
-	cout << "vorticity: max " << maxV << " / mean " << meanV/mesh.numTris() << endl;
-}
-
-PYTHON void smoothVorticity(VortexSheetMesh& mesh, int iter=1, Real sigma=0.2, Real alpha=0.8)
-{
-	const Real mult = -0.5 / sigma / sigma;
-	
-	// pre-calculate positions and weights
-	vector<Vec3> vort(mesh.numTris()), pos(mesh.numTris());    
-	vector<Real> weights(3*mesh.numTris());
-	vector<int> index(3*mesh.numTris());
-	for(int i=0; i<mesh.numTris(); i++) {
-		pos[i] = mesh.getFaceCenter(i);
-		mesh.sheet(i).vorticitySmoothed = mesh.sheet(i).vorticity;
-	}
-	for(int i=0; i<mesh.numTris(); i++) {
-		for (int c=0; c<3; c++) {
-			int oc = mesh.corners(i,c).opposite;
-			if (oc>=0) {
-				int t = mesh.corners(oc).tri;
-				weights[3*i+c] = exp(normSquare(pos[t]-pos[i])*mult);
-				index[3*i+c] = t;
-			}
-			else {
-				weights[3*i+c] = 0;
-				index[3*i+c] = 0;
-			}
-		}        
-	}
-		
-	for (int it=0; it<iter; ++it) {
-		// first, preload
-		for(int i=0; i<mesh.numTris(); i++) vort[i] = mesh.sheet(i).vorticitySmoothed;
-			
-		for(int i=0,idx=0; i<mesh.numTris(); i++) {            
-			// loop over adjacent tris
-			Real sum=1.0f;
-			Vec3 v=vort[i];
-			for (int c=0;c<3;c++,idx++) {
-				Real w = weights[index[idx]];
-				v += w*vort[index[idx]];
-				sum += w;
-			}
-			mesh.sheet(i).vorticitySmoothed = v/sum;
-		}
-	}
-	for(int i=0; i<mesh.numTris(); i++) mesh.sheet(i).vorticitySmoothed *= alpha;
-}
-
-//! Seed Vortex Particles inside shape with K41 characteristics
-PYTHON void VPseedK41(VortexParticleSystem& system, Shape* shape, Real strength=0, Real sigma0=0.2, Real sigma1=1.0, Real probability=1.0, Real N=3.0) {
-	Grid<Real> temp(system.getParent());
-	const Real dt = system.getParent()->getDt();
-	static RandomStream rand(3489572);
-	Real s0 = pow( (Real)sigma0, (Real)(-N+1.0) );
-	Real s1 = pow( (Real)sigma1, (Real)(-N+1.0) );
-	
-	FOR_IJK(temp) {
-		if (shape->isInsideGrid(i,j,k)) {
-			if (rand.getReal() < probability*dt) {
-				Real p = rand.getReal();
-				Real sigma = pow( (1.0-p)*s0 + p*s1, 1./(-N+1.0) );
-				Vec3 randDir (rand.getReal(), rand.getReal(), rand.getReal());
-				Vec3 posUpd (i+rand.getReal(), j+rand.getReal(), k+rand.getReal());
-				normalize(randDir);
-				Vec3 vorticity = randDir * strength * pow( (Real)sigma, (Real)(-10./6.+N/2.0) );
-				system.add(VortexParticleData(posUpd, vorticity, sigma)); 
-			}
-		}
-	}    
-}
-		
-//! Vortex-in-cell integration
-PYTHON void VICintegration(VortexSheetMesh& mesh, Real sigma, Grid<Vec3>& vel, FlagGrid& flags,
-					  Grid<Vec3>* vorticity=NULL, Real cgMaxIterFac=1.5, Real cgAccuracy=1e-3, Real scale = 0.01, int precondition=0) {
-	
-	MuTime t0;
-	const Real fac = 16.0; // experimental factor to balance out regularization
-	
-	// if no vort grid is given, use a temporary one
-	Grid<Vec3> vortTemp(mesh.getParent());    
-	Grid<Vec3>& vort = (vorticity) ? (*vorticity) : (vortTemp);
-	vort.clear();
-	
-	// map vorticity to grid using Peskin kernel
-	int sgi = ceil(sigma);
-	Real pkfac=M_PI/sigma;
-	const int numTris = mesh.numTris();
-	for (int t=0; t<numTris; t++) {
-		Vec3 pos = mesh.getFaceCenter(t);
-		Vec3 v = mesh.sheet(t).vorticity * mesh.getFaceArea(t) * fac;
-					
-		// inner kernel
-		// first, summate 
-		Real sum=0;
-		for (int i=-sgi; i<sgi; i++) {
-			if (pos.x+i < 0 || (int)pos.x+i >= vort.getSizeX()) continue;
-			for (int j=-sgi; j<sgi; j++) {
-				if (pos.y+j < 0 || (int)pos.y+j >= vort.getSizeY()) continue;            
-				for (int k=-sgi; k<sgi; k++) {
-					if (pos.z+k < 0 || (int)pos.z+k >= vort.getSizeZ()) continue;                                
-					Vec3i cell(pos.x+i, pos.y+j, pos.z+k);
-					if (!flags.isFluid(cell)) continue;
-					Vec3 d = pos - Vec3(i+0.5+floor(pos.x), j+0.5+floor(pos.y), k+0.5+floor(pos.z));
-					Real dl = norm(d);
-					if (dl > sigma) continue;
-					// precalc Peskin kernel
-					sum += 1.0 + cos(dl * pkfac);
-				}
-			}
-		}
-		// then, apply normalized kernel
-		Real wnorm = 1.0/sum;
-		for (int i=-sgi; i<sgi; i++) {
-			if (pos.x+i < 0 || (int)pos.x+i >= vort.getSizeX()) continue;
-			for (int j=-sgi; j<sgi; j++) {
-				if (pos.y+j < 0 || (int)pos.y+j >= vort.getSizeY()) continue;            
-				for (int k=-sgi; k<sgi; k++) {
-					if (pos.z+k < 0 || (int)pos.z+k >= vort.getSizeZ()) continue;                                
-					Vec3i cell(pos.x+i, pos.y+j, pos.z+k);  
-					if (!flags.isFluid(cell)) continue;                    
-					Vec3 d = pos - Vec3(i+0.5+floor(pos.x), j+0.5+floor(pos.y), k+0.5+floor(pos.z));
-					Real dl = norm(d);
-					if (dl > sigma) continue;
-					Real w = (1.0 + cos(dl * pkfac))*wnorm;
-					vort(cell) += v * w;
-				}
-			}
-		}
-	}
-	
-	// Prepare grids for poisson solve
-	Grid<Vec3> vortexCurl(mesh.getParent());
-	Grid<Real> rhs(mesh.getParent());
-	Grid<Real> solution(mesh.getParent());
-	Grid<Real> residual(mesh.getParent());
-	Grid<Real> search(mesh.getParent());
-	Grid<Real> temp1(mesh.getParent());
-	Grid<Real> A0(mesh.getParent());
-	Grid<Real> Ai(mesh.getParent());
-	Grid<Real> Aj(mesh.getParent());
-	Grid<Real> Ak(mesh.getParent());
-	Grid<Real> pca0(mesh.getParent());
-	Grid<Real> pca1(mesh.getParent());
-	Grid<Real> pca2(mesh.getParent());
-	Grid<Real> pca3(mesh.getParent());
-	
-	MakeLaplaceMatrix (flags, A0, Ai, Aj, Ak);    
-	CurlOp(vort, vortexCurl);    
-	
-	// Solve vector poisson equation
-	for (int c=0; c<3; c++) {
-		// construct rhs    
-		if (vel.getType() & GridBase::TypeMAC)
-			GetShiftedComponent(vortexCurl, rhs, c);
-		else
-			GetComponent(vortexCurl, rhs, c);
-				
-		// prepare CG solver
-		const int maxIter = (int)(cgMaxIterFac * vel.getSize().max());    
-		GridCgInterface *gcg = new GridCg<ApplyMatrix>(solution, rhs, residual, search, flags, temp1, &A0, &Ai, &Aj, &Ak );
-		gcg->setAccuracy(cgAccuracy); 
-		gcg->setUseResNorm(true);
-		gcg->setPreconditioner( (GridCgInterface::PreconditionType)precondition, &pca0, &pca1, &pca2, &pca3); 
-		
-		// iterations
-		for (int iter=0; iter<maxIter; iter++) {
-			if (!gcg->iterate()) iter=maxIter;
-		} 
-		debMsg("VICintegration CG iterations:"<<gcg->getIterations()<<", res:"<<gcg->getSigma(), 1);
-		delete gcg;
-	
-		// copy back
-		solution *= scale;
-		SetComponent(vel, solution, c);
-	}
-}
-
-//! Obtain density field from levelset with linear gradient of size sigma over the interface
-PYTHON void densityFromLevelset(LevelsetGrid& phi, Grid<Real>& density, Real value=1.0, Real sigma=1.0) {
-	FOR_IJK(phi) {
-		// remove boundary
-		if (i<2 || j<2 || k<2 || i>=phi.getSizeX()-2 || j>=phi.getSizeY()-2 || k>=phi.getSizeZ()-2)
-			density(i,j,k) = 0;
-		else if (phi(i,j,k) < -sigma)
-			density(i,j,k) = value;
-		else if (phi(i,j,k) > sigma)
-			density(i,j,k) = 0;
-		else
-			density(i,j,k) = clamp((Real)(0.5*value/sigma*(1.0-phi(i,j,k))), _0, value);
-	}    
-}
-
-} // namespace
\ No newline at end of file
diff --git a/source/blender/python/manta_pp/source/waveletturbulence.cpp b/source/blender/python/manta_pp/source/waveletturbulence.cpp
deleted file mode 100644
index 2b496da717c..00000000000
--- a/source/blender/python/manta_pp/source/waveletturbulence.cpp
+++ /dev/null
@@ -1,296 +0,0 @@
-/******************************************************************************
- *
- * MantaFlow fluid solver framework
- * Copyright 2011 Tobias Pfaff, Nils Thuerey 
- *
- * This program is free software, distributed under the terms of the
- * GNU General Public License (GPL) 
- * http://www.gnu.org/licenses
- *
- * Functions for calculating wavelet turbulence,
- * plus helpers to compute vorticity, and strain rate magnitude
- *
- ******************************************************************************/
-
-#include "vectorbase.h"
-#include "shapes.h"
-#include "commonkernels.h"
-#include "noisefield.h"
-
-using namespace std;
-
-namespace Manta {
-
-
-//! Apply vector noise to grid, this is a simplified version - no position scaling or UVs
-KERNEL 
-void knApplySimpleNoiseVec(FlagGrid& flags, Grid<Vec3>& target, WaveletNoiseField& noise, 
-					  Real scale, Grid<Real>* weight ) 
-{
-	if ( !flags.isFluid(i,j,k) ) return; 
-	Real factor = 1;
-	if(weight) factor = (*weight)(i,j,k);
-	target(i,j,k) += noise.evaluateCurl( Vec3(i,j,k) ) * scale * factor;
-}
-PYTHON void applySimpleNoiseVec3(FlagGrid& flags, Grid<Vec3>& target, WaveletNoiseField& noise, 
-							Real scale=1.0 , Grid<Real>* weight=NULL )
-{
-	// note - passing a MAC grid here is slightly inaccurate, we should evaluate each component separately
-	knApplySimpleNoiseVec(flags, target, noise, scale , weight );
-}
-
-
-//! Simple noise for a real grid , follows applySimpleNoiseVec3
-KERNEL 
-void knApplySimpleNoiseReal(FlagGrid& flags, Grid<Real>& target, WaveletNoiseField& noise, 
-					  Real scale, Grid<Real>* weight ) 
-{
-	if ( !flags.isFluid(i,j,k) ) return; 
-	Real factor = 1;
-	if(weight) factor = (*weight)(i,j,k);
-	target(i,j,k) += noise.evaluate( Vec3(i,j,k) ) * scale * factor;
-}
-PYTHON void applySimpleNoiseReal(FlagGrid& flags, Grid<Real>& target, WaveletNoiseField& noise, 
-							Real scale=1.0 , Grid<Real>* weight=NULL )
-{
-	knApplySimpleNoiseReal(flags, target, noise, scale , weight );
-}
-
-
-
-//! Apply vector-based wavelet noise to target grid
-//! This is the version with more functionality - supports uv grids, and on-the-fly interpolation
-//! of input grids.
-KERNEL 
-void knApplyNoiseVec(FlagGrid& flags, Grid<Vec3>& target, WaveletNoiseField& noise, 
-					  Real scale, Real scaleSpatial, Grid<Real>* weight, Grid<Vec3>* uv, bool uvInterpol, const Vec3& sourceFactor ) 
-{
-	if ( !flags.isFluid(i,j,k) ) return;
-
-	// get weighting, interpolate if necessary
-	Real w = 1;
-	if(weight) {
-		if(!uvInterpol) {
-			w = (*weight)(i,j,k);
-		} else {
-			w = weight->getInterpolated( Vec3(i,j,k) * sourceFactor );
-		}
-	}
-
-	// compute position where to evaluate the noise
-	Vec3 pos = Vec3(i,j,k);
-	if(uv) {
-		if(!uvInterpol) {
-			pos = (*uv)(i,j,k);
-		} else {
-			pos = uv->getInterpolated( Vec3(i,j,k) * sourceFactor );
-			// uv coordinates are in local space - so we need to adjust the values of the positions
-			pos /= sourceFactor;
-		}
-	}
-	pos *= scaleSpatial;
-
-	Vec3 noiseVec3 = noise.evaluateCurl( pos ) * scale * w; 
-	//noiseVec3=pos; // debug , show interpolated positions
-	target(i,j,k) += noiseVec3;
-} 
-PYTHON void applyNoiseVec3(FlagGrid& flags, Grid<Vec3>& target, WaveletNoiseField& noise, 
-							Real scale=1.0 , Real scaleSpatial=1.0 , Grid<Real>* weight=NULL , Grid<Vec3>* uv=NULL )
-{
-	// check whether the uv grid has a different resolution
-	bool uvInterpol = false; 
-	// and pre-compute conversion (only used if uvInterpol==true)
-	// used for both uv and weight grid...
-	Vec3 sourceFactor = Vec3(1.);
-	if(uv) {
-		uvInterpol = (target.getSize() != uv->getSize());
-		sourceFactor = calcGridSizeFactor( uv->getSize(), target.getSize() );
-	} else if(weight) {
-		uvInterpol = (target.getSize() != weight->getSize());
-		sourceFactor = calcGridSizeFactor( weight->getSize(), target.getSize() );
-	}
-	if(uv && weight) assertMsg( uv->getSize() == weight->getSize(), "UV and weight grid have to match!");
-
-	// note - passing a MAC grid here is slightly inaccurate, we should evaluate each component separately
-	knApplyNoiseVec(flags, target, noise, scale, scaleSpatial, weight , uv,uvInterpol,sourceFactor );
-}
-
-
-
-//! Compute energy of a staggered velocity field (at cell center)
-KERNEL 
-void KnApplyComputeEnergy( FlagGrid& flags, MACGrid& vel, Grid<Real>& energy ) 
-{
-	Real e = 0.f;
-	if ( flags.isFluid(i,j,k) ) {
-		Vec3 v = vel.getCentered(i,j,k);
-		e = 0.5 * v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
-	}
-	energy(i,j,k) = e;
-}
-
-PYTHON void computeEnergy( FlagGrid& flags, MACGrid& vel, Grid<Real>& energy )
-{
-	KnApplyComputeEnergy( flags, vel, energy );
-}
-
-
-
-//!interpolate grid from one size to another size
-KERNEL 
-void KnInterpolateGrid(Grid<Real>& target, Grid<Real>& source, const Vec3& sourceFactor)
-{
-	Vec3 pos = Vec3(i,j,k) * sourceFactor;
-	if(!source.is3D()) pos[2] = 0; // allow 2d -> 3d
-	target(i,j,k) = source.getInterpolated(pos);
-}
-
-PYTHON void interpolateGrid( Grid<Real>& target, Grid<Real>& source )
-{
-	Vec3 sourceFactor = calcGridSizeFactor( source.getSize(), target.getSize() );
-
-	// a brief note on a mantaflow specialty: the target grid has to be the first argument here!
-	// the parent fluidsolver object is taken from the first grid, and it determines the size of the
-	// loop for the kernel call. as we're writing into target, it's important to loop exactly over
-	// all cells of the target grid... (note, when calling the plugin in python, it doesnt matter anymore).
-
-	KnInterpolateGrid(target, source, sourceFactor);
-}
-
-
-//!interpolate a mac velocity grid from one size to another size
-KERNEL 
-void KnInterpolateMACGrid(MACGrid& target, MACGrid& source, const Vec3& sourceFactor)
-{
-	Vec3 pos = Vec3(i,j,k) * sourceFactor;
-
-	Real vx = source.getInterpolated(pos - Vec3(0.5,0,0))[0];
-	Real vy = source.getInterpolated(pos - Vec3(0,0.5,0))[1];
-	Real vz = 0.f;
-	if(source.is3D()) vz = source.getInterpolated(pos - Vec3(0,0,0.5))[2];
-
-	target(i,j,k) = Vec3(vx,vy,vz);
-}
-
-PYTHON void interpolateMACGrid(MACGrid& target, MACGrid& source)
-{
-	Vec3 sourceFactor = calcGridSizeFactor( source.getSize(), target.getSize() );
-
-	// see interpolateGrid for why the target grid needs to come first in the parameters!
-
-	KnInterpolateMACGrid(target, source, sourceFactor);
-}
-
-PYTHON void computeWaveletCoeffs(Grid<Real>& input)
-{
-	Grid<Real> temp1(input.getParent()), temp2(input.getParent());
-	WaveletNoiseField::computeCoefficients(input, temp1, temp2);
-}
-
-// note - alomst the same as for vorticity confinement
-PYTHON void computeVorticity(MACGrid& vel, Grid<Vec3>& vorticity, Grid<Real>* norm) {
-	Grid<Vec3> velCenter(vel.getParent());
-	GetCentered(velCenter, vel);
-	CurlOp(velCenter, vorticity);
-	if(norm) GridNorm( *norm, vorticity);
-}
-
-// note - very similar to KnComputeProductionStrain, but for use as wavelet turb weighting
-KERNEL(bnd=1) 
-void KnComputeStrainRateMag(const MACGrid& vel, const Grid<Vec3>& velCenter, Grid<Real>& prod ) 
-{
-	// compute Sij = 1/2 * (dU_i/dx_j + dU_j/dx_i)
-	Vec3 diag = Vec3(vel(i+1,j,k).x, vel(i,j+1,k).y, 0. ) - vel(i,j,k);
-	if(vel.is3D()) diag[2] += vel(i,j,k+1).z;
-	else           diag[2]  = 0.;
-
-	Vec3 ux =         0.5*(velCenter(i+1,j,k)-velCenter(i-1,j,k));
-	Vec3 uy =         0.5*(velCenter(i,j+1,k)-velCenter(i,j-1,k));
-	Vec3 uz;
-	if(vel.is3D()) uz=0.5*(velCenter(i,j,k+1)-velCenter(i,j,k-1));
-
-	Real S12 = 0.5*(ux.y+uy.x);
-	Real S13 = 0.5*(ux.z+uz.x);
-	Real S23 = 0.5*(uy.z+uz.y);
-	Real S2 = square(diag.x) + square(diag.y) + square(diag.z) +
-		2.0*square(S12) + 2.0*square(S13) + 2.0*square(S23);
-	prod(i,j,k) = S2;
-}
-PYTHON void computeStrainRateMag(MACGrid& vel, Grid<Real>& mag) {
-	Grid<Vec3> velCenter(vel.getParent());
-	GetCentered(velCenter, vel);
-	KnComputeStrainRateMag(vel, velCenter, mag);
-}
-
-
-// extrapolate a real grid into a flagged region (based on initial flags)
-// by default extrapolates from fluid to obstacle cells
-template<class T> 
-void extrapolSimpleFlagsHelper (FlagGrid& flags, Grid<T>& val, int distance = 4, 
-									int flagFrom=FlagGrid::TypeFluid, int flagTo=FlagGrid::TypeObstacle ) 
-{
-	Grid<int> tmp( flags.getParent() );
-	int dim = (flags.is3D() ? 3:2);
-	const Vec3i nb[6] = { 
-		Vec3i(1 ,0,0), Vec3i(-1,0,0),
-		Vec3i(0,1 ,0), Vec3i(0,-1,0),
-		Vec3i(0,0,1 ), Vec3i(0,0,-1) };
-
-	// remove all fluid cells (set to 1)
-	tmp.clear();
-	bool foundTarget = false;
-	FOR_IJK_BND(flags,0) {
-		if (flags(i,j,k) & flagFrom) 
-			tmp( Vec3i(i,j,k) ) = 1;
-		if (!foundTarget && (flags(i,j,k) & flagTo)) foundTarget=true;
-	}
-	// optimization, skip extrapolation if we dont have any cells to extrapolate to
-	if(!foundTarget) {
-		debMsg("No target cells found, skipping extrapolation", 1);
-		return;
-	}
-
-	// extrapolate for given distance
-	for(int d=1; d<1+distance; ++d) {
-
-		// TODO, parallelize
-		FOR_IJK_BND(flags,1) {
-			if (tmp(i,j,k) != 0)          continue;
-			if (!(flags(i,j,k) & flagTo)) continue;
-
-			// copy from initialized neighbors
-			Vec3i p(i,j,k);
-			int nbs = 0;
-			T avgVal = 0.;
-			for (int n=0; n<2*dim; ++n) {
-				if (tmp(p+nb[n]) == d) {
-					avgVal += val(p+nb[n]);
-					nbs++;
-				}
-			}
-
-			if(nbs>0) {
-				tmp(p) = d+1;
-				val(p) = avgVal / nbs;
-			}
-		}
-
-	} // distance 
-}
-PYTHON void extrapolateSimpleFlags (FlagGrid& flags, GridBase* val, int distance = 4, 
-									int flagFrom=FlagGrid::TypeFluid, int flagTo=FlagGrid::TypeObstacle ) 
-{
-	if (val->getType() & GridBase::TypeReal) {
-		extrapolSimpleFlagsHelper<Real>(flags,*((Grid<Real>*) val),distance,flagFrom,flagTo);
-	}
-	else if (val->getType() & GridBase::TypeInt) {    
-		extrapolSimpleFlagsHelper<int >(flags,*((Grid<int >*) val),distance,flagFrom,flagTo);
-	}
-	else if (val->getType() & GridBase::TypeVec3) {    
-		extrapolSimpleFlagsHelper<Vec3>(flags,*((Grid<Vec3>*) val),distance,flagFrom,flagTo);
-	}
-	else
-		errMsg("extrapolateSimpleFlags: Grid Type is not supported (only int, Real, Vec3)");    
-}
-
-} // namespace