integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G11, new interparticlePotential::PsiShanChen93<T>(rho0) ), lattice.getBoundingBox(), latticeOne, processorLevel );

integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G22, new interparticlePotential::PsiShanChen93<T>(rho0) ),lattice.getBoundingBox(), latticeTwo, processorLevel );

integrateProcessingFunctional (new ShanChenMultiComponentProcessor2D<T,DESCRIPTOR>(G12,constOmegaValues), Box2D(1,nx-2,0,ny-1),
            blockLattices, processorLevel );

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");
.........

     // Main loop over time iterations.
    for (int iT=0; iT<maxIter; ++iT) {

     if (iT > 20){
       Array<T,2> externalforce (0.,0.001);
       setExternalField (insideFluid, Box2D(0, nx-1, 0, ny-1), 1, &externalforce[0]);
     }

void Dynamics<T,Descriptor>::setExternalField (
        Cell<T,Descriptor>& cell, plint pos, plint size, const T* ext)

template<typename T, template<typename U> class Descriptor>
class MyFunctional : public BoxProcessingFunctional3D_S<T>  {
  private:
    plint z;
    plint xmax, ymax;
    plint area;
  public:
    MyFunctional (plint z_, plint xmax_, plint ymax_) : z(z_), xmax(xmax_), ymax(ymax_), area(0) {
      }
      
    virtual void process (MultiScalarField3D<bool>& matrix) {
      for (plint i=0; i<xmax; i++) {
	for (plint j=0; j<ymax; j++) {
	  if (matrix.get(i, j, z)) area++;
	}
      }
    }
    
    virtual MyFunctional<T, Descriptor>* clone() const  {
      return new MyFunctional<T, Descriptor> (*this);
    }
};


void CreateMyFunctional (MultiScalarField3D<bool> &matrix, plint z, plint xmax, plint ymax)	{
  Box3D foo(0,xmax,0,ymax,z,z);
  applyProcessingFunctional (new MyFunctional<T, DESCRIPTOR> (z, xmax, ymax), foo, matrix);
  
}

      PROGRAM COUETTE
      DIMENSION U(101)
      COMMON/TDMA/A(101),B(101),C(101),D(101)
C-----IMAX IS NUMBER OF GRID LINES
      IMAX=20
      WRITE(6,*) 'HERE 1'
C-----ISKIP IS THE LINES SKIPPED BETWEEN PRINTOUT
       ISKIP=1
       IM1=IMAX-1
C-----BOUNDRARY CONDITIONS
       U(1)=0.0
       U(IMAX)=1.0
C-----EVALUATION OF TDMA COEFFICIENTS
       DY=1./FLOAT(IM1)
       DO 10 I=2,IM1
          AI=FLOAT(I)
          Y=(AI-1.5)*DY
          A(I)=EXP(1.92*(1.-Y))
          Y=(AI-.5)*DY
          C(I)=EXP(1.92*(1.-Y))
          B(I)=-A(I)-C(I)
  10      D(I)=0.0
      WRITE(6,*) 'HERE 1'
C-----CALL THE TRIDIAGONAL ALGORITHM TO OBTAIN VELOCITIES
      CALL TRIDIA(1,IMAX,U)
C-----PRINTOUT OF VALUES
       WRITE(6,1000)
       WRITE(6,1001) (I,U(I),I=1,IMAX,ISKIP)
       PAUSE
1000   FORMAT('VELOCITY DISTIBUTION BETWEEN PLATES *****'
     1/'       I         U(I)')
1001   FORMAT(I4,5X,F10.5)
       END
C
C*****TRIDIAGONAL SUBSUBROUTINE WHERE 'IT' IS THE LAST LINE, 'IB' IS
C     THE FIRST LINE AND 'PHI' IS THE DEPENDENT VARIABLE THE
C     NOMENCLATURE CORRESPONDS TO EQNS 16-48 THRU 16-53 IN THE C NEXT
      SUBROUTINE TRIDIA (IB,IT,PHI)
      DIMENSION E(101),F(101),PHI(IT)
      COMMON/TDMA/A(101),B(101),C(101),D(101)
       IM1=IT-1
       IP1=IB+1
       IP2=IB+2
       WRITE(6,*) 'HERE 1'
C-----SETTING UP BOUNDARY VALUES
       D(IP1)=D(IP1)-A(IP1)*PHI(IB)
       A(IP1)=0.0
       D(IM1)=D(IM1)-C(IM1)*PHI(IT)
       C(IM1)=0.0
       E(IP1+1)=-C(IP1)/B(IP1)
       F(IP1+1)=D(IP1)/B(IP1)
C-----EVALUATION OF E'S AND F'S
       DO 20 K=IP2,IM1
          FAC=1./(A(K)*E(K)+B(K))
          E(K+1)=-C(K)*FAC
   20     F(K+1)=(D(K)-A(K)*F(K))*FAC
C-----USING E'S AND F'S TO CALCULATE PHI
       PHI(IM1)=F(IT)
       DO 30 L=IP2,IM1
          LL=IT+IP1-L
   30     PHI(LL-1)=E(LL)*PHI(LL)+F(LL)
       RETURN
       END

This file is part of the Palabos library.
* Copyright (C) 2009 Jonas Latt
* E-mail contact: jonas@lbmethod.org
* The most recent release of Palabos can be downloaded at 
* [www.lbmethod.org]
*
* The library Palabos is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* The library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see [www.gnu.org].
*/

#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>
#include <iomanip>

using namespace plb;
using namespace std;

typedef double T;


#define DESCRIPTOR descriptors::ShanChenD2Q9Descriptor



void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
			  plint iter)
{
	//T dx = parameters.getDeltaX();
	// T dt = parameters.getDeltaT();
	VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 6), 1.0);
	vtkOut.writeData<float>(*computeDensity(lattice), "Density", 1.0);
}
int main(int argc, char *argv[])
{
	plbInit(&argc, &argv);
	global::timer("simTime").start();
	global::directories().setOutputDir("./tmp/");

	OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition = createEquilibriumBoundaryCondition2D<T,DESCRIPTOR>();


	// For the choice of the parameters G, rho0, and psi0, we refer to the book
	//   Michael C. Sukop and Daniel T. Thorne (2006), 
	//   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
	//   Springer-Verlag Berlin/Heidelberg.


	T tIni = global::timer("simTime").stop();
	global::timer("simTime").start();

	const T omega = 1.0;
	const int nx   = 272;
	const int ny   = 272;
	//const int nz   = 32;
	const T G      = -92.0; //setting the force to zero
	const int vapour_inlet=10;

	const int maxIter  = 100001;//100001;
	const int saveIter = 10;
	const int statIter = 10;
	plint evalTime =100;


	const T liquid_fluid = 200.0;
	const T psi0 = 4.0;
	const T vel=0.05;
	const T virt_density=liquid_fluid+5.0;
	const T vapour_fluid=40.0;
	const int inlet=10;


	pcout << "Reading the geometry file." << endl;
	std::string fNameIn  = "SOslice.dat";
	MultiScalarField2D<T> geometry(nx,ny);
	plb_ifstream geometryFile(fNameIn.c_str());
	if (!geometryFile.is_open()) {
		pcout << "Error: could not open geometry file " << fNameIn << endl;
		return -1;
	}
	geometryFile >> geometry;
	//reset geometry file everything is a pore
	//geometry.reset();

	MultiBlockLattice2D<T, DESCRIPTOR> lattice ( nx,ny, new ExternalMomentRegularizedBGKdynamics <T, DESCRIPTOR>(omega) );

	lattice.periodicity().toggle(0, false);	
	lattice.periodicity().toggle(1, false);

	// Add the data processor which implements the Shan/Chen interaction potential.
	plint processorLevel = 1;
	integrateProcessingFunctional (
		new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
		G, new interparticlePotential::PsiShanChen94<T>(psi0,liquid_fluid) ),
		lattice.getBoundingBox(),
		lattice,
		processorLevel );

	initializeAtEquilibrium (
		lattice, Box2D(0, inlet, 0, ny-1),liquid_fluid,
		Array<T,2>(0.0,0.0) );
	initializeAtEquilibrium (
		lattice, Box2D(inlet+1, nx-1, 0,ny-1),vapour_fluid,
		Array<T,2>(0.0, 0.0) );

		//OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition = createEquilibriumBoundaryCondition2D<T,DESCRIPTOR>();


	//y boundaries
	setBoundaryVelocity(lattice, Box2D(0, inlet, ny-1, ny-1), Array<T,2>(0.0,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(0, inlet, ny-1, ny-2),liquid_fluid,
		Array<T,2>(0.0,0.0) );
	setBoundaryVelocity(lattice, Box2D(inlet+1, nx-1, ny-1, ny-1), Array<T,2>(0.0,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(inlet+1, nx-1, ny-1, ny-2),vapour_fluid,
		Array<T,2>(0.0,0.0) );


	setBoundaryVelocity(lattice, Box2D(0, inlet, 0, 0), Array<T,2>(0.0,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(0, inlet, 0, 0),liquid_fluid,
		Array<T,2>(0.0, 0.0) );
	setBoundaryVelocity(lattice, Box2D(inlet+1, nx-1, 0, 0), Array<T,2>(0.0,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(inlet+1, nx-1, 0, 0),vapour_fluid,
		Array<T,2>(0.0, 0.0) );

	//x boundaries
	setBoundaryVelocity(lattice, Box2D(nx-1, nx-1, 0, ny-1), Array<T,2>(vel,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(nx-1, nx-1, 1, ny-2),vapour_fluid,
		Array<T,2>(vel,0.0) );//0.00005
	setBoundaryVelocity(lattice, Box2D(0, 0, 0,ny-1), Array<T,2>(vel,0.) );
	initializeAtEquilibrium (
		lattice, Box2D(0, 0, 1,ny-2),liquid_fluid,
		Array<T,2>(vel, 0.0) );

	// Where "geometry" evaluates to 1, use bounce-back.//0.5 adhesive force
	defineDynamics(lattice, geometry, new BounceBack<T, DESCRIPTOR>(virt_density), true);
	defineDynamics(lattice, Box2D(Box2D(0, nx-1, 0, 0)), new BounceBack<T, DESCRIPTOR>(virt_density));
	defineDynamics(lattice, Box2D(Box2D(0, nx-1, ny-1, ny-1)), new BounceBack<T, DESCRIPTOR>(virt_density));



	lattice.initialize();

	pcout << "Starting simulation" << endl;
	// The value-tracer is used to stop the simulation once is has converged.
	// 1st parameter:velocity
	// 2nd parameter:size
	// 3rd parameters:threshold
	// 1st and second parameters ae used for the length of the time average (size/velocity)
	util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);
	int iT=1;

	for (int iT; iT<maxIter; ++iT) {


	//initializeAtEquilibrium (
		//lattice, Box2D(0, 10, 0, ny-1),200.0,
		//Array<T,2>(0.05,0.0) );
		//lattice, Box2D(0, 20, 0, ny-1),200.0 );
		if (iT%statIter==0) {
			auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
			converge.takeValue(getStoredAverageEnergy(lattice),true);
			pcout << "iteration "  << iT << endl;
			pcout << "Average density: " << computeAverage(*computeDensity(lattice, Box2D(0, nx-1, 0, ny-1 ))) << endl;
			pcout << "Minimum density: " << computeMin(*rho) << endl;
			pcout << "Maximum density: " << computeMax(*rho) << endl;

		}
		if (iT %saveIter==0) {
			ImageWriter<T>("leeloo.map").writeGif (
			createFileName("rho", iT, 6),*computeDensity(lattice, Box2D(0, nx-1, 0, ny-1 )), vapour_fluid, virt_density );			               //writeVTK(lattice,iT);
			char path[200]="density";//createFileName("density", iT, 6)
			char iter[100]="";
			sprintf(iter, "%d", iT);
			strcat(path, iter);			
			strcat(path, ".dat");
			pcout << path << endl;
			plb_ofstream ofile(path);
			ofile << setprecision(10) << *computeDensity(lattice) << endl;

		}
		lattice.collideAndStream();
		converge.takeValue(getStoredAverageEnergy(lattice),true);
		if (converge.hasConverged())
		{
			// break;
		}
	}
	//convert(convert -delay 5 rho*.gif animation.gif);
	T tEnd = global::timer("simTime").stop();
	T totalTime = tEnd-tIni;
	pcout << "number of processors: " << global::mpi().getSize() << endl;
	pcout << "simulation time: " << totalTime << endl;

}

tutorial_2_2.cpp: In function ‘int main(int, char**)’:
tutorial_2_2.cpp:128: error: ‘MultiBlockDistribution2D’ was not declared in this scope
tutorial_2_2.cpp:128: error: expected `;' before ‘blockDistribution’
tutorial_2_2.cpp:129: error: ‘blockDistribution’ was not declared in this scope
tutorial_2_2.cpp:139: error: expected primary-expression before ‘>’ token
tutorial_2_2.cpp:139: error: expected primary-expression before ‘)’ token
tutorial_2_2.cpp:140: error: expected primary-expression before ‘>’ token
tutorial_2_2.cpp:140: error: expected primary-expression before ‘)’ token

python ../../../scons/scons.py -f ../../../SConstruct palabosRoot=../../.. projectFiles="tutorial_2_2.cpp" precompiled=false optimize=true debug=false profile=false MPIparallel=false SMPparallel=false usePOSIX=true useCVMLCPP=false serialCXX=g++ parallelCXX=mpicxx compileFlags="" linkFlags="" optimFlags="-O3" debugFlags="-g" profileFlags="-pg" libraryPaths="" includePaths="" libraries=""
scons: Reading SConscript files ...
scons: done reading SConscript files.
scons: Building targets ...
g++ -o tutorial_2_2.o -c -O3 -DPLB_USE_POSIX -~/konstantin/palabos-v0.7r0/src -~/konstantin/palabos-v0.7r0/externalLibraries tutorial_2_2.cpp
tutorial_2_2.cpp: In function ‘int main(int, char**)’:
tutorial_2_2.cpp:128: error: ‘MultiBlockDistribution2D’ was not declared in this scope
tutorial_2_2.cpp:128: error: expected `;' before ‘blockDistribution’
tutorial_2_2.cpp:129: error: ‘blockDistribution’ was not declared in this scope
tutorial_2_2.cpp:139: error: expected primary-expression before ‘>’ token
tutorial_2_2.cpp:139: error: expected primary-expression before ‘)’ token
tutorial_2_2.cpp:140: error: expected primary-expression before ‘>’ token
tutorial_2_2.cpp:140: error: expected primary-expression before ‘)’ token
scons: *** [tutorial_2_2.o] Error 1
scons: building terminated because of errors.
make: *** [compile] Error 2

# HG changeset patch
# User arne@brand.medma
# Date 1276505021 -7200
# Branch babtests
# Node ID 87a577c922aa66a93d282b6f7c3e2b8e0363fc00
# Parent  4fe639a63a13fb937417a5eb228db924a5905e56
Added greyscale output

diff -r 4fe639a63a13 -r 87a577c922aa src/io/colormaps.cpp
--- a/src/io/colormaps.cpp	Do Jun 10 17:43:00 2010 +0200
+++ b/src/io/colormaps.cpp	Mo Jun 14 10:43:41 2010 +0200
@@ -324,6 +324,25 @@
     return leeLooBlue;
 }
 
+PiecewiseFunction generateGreyRGB() {
+    double p0  =  0.;
+    double p2  =  3./8.;
+    double p3  =  5./8.;
+    double p4  =  7./8.;
+    double p5  =  1.;
+    double p6  =  9./8.;
+
+
+    PiecewiseFunction greyRGB;
+    greyRGB.addPiece (Piece(p0, p2), new LinearFunction(p0, p2, 1., 0.3) );
+    greyRGB.addPiece (Piece(p2, p3), new LinearFunction(p2, p3, 0.3, 0.2) );
+    greyRGB.addPiece (Piece(p3, p4), new LinearFunction(p3, p4, 0.2, 0.1) );
+    greyRGB.addPiece (Piece(p4, p5), new LinearFunction(p4, p6, 0.1, 0.) );
+
+    return greyRGB;
+}
+
+
 ColorMap generateMap(std::string mapName) {
     if (mapName == "earth") {
         return ColorMap (
@@ -355,6 +374,12 @@
                 generateLeeLooGreen(),
                 generateLeeLooBlue() );
     }
+    else if (mapName == "grey") {
+        return ColorMap (
+                generateGreyRGB(),
+                generateGreyRGB(),
+                generateGreyRGB() );
+    }
     return ColorMap (
             generateLeeLooRed(),
             generateLeeLooGreen(),

[PC183:14967] *** Process received signal ***
[PC183:14967] Signal: Segmentation fault (11)
[PC183:14967] Signal code: Address not mapped (1)
[PC183:14967] Failing at address: 0xbf0bb60c
[PC183:14967] [ 0] [0xb77da410]
[PC183:14967] [ 1] finderror02(_ZN3plb14BlockLattice3DIdNS_11descriptors21ForcedD3Q19DescriptorEE16collideAndStreamENS_5Box3DE+0x30e) [0x80ee19e]
[PC183:14967] [ 2] finderror02(_ZN3plb19MultiBlockLattice3DIdNS_11descriptors21ForcedD3Q19DescriptorEE16collideAndStreamEv+0xe2) [0x80e64a2]