Updated manual pages, added working example (tutorial).
authorAndras Laszlo <laszlo.andras@wigner.mta.hu>
Fri, 28 Aug 2015 17:05:57 +0000 (19:05 +0200)
committerAndras Laszlo <laszlo.andras@wigner.mta.hu>
Fri, 28 Aug 2015 17:05:57 +0000 (19:05 +0200)
19 files changed:
include/gridripper/phys/example/kgminkowski/Factory.h [new file with mode: 0644]
include/gridripper/phys/example/kgminkowski/Hunch.h [new file with mode: 0644]
include/gridripper/phys/example/kgminkowski/KGMinkowski.h [new file with mode: 0644]
include/gridripper/phys/example/kgminkowski/packages.dox [new file with mode: 0644]
include/gridripper/phys/example/package.dox [new file with mode: 0644]
inputs/example/kgminkowski/analysis/Makefile [new file with mode: 0644]
inputs/example/kgminkowski/analysis/bin/arg.list [new file with mode: 0644]
inputs/example/kgminkowski/analysis/bin/go.sh [new file with mode: 0755]
inputs/example/kgminkowski/analysis/doc/README [new file with mode: 0644]
inputs/example/kgminkowski/analysis/inc/config.h [new file with mode: 0644]
inputs/example/kgminkowski/analysis/inc/pstream.h [new file with mode: 0644]
inputs/example/kgminkowski/analysis/inp/intq.conf [new file with mode: 0644]
inputs/example/kgminkowski/analysis/lib/config.cc [new file with mode: 0644]
inputs/example/kgminkowski/analysis/src/intq.cc [new file with mode: 0644]
inputs/example/kgminkowski/kgminkowski.in [new file with mode: 0644]
man/gridripperapp.man
src/phys/example/kgminkowski/Factory.cxx [new file with mode: 0644]
src/phys/example/kgminkowski/Hunch.cxx [new file with mode: 0644]
src/phys/example/kgminkowski/KGMinkowski.cxx [new file with mode: 0644]

diff --git a/include/gridripper/phys/example/kgminkowski/Factory.h b/include/gridripper/phys/example/kgminkowski/Factory.h
new file mode 100644 (file)
index 0000000..7806553
--- /dev/null
@@ -0,0 +1,19 @@
+#ifndef gridripper_phys_example_kgminkowski_Factory_h
+#define gridripper_phys_example_kgminkowski_Factory_h
+#include <gridripper/factory_inc.h>
+
+/**
+ * Factory for wave equation over Minkowski spacetime 
+ * (multipole expansion is used).
+ *
+ * @version 1.1.1, 08/28/2015
+ * @since   GridRipper 1.1.1, 08/28/2015
+ * @author  Andras Laszlo
+ */
+namespace gridripper { namespace phys { namespace example {
+namespace kgminkowski {
+#include <gridripper/factory_decl.h>
+} } } } // namespace gridripper::phys::example::kgminkowski;
+
+
+#endif /* gridripper_phys_example_kgminkowski_Factory_h */
diff --git a/include/gridripper/phys/example/kgminkowski/Hunch.h b/include/gridripper/phys/example/kgminkowski/Hunch.h
new file mode 100644 (file)
index 0000000..14acd3b
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef gridripper_phys_example_kgminkowski_Hunch_h
+#define gridripper_phys_example_kgminkowski_Hunch_h
+
+
+#include <gridripper/Parameters.h>
+#include <gridripper/amr1d/FuncInitCond.h>
+#include <gridripper/amr1d/PDE.h>
+#include <gridripper/amr1d/FuncInitCond.h>
+
+
+namespace gridripper { namespace phys { namespace example {
+namespace kgminkowski {
+
+
+using namespace gridripper;
+using namespace gridripper::util;
+using namespace gridripper::amr1d;
+using namespace std;
+
+
+/**
+ * Wave packet initial condition for wave equation with 
+ * multipole expansion method over Minkowski spacetime.
+ *
+ * @version 1.1.1, 08/28/2015
+ * @since   1.1.1, 08/28/2015
+ * @author  Andras Laszlo
+ */
+class Hunch: public FuncInitCond
+{
+private:
+    /** Maximal multipole order. */
+    int maxOrder;
+    /** Array size of multipole vector. */
+    int arraySize;
+
+    // The evolved field is f=Phi*r.
+
+    /** Wave packet parameters. */
+    GReal_t amplitudeF, centerF, widthF, rampF;
+    /** Angular momentum quantum numbers. */
+    int lF, mF;
+    /** Wave packet parameters. */
+    GReal_t amplitudeFt, centerFt, widthFt, rampFt;
+    /** Angular momentum quantum numbers. */
+    int lFt, mFt;
+    /** Parameter omega telling t phase. */
+    GReal_t omega;
+
+    /** Minimal and maximal spatial coordinates. */
+    GReal_t minr;
+    GReal_t maxr;
+
+    /** Initial time. */
+    GReal_t t0;
+
+    /** Message string. */
+    string message;
+
+public:
+    Hunch(string& args, const Parameters* p, const PDE& pde)
+         throw(InitCond::Exception&, IllegalArgumentException&);
+    ~Hunch();
+    GReal_t getTime() const;
+
+    void init(PDE* pde, GReal_t r, FieldWrapper& v)
+       throw(InitCond::Exception&);
+
+    string getMessage() const;
+};
+
+
+} } } } // namespace gridripper::phys::example::kgminkowski
+
+
+#endif /* gridripper_phys_example_kgminkowski_Hunch_h */
diff --git a/include/gridripper/phys/example/kgminkowski/KGMinkowski.h b/include/gridripper/phys/example/kgminkowski/KGMinkowski.h
new file mode 100644 (file)
index 0000000..f556e74
--- /dev/null
@@ -0,0 +1,206 @@
+#ifndef gridripper_phys_example_kgminkowski_KGMinkowski_h
+#define gridripper_phys_example_kgminkowski_KGMinkowski_h
+
+
+#include <gridripper/Parameters.h>
+#include <gridripper/amr1d/PDE.h>
+#include <gridripper/lang/IllegalArgumentException.h>
+#include <gridripper/multipole/ScalarFieldMP.h>
+
+#include <time.h>
+
+namespace gridripper { namespace phys { namespace example {
+namespace kgminkowski {
+
+
+using namespace gridripper;
+using namespace gridripper::multipole;
+using namespace std;
+
+
+/**
+ * Wave equation, using multipole expansion over Minkowski spacetime.
+ * Coordinates \f$t\f$, \f$r\f$, \f$\vartheta\f$, \f$\varphi\f$ 
+ * (polar). Scalar field: \f$\Phi\f$.
+ * 
+ * Lagrange form (with 1,-1,-1,-1 metric sign convention):
+ * \f$\left(
+ *     g^{ab}\overline{\mathrm{i}D_{a}(\Phi)}\mathrm{i}D_{b}(\Phi) 
+ *     -V(\Phi)\right)\mathrm{d}v_{g}\f$.
+ * Here \f$V(\Phi):=0\f$ and 
+ * \f$D=\nabla\f$.
+ * Field equation:
+ * \f$g^{ab}\mathrm{i}D_{a}\mathrm{i}D_{b}\Phi=V'(\Phi)\f$.
+ * Here \f$V'(\Phi)=0\f$.
+ * \f$T_{ab}=\mathrm{Re}(\overline{\mathrm{i}D_{a}(\Phi)}\mathrm{i}D_{b}(\Phi))
+ *     -\frac{1}{2}g_{ab}\left(g^{cd}
+ *     \bar{\mathrm{i}D_{c}(\Phi)}\mathrm{i}D_{d}(\Phi)-V(\Phi)
+ *     \right)\f$.
+ *
+ * \f$g=1 dt\vee dt
+ *     +0 dt\vee d\varphi
+ *     -0 dr\vee dr
+ *     -r^2 d\vartheta\vee d\vartheta
+ *     -r^2\sin^{2}{\vartheta} d\varphi\vee d\varphi\f$.
+ *
+ * The actually evolved field is \f$r\Phi\f$ in the according 
+ * transformed field equation.
+ *
+ * @version 06/22/2009
+ * @since   GridRipper 0.5, 05/01/2007
+ * @author  Andras Laszlo
+ */
+class KGMinkowski: public PDE
+{
+public:
+
+    /** Error tolerance for checking equality of floating point numbers. */
+    static const GReal_t EPSILON;
+
+    /** Field component indices. */
+    enum { ind_f=0, ind_ft=1, ind_fr=2 };
+
+    /** Maximal multipole order (>=0), calculated by the Factory. */
+    const int maxOrder;
+    /**
+     * Size of a multipole array, 
+     * (maxOrder+1)*(maxOrder+1)*sizeof(GComplex_t)/sizeof(GReal_t).
+     * Calculated by the Factory.
+     */
+    const int arraySize;
+
+    /** Indices of field component sectors (calculated by the constructor). */
+    const int I_f;     // <- ind_f*arraySize
+    const int I_ft;    // <- ind_ft*arraySize
+    const int I_fr;    // <- ind_fr*arraySize
+    /** Number of independent components. */
+    static const int NUM_INDEP_COMPS=2;
+    /** Number of constraints. */
+    static const int NUM_CONSTRAINTS=1;
+    /** Number of all components (independent+constraints). */
+    static const int NUM_COMPS=NUM_INDEP_COMPS+NUM_CONSTRAINTS;
+
+    /** Component and constraint names, and coordinate labels (see code). */
+    static const string COMP_NAMES[NUM_COMPS];
+    static const string CONSTRAINT_NAMES[NUM_CONSTRAINTS];
+    static const string COORD_LABELS[2];
+
+private:
+
+    /** Result vector (time derivative). */
+    tvalarray<GReal_t> df;
+
+    /** Buffers for regularization trick around the origin. */
+    static const int NUM_AUX_POINTS=4;
+    tvector< FieldComponents<GReal_t> > fieldat;
+
+    // The field is f=Phi*r.
+
+    /** Minimal and maximal spatial coordinate. */
+    GReal_t minr;
+    GReal_t maxr;
+
+    /** Is time reversed? (Backward evolution?) */
+    bool isTimeReversed;
+
+    /* Start time. */
+    GReal_t t0;
+
+    /** Bits to store the exclude flags (to forbid dissipation term). */
+    int exclude;
+
+public:
+
+    KGMinkowski(const Parameters* p, const int maxorder, const int arraysize) throw(IllegalArgumentException&);
+
+    ~KGMinkowski();
+
+    GReal_t getMinX() const;
+    GReal_t getMaxX() const;
+    GReal_t getPhysicalMinX() const;
+    GReal_t getPhysicalMaxX() const;
+    bool isXPeriodic() const;
+
+    /**
+     * Evaluate the field equation. This is the key algorithm.
+     * 
+     * @param  dF   : result array pointer
+     * @param  t    : temporal coordinate
+     * @param  r    : spatial coordinate
+     * @param  F    : the field array at the given point
+     * @param  F_r  : the field spatial derivative at the given point
+     * @param  G    : the grid
+     * @param  i    : the grid index of the given spatial point
+     * @param  d    : the spatial differentiation scheme operator
+     * @param  dt   : the actual temporal step
+     * @return exclusion bitmask (to kill dissipation term for some components)
+     */
+    int eval(GReal_t* dF, GReal_t t, GReal_t r,
+            const GReal_t* F, const GReal_t* F_r,
+            const Grid& G, int i, Grad& d, GReal_t dt);
+
+    /**
+     * Evaluate the constrained components. This is also a key algorithm.
+     * Only executed for the initial data! I.e. to be used only 
+     * for conserved constraints!
+     * 
+     * @param  G : the grid
+     * @param  i : the position index
+     * @param  v : result array
+     */
+    void evalConstrainedComponents(Grid& G, int i, FieldWrapper& v);
+
+     /**
+     * Sets nonevolving field components.
+     * Called after the initial data generation and after each time step.
+     * Can be used for non-conserved constraints.
+     * 
+     * @param  G : the grid 
+     * @param  d : differentiation scheme
+     */
+    void setNonevolvingComponents(Grid& G, Grad& d);
+
+    FieldWrapper* createField(int level) const;
+
+    /**
+     * Can the field be extended beyond the specified bound?
+     * 
+     * @param  where : the grid index
+     * @return <code>true</code> : if can be extended at the point
+     */
+    bool canBeExtended(int where) const;
+
+    /** Extend left side. */
+    void mirrorLeft(GReal_t r, FieldWrapper& f) const;
+
+    /** Extend right side. */
+    void mirrorRight(GReal_t r, FieldWrapper& f) const;
+
+    /**
+     * Tests whether refined field data exists at a given point.
+     * 
+     * @param  G : the grid
+     * @param  i : refined grid point index
+     * @return <code>true</code> : if the point contains refined data
+     */
+    bool hasExtendedRefinedFieldData(const Grid& G, int i) const;
+
+    /** Get the density quantities. */
+    PhysicalQuantityArray getPhysicalQuantities() const;
+
+    /** Get the marker quantities. */
+    MarkerQuantityArray getMarkerQuantities() const;
+
+private:
+
+    /** Some helper function for labeling. */
+    static tvector<string> componentNames(int maxOrder);
+    static tvector<string> constraintNames(int maxOrder);
+
+};
+
+
+} } } } // namespace gridripper::phys::example::kgminkowski
+
+
+#endif /* gridripper_phys_example_kgminkowski_KGMinkowski_h */
diff --git a/include/gridripper/phys/example/kgminkowski/packages.dox b/include/gridripper/phys/example/kgminkowski/packages.dox
new file mode 100644 (file)
index 0000000..c9d34cc
--- /dev/null
@@ -0,0 +1,7 @@
+/** \dir gridripper/phys/example/kgminkowski
+Wave equation over 3+1 dimensional Minkowski spacetime.
+*/
+
+/** \namespace gridripper::phys::example::kgminkowski
+Wave equation over 3+1 dimensional Minkowski spacetime.
+*/
diff --git a/include/gridripper/phys/example/package.dox b/include/gridripper/phys/example/package.dox
new file mode 100644 (file)
index 0000000..e98f9a6
--- /dev/null
@@ -0,0 +1,7 @@
+/** \dir gridripper/phys/example
+PDE classes for educative examples.
+*/
+
+/** \namespace gridripper::phys::example
+PDE classes for educative examples.
+*/
diff --git a/inputs/example/kgminkowski/analysis/Makefile b/inputs/example/kgminkowski/analysis/Makefile
new file mode 100644 (file)
index 0000000..6c50a40
--- /dev/null
@@ -0,0 +1,183 @@
+########################################################################
+########################################################################
+#                                                                      #
+# This is an automatic 'Makefile'. One generally does not have to edit #
+# this file. Just put the program sources (source files with the       #
+# program body) into the 'src' directory, put the library sources      #
+# (source files to be linked to the programs, or the compiled          #
+# libraries or shared libraries themselves) into the 'lib'             #
+# directory, put the header sources (source files to be included in)   #
+# into the 'inc' directory. You can directly refer to include files in #
+# 'inc' directory, as they are added to the include search path.       #
+# Then, by saying 'make' (or 'make bin/targetfile', where              #
+# 'targetfile' is substituted by an appropriate name) source is        #
+# compiled.                                                            #
+#                                                                      #
+# File naming conventions:                                             #
+# - header files: '*.h' for C, '*.h', '*.H' for C++,                   #
+#   '*.inc' for FORTRAN,                                               #
+# - source files: '*.c' for C, '*.cc', '*.cpp', '*.cxx', '*.c++',      #
+#   '*.C' for C++, '*.f', '*.F' for FORTRAN.                           #
+# - library files: 'lib*.a' for static libraries, 'lib*.so' for shared #
+#   libraries.                                                         #
+#                                                                      #
+# The compilers are specified by the 'C', 'CPP' and 'F' variables for  #
+# C, C++ and FORTRAN. The linker is specified by 'L' variable.         #
+#                                                                      #
+# If you want to pass some global switches while compilation and       #
+# linking, edit the variable 'ARGS'. If you want to change             #
+# optimization settings, edit the variable 'OPT'.                      #
+#                                                                      #
+# If you want to pass arguments to the linker, then edit 'LARGS'       #
+# variable.                                                            #
+# If you want to pass arguments to the C, C++ or FORTRAN compilers,    #
+# respectively, then edit 'CARGS', 'CPPARGS', 'FARGS'.                 #
+#                                                                      #
+# To compile in debug mode, say 'make debug' instead of 'make'.        #
+#                                                                      #
+########################################################################
+########################################################################
+
+
+####################### User definitions section: ######################
+
+
+### C compiler:
+ifeq "$(strip $(USEDC))" ""
+ USEDC         = gcc
+endif
+### C++ compiler:
+ifeq "$(strip $(USEDCPP))" ""
+ USEDCPP       = g++
+endif
+### FORTRAN compiler:
+ifeq "$(strip $(USEDF))" ""
+ USEDF         = g77
+endif
+### Linker:
+USEDL          = $(USEDCPP)
+
+### Global arguments:
+ARGS           = 
+### Global optimization arguments:
+OPT            = -O3
+#OPT           = -O0
+#OPT           = -O2
+#OPT           = -O3 -fno-fast-math -ffloat-store
+
+### Arguments to be passed when C compiling. E.g. -Idir switches:
+CARGS          = -Wall
+### Arguments to be passed when C++ compiling. E.g. -Idir switches:
+CPPARGS                = -Wall `pkg-config gridripperapp --cflags` `pkg-config gridripper --cflags` `blop-config --cflags`
+### Arguments to be passed when FORTRAN compiling. E.g. -Idir switches:
+FARGS          = -Wall
+### Arguments to be passed while linking. E.g. -Llibdir and -llib 
+### switches:
+LARGS          = -Wall `pkg-config gridripperapp --libs` `pkg-config gridripperapp --libs` -L/usr/lib64 `blop-config --libs`
+
+### Debug flags if debug mode is switched on:
+ifeq "$(strip $(DEBUG))" "TRUE"
+ CARGS += -g
+ CPPARGS += -g
+ FARGS += -g
+ LARGS += -g
+endif
+
+### Static flags if static mode is switched on:
+ifeq "$(strip $(STATIC))" "TRUE"
+ LARGS += -static
+endif
+
+### Other stuff:
+
+
+#################### End of user definitions section. ##################
+
+
+### Extract source file names etc. from directories:
+CHEADERS       = $(filter %.h,$(wildcard inc/*))
+CPPHEADERS     = $(filter %.h %.H,$(wildcard inc/*))
+FHEADERS       = $(filter %.inc,$(wildcard inc/*))
+CSOURCES       = $(filter %.c,$(wildcard src/*))
+CPPSOURCES     = $(filter %.cc %.cpp %.cxx %.c++ %.C,$(wildcard src/*))
+FSOURCES       = $(filter %.f %.F,$(wildcard src/*))
+CLIBSOURCES    = $(filter %.c,$(wildcard lib/*))
+CPPLIBSOURCES  = $(filter %.cc %.cpp %.cxx %.c++ %.C,$(wildcard lib/*))
+FLIBSOURCES    = $(filter %.f %.F,$(wildcard lib/*))
+OBJECTS                = $(patsubst %.c,%.o,$(CSOURCES)) $(patsubst %.cc,%.o,$(patsubst %.cpp,%.o,$(patsubst %.cxx,%.o,$(patsubst %.c++,%.o,$(patsubst %.C,%.o,$(CPPSOURCES)))))) $(patsubst %.f,%.o,$(patsubst %.F,%.o,$(FSOURCES)))
+LIBOBJECTS     = $(patsubst %.c,%.o,$(CLIBSOURCES)) $(patsubst %.cc,%.o,$(patsubst %.cpp,%.o,$(patsubst %.cxx,%.o,$(patsubst %.c++,%.o,$(patsubst %.C,%.o,$(CPPLIBSOURCES)))))) $(patsubst %.f,%.o,$(patsubst %.F,%.o,$(FLIBSOURCES)))
+LIBS           = $(filter lib/lib%.a,$(wildcard lib/*))
+BINARIES       = $(patsubst src/%.o,bin/%,$(OBJECTS))
+CINCFLAGS      = -Iinc
+CPPINCFLAGS    = -Iinc
+FINCFLAGS      = -Iinc
+LIBFLAGS       = -Llib $(patsubst lib/lib%.a,-l%,$(LIBS))
+
+### These entries do not correspond to file names:
+.PHONY         : all debug clean Clean
+
+### These are the secondary files:
+.SECONDARY     : $(OBJECTS) $(LIBOBJECTS)
+
+### Make everything. This is the default:
+all            : $(BINARIES)
+
+### Make everything, with debug mode:
+debug          :
+       $(MAKE) DEBUG=TRUE
+
+### Make everything, dynamic linking:
+dynamic                :
+       rm -f $(BINARIES)
+       $(MAKE) STATIC=FALSE
+
+### Make everything, static linking:
+static         :
+       rm -f $(BINARIES)
+       $(MAKE) STATIC=TRUE
+
+### Link objects to create binaries:
+bin/%          : src/%.o $(LIBOBJECTS) $(LIBS) Makefile
+       $(USEDL) $(ARGS) $(OPT) -o $@ $< $(LIBOBJECTS) $(LIBFLAGS) $(LARGS)
+
+### Compile library sources:
+lib/%.o                : lib/%.c $(CHEADERS) Makefile
+       $(USEDC) $(ARGS) $(OPT) -c -o $@ $< $(CINCFLAGS) $(CARGS)
+lib/%.o                : lib/%.cc $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+lib/%.o                : lib/%.cpp $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+lib/%.o                : lib/%.cxx $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+lib/%.o                : lib/%.c++ $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+lib/%.o                : lib/%.C $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+lib/%.o                : lib/%.f $(FHEADERS) Makefile
+       $(USEDF) $(ARGS) $(OPT) -c -o $@ $< $(FINCFLAGS) $(FARGS)
+lib/%.o                : lib/%.F $(FHEADERS) Makefile
+       $(USEDF) $(ARGS) $(OPT) -c -o $@ $< $(FINCFLAGS) $(FARGS)
+
+### Compile program sources:
+src/%.o                : src/%.c $(CHEADERS) Makefile
+       $(USEDC) $(ARGS) $(OPT) -c -o $@ $< $(CINCFLAGS) $(CARGS)
+src/%.o                : src/%.cc $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+src/%.o                : src/%.cpp $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+src/%.o                : src/%.cxx $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+src/%.o                : src/%.c++ $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+src/%.o                : src/%.C $(CPPHEADERS) Makefile
+       $(USEDCPP) $(ARGS) $(OPT) -c -o $@ $< $(CPPINCFLAGS) $(CPPARGS)
+src/%.o                : src/%.f $(FHEADERS) Makefile
+       $(USEDF) $(ARGS) $(OPT) -c -o $@ $< $(FINCFLAGS) $(FARGS)
+src/%.o                : src/%.F $(FHEADERS) Makefile
+       $(USEDF) $(ARGS) $(OPT) -c -o $@ $< $(FINCFLAGS) $(FARGS)
+
+### Clean up sources:
+clean          :
+       rm -f $(OBJECTS) $(BINARIES)
+Clean          :
+       rm -f $(OBJECTS) $(BINARIES) $(LIBOBJECTS)
diff --git a/inputs/example/kgminkowski/analysis/bin/arg.list b/inputs/example/kgminkowski/analysis/bin/arg.list
new file mode 100644 (file)
index 0000000..65a33a7
--- /dev/null
@@ -0,0 +1,4 @@
+#################### Config file to determine program sequence #########
+
+### Run intq:
+intq inp/intq.conf
diff --git a/inputs/example/kgminkowski/analysis/bin/go.sh b/inputs/example/kgminkowski/analysis/bin/go.sh
new file mode 100755 (executable)
index 0000000..93d867f
--- /dev/null
@@ -0,0 +1,133 @@
+#!/bin/bash
+########################################################################
+########################################################################
+#                                                                      #
+# This is an automatic 'go.sh' script. One generally does not have to  #
+# edit this file. Just put the program sources (source files with the  #
+# program body) into the 'src' directory, put the library sources      #
+# (source files to be linked to the programs, or the compiled          #
+# libraries or shared libraries themselves) into the 'lib'             #
+# directory, put the header sources (source files to be included in)   #
+# into the 'inc' directory. You can directly refer to include files in #
+# 'inc' directory, as they are added to the include search path.       #
+#                                                                      #
+# By running this 'go.sh' shell script in 'bin' directory (from        #
+# anywhere), everything is done. This will compile source (by invoking #
+# automatic 'Makefile' in the previous directory), remove old outputs, #
+# and run the program binaries.                                        #
+#                                                                      #
+# If you want to pass arguments to the program binaries, edit the      #
+# 'PARAMS' variable. Then every binary will be run once and only once  #
+# in alphabetic order, with these parameters. If you want to run the   #
+# programs in a more complex way, then create a file in directory      #
+# 'bin' named 'arg.list', containing lines in the following format:    #
+# binary_name_0 some_arguments_0                                       #
+# binary_name_1 some_arguments_1                                       #
+# ...                                                                  #
+# The above scheme should be self-explanatory. Empty lines or lines    #
+# beginning with '#' are skipped. Each argument line is appended with  #
+# the content of 'PARAMS' environmental variable.                      #
+#                                                                      #
+# Refer to 'Makefile' for compilation issues.                          #
+#                                                                      #
+########################################################################
+########################################################################
+
+
+####################### User definitions section: ######################
+
+
+### Pass these parameters to the programs:
+PARAMS=
+
+### Output directories to be cleaned up before running:
+OUTDIRS="out/bplots out/data out/hist out/plots"
+
+
+#################### End of user definitions section. ##################
+
+### Get directory, where the programs are located:
+DIR=$PWD
+DIROLD=$DIR
+ROUTE=$0
+NF=$((`echo -n $ROUTE | tr -c -d '/' | wc -c`))
+if [[ $NF > 0 ]] ; then 
+    ROUTE=`echo -n $ROUTE | cut --delimiter='/' --fields=1-$NF`
+else
+    ROUTE=""
+fi
+DIR=$DIR/$ROUTE/..
+
+### Change to directory, where the programs are located:
+cd $DIR
+
+### Get binary names:
+BINARIESORIG=`cd src ; ls *.c *.cc *.cpp *.cxx *.c++ *.C *.f *.F 2>/dev/null ; cd ..`
+BINARIES=""
+for B in $BINARIESORIG ; do
+    NF=$((`echo -n $B | tr -c -d '.' | wc -c`))
+    BINARIES="$BINARIES `echo -n $B | cut --delimiter='.' --fields=1-$NF`"
+done
+
+### Extract programs from file 'arg.list' in directory 'bin', if exists:
+COMPBINARIES=""
+NP=0
+if [[ -e bin/arg.list ]] ; then
+    for LINE in `cat bin/arg.list | tr ' ' '^' | tr '\t' '^'` ; do
+       LINE=`echo -n "$LINE" | tr '^' ' '`
+       NW=0
+       for WORD in $LINE ; do
+           if (( $NW == 0 )) && [[ "`echo -n $WORD | cut -c 1`" == "#" ]] ; then break ; fi
+           if (( $NW == 0 )) && ! ( echo -n "$BINARIES" | grep "$WORD" >&/dev/null ) ; then
+               echo "Cannot find $WORD among binary names! Nothing done. Exiting."
+               cd $DIROLD
+               exit 1
+           fi
+           if (( $NW == 0 )) ; then
+               PROGRAM[$NP]=$WORD
+               if ! ( echo -n "$COMPBINARIES" | grep "$WORD" >&/dev/null ) ; then COMPBINARIES="$COMPBINARIES $WORD" ; fi
+           else
+               PROGRAM[$NP]="${PROGRAM[$NP]} $WORD"
+           fi
+           NW=$(($NW+1))
+       done
+       if (( $NW == 0 )) ; then continue ; fi
+       PROGRAM[$NP]="${PROGRAM[$NP]} $PARAMS"
+       NP=$(($NP+1))
+    done
+fi
+
+### As a fallback take all the binaries:
+if (( $NP == 0 )) ; then
+    for WORD in $BINARIES ; do
+       PROGRAM[$NP]="$WORD $PARAMS"
+       NP=$(($NP+1))
+    done
+    COMPBINARIES=$BINARIES
+fi
+
+### Compile and link sources to create binaries:
+for B in $COMPBINARIES ; do
+    echo "Compliling $B ..."
+    make -s bin/$B
+    if (( $? != 0 )) ; then echo -e "go.sh:\tCompilation of $B failed!" ; cd $DIROLD ; exit 1 ; fi
+done
+
+### Remove old outputs:
+echo "Removing old outputs ..."
+for F in $OUTDIRS ; do
+    mkdir -p $F
+    rm -rf $F/*
+done
+
+### Run programs:
+N=0
+while (( $N < $NP )) ; do
+    echo "Starting ${PROGRAM[$N]} ..."
+    nice -+19 ./bin/${PROGRAM[$N]}
+    if (( $? != 0 )) ; then echo -e "go.sh:\tRunning of ${PROGRAM[$N]} failed!" ; cd $DIROLD ; exit 1 ; fi
+    N=$(($N+1))
+done
+
+### Change to original directory:
+cd $DIROLD
diff --git a/inputs/example/kgminkowski/analysis/doc/README b/inputs/example/kgminkowski/analysis/doc/README
new file mode 100644 (file)
index 0000000..f55361e
--- /dev/null
@@ -0,0 +1,28 @@
+
+This is a program package to analyse the BData files generated by 
+GridRipper for the KerrHiggs PDE. Run the program without arguments 
+to get description of program arguments.
+
+Compilation and running is described below.
+
+This is a framework for simple jobs ran on simple computers. It contains 
+an automated 'Makefile' and and automated starter shell script 'bin/go.sh'. 
+
+The 'Makefile' is automated in the sense that one only have to specify 
+libraries, switches etc. in the 'User definitions section'. Then, the sources 
+and library sources are autodetected, compiled and linked. It is capable 
+to compile and link C, C++ and FORTRAN sources and libraries.
+
+The starter shell script 'bin/go.sh' autodetects the sources, 
+compiles and links them with the above mentioned 'Makefile', and runs 
+them with appropriate arguments. There are two ways of operating it.
+1) Passing parameters to the variable 'PARAMS' in the 
+'User definitions section' of 'bin/go.sh'. In this case each binary 
+will be ran in alphabetic order after each other, with the parameters 
+'$PARAMS'.
+2) Prepare a file called 'bin/arg.list'. Each of its line can be of the 
+format:
+binary_name arguments
+In this way one can run arbitrary program binaries with arbitrary arguments, 
+in arbitrary order. The string '$PARAMS' will be appended to each argument 
+list. Empty lines and lines beginning with '#' are skipped in 'arg.list'.
diff --git a/inputs/example/kgminkowski/analysis/inc/config.h b/inputs/example/kgminkowski/analysis/inc/config.h
new file mode 100644 (file)
index 0000000..0cd4061
--- /dev/null
@@ -0,0 +1,344 @@
+/**
+ * config.h  Declares tools for reading config files, consisting of
+ *           lines of the format:
+ *           token      value
+ *           pairs. Empty lines, lines containing only delimiters
+ *           (tabs and spaces), and lines beginning with '#' are
+ *           ignored.
+ *           A configfile is read by:
+ *           // Declare config variable:
+ *           config conf;
+ *           // Read the content of a configfile into conf:
+ *           conf.append("configfile");
+ *           // Get the value associated to "token", interpreted
+ *           // explicitely as of type T:
+ *           T t=getconfig<T>(conf, "token");
+ *           // (In case of failure, automatic warning is displayed.)
+ *           // One may omit the template argument, to interpret the
+ *           // value automatically (implicitly) of the appropriate
+ *           // type:
+ *           T t=getconfig(conf, "token");
+ *           // Clear the content of the config variable:
+ *           conf.clear();
+ *           There are also useful system resource usage functions.
+ */
+
+
+#ifndef __CONFIG_H
+#define __CONFIG_H
+
+
+// Uncomment this if you don not want to launch memory watcher 
+// satellite program:
+#define __NO_MEMORY_WATCHER
+
+
+// Uncomment this if you don not want automatic logging 
+// when program finishes:
+//#define __NO_AUTO_LOGGING
+
+
+#include <exception>
+#include <string>
+#include <map>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <ctime>
+#include <sys/types.h>
+#include <sys/resource.h>
+#include <unistd.h>
+#include <pstream.h>
+
+
+// Forward declaration of class config_entry. This stores a config file 
+// entry.
+class config_entry;
+
+
+// Forward declaration of class config. This stores a config file 
+// content (an array of config_entry variables).
+class config;
+
+
+// Forward declaration of template function getconfig. This extracts an 
+// entry from a config variable, indexed by a token, explicitly 
+// interpreted as of a given specified type.
+template <typename T>
+T getconfig(const config&, const std::string&);
+// The same, but reads from file, specified by an input stream.
+template <typename T>
+T getconfig(std::istream&, const std::string&);
+config_entry getconfig(std::istream&, const std::string&);
+// The same, but reads from file, specified by a file name.
+template <typename T>
+T getconfig(const std::string&, const std::string&);
+
+// Forward declaration of function getconfig. This extracts an entry 
+// from a config variable, indexed by a token, but does not interpret 
+// it in a forced way (automatic conversion will interpret it as some 
+// type later, implicitly).
+extern config_entry getconfig(const config&, const std::string&);
+// The same, but reads from file, specified by an input stream:
+extern config_entry getconfig(std::istream&, const std::string&);
+// The same, but reads from file, specified by a file name:
+extern config_entry getconfig(const std::string&, const std::string&);
+
+
+// Forward declaration of class config_exception. This is an exception 
+// thrown when problem occurs during parsing of a configuration content.
+class config_exception;
+
+
+/**
+ * Stores a config file entry, which can either contain a string
+ * or a number (long double, int, unsigned int, bool), or a char.
+ */
+class config_entry
+{
+    protected:
+       /// Original string value:
+       std::string string_;
+       /// string_ interpreted as long double:
+       long double longdouble_;
+       /// string -> long double conversion successfullness flag:
+       bool success_;
+    public:
+       /// Default constructor:
+       config_entry();
+       /// Copy constructor:
+       config_entry(const config_entry&);
+       /// Reinitialization from an std::string:
+       config_entry& reinit(const std::string&);
+       /// Constructor from an std::string:
+       config_entry(const std::string&);
+       /// Destructor:
+       ~config_entry();
+       /// Implicit conversion to std::string:
+       operator std::string() const;
+       /// Implicit conversion to float:
+       operator float() const;
+       /// Implicit conversion to double:
+       operator double() const;
+       /// Implicit conversion to long double:
+       operator long double() const;
+       /// Implicit conversion to int:
+       operator int() const;
+       /// Implicit conversion to unsigned int:
+       operator unsigned int() const;
+       /// Implicit conversion to bool:
+       operator bool() const;
+       /// Implicit conversion to char:
+       operator char() const;
+       /// Return the string data field:
+       std::string string() const;
+       /// Return the long double data field:
+       long double ldbl() const;
+       /// Return a flag, which tells if initial string -> long double conversion was successful:
+       bool success() const;
+       /// Assignment:
+       const config_entry& operator=(const config_entry&);
+       /// Friend function getconfig:
+       friend config_entry getconfig(const config&, const std::string&);
+};
+
+
+/**
+ * Declaration of a config file content container (class config).
+ */
+class config : protected std::map< std::string, config_entry >
+{
+    public:
+       /// Default constructor:
+       config();
+       /// Copy constructor:
+       config(const config&);
+       /// Destructor:
+       ~config();
+       /// Append the content of a config file:
+       config& append(std::istream&);
+       /// Append the content of a config file, specified by the name:
+       config& append(const std::string&);
+       /// Constructor with the content of a config file:
+       config(std::istream&);
+       /// Constructor with the content of a config file, specified by the name:
+       config(const std::string&);
+       /// Clear content:
+       config& clear();
+       /// Friend function getconfig:
+       friend config_entry getconfig(const config&, const std::string&);
+};
+
+
+/**
+ * Implementation of getconfig<> template functions.
+ */
+template <typename T>
+T getconfig(const config &conf, const std::string &token)
+{
+    return (T)getconfig(conf, token);
+}
+
+template <typename T>
+T getconfig(std::istream &in, const std::string &token)
+{
+    return (T)getconfig(in, token);
+}
+
+template <typename T>
+T getconfig(const std::string &filename, const std::string &token)
+{
+    return (T)getconfig(filename, token);
+}
+
+// See the implementation of getconfig functions without template 
+// arguments in config.cc.
+
+
+/**
+ * File numbering function. Extends the front digits with appropriate 
+ * number of zeros. (Arguments: index, maximal index):
+ */
+extern std::string filenumbering(const unsigned int, const unsigned int);
+
+
+/**
+ * Get the revision string.
+ */
+extern std::string getrevision();
+
+
+/**
+ * Get the elapsed time in seconds since the Epoch 
+ * (00:00:00.00 UTC, 1 Jauary 1970):
+ */
+extern unsigned int gettime();
+
+
+/**
+ * Get the actual time and date in a human-readable string format:
+ */
+extern std::string gettime_hr();
+
+
+/**
+ * Get the user time of the in seconds, used up so far by the program:
+ */
+extern unsigned int getusertime();
+
+
+/**
+ * Get the actual CPU time in seconds, used up so far by the program:
+ */
+extern unsigned int getcputime();
+
+
+/**
+ * Get currently used memory in MegaBytes:
+ */
+extern unsigned int getmem();
+
+
+/**
+ * Write logging data into file:
+ */
+extern void putlogdata(const std::string&);
+
+
+/**
+ * Exception class for config exception handling.
+ */
+class config_exception : public std::exception
+{
+    protected:
+       const std::string function_name_;
+       const std::string message_;
+    public:
+       config_exception(const std::string& function_name = std::string(), const std::string& message = std::string());
+       const std::string& function_name() const;
+       const std::string& message() const;
+       virtual const char* what() const throw();
+       virtual ~config_exception() throw();
+       static void print_error(const std::string& function_name = std::string(), const std::string& message = std::string());
+};
+
+
+/**
+ * This class stores program summary information (starting time etc.):
+ */
+class summaryinfo
+{
+    // Private, not protected, so user cannot use these, 
+    // even with inherited class:
+    private:
+       /// Starting time in seconds from the Epoch (00:00:00.00 UTC, 1 Jauary 1970):
+       static const unsigned int starttime_;
+       /// Starting time in human-readable format:
+       static const std::string starttime_hr_;
+       /// Current maximal used memory in MegaBytes:
+       static unsigned int maxmemory_;
+#ifndef __NO_MEMORY_WATCHER
+       /// An auxiliary stream for communicating with the below stream:
+       static std::ipstream& auxmaxmemstream_;
+       /// Maximal memory reading stream:
+       static std::ipstream& maxmemstream_;
+#endif
+       /// Logfile stream:
+       static std::ostream &logfile_;
+       /// Logfile name:
+       static std::string logfilename_;
+       /// Determine if this is the first copy of this class:
+       const bool firstcopy_;
+       static bool firstcopydeclared_;
+       /// Determine if already written:
+       static bool iswritten_;
+       /// Write/not write logfile:
+       static bool writelogfile_;
+    public:
+       /// Default constructor:
+       summaryinfo();
+    private:
+       /// Copy constructor (so that user cannot call it):
+       summaryinfo(const summaryinfo&);
+    public:
+       /// Destructor (writes summary logfile, if requested):
+       ~summaryinfo();
+    private:
+       /// Assignment operator (so that user cannot call it):
+       const summaryinfo& operator=(const summaryinfo&);
+    public:
+       /// Change the logfile stream (cerr by default):
+       static void logfile(const std::ostream&);
+       /// Get the logfile stream:
+       static std::ostream& logfile();
+       /// Change the name of logfile ("", directed to cerr by default):
+       static void logfilename(const std::string&);
+       /// Get the name of logfile:
+       static std::string logfilename();
+       /// Require/not require logfile writing (require by default):
+       static void write_logfile(const bool);
+       /// Determine if logfile is required:
+       static bool write_logfile();
+       /// Get start time:
+       static unsigned int getstarttime();
+       /// Get start time in human-readable format:
+       static std::string getstarttime_hr();
+       /// Acquire maximal memory:
+       static void acquiremaxmem();
+       /// Show used maximal memory:
+       static unsigned int showmaxmem();
+};
+
+
+#ifndef __NO_AUTO_LOGGING
+/**
+ * Summary info variable. The constructor of this will record the 
+ * starting time of the program. The destructor of this will write a 
+ * summary logfile.
+ */
+extern const summaryinfo __summaryinfo;
+#endif
+
+
+#endif /* __CONFIG_H */
diff --git a/inputs/example/kgminkowski/analysis/inc/pstream.h b/inputs/example/kgminkowski/analysis/inc/pstream.h
new file mode 100644 (file)
index 0000000..933858a
--- /dev/null
@@ -0,0 +1,2122 @@
+/* $Id: pstream.h,v 1.90 2005/06/11 09:25:06 redi Exp $
+PStreams - POSIX Process I/O for C++
+Copyright (C) 2001,2002,2003,2004 Jonathan Wakely
+
+This file is part of PStreams.
+
+PStreams is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as
+published by the Free Software Foundation; either version 2.1 of
+the License, or (at your option) any later version.
+
+PStreams 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with PStreams; if not, write to the Free Software Foundation, Inc.,
+59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+/**
+ * @file pstream.h
+ * @brief Declares all PStreams classes.
+ * @author Jonathan Wakely
+ *
+ * Defines classes redi::ipstream, redi::opstream, redi::pstream
+ * and redi::rpstream.
+ */
+
+///// I changed 'redi' to 'std' !!! /////
+///// I added '#include <string>' !!! /////
+///// I added '#include <sstream>' !!! /////
+///// I added 'isinpipe' and 'isoutpipe' functions!!! /////
+///// I added 'openin' and 'openout' functions!!! /////
+///// I added 'pstream_common::kill' function!!! /////
+
+#ifndef REDI_PSTREAM_H_SEEN
+#define REDI_PSTREAM_H_SEEN
+
+#include <ios>
+#include <streambuf>
+#include <istream>
+#include <ostream>
+#include <string>
+#include <vector>
+#include <algorithm>    // for min()
+#include <cstring>      // for memcpy(), memmove() etc.
+#include <cerrno>       // for errno
+#include <cstddef>      // for size_t
+#include <cstdlib>      // for exit()
+#include <sys/types.h>  // for pid_t
+#include <sys/wait.h>   // for waitpid()
+#include <sys/ioctl.h>  // for ioctl() and FIONREAD
+#if defined(__sun)
+# include <sys/filio.h> // for FIONREAD on Solaris 2.5
+#endif
+#include <unistd.h>     // for pipe() fork() exec() and filedes functions
+#include <signal.h>     // for kill()
+#include <fcntl.h>      // for fcntl()
+#if REDI_EVISCERATE_PSTREAMS
+# include <stdio.h>     // for FILE, fdopen()
+#endif
+#include <string>
+#include <sstream>
+
+
+/// The library version.
+#define PSTREAMS_VERSION 0x0052   // 0.5.2
+
+/**
+ *  @namespace redi
+ *  @brief  All PStreams classes are declared in namespace redi.
+ *
+ *  Like the standard IOStreams, PStreams is a set of class templates,
+ *  taking a character type and traits type. As with the standard streams
+ *  they are most likely to be used with @c char and the default
+ *  traits type, so typedefs for this most common case are provided.
+ *
+ *  The @c pstream_common class template is not intended to be used directly,
+ *  it is used internally to provide the common functionality for the
+ *  other stream classes.
+ */
+//namespace redi
+namespace std
+{
+  /// Common base class providing constants and typenames.
+  struct pstreams
+  {
+    /// Type used to specify how to connect to the process.
+    typedef std::ios_base::openmode           pmode;
+
+    /// Type used to hold the arguments for a command.
+    typedef std::vector<std::string>          argv_type;
+
+    /// Type used for file descriptors.
+    typedef int                               fd_type;
+
+    static const pmode pstdin  = std::ios_base::out; ///< Write to stdin
+    static const pmode pstdout = std::ios_base::in;  ///< Read from stdout
+    static const pmode pstderr = std::ios_base::app; ///< Read from stderr
+
+  protected:
+    enum { bufsz = 32 };  ///< Size of pstreambuf buffers.
+    enum { pbsz  = 2 };   ///< Number of putback characters kept.
+  };
+
+  /// Class template for stream buffer.
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class basic_pstreambuf
+    : public std::basic_streambuf<CharT, Traits>
+    , public pstreams
+    {
+    public:
+      // Type definitions for dependent types
+      typedef CharT                             char_type;
+      typedef Traits                            traits_type;
+      typedef typename traits_type::int_type    int_type;
+      typedef typename traits_type::off_type    off_type;
+      typedef typename traits_type::pos_type    pos_type;
+      /** @deprecated use pstreams::fd_type instead. */
+      typedef fd_type                           fd_t;
+
+      /// Default constructor.
+      basic_pstreambuf();
+
+      /// Constructor that initialises the buffer with @a command.
+      basic_pstreambuf(const std::string& command, pmode mode);
+
+      /// Constructor that initialises the buffer with @a file and @a argv.
+      basic_pstreambuf( const std::string& file,
+                        const argv_type& argv,
+                        pmode mode );
+
+      /// Destructor.
+      ~basic_pstreambuf();
+
+      /// Initialise the stream buffer with @a command.
+      basic_pstreambuf*
+      open(const std::string& command, pmode mode);
+
+      /// Initialise the stream buffer with @a file and @a argv.
+      basic_pstreambuf*
+      open(const std::string& file, const argv_type& argv, pmode mode);
+
+      /// Close the stream buffer and wait for the process to exit.
+      basic_pstreambuf*
+      close();
+
+      /// Send a signal to the process.
+      basic_pstreambuf*
+      kill(int signal = SIGTERM);
+
+      /// Close the pipe connected to the process' stdin.
+      void
+      peof();
+
+      /// Change active input source.
+      bool
+      read_err(bool readerr = true);
+
+      /// Report whether the stream buffer has been initialised.
+      bool
+      is_open() const;
+
+      /// Report whether the process has exited.
+      bool
+      exited();
+
+#if REDI_EVISCERATE_PSTREAMS
+      /// Obtain FILE pointers for each of the process' standard streams.
+      std::size_t
+      fopen(std::FILE*& in, std::FILE*& out, std::FILE*& err);
+#endif
+
+      /// Return the exit status of the process.
+      int
+      status() const;
+
+      /// Return the error number for the most recent failed operation.
+      int
+      error() const;
+
+    protected:
+      /// Transfer characters to the pipe when character buffer overflows.
+      int_type
+      overflow(int_type c);
+
+      /// Transfer characters from the pipe when the character buffer is empty.
+      int_type
+      underflow();
+
+      /// Make a character available to be returned by the next extraction.
+      int_type
+      pbackfail(int_type c = traits_type::eof());
+
+      /// Write any buffered characters to the stream.
+      int
+      sync();
+
+      /// Insert multiple characters into the pipe.
+      std::streamsize
+      xsputn(const char_type* s, std::streamsize n);
+
+      /// Insert a sequence of characters into the pipe.
+      std::streamsize
+      write(const char_type* s, std::streamsize n);
+
+      /// Extract a sequence of characters from the pipe.
+      std::streamsize
+      read(char_type* s, std::streamsize n);
+
+      /// Report how many characters can be read from active input without blocking.
+      std::streamsize
+      showmanyc();
+
+    protected:
+      /// Enumerated type to indicate whether stdout or stderr is to be read.
+      enum buf_read_src { rsrc_out = 0, rsrc_err = 1 };
+
+      /// Initialise pipes and fork process.
+      pid_t
+      fork(pmode mode);
+
+      /// Wait for the child process to exit.
+      int
+      wait(bool nohang = false);
+
+      /// Return the file descriptor for the output pipe.
+      fd_type&
+      wpipe();
+
+      /// Return the file descriptor for the active input pipe.
+      fd_type&
+      rpipe();
+
+      /// Return the file descriptor for the specified input pipe.
+      fd_type&
+      rpipe(buf_read_src which);
+
+      void
+      create_buffers(pmode mode);
+
+      void
+      destroy_buffers(pmode mode);
+
+      /// Writes buffered characters to the process' stdin pipe.
+      bool
+      empty_buffer();
+
+      bool
+      fill_buffer();
+
+      /// Return the active input buffer.
+      char_type*
+      rbuffer();
+
+      buf_read_src
+      switch_read_buffer(buf_read_src);
+
+    private:
+      basic_pstreambuf(const basic_pstreambuf&);
+      basic_pstreambuf& operator=(const basic_pstreambuf&);
+
+      void
+      init_rbuffers();
+
+      pid_t         ppid_;        // pid of process
+      fd_type       wpipe_;       // pipe used to write to process' stdin
+      fd_type       rpipe_[2];    // two pipes to read from, stdout and stderr
+      char_type*    wbuffer_;
+      char_type*    rbuffer_[2];
+      char_type*    rbufstate_[3];
+      /// Index into rpipe_[] to indicate active source for read operations.
+      buf_read_src  rsrc_;
+      int           status_;      // hold exit status of child process
+      int           error_;       // hold errno if fork() or exec() fails
+    };
+
+  /// Class template for common base class.
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class pstream_common
+    : virtual public std::basic_ios<CharT, Traits>
+    , virtual public pstreams
+    {
+    protected:
+      typedef basic_pstreambuf<CharT, Traits>       streambuf_type;
+
+      /// Default constructor.
+      pstream_common();
+
+      /// Constructor that initialises the stream by starting a process.
+      pstream_common(const std::string& command, pmode mode);
+
+      /// Constructor that initialises the stream by starting a process.
+      pstream_common(const std::string& file, const argv_type& argv, pmode mode);
+
+      /// Pure virtual destructor.
+      virtual
+      ~pstream_common() = 0;
+
+      /// Start a process.
+      void
+      do_open(const std::string& command, pmode mode);
+
+      /// Start a process.
+      void
+      do_open(const std::string& file, const argv_type& argv, pmode mode);
+
+    public:
+      /// Close the pipe.
+      void
+      close();
+
+      /// Kill the pipe.
+      void
+      kill(int signal=SIGTERM);
+
+      /// Report whether the stream's buffer has been initialised.
+      bool
+      is_open() const;
+
+      /// Return the command used to initialise the stream.
+      const std::string&
+      command() const;
+
+      /// Return a pointer to the stream buffer.
+      streambuf_type*
+      rdbuf() const;
+
+#if REDI_EVISCERATE_PSTREAMS
+      /// Obtain FILE pointers for each of the process' standard streams.
+      std::size_t
+      fopen(std::FILE*& in, std::FILE*& out, std::FILE*& err);
+#endif
+
+    protected:
+      std::string       command_; ///< The command used to start the process.
+      streambuf_type    buf_;     ///< The stream buffer.
+    };
+
+
+  /**
+   * @class basic_ipstream
+   * @brief Class template for Input PStreams.
+   *
+   * Reading from an ipstream reads the command's standard output and/or
+   * standard error (depending on how the ipstream is opened)
+   * and the command's standard input is the same as that of the process
+   * that created the object, unless altered by the command itself.
+   */
+
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class basic_ipstream
+    : public std::basic_istream<CharT, Traits>
+    , public pstream_common<CharT, Traits>
+    , virtual public pstreams
+    {
+      typedef std::basic_istream<CharT, Traits>     istream_type;
+      typedef pstream_common<CharT, Traits>         pbase_type;
+
+      using pbase_type::buf_;  // declare name in this scope
+
+    public:
+      /// Type used to specify how to connect to the process.
+      typedef typename pbase_type::pmode            pmode;
+
+      /// Type used to hold the arguments for a command.
+      typedef typename pbase_type::argv_type        argv_type;
+
+      /// Default constructor, creates an uninitialised stream.
+      basic_ipstream()
+      : istream_type(NULL), pbase_type()
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      basic_ipstream(const std::string& command, pmode mode = pstdout)
+      : istream_type(NULL), pbase_type(command, mode|pstdout)
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      basic_ipstream( const std::string& file,
+                      const argv_type& argv,
+                      pmode mode = pstdout )
+      : istream_type(NULL), pbase_type(file, argv, mode|pstdout)
+      { }
+
+      /**
+       * @brief Destructor.
+       *
+       * Closes the stream and waits for the child to exit.
+       */
+      ~basic_ipstream()
+      { }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a %command , @a mode|pstdout ).
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      void
+      open(const std::string& command, pmode mode = pstdout)
+      {
+        this->do_open(command, mode|pstdout);
+      }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a file , @a argv , @a mode|pstdout ).
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      void
+      open( const std::string& file,
+            const argv_type& argv,
+            pmode mode = pstdout )
+      {
+        this->do_open(file, argv, mode|pstdout);
+      }
+
+      /**
+       * @brief Set streambuf to read from process' @c stdout.
+       * @return  @c *this
+       */
+      basic_ipstream&
+      out()
+      {
+        this->buf_.read_err(false);
+        return *this;
+      }
+
+      /**
+       * @brief Set streambuf to read from process' @c stderr.
+       * @return  @c *this
+       */
+      basic_ipstream&
+      err()
+      {
+        this->buf_.read_err(true);
+        return *this;
+      }
+    };
+
+
+  /**
+   * @class basic_opstream
+   * @brief Class template for Output PStreams.
+   *
+   * Writing to an open opstream writes to the standard input of the command;
+   * the command's standard output is the same as that of the process that
+   * created the pstream object, unless altered by the command itself.
+   */
+
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class basic_opstream
+    : public std::basic_ostream<CharT, Traits>
+    , public pstream_common<CharT, Traits>
+    , virtual public pstreams
+    {
+      typedef std::basic_ostream<CharT, Traits>     ostream_type;
+      typedef pstream_common<CharT, Traits>         pbase_type;
+
+      using pbase_type::buf_;  // declare name in this scope
+
+    public:
+      /// Type used to specify how to connect to the process.
+      typedef typename pbase_type::pmode            pmode;
+
+      /// Type used to hold the arguments for a command.
+      typedef typename pbase_type::argv_type        argv_type;
+
+      /// Default constructor, creates an uninitialised stream.
+      basic_opstream()
+      : ostream_type(NULL), pbase_type()
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      basic_opstream(const std::string& command, pmode mode = pstdin)
+      : ostream_type(NULL), pbase_type(command, mode|pstdin)
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      basic_opstream( const std::string& file,
+                      const argv_type& argv,
+                      pmode mode = pstdin )
+      : ostream_type(NULL), pbase_type(file, argv, mode|pstdin)
+      { }
+
+      /**
+       * @brief Destructor
+       *
+       * Closes the stream and waits for the child to exit.
+       */
+      ~basic_opstream() { }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a %command , @a mode|pstdin ).
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      void
+      open(const std::string& command, pmode mode = pstdin)
+      {
+        this->do_open(command, mode|pstdin);
+      }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a file , @a argv , @a mode|pstdin ).
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      void
+      open( const std::string& file,
+            const argv_type& argv,
+            pmode mode = pstdin)
+      {
+        this->do_open(file, argv, mode|pstdin);
+      }
+    };
+
+
+  /**
+   * @class basic_pstream
+   * @brief Class template for Bidirectional PStreams.
+   *
+   * Writing to a pstream opened with @c pmode @c pstdin writes to the
+   * standard input of the command.
+   * Reading from a pstream opened with @c pmode @c pstdout and/or @c pstderr
+   * reads the command's standard output and/or standard error.
+   * Any of the process' @c stdin, @c stdout or @c stderr that is not
+   * connected to the pstream (as specified by the @c pmode)
+   * will be the same as the process that created the pstream object,
+   * unless altered by the command itself.
+   */
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class basic_pstream
+    : public std::basic_iostream<CharT, Traits>
+    , public pstream_common<CharT, Traits>
+    , virtual public pstreams
+    {
+      typedef std::basic_iostream<CharT, Traits>    iostream_type;
+      typedef pstream_common<CharT, Traits>         pbase_type;
+
+      using pbase_type::buf_;  // declare name in this scope
+
+    public:
+      /// Type used to specify how to connect to the process.
+      typedef typename pbase_type::pmode            pmode;
+
+      /// Type used to hold the arguments for a command.
+      typedef typename pbase_type::argv_type        argv_type;
+
+      /// Default constructor, creates an uninitialised stream.
+      basic_pstream()
+      : iostream_type(NULL), pbase_type()
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      basic_pstream(const std::string& command, pmode mode = pstdout|pstdin)
+      : iostream_type(NULL), pbase_type(command, mode)
+      { }
+
+      /**
+       * @brief Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      basic_pstream( const std::string& file,
+                     const argv_type& argv,
+                     pmode mode = pstdout|pstdin )
+      : iostream_type(NULL), pbase_type(file, argv, mode)
+      { }
+
+      /**
+       * @brief Destructor
+       *
+       * Closes the stream and waits for the child to exit.
+       */
+      ~basic_pstream() { }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a %command , @a mode ).
+       *
+       * @param command  a string containing a shell command.
+       * @param mode     the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      void
+      open(const std::string& command, pmode mode = pstdout|pstdin)
+      {
+        this->do_open(command, mode);
+      }
+
+      /**
+       * @brief Start a process.
+       *
+       * Calls do_open( @a file , @a argv , @a mode ).
+       *
+       * @param file  a string containing the pathname of a program to execute.
+       * @param argv  a vector of argument strings passed to the new program.
+       * @param mode  the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      void
+      open( const std::string& file,
+            const argv_type& argv,
+            pmode mode = pstdout|pstdin )
+      {
+        this->do_open(file, argv, mode);
+      }
+
+      /**
+       * @brief Set streambuf to read from process' @c stdout.
+       * @return  @c *this
+       */
+      basic_pstream&
+      out()
+      {
+        this->buf_.read_err(false);
+        return *this;
+      }
+
+      /**
+       * @brief Set streambuf to read from process' @c stderr.
+       * @return  @c *this
+       */
+      basic_pstream&
+      err()
+      {
+        this->buf_.read_err(true);
+        return *this;
+      }
+    };
+
+
+  /**
+   * @class basic_rpstream
+   * @brief template for Restricted PStreams.
+   *
+   * Writing to an rpstream opened with @c pmode @c pstdin writes to the
+   * standard input of the command.
+   * It is not possible to read directly from an rpstream object, to use
+   * an rpstream as in istream you must call either basic_rpstream::out()
+   * or basic_rpstream::err(). This is to prevent accidental reads from
+   * the wrong input source. If the rpstream was not opened with @c pmode
+   * @c pstderr then the class cannot read the process' @c stderr, and
+   * basic_rpstream::err() will return an istream that reads from the
+   * process' @c stdout, and vice versa.
+   * Reading from an rpstream opened with @c pmode @c pstdout and/or
+   * @c pstderr reads the command's standard output and/or standard error.
+   * Any of the process' @c stdin, @c stdout or @c stderr that is not
+   * connected to the pstream (as specified by the @c pmode)
+   * will be the same as the process that created the pstream object,
+   * unless altered by the command itself.
+   */
+
+  template <typename CharT, typename Traits = std::char_traits<CharT> >
+    class basic_rpstream
+    : public std::basic_ostream<CharT, Traits>
+    , private std::basic_istream<CharT, Traits>
+    , private pstream_common<CharT, Traits>
+    , virtual public pstreams
+    {
+      typedef std::basic_ostream<CharT, Traits>     ostream_type;
+      typedef std::basic_istream<CharT, Traits>     istream_type;
+      typedef pstream_common<CharT, Traits>         pbase_type;
+
+      using pbase_type::buf_;  // declare name in this scope
+
+    public:
+      /// Type used to specify how to connect to the process.
+      typedef typename pbase_type::pmode            pmode;
+
+      /// Type used to hold the arguments for a command.
+      typedef typename pbase_type::argv_type        argv_type;
+
+      /// Default constructor, creates an uninitialised stream.
+      basic_rpstream()
+      : ostream_type(NULL), istream_type(NULL), pbase_type()
+      { }
+
+      /**
+       * @brief  Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param command a string containing a shell command.
+       * @param mode the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      basic_rpstream(const std::string& command, pmode mode = pstdout|pstdin)
+      : ostream_type(NULL) , istream_type(NULL) , pbase_type(command, mode)
+      { }
+
+      /**
+       * @brief  Constructor that initialises the stream by starting a process.
+       *
+       * Initialises the stream buffer by calling do_open() with the supplied
+       * arguments.
+       *
+       * @param file a string containing the pathname of a program to execute.
+       * @param argv a vector of argument strings passed to the new program.
+       * @param mode the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      basic_rpstream( const std::string& file,
+                      const argv_type& argv,
+                      pmode mode = pstdout|pstdin )
+      : ostream_type(NULL), istream_type(NULL), pbase_type(file, argv, mode)
+      { }
+
+      /// Destructor
+      ~basic_rpstream() { }
+
+      /**
+       * @brief  Start a process.
+       *
+       * Calls do_open( @a %command , @a mode ).
+       *
+       * @param command a string containing a shell command.
+       * @param mode the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, pmode)
+       */
+      void
+      open(const std::string& command, pmode mode = pstdout|pstdin)
+      {
+        this->do_open(command, mode);
+      }
+
+      /**
+       * @brief  Start a process.
+       *
+       * Calls do_open( @a file , @a argv , @a mode ).
+       *
+       * @param file a string containing the pathname of a program to execute.
+       * @param argv a vector of argument strings passed to the new program.
+       * @param mode the I/O mode to use when opening the pipe.
+       * @see   do_open(const std::string&, const argv_type&, pmode)
+       */
+      void
+      open( const std::string& file,
+            const argv_type& argv,
+            pmode mode = pstdout|pstdin )
+      {
+        this->do_open(file, argv, mode);
+      }
+
+      /**
+       * @brief  Obtain a reference to the istream that reads
+       *         the process' @c stdout.
+       * @return @c *this
+       */
+      istream_type&
+      out()
+      {
+        this->buf_.read_err(false);
+        return *this;
+      }
+
+      /**
+       * @brief  Obtain a reference to the istream that reads
+       *         the process' @c stderr.
+       * @return @c *this
+       */
+      istream_type&
+      err()
+      {
+        this->buf_.read_err(true);
+        return *this;
+      }
+    };
+
+
+  /// Type definition for common template specialisation.
+  typedef basic_pstreambuf<char> pstreambuf;
+  /// Type definition for common template specialisation.
+  typedef basic_ipstream<char> ipstream;
+  /// Type definition for common template specialisation.
+  typedef basic_opstream<char> opstream;
+  /// Type definition for common template specialisation.
+  typedef basic_pstream<char> pstream;
+  /// Type definition for common template specialisation.
+  typedef basic_rpstream<char> rpstream;
+
+
+  /**
+   * When inserted into an output pstream the manipulator calls
+   * basic_pstreambuf<C,T>::peof() to close the output pipe,
+   * causing the child process to receive the end-of-file indicator
+   * on subsequent reads from its @c stdin stream.
+   *
+   * @brief   Manipulator to close the pipe connected to the process' stdin.
+   * @param   s  An output PStream class.
+   * @return  The stream object the manipulator was invoked on.
+   * @see     basic_pstreambuf<C,T>::peof()
+   * @relates basic_opstream basic_pstream basic_rpstream
+   */
+  template <typename C, typename T>
+    inline std::basic_ostream<C,T>&
+    peof(std::basic_ostream<C,T>& s)
+    {
+      typedef basic_pstreambuf<C,T> pstreambuf;
+      if (pstreambuf* p = dynamic_cast<pstreambuf*>(s.rdbuf()))
+        p->peof();
+      return s;
+    }
+
+
+  /*
+   * member definitions for pstreambuf
+   */
+
+
+  /**
+   * @class basic_pstreambuf
+   * Provides underlying streambuf functionality for the PStreams classes.
+   */
+
+  /** Creates an uninitialised stream buffer. */
+  template <typename C, typename T>
+    inline
+    basic_pstreambuf<C,T>::basic_pstreambuf()
+    : ppid_(-1)   // initialise to -1 to indicate no process run yet.
+    , wpipe_(-1)
+    , wbuffer_(NULL)
+    , rsrc_(rsrc_out)
+    , status_(-1)
+    , error_(0)
+    {
+      init_rbuffers();
+    }
+
+  /**
+   * Initialises the stream buffer by calling open() with the supplied
+   * arguments.
+   *
+   * @param command a string containing a shell command.
+   * @param mode    the I/O mode to use when opening the pipe.
+   * @see   open()
+   */
+  template <typename C, typename T>
+    inline
+    basic_pstreambuf<C,T>::basic_pstreambuf(const std::string& command, pmode mode)
+    : ppid_(-1)   // initialise to -1 to indicate no process run yet.
+    , wpipe_(-1)
+    , wbuffer_(NULL)
+    , rsrc_(rsrc_out)
+    , status_(-1)
+    , error_(0)
+    {
+      init_rbuffers();
+      open(command, mode);
+    }
+
+  /**
+   * Initialises the stream buffer by calling open() with the supplied
+   * arguments.
+   *
+   * @param file  a string containing the name of a program to execute.
+   * @param argv  a vector of argument strings passsed to the new program.
+   * @param mode  the I/O mode to use when opening the pipe.
+   * @see   open()
+   */
+  template <typename C, typename T>
+    inline
+    basic_pstreambuf<C,T>::basic_pstreambuf( const std::string& file,
+                                             const argv_type& argv,
+                                             pmode mode )
+    : ppid_(-1)   // initialise to -1 to indicate no process run yet.
+    , wpipe_(-1)
+    , wbuffer_(NULL)
+    , rsrc_(rsrc_out)
+    , status_(-1)
+    , error_(0)
+    {
+      init_rbuffers();
+      open(file, argv, mode);
+    }
+
+  /**
+   * Closes the stream by calling close().
+   * @see close()
+   */
+  template <typename C, typename T>
+    inline
+    basic_pstreambuf<C,T>::~basic_pstreambuf()
+    {
+      close();
+    }
+
+  /**
+   * Starts a new process by passing @a command to the shell
+   * and opens pipes to the process with the specified @a mode.
+   *
+   * Will duplicate the actions of  the  shell  in searching for an
+   * executable file if the specified file name does not contain a slash (/)
+   * character.
+   *
+   * There is no way to tell whether the shell command succeeded, this
+   * function will always succeed unless resource limits (such as
+   * memory usage, or number of processes or open files) are exceeded.
+   * This means is_open() will return true even if @a command cannot
+   * be executed.
+   *
+   * @param   command  a string containing a shell command.
+   * @param   mode     a bitwise OR of one or more of @c out, @c in, @c err.
+   * @return  NULL if the shell could not be started or the
+   *          pipes could not be opened, @c this otherwise.
+   * @see     <b>execlp</b>(3)
+   */
+  template <typename C, typename T>
+    basic_pstreambuf<C,T>*
+    basic_pstreambuf<C,T>::open(const std::string& command, pmode mode)
+    {
+#if 0
+      const std::string argv[] = { "sh", "-c", command };
+      return this->open("sh", std::vector<std::string>(argv, argv+3), mode);
+#else
+      basic_pstreambuf<C,T>* ret = NULL;
+
+      if (!is_open())
+      {
+        switch(fork(mode))
+        {
+        case 0 :
+          // this is the new process, exec command
+          ::execlp("sh", "sh", "-c", command.c_str(), (void*)NULL);
+
+          // can only reach this point if exec() failed
+
+          // parent can get exit code from waitpid()
+          ::_exit(errno);
+          // using std::exit() would make static dtors run twice
+
+        case -1 :
+          // couldn't fork, error already handled in pstreambuf::fork()
+          break;
+
+        default :
+          // this is the parent process
+          // activate buffers
+          create_buffers(mode);
+          ret = this;
+        }
+      }
+      return ret;
+#endif
+    }
+
+  /**
+   * @brief  Helper function to close a file descriptor.
+   *
+   * Inspects @a fd and calls <b>close</b>(3) if it has a non-negative value.
+   *
+   * @param   fd  a file descriptor.
+   * @relates basic_pstreambuf
+   */
+  inline void
+  close_fd(pstreams::fd_type& fd)
+  {
+    if (fd >= 0 && ::close(fd) == 0)
+      fd = -1;
+  }
+
+  /**
+   * @brief  Helper function to close an array of file descriptors.
+   *
+   * Calls @c close_fd() on each member of the array.
+   * The length of the array is determined automatically by
+   * template argument deduction to avoid errors.
+   *
+   * @param   fds  an array of file descriptors.
+   * @relates basic_pstreambuf
+   */
+  template <int N>
+    inline void
+    close_fd_array(pstreams::fd_type (&fds)[N])
+    {
+      for (std::size_t i = 0; i < N; ++i)
+        close_fd(fds[i]);
+    }
+
+  /**
+   * Starts a new process by executing @a file with the arguments in
+   * @a argv and opens pipes to the process with the specified @a mode.
+   *
+   * By convention @c argv[0] should be the file name of the file being
+   * executed.
+   * Will duplicate the actions of  the  shell  in searching for an
+   * executable file if the specified file name does not contain a slash (/)
+   * character.
+   *
+   * Iff @a file is successfully executed then is_open() will return true.
+   * Note that exited() will return true if file cannot be executed, since
+   * the child process will have exited.
+   *
+   * @param   file  a string containing the pathname of a program to execute.
+   * @param   argv  a vector of argument strings passed to the new program.
+   * @param   mode  a bitwise OR of one or more of @c out, @c in and @c err.
+   * @return  NULL if a pipe could not be opened or if the program could
+   *          not be executed, @c this otherwise.
+   * @see     <b>execvp</b>(3)
+   */
+  template <typename C, typename T>
+    basic_pstreambuf<C,T>*
+    basic_pstreambuf<C,T>::open( const std::string& file,
+                                 const argv_type& argv,
+                                 pmode mode )
+    {
+      basic_pstreambuf<C,T>* ret = NULL;
+
+      if (!is_open())
+      {
+        // constants for read/write ends of pipe
+        enum { RD, WR };
+
+        // open another pipe and set close-on-exec
+        fd_type ck_exec[] = { -1, -1 };
+        if (-1 == ::pipe(ck_exec)
+            || -1 == ::fcntl(ck_exec[RD], F_SETFD, FD_CLOEXEC)
+            || -1 == ::fcntl(ck_exec[WR], F_SETFD, FD_CLOEXEC))
+        {
+          error_ = errno;
+          close_fd_array(ck_exec);
+        }
+        else
+        {
+          switch(fork(mode))
+          {
+          case 0 :
+            // this is the new process, exec command
+            {
+              char** arg_v = new char*[argv.size()+1];
+              for (std::size_t i = 0; i < argv.size(); ++i)
+              {
+                const std::string& src = argv[i];
+                char*& dest = arg_v[i];
+                dest = new char[src.size()+1];
+                dest[ src.copy(dest, src.size()) ] = '\0';
+              }
+              arg_v[argv.size()] = NULL;
+
+              ::execvp(file.c_str(), arg_v);
+
+              // can only reach this point if exec() failed
+
+              // parent can get error code from ck_exec pipe
+              error_ = errno;
+
+              ::write(ck_exec[WR], &error_, sizeof(error_));
+              ::close(ck_exec[WR]);
+              ::close(ck_exec[RD]);
+
+              ::_exit(error_);
+              // using std::exit() would make static dtors run twice
+            }
+
+          case -1 :
+            // couldn't fork, error already handled in pstreambuf::fork()
+            close_fd_array(ck_exec);
+            break;
+
+          default :
+            // this is the parent process
+
+            // check child called exec() successfully
+            ::close(ck_exec[WR]);
+            switch (::read(ck_exec[RD], &error_, sizeof(error_)))
+            {
+            case 0:
+              // activate buffers
+              create_buffers(mode);
+              ret = this;
+              break;
+            case -1:
+              error_ = errno;
+              break;
+            default:
+              // error_ contains error code from child
+              // call wait() to clean up and set ppid_ to 0
+              this->wait();
+              break;
+            }
+            ::close(ck_exec[RD]);
+          }
+        }
+      }
+      return ret;
+    }
+
+  /**
+   * Creates pipes as specified by @a mode and calls @c fork() to create
+   * a new process. If the fork is successful the parent process stores
+   * the child's PID and the opened pipes and the child process replaces
+   * its standard streams with the opened pipes.
+   *
+   * If an error occurs the error code will be set to one of the possile
+   * errors for @c pipe() or @c fork().
+   * See your system's documentation for these error codes.
+   *
+   * @param   mode  an OR of pmodes specifying which of the child's
+   *                standard streams to connect to.
+   * @return  On success the PID of the child is returned in the parent's
+   *          context and zero is returned in the child's context.
+   *          On error -1 is returned and the error code is set appropriately.
+   */
+  template <typename C, typename T>
+    pid_t
+    basic_pstreambuf<C,T>::fork(pmode mode)
+    {
+      pid_t pid = -1;
+
+      // Three pairs of file descriptors, for pipes connected to the
+      // process' stdin, stdout and stderr
+      // (stored in a single array so close_fd_array() can close all at once)
+      fd_type fd[] = { -1, -1, -1, -1, -1, -1 };
+      fd_type* const pin = fd;
+      fd_type* const pout = fd+2;
+      fd_type* const perr = fd+4;
+
+      // constants for read/write ends of pipe
+      enum { RD, WR };
+
+      // N.B.
+      // For the pstreambuf pin is an output stream and
+      // pout and perr are input streams.
+
+      if (!error_ && mode&pstdin && ::pipe(pin))
+        error_ = errno;
+
+      if (!error_ && mode&pstdout && ::pipe(pout))
+        error_ = errno;
+
+      if (!error_ && mode&pstderr && ::pipe(perr))
+        error_ = errno;
+
+      if (!error_)
+      {
+        pid = ::fork();
+        switch (pid)
+        {
+          case 0 :
+          {
+            // this is the new process
+
+            // for each open pipe close one end and redirect the
+            // respective standard stream to the other end
+
+            if (*pin >= 0)
+            {
+              ::close(pin[WR]);
+              ::dup2(pin[RD], STDIN_FILENO);
+              ::close(pin[RD]);
+            }
+            if (*pout >= 0)
+            {
+              ::close(pout[RD]);
+              ::dup2(pout[WR], STDOUT_FILENO);
+              ::close(pout[WR]);
+            }
+            if (*perr >= 0)
+            {
+              ::close(perr[RD]);
+              ::dup2(perr[WR], STDERR_FILENO);
+              ::close(perr[WR]);
+            }
+            break;
+          }
+          case -1 :
+          {
+            // couldn't fork for some reason
+            error_ = errno;
+            // close any open pipes
+            close_fd_array(fd);
+            break;
+          }
+          default :
+          {
+            // this is the parent process, store process' pid
+            ppid_ = pid;
+
+            // store one end of open pipes and close other end
+            if (*pin >= 0)
+            {
+              wpipe_ = pin[WR];
+              ::close(pin[RD]);
+            }
+            if (*pout >= 0)
+            {
+              rpipe_[rsrc_out] = pout[RD];
+              ::close(pout[WR]);
+            }
+            if (*perr >= 0)
+            {
+              rpipe_[rsrc_err] = perr[RD];
+              ::close(perr[WR]);
+            }
+
+            if (rpipe_[rsrc_out] == -1 && rpipe_[rsrc_err] >= 0)
+            {
+              // reading stderr but not stdout, so use stderr for all reads
+              read_err(true);
+            }
+          }
+        }
+      }
+      else
+      {
+        // close any pipes we opened before failure
+        close_fd_array(fd);
+      }
+      return pid;
+    }
+
+  /**
+   * Closes all pipes and calls wait() to wait for the process to finish.
+   * If an error occurs the error code will be set to one of the possible
+   * errors for @c waitpid().
+   * See your system's documentation for these errors.
+   *
+   * @return  @c this on successful close or @c NULL if there is no
+   *          process to close or if an error occurs.
+   */
+  template <typename C, typename T>
+    basic_pstreambuf<C,T>*
+    basic_pstreambuf<C,T>::close()
+    {
+      basic_pstreambuf<C,T>* ret = NULL;
+      if (is_open())
+      {
+        sync();
+
+        destroy_buffers(pstdin|pstdout|pstderr);
+
+        // close pipes before wait() so child gets EOF/SIGPIPE
+        close_fd(wpipe_);
+        close_fd_array(rpipe_);
+
+        if (wait() == 1)
+        {
+          ret = this;
+        }
+      }
+      return ret;
+    }
+
+  /**
+   *  Called on construction to initialise the arrays used for reading.
+   */
+  template <typename C, typename T>
+    inline void
+    basic_pstreambuf<C,T>::init_rbuffers()
+    {
+      rpipe_[rsrc_out] = rpipe_[rsrc_err] = -1;
+      rbuffer_[rsrc_out] = rbuffer_[rsrc_err] = NULL;
+      rbufstate_[0] = rbufstate_[1] = rbufstate_[2] = NULL;
+    }
+
+  template <typename C, typename T>
+    void
+    basic_pstreambuf<C,T>::create_buffers(pmode mode)
+    {
+      if (mode & pstdin)
+      {
+        delete[] wbuffer_;
+        wbuffer_ = new char_type[bufsz];
+        this->setp(wbuffer_, wbuffer_ + bufsz);
+      }
+      if (mode & pstdout)
+      {
+        delete[] rbuffer_[rsrc_out];
+        rbuffer_[rsrc_out] = new char_type[bufsz];
+        if (rsrc_ == rsrc_out)
+          this->setg(rbuffer_[rsrc_out] + pbsz, rbuffer_[rsrc_out] + pbsz,
+              rbuffer_[rsrc_out] + pbsz);
+      }
+      if (mode & pstderr)
+      {
+        delete[] rbuffer_[rsrc_err];
+        rbuffer_[rsrc_err] = new char_type[bufsz];
+        if (rsrc_ == rsrc_err)
+          this->setg(rbuffer_[rsrc_err] + pbsz, rbuffer_[rsrc_err] + pbsz,
+              rbuffer_[rsrc_err] + pbsz);
+      }
+    }
+
+  template <typename C, typename T>
+    void
+    basic_pstreambuf<C,T>::destroy_buffers(pmode mode)
+    {
+      if (mode & pstdin)
+      {
+        this->setp(NULL, NULL);
+        delete[] wbuffer_;
+        wbuffer_ = NULL;
+      }
+      if (mode & pstdout)
+      {
+        if (rsrc_ == rsrc_out)
+          this->setg(NULL, NULL, NULL);
+        delete[] rbuffer_[rsrc_out];
+        rbuffer_[rsrc_out] = NULL;
+      }
+      if (mode & pstderr)
+      {
+        if (rsrc_ == rsrc_err)
+          this->setg(NULL, NULL, NULL);
+        delete[] rbuffer_[rsrc_err];
+        rbuffer_[rsrc_err] = NULL;
+      }
+    }
+
+  template <typename C, typename T>
+    typename basic_pstreambuf<C,T>::buf_read_src
+    basic_pstreambuf<C,T>::switch_read_buffer(buf_read_src src)
+    {
+      if (rsrc_ != src)
+      {
+        char_type* tmpbufstate[] = {this->eback(), this->gptr(), this->egptr()};
+        this->setg(rbufstate_[0], rbufstate_[1], rbufstate_[2]);
+        for (std::size_t i = 0; i < 3; ++i)
+          rbufstate_[i] = tmpbufstate[i];
+        rsrc_ = src;
+      }
+      return rsrc_;
+    }
+
+  /**
+   * Suspends execution and waits for the associated process to exit, or
+   * until a signal is delivered whose action is to terminate the current
+   * process or to call a signal handling function. If the process has
+   * already exited wait() returns immediately.
+   *
+   * @param   nohang  true to return immediately if the process has not exited.
+   * @return  1 if the process has exited.
+   *          0 if @a nohang is true and the process has not exited yet.
+   *          -1 if no process has been started or if an error occurs,
+   *          in which case the error can be found using error().
+   */
+  template <typename C, typename T>
+    int
+    basic_pstreambuf<C,T>::wait(bool nohang)
+    {
+      int exited = -1;
+      if (is_open())
+      {
+        int status;
+        switch(::waitpid(ppid_, &status, nohang ? WNOHANG : 0))
+        {
+          case 0 :
+            // nohang was true and process has not exited
+            exited = 0;
+            break;
+          case -1 :
+            error_ = errno;
+            break;
+          default :
+            // process has exited
+            ppid_ = 0;
+            status_ = status;
+            exited = 1;
+            destroy_buffers(pstdin|pstdout|pstderr);
+            close_fd(wpipe_);
+            close_fd_array(rpipe_);
+            break;
+        }
+      }
+      return exited;
+    }
+
+  /**
+   * Sends the specified signal to the process.  A signal can be used to
+   * terminate a child process that would not exit otherwise.
+   *
+   * If an error occurs the error code will be set to one of the possible
+   * errors for @c kill().  See your system's documentation for these errors.
+   *
+   * @param   signal  A signal to send to the child process.
+   * @return  @c this or @c NULL if @c kill() fails.
+   */
+  template <typename C, typename T>
+    inline basic_pstreambuf<C,T>*
+    basic_pstreambuf<C,T>::kill(int signal)
+    {
+      basic_pstreambuf<C,T>* ret = NULL;
+      if (is_open())
+      {
+        if (::kill(ppid_, signal))
+          error_ = errno;
+        else
+        {
+          // TODO call exited() to check for exit and clean up? leave to user?
+          ret = this;
+        }
+      }
+      return ret;
+    }
+
+  /**
+   *  @return  True if the associated process has exited, false otherwise.
+   *  @see     basic_pstreambuf<C,T>::close()
+   */
+  template <typename C, typename T>
+    inline bool
+    basic_pstreambuf<C,T>::exited()
+    {
+      return ppid_ == 0 || wait(true)==1;
+    }
+
+
+  /**
+   *  @return  The exit status of the child process, or -1 if close()
+   *           has not yet been called to wait for the child to exit.
+   *  @see     basic_pstreambuf<C,T>::close()
+   */
+  template <typename C, typename T>
+    inline int
+    basic_pstreambuf<C,T>::status() const
+    {
+      return status_;
+    }
+
+  /**
+   *  @return  The error code of the most recently failed operation, or zero.
+   */
+  template <typename C, typename T>
+    inline int
+    basic_pstreambuf<C,T>::error() const
+    {
+      return error_;
+    }
+
+  /**
+   *  Closes the output pipe, causing the child process to receive the
+   *  end-of-file indicator on subsequent reads from its @c stdin stream.
+   */
+  template <typename C, typename T>
+    inline void
+    basic_pstreambuf<C,T>::peof()
+    {
+      sync();
+      destroy_buffers(pstdin);
+      close_fd(wpipe_);
+    }
+
+  /**
+   * @return  true if a previous call to open() succeeded and wait() has
+   *          not been called and determined that the process has exited,
+   *          false otherwise.
+   * @warning This function can not be used to determine whether the
+   *          command used to initialise the buffer was successfully
+   *          executed or not. If the shell command failed this function
+   *          will still return true.
+   *          You can use exited() to see if it's still open.
+   */
+  template <typename C, typename T>
+    inline bool
+    basic_pstreambuf<C,T>::is_open() const
+    {
+      return ppid_ > 0;
+    }
+
+  /**
+   * Toggle the stream used for reading. If @a readerr is @c true then the
+   * process' @c stderr output will be used for subsequent extractions, if
+   * @a readerr is false the the process' stdout will be used.
+   * @param   readerr  @c true to read @c stderr, @c false to read @c stdout.
+   * @return  @c true if the requested stream is open and will be used for
+   *          subsequent extractions, @c false otherwise.
+   */
+  template <typename C, typename T>
+    inline bool
+    basic_pstreambuf<C,T>::read_err(bool readerr)
+    {
+      buf_read_src src = readerr ? rsrc_err : rsrc_out;
+      if (rpipe_[src]>=0)
+      {
+        switch_read_buffer(src);
+        return true;
+      }
+      return false;
+    }
+
+  /**
+   * Called when the internal character buffer is not present or is full,
+   * to transfer the buffer contents to the pipe.
+   *
+   * @param   c  a character to be written to the pipe.
+   * @return  @c traits_type::not_eof(c) if @a c is equal to @c
+   *          traits_type::eof(). Otherwise returns @a c if @a c can be
+   *          written to the pipe, or @c traits_type::eof() if not.
+   */
+  template <typename C, typename T>
+    typename basic_pstreambuf<C,T>::int_type
+    basic_pstreambuf<C,T>::overflow(int_type c)
+    {
+      if (!empty_buffer())
+        return traits_type::eof();
+      else if (!traits_type::eq_int_type(c, traits_type::eof()))
+        return this->sputc(c);
+      else
+        return traits_type::not_eof(c);
+    }
+
+
+  template <typename C, typename T>
+    int
+    basic_pstreambuf<C,T>::sync()
+    {
+      return !exited() && empty_buffer() ? 0 : -1;
+    }
+
+  /**
+   * @param   s  character buffer.
+   * @param   n  buffer length.
+   * @return  the number of characters written.
+   */
+  template <typename C, typename T>
+    std::streamsize
+    basic_pstreambuf<C,T>::xsputn(const char_type* s, std::streamsize n)
+    {
+      if (n < this->epptr() - this->pptr())
+      {
+        std::memcpy(this->pptr(), s, n * sizeof(char_type));
+        this->pbump(n);
+        return n;
+      }
+      else
+      {
+        for (std::streamsize i = 0; i < n; ++i)
+        {
+          if (traits_type::eq_int_type(this->sputc(s[i]), traits_type::eof()))
+            return i;
+        }
+        return n;
+      }
+    }
+
+  /**
+   * @return  true if the buffer was emptied, false otherwise.
+   */
+  template <typename C, typename T>
+    bool
+    basic_pstreambuf<C,T>::empty_buffer()
+    {
+      const std::streamsize count = this->pptr() - this->pbase();
+      const std::streamsize written = this->write(this->wbuffer_, count);
+      if (count > 0 && written == count)
+      {
+        this->pbump(-written);
+        return true;
+      }
+      return false;
+    }
+
+  /**
+   * Called when the internal character buffer is is empty, to re-fill it
+   * from the pipe.
+   *
+   * @return The first available character in the buffer,
+   * or @c traits_type::eof() in case of failure.
+   */
+  template <typename C, typename T>
+    typename basic_pstreambuf<C,T>::int_type
+    basic_pstreambuf<C,T>::underflow()
+    {
+      if (this->gptr() < this->egptr() || fill_buffer())
+        return traits_type::to_int_type(*this->gptr());
+      else
+        return traits_type::eof();
+    }
+
+  /**
+   * Attempts to make @a c available as the next character to be read by
+   * @c sgetc().
+   *
+   * @param   c   a character to make available for extraction.
+   * @return  @a c if the character can be made available,
+   *          @c traits_type::eof() otherwise.
+   */
+  template <typename C, typename T>
+    typename basic_pstreambuf<C,T>::int_type
+    basic_pstreambuf<C,T>::pbackfail(int_type c)
+    {
+      if (this->gptr() != this->eback())
+      {
+        this->gbump(-1);
+        if (!traits_type::eq_int_type(c, traits_type::eof()))
+          *this->gptr() = traits_type::to_char_type(c);
+        return traits_type::not_eof(c);
+      }
+      else
+         return traits_type::eof();
+    }
+
+  template <typename C, typename T>
+    std::streamsize
+    basic_pstreambuf<C,T>::showmanyc()
+    {
+      int avail = 0;
+#ifdef FIONREAD
+      if (ioctl(rpipe(), FIONREAD, &avail) == -1)
+        avail = -1;
+      else
+#endif
+      if (const std::ptrdiff_t buflen = this->gptr() - this->eback())
+        avail += buflen;
+      return std::streamsize(avail);
+    }
+
+  /**
+   * @return  true if the buffer was filled, false otherwise.
+   */
+  template <typename C, typename T>
+    bool
+    basic_pstreambuf<C,T>::fill_buffer()
+    {
+      const std::streamsize pb1 = this->gptr() - this->eback();
+      const std::streamsize pb2 = pbsz;
+      const std::streamsize npb = std::min(pb1, pb2);
+
+      std::memmove( rbuffer() + pbsz - npb,
+                    this->gptr() - npb,
+                    npb * sizeof(char_type) );
+
+      const std::streamsize rc = read(rbuffer() + pbsz, bufsz - pbsz);
+
+      if (rc > 0)
+      {
+        this->setg( rbuffer() + pbsz - npb,
+                    rbuffer() + pbsz,
+                    rbuffer() + pbsz + rc );
+        return true;
+      }
+      else
+      {
+        this->setg(NULL, NULL, NULL);
+        return false;
+      }
+    }
+
+  /**
+   * Writes up to @a n characters to the pipe from the buffer @a s.
+   * This currently only works for fixed width character encodings where
+   * each character uses @c sizeof(char_type) bytes.
+   *
+   * @param   s  character buffer.
+   * @param   n  buffer length.
+   * @return  the number of characters written.
+   */
+  template <typename C, typename T>
+    inline std::streamsize
+    basic_pstreambuf<C,T>::write(const char_type* s, std::streamsize n)
+    {
+      return wpipe() >= 0 ? ::write(wpipe(), s, n * sizeof(char_type)) : 0;
+    }
+
+  /**
+   * Reads up to @a n characters from the pipe to the buffer @a s.
+   * This currently only works for fixed width character encodings where
+   * each character uses @c sizeof(char_type) bytes.
+   *
+   * @param   s  character buffer.
+   * @param   n  buffer length.
+   * @return  the number of characters read.
+   */
+  template <typename C, typename T>
+    inline std::streamsize
+    basic_pstreambuf<C,T>::read(char_type* s, std::streamsize n)
+    {
+      return rpipe() >= 0 ? ::read(rpipe(), s, n * sizeof(char_type)) : 0;
+    }
+
+  /** @return a reference to the output file descriptor */
+  template <typename C, typename T>
+    inline typename basic_pstreambuf<C,T>::fd_type&
+    basic_pstreambuf<C,T>::wpipe()
+    {
+      return wpipe_;
+    }
+
+  /** @return a reference to the active input file descriptor */
+  template <typename C, typename T>
+    inline typename basic_pstreambuf<C,T>::fd_type&
+    basic_pstreambuf<C,T>::rpipe()
+    {
+      return rpipe_[rsrc_];
+    }
+
+  /** @return a reference to the specified input file descriptor */
+  template <typename C, typename T>
+    inline typename basic_pstreambuf<C,T>::fd_type&
+    basic_pstreambuf<C,T>::rpipe(buf_read_src which)
+    {
+      return rpipe_[which];
+    }
+
+  /** @return a pointer to the start of the active input buffer area. */
+  template <typename C, typename T>
+    inline typename basic_pstreambuf<C,T>::char_type*
+    basic_pstreambuf<C,T>::rbuffer()
+    {
+      return rbuffer_[rsrc_];
+    }
+
+
+  /*
+   * member definitions for pstream_common
+   */
+
+  /**
+   * @class pstream_common
+   * Abstract Base Class providing common functionality for basic_ipstream,
+   * basic_opstream and basic_pstream.
+   * pstream_common manages the basic_pstreambuf stream buffer that is used
+   * by the derived classes to initialise an IOStream class.
+   */
+
+  /** Creates an uninitialised stream. */
+  template <typename C, typename T>
+    inline
+    pstream_common<C,T>::pstream_common()
+    : std::basic_ios<C,T>(NULL)
+    , command_()
+    , buf_()
+    {
+      this->init(&buf_);
+    }
+
+  /**
+   * Initialises the stream buffer by calling
+   * do_open( @a command , @a mode )
+   *
+   * @param command a string containing a shell command.
+   * @param mode    the I/O mode to use when opening the pipe.
+   * @see   do_open(const std::string&, pmode)
+   */
+  template <typename C, typename T>
+    inline
+    pstream_common<C,T>::pstream_common(const std::string& command, pmode mode)
+    : std::basic_ios<C,T>(NULL)
+    , command_(command)
+    , buf_()
+    {
+      this->init(&buf_);
+      do_open(command, mode);
+    }
+
+  /**
+   * Initialises the stream buffer by calling
+   * do_open( @a file , @a argv , @a mode )
+   *
+   * @param file  a string containing the pathname of a program to execute.
+   * @param argv  a vector of argument strings passed to the new program.
+   * @param mode  the I/O mode to use when opening the pipe.
+   * @see do_open(const std::string&, const argv_type&, pmode)
+   */
+  template <typename C, typename T>
+    inline
+    pstream_common<C,T>::pstream_common( const std::string& file,
+                                         const argv_type& argv,
+                                         pmode mode )
+    : std::basic_ios<C,T>(NULL)
+    , command_(file)
+    , buf_()
+    {
+      this->init(&buf_);
+      do_open(file, argv, mode);
+    }
+
+  /**
+   * This is a pure virtual function to make @c pstream_common abstract.
+   * Because it is the destructor it will be called by derived classes
+   * and so must be defined.  It is also protected, to discourage use of
+   * the PStreams classes through pointers or references to the base class.
+   *
+   * @sa If defining a pure virtual seems odd you should read
+   * http://www.gotw.ca/gotw/031.htm (and the rest of the site as well!)
+   */
+  template <typename C, typename T>
+    inline
+    pstream_common<C,T>::~pstream_common()
+    {
+    }
+
+  /**
+   * Calls rdbuf()->open( @a command , @a mode )
+   * and sets @c failbit on error.
+   *
+   * @param command a string containing a shell command.
+   * @param mode    the I/O mode to use when opening the pipe.
+   * @see   basic_pstreambuf::open(const std::string&, pmode)
+   */
+  template <typename C, typename T>
+    inline void
+    pstream_common<C,T>::do_open(const std::string& command, pmode mode)
+    {
+      if (!buf_.open((command_=command), mode))
+        this->setstate(std::ios_base::failbit);
+    }
+
+  /**
+   * Calls rdbuf()->open( @a file, @a  argv, @a mode )
+   * and sets @c failbit on error.
+   *
+   * @param file  a string containing the pathname of a program to execute.
+   * @param argv  a vector of argument strings passed to the new program.
+   * @param mode  the I/O mode to use when opening the pipe.
+   * @see   basic_pstreambuf::open(const std::string&, const argv_type&, pmode)
+   */
+  template <typename C, typename T>
+    inline void
+    pstream_common<C,T>::do_open( const std::string& file,
+                                  const argv_type& argv,
+                                  pmode mode )
+    {
+      if (!buf_.open((command_=file), argv, mode))
+        this->setstate(std::ios_base::failbit);
+    }
+
+  /** Calls rdbuf->close() and sets @c failbit on error. */
+  template <typename C, typename T>
+    inline void
+    pstream_common<C,T>::close()
+    {
+      if (!buf_.close())
+        this->setstate(std::ios_base::failbit);
+    }
+
+  /** Calls rdbuf->kill(int) and sets @c failbit on error. */
+  template <typename C, typename T>
+    inline void
+    pstream_common<C,T>::kill(int signal)
+    {
+      if (!buf_.kill(signal))
+        this->setstate(std::ios_base::failbit);
+    }
+
+  /**
+   * @return  rdbuf()->is_open().
+   * @see     basic_pstreambuf::is_open()
+   */
+  template <typename C, typename T>
+    inline bool
+    pstream_common<C,T>::is_open() const
+    {
+      return buf_.is_open();
+    }
+
+  /** @return a string containing the command used to initialise the stream. */
+  template <typename C, typename T>
+    inline const std::string&
+    pstream_common<C,T>::command() const
+    {
+      return command_;
+    }
+
+  /** @return a pointer to the private stream buffer member. */
+  // TODO  document behaviour if buffer replaced.
+  template <typename C, typename T>
+    inline typename pstream_common<C,T>::streambuf_type*
+    pstream_common<C,T>::rdbuf() const
+    {
+      return const_cast<streambuf_type*>(&buf_);
+    }
+
+
+#if REDI_EVISCERATE_PSTREAMS
+  /**
+   * @def REDI_EVISCERATE_PSTREAMS
+   * If this macro has a non-zero value then certain internals of the
+   * @c basic_pstreambuf template class are exposed. In general this is
+   * a Bad Thing, as the internal implementation is largely undocumented
+   * and may be subject to change at any time, so this feature is only
+   * provided because it might make PStreams useful in situations where
+   * it is necessary to do Bad Things.
+   */
+
+  /**
+   * @warning  This function exposes the internals of the stream buffer and
+   *           should be used with caution. It is the caller's responsibility
+   *           to flush streams etc. in order to clear any buffered data.
+   *           The POSIX.1 function <b>fdopen</b>(3) is used to obtain the
+   *           @c FILE pointers from the streambuf's private file descriptor
+   *           members so consult your system's documentation for
+   *           <b>fdopen</b>(3).
+   *
+   * @param   in    A FILE* that will refer to the process' stdin.
+   * @param   out   A FILE* that will refer to the process' stdout.
+   * @param   err   A FILE* that will refer to the process' stderr.
+   * @return  An OR of zero or more of @c pstdin, @c pstdout, @c pstderr.
+   *
+   * For each open stream shared with the child process a @c FILE* is
+   * obtained and assigned to the corresponding parameter. For closed
+   * streams @c NULL is assigned to the parameter.
+   * The return value can be tested to see which parameters should be
+   * @c !NULL by masking with the corresponding @c pmode value.
+   *
+   * @see <b>fdopen</b>(3)
+   */
+  template <typename C, typename T>
+    std::size_t
+    basic_pstreambuf<C,T>::fopen(std::FILE*& in, std::FILE*& out, std::FILE*& err)
+    {
+      in = out = err = NULL;
+      std::size_t open_files = 0;
+      if (wpipe() > -1)
+      {
+        if ((in = ::fdopen(wpipe(), "w")))
+        {
+            open_files |= pstdin;
+        }
+      }
+      if (rpipe(rsrc_out) > -1)
+      {
+        if ((out = ::fdopen(rpipe(rsrc_out), "r")))
+        {
+            open_files |= pstdout;
+        }
+      }
+      if (rpipe(rsrc_err) > -1)
+      {
+        if ((err = ::fdopen(rpipe(rsrc_err), "r")))
+        {
+            open_files |= pstderr;
+        }
+      }
+      return open_files;
+    }
+
+  /**
+   *  @warning This function exposes the internals of the stream buffer and
+   *  should be used with caution.
+   *
+   *  @param  in   A FILE* that will refer to the process' stdin.
+   *  @param  out  A FILE* that will refer to the process' stdout.
+   *  @param  err  A FILE* that will refer to the process' stderr.
+   *  @return A bitwise-or of zero or more of @c pstdin, @c pstdout, @c pstderr.
+   *  @see    basic_pstreambuf::fopen()
+   */
+  template <typename C, typename T>
+    inline std::size_t
+    pstream_common<C,T>::fopen(std::FILE*& in, std::FILE*& out, std::FILE*& err)
+    {
+      return buf_.fopen(in, out, err);
+    }
+
+#endif // REDI_EVISCERATE_PSTREAMS
+
+
+// Returns 0 for ifstream, 1 for ipstream, 2 for istringstream:
+inline int isinpipe(const string &commandin, string &commandout)
+{
+    if ( commandin.length()>=1 )
+    {
+       if ( commandin[commandin.length()-1]=='|' )
+       {
+           commandout=commandin.substr(0, commandin.length()-1);
+           return 1;
+       }
+    }
+    if ( commandin.length()>=2 )
+    {
+       if ( commandin.substr(0, 2)=="<<" )
+       {
+           commandout=commandin.substr(2, commandin.length()-2);
+           return 2;
+       }
+    }
+    if ( commandin.length()>=1 )
+    {
+       if ( commandin[0]=='<' )
+       {
+           commandout=commandin.substr(1, commandin.length()-1);
+           return 0;
+       }
+    }
+    commandout=commandin;
+    return 0;
+}
+
+
+// Returns 0 for ofstream, 1 for opstream, 2 for ofstream in append mode:
+inline int isoutpipe(const string &commandin, string &commandout)
+{
+    if ( commandin.length()>=1 )
+    {
+       if ( commandin[0]=='|' )
+       {
+           commandout=commandin.substr(1, commandin.length()-1);
+           return 1;
+       }
+       if ( commandin[0]=='>' )
+       {
+           commandout=commandin.substr(1, commandin.length()-1);
+           return 0;
+       }
+       if ( commandin.length()>=2 )
+       {
+           if ( commandin.substr(0, 2)==">>" )
+           {
+               commandout=commandin.substr(2, commandin.length()-2);
+               return 2;
+           }
+       }
+    }
+    commandout=commandin;
+    return 0;
+}
+
+
+inline istream* openin(const string &commandin, const ios::openmode &mode=ios::in)
+{
+    istream *result=0;
+    string commandout;
+    int flag=isinpipe(commandin, commandout);
+    if ( flag==0 ) result=new ifstream(commandout.c_str(), mode);
+    else if ( flag==1 ) result=new ipstream(commandout.c_str(), mode);
+    else if ( flag==2 ) result=new istringstream(commandout.c_str(), mode);
+    else result=new ifstream(commandout.c_str(), mode);
+    return result;
+}
+
+
+inline ostream* openout(const string &commandin, const ios::openmode &mode=ios::out)
+{
+    ostream *result=0;
+    string commandout;
+    int flag=isoutpipe(commandin, commandout);
+    if ( flag==0 ) result=new ofstream(commandout.c_str(), mode);
+    else if ( flag==1 ) result=new opstream(commandout.c_str(), mode);
+    else if ( flag==2 ) result=new ofstream(commandout.c_str(), ios::app|mode);
+    else result=new ofstream(commandout.c_str(), mode);
+    return result;
+}
+
+
+} // namespace redi
+
+
+/**
+ * @mainpage PStreams Reference
+ * @htmlinclude mainpage.html
+ */
+
+#endif  // REDI_PSTREAM_H_SEEN
+
+// vim: ts=2 sw=2 expandtab
+
diff --git a/inputs/example/kgminkowski/analysis/inp/intq.conf b/inputs/example/kgminkowski/analysis/inp/intq.conf
new file mode 100644 (file)
index 0000000..6a17988
--- /dev/null
@@ -0,0 +1,16 @@
+#################### Configuration file for intq #######################
+
+### Input file name:
+inputFile      /home/laszloa/GIT/GR/gridripperprod/inputs/example/kgminkowski/kgminkowski.in
+
+### Output file name, containing the itegrated data:
+intqFile       out/data/intq.dat
+
+### Output file name, containing density data:
+densityFile    out/data/density.dat
+
+### Logfile name:
+logFile                out/data/intq.log
+
+### EPSILON:
+EPSILON                1e-8
diff --git a/inputs/example/kgminkowski/analysis/lib/config.cc b/inputs/example/kgminkowski/analysis/lib/config.cc
new file mode 100644 (file)
index 0000000..18bbed6
--- /dev/null
@@ -0,0 +1,722 @@
+/**
+ * config.cc  Implements tools for reading config files.
+ */
+
+
+#include "config.h"
+
+
+//////////////////// Implemetation of class config_exception ///////////
+
+
+config_exception::config_exception(const std::string& function_name, const std::string& message)
+ : function_name_(function_name), message_(message)
+{
+}
+
+
+const std::string& config_exception::function_name() const
+{
+    return function_name_;
+}
+
+
+const std::string& config_exception::message() const
+{
+    return message_;
+}
+
+
+const char* config_exception::what() const throw()
+{
+    return (std::string("[config] ")+message_+"\n[config]\t"+function_name_+"\n").c_str();
+}
+
+
+config_exception::~config_exception() throw()
+{
+}
+
+
+void config_exception::print_error(const std::string& function_name, const std::string& message)
+{
+    std::cerr<<"[config] "<<message<<"\n[config]\t"<<function_name<<std::endl;
+}
+
+
+//////////////////// Implemetation of class config_entry ///////////////
+
+
+config_entry::config_entry()
+ : string_(), longdouble_(0.0), success_(false)
+{
+}
+
+
+config_entry::config_entry(const config_entry &other)
+ : string_(other.string_), longdouble_(other.longdouble_), success_(other.success_)
+{
+}
+
+
+config_entry& config_entry::reinit(const std::string &str)
+{
+    string_=str;
+    std::istringstream iss(string_);
+    long double zero=0.0;
+    long double minusone=-1.0;
+    if ( !(iss>>longdouble_) )
+    {
+       if ( string_=="inf" ) { longdouble_=1.0/zero; success_=true; }
+       else if ( string_=="-inf" ) { longdouble_=-1.0/zero; success_=true; }
+       else if ( string_=="nan" ) { longdouble_=::sqrtl(minusone); success_=true; }
+       else { longdouble_=0.0; success_=false; }
+    }
+    else { success_=true; }
+    return *this;
+}
+
+
+config_entry::config_entry(const std::string &str)
+ : string_(), longdouble_(0.0), success_(false)
+{
+    reinit(str);
+}
+
+
+config_entry::~config_entry()
+{
+}
+
+
+config_entry::operator std::string() const
+{
+    return string_;
+}
+
+
+config_entry::operator float() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as float.";
+       const std::string func = "config_entry::operator float() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return (float)longdouble_;
+}
+
+
+config_entry::operator double() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as double.";
+       const std::string func = "config_entry::operator double() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return (double)longdouble_;
+}
+
+
+config_entry::operator long double() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as long double.";
+       const std::string func = "config_entry::operator long double() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return longdouble_;
+}
+
+
+config_entry::operator int() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as int.";
+       const std::string func = "config_entry::operator int() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    const int result=(int)longdouble_;
+    if ( (long double)result!=longdouble_ )
+    {
+       const std::string msg = std::string("Problems while interpreting entry \"") + string_ + "\" as int.";
+       const std::string func = "config_entry::operator int() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return result;
+}
+
+
+config_entry::operator unsigned int() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as unsigned int.";
+       const std::string func = "config_entry::operator unsigned int() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    const unsigned int result=(unsigned int)longdouble_;
+    if ( (double)result!=longdouble_ )
+    {
+       const std::string msg = std::string("Problems while interpreting entry \"") + string_ + "\" as unsigned int.";
+       const std::string func = "config_entry::operator unsigned int() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return result;
+}
+
+
+config_entry::operator bool() const
+{
+    if ( !success_ )
+    {
+       const std::string msg = std::string("Could not interpret entry \"") + string_ + "\" as bool.";
+       const std::string func = "config_entry::operator bool() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    const bool result=(bool)longdouble_;
+    if ( (long double)result!=longdouble_ )
+    {
+       const std::string msg = std::string("Problems while interpreting entry \"") + string_ + "\" as bool.";
+       const std::string func = "config_entry::operator bool() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return result;
+}
+
+
+config_entry::operator char() const
+{
+    if ( string_.empty() ) return '\0';
+    if ( string_.length()>1 )
+    {
+       const std::string msg = std::string("String \"") + string_ + "\" consists of more than 1 characters.";
+       const std::string func = "config_entry::operator char() const";
+       config_exception::print_error(func, msg);
+       throw config_exception(func, msg);
+    }
+    return string_[0];
+}
+
+
+std::string config_entry::string() const
+{
+    return string_;
+}
+
+
+long double config_entry::ldbl() const
+{
+    return longdouble_;
+}
+
+
+bool config_entry::success() const
+{
+    return success_;
+}
+
+
+const config_entry& config_entry::operator=(const config_entry &other)
+{
+    string_=other.string_;
+    longdouble_=other.longdouble_;
+    success_=other.success_;
+    return *this;
+}
+
+
+//////////////////// Implementation of class config ////////////////////
+
+config::config()
+ : std::map< std::string, config_entry >()
+{
+}
+
+
+config::config(const config &other)
+ : std::map< std::string, config_entry >(other)
+{
+}
+
+
+config::~config()
+{
+}
+
+
+config& config::append(std::istream &file)
+{
+    std::string linebuff;
+    std::istringstream iss;
+    std::string token, value_str1, value_str, command, buff;
+    while ( std::getline(file, linebuff) )
+    {
+       if ( linebuff.length()==0 ) continue;
+       if ( linebuff[0]=='#' ) continue;
+       unsigned int n = 0;
+       for ( unsigned int i=0 ; i<linebuff.length() ; ++i )
+       {
+           if ( linebuff[i]==' ' || linebuff[i]=='\t' || linebuff[i]=='=' ) break;
+           ++n;
+       }
+       linebuff[n] = ' ';
+       iss.str(linebuff);
+       if ( !(iss>>token) ) { iss.str(""); iss.clear(); continue; }
+       if ( !(iss>>value_str1) )
+       {
+           const std::string msg = std::string("Could not extract config value from line:\n[config] ")+linebuff;
+           const std::string func = "config& config::append(std::ifstream&)";
+           config_exception::print_error(func, msg);
+           throw(config_exception(func, msg));
+           iss.str("");
+           iss.clear();
+           continue;
+       }
+       iss.str("");
+       iss.clear();
+       const unsigned int k = linebuff.find(token);
+       buff = linebuff.substr(k+token.length()+1, std::string::npos);
+       const unsigned int m = buff.find(value_str1);
+       value_str = buff.substr(m, std::string::npos);
+       if ( !value_str.empty() && value_str[0]=='`' && value_str[value_str.length()-1]=='`' )
+       {
+           command = value_str.substr(1, value_str.length()-2);
+           std::ipstream* ipstr = new std::ipstream(command.c_str());
+           if ( !getline(*ipstr, buff) )
+           {
+               const std::string msg = std::string("Could not read command pipe:\n[config] ")+command;
+               const std::string func = "config& config::append(std::ifstream&)";
+               config_exception::print_error(func, msg);
+               delete ipstr;
+               throw(config_exception(func, msg));
+           }
+           delete ipstr;
+           value_str = buff;
+       }
+       (*this)[token]=config_entry(value_str);
+    }
+    return *this;
+}
+
+
+config& config::append(const std::string &filename)
+{
+    std::istream *file=std::openin(filename);
+    if ( !(*file) )
+    {
+       const std::string msg = std::string("Could not open config file ")+filename+" !";
+       const std::string func = "config& config::append(const std::string&)";
+       config_exception::print_error(func, msg);
+       delete file;
+       throw(config_exception(func, msg));
+       return *this;
+    }
+    append(*file);
+    delete file;
+    return *this;
+}
+
+
+config::config(std::istream &file)
+ : std::map< std::string, config_entry >()
+{
+    append(file);
+}
+
+
+config::config(const std::string &filename)
+ : std::map< std::string, config_entry >()
+{
+    append(filename);
+}
+
+
+config& config::clear()
+{
+    std::map< std::string, config_entry >::clear();
+    return *this;
+}
+
+
+//////////////////// Implementation of getconfig functions /////////////
+
+
+config_entry getconfig(const config &conf, const std::string &token)
+{
+    const config::const_iterator it = conf.find(token);
+    if ( it == conf.end() )
+    {
+       const std::string msg = token+" is not specified in configfile!";
+       const std::string func = "config_entry getconfig(config&, const std::string&)";
+       config_exception::print_error(func, msg);
+       throw(config_exception(func, msg));
+       return config_entry();
+    }
+    return (it->second);
+}
+
+
+config_entry getconfig(std::istream &in, const std::string &token)
+{
+    const config conf(in);
+    return getconfig(conf, token);
+}
+
+
+config_entry getconfig(const std::string &filename, const std::string &token)
+{
+    const config conf(filename);
+    return getconfig(conf, token);
+}
+
+
+//////////////////// Implementation of filenumbering function //////////
+
+
+std::string filenumbering(const unsigned int no, const unsigned int maxno)
+{
+    const int ndigits=(no!=0 ? (int)(::log10((double)no))+1 : 1);
+    const int ndigitsmax=(maxno!=0 ? (int)(::log10((double)maxno))+1 : 1);
+    std::ostringstream oss;
+    for ( int i=0 ; i<(int)ndigitsmax-(int)ndigits ; ++i ) oss<<"0";
+    oss<<no;
+    return oss.str();
+}
+
+
+//////////////////// Implementation of getrevision function ////////////
+
+
+#ifndef REVISIONSTRING
+#define REVISIONSTRING "unknown"
+#endif
+std::string getrevision()
+{
+    return std::string(REVISIONSTRING);
+}
+
+
+//////////////////// Implementation of gettime function ////////////////
+
+
+unsigned int gettime()
+{
+    time_t TIME;
+    time(&TIME);
+    return (unsigned int)TIME;
+}
+
+
+//////////////////// Implementation of gettime_hr function /////////////
+
+
+std::string gettime_hr()
+{
+    time_t TIME;
+    char chartime[32];
+    time(&TIME);
+    ctime_r(&TIME, chartime);
+    return std::string(chartime);
+}
+
+
+//////////////////// Implementation of getusertime function ////////////
+
+
+unsigned int getusertime()
+{
+    rusage usage;
+    getrusage(RUSAGE_SELF, &usage);
+    return usage.ru_utime.tv_sec;
+}
+
+
+//////////////////// Implementation of getcputime function /////////////
+
+
+unsigned int getcputime()
+{
+    rusage usage;
+    getrusage(RUSAGE_SELF, &usage);
+    return usage.ru_stime.tv_sec;
+}
+
+
+//////////////////// Implementation of getmem function ////////////////
+
+
+unsigned int getmem()
+{
+    const pid_t pid=getpid();
+    std::ostringstream str;
+    str<<"/proc/"<<pid<<"/stat";
+    std::ifstream file;
+    file.open(str.str().c_str());
+    std::string linebuff;
+    unsigned int i=0;
+    while ( std::getline(file, linebuff, ' ') )
+    {
+       ++i;
+       if ( i>=23 ) break;
+    }
+    file.close();
+    file.clear();
+    std::istringstream iss(linebuff);
+    unsigned int mem=0;
+    if ( i<23 || !(iss>>mem) )
+    {
+       const std::string msg = "Could not get memory!";
+       const std::string func = "unsigned int getmem()";
+       config_exception::print_error(func, msg);
+       return 0;
+    }
+    return mem/1048576;
+}
+
+
+//////////////////// Implementation of putlogdata function /////////////
+
+
+void putlogdata(const std::string &filename)
+{
+    std::ostream *logfile=std::openout(filename);
+    if ( !(*logfile) )
+    {
+       const std::string msg = std::string("Could not open logfile ")+filename+" !";
+       const std::string func = "void putlogdata(const std::string&, const unsigned int)";
+       config_exception::print_error(func, msg);
+       delete logfile;
+       return;
+    }
+    if ( !((*logfile)<<"Revision: "<<getrevision()<<"\n"<<"Started at: "<<summaryinfo::getstarttime_hr()<<"Running time: "<<getusertime()<<" seconds\n"<<"CPU time: "<<getcputime()<<" seconds\n"<<"Memory usage: "<<getmem()<<" MegaBytes\n"<<"Logged at: "<<gettime_hr()<<std::endl) )
+    {
+       const std::string msg = std::string("Could not write into logfile ")+filename+" !";
+       const std::string func = "void putlogdata(const std::string&, const unsigned int)";
+       config_exception::print_error(func, msg);
+       delete logfile;
+       return;
+    }
+    delete logfile;
+}
+
+
+//////////////////// Implementation of class summaryinfo ///////////////
+
+
+#ifndef __NO_MEMORY_WATCHER
+static std::ipstream& mkauxmaxmemstream()
+{
+    const pid_t pid=getpid();
+    if ( pid==0 )
+    {
+       const std::string msg = "Zero pid!";
+       const std::string func = "static std::ipstream& mkauxmaxmemstream()";
+       config_exception::print_error(func, msg);
+    }
+    std::ostringstream strpid;
+    strpid<<pid;
+    std::ipstream *auxmaxmemstream=new std::ipstream(((const std::string)"while (( "+strpid.str()+" )) ; do sleep 1h ; done").c_str());
+    return *auxmaxmemstream;
+}
+
+
+static std::ipstream& mkmaxmemstream()
+{
+    const pid_t pid=getpid();
+    std::ostringstream strpid;
+    strpid<<pid;
+    std::ipstream *maxmemstream=new std::ipstream(((const std::string)"__MEM=0 ; __MAXMEM=0 ; while ( ps aux | grep -v \"grep\" | grep \"while (( "+strpid.str()+" ))\" >&/dev/null ) ; do __MEM=`cat /proc/"+strpid.str()+"/stat 2>/dev/null | cut -f 23 -d ' '` ; if (( $__MEM > $__MAXMEM )) ; then __MAXMEM=$__MEM ; fi ; sleep 1s ; done ; echo \"$__MAXMEM\" ; unset __MEM ; unset __MAXMEM").c_str());
+    return *maxmemstream;
+}
+#endif
+
+
+
+
+const unsigned int summaryinfo::starttime_=gettime();
+
+const std::string summaryinfo::starttime_hr_=gettime_hr();
+
+unsigned int summaryinfo::maxmemory_=getmem();
+
+#ifndef __NO_MEMORY_WATCHER
+std::ipstream& summaryinfo::auxmaxmemstream_=mkauxmaxmemstream();
+
+std::ipstream& summaryinfo::maxmemstream_=mkmaxmemstream();
+#endif
+
+std::ostream& summaryinfo::logfile_=std::cerr;
+
+std::string summaryinfo::logfilename_="";
+
+bool summaryinfo::firstcopydeclared_=false;
+
+bool summaryinfo::iswritten_=false;
+
+bool summaryinfo::writelogfile_=true;
+
+
+summaryinfo::summaryinfo()
+ : firstcopy_(!firstcopydeclared_)
+{
+    firstcopydeclared_=true;
+}
+
+
+summaryinfo::summaryinfo(const summaryinfo &other)
+ : firstcopy_(!firstcopydeclared_)
+{
+    firstcopydeclared_=true;
+}
+
+
+summaryinfo::~summaryinfo()
+{
+    if ( firstcopy_==true )
+    {
+       acquiremaxmem();
+#ifndef __NO_MEMORY_WATCHER
+       auxmaxmemstream_.kill();
+       auxmaxmemstream_.close();
+       auxmaxmemstream_.clear();
+       unsigned int membuff=0;
+       if ( !(maxmemstream_>>membuff) )
+       {
+           const std::string msg = "Could not read maximal memory stream!";
+           const std::string func = "summaryinfo::~summaryinfo()";
+           config_exception::print_error(func, msg);
+           membuff=0;
+       }
+       maxmemstream_.close();
+       maxmemstream_.clear();
+       membuff/=1048576;
+       if ( membuff>maxmemory_ ) maxmemory_=membuff;
+#endif
+    }
+    if ( firstcopy_==true && iswritten_==false && writelogfile_==true )
+    {
+       if ( !logfilename_.empty() )
+       {
+           std::ostream *logfile=std::openout(logfilename_);
+           if ( !(*logfile) )
+           {
+               const std::string msg = std::string("Could not open logfile ")+logfilename_+" !";
+               const std::string func = "summaryinfo::~summaryinfo()";
+               config_exception::print_error(func, msg);
+               delete logfile;
+               iswritten_=true;
+               return;
+           }
+           if ( !((*logfile)<<"Started at: "<<getstarttime_hr()<<"Running time: "<<getusertime()<<" seconds\n"<<"CPU time: "<<getcputime()<<" seconds\n"<<"Memory usage: "<<maxmemory_<<" MegaBytes\n"<<"Ended at: "<<gettime_hr()<<std::endl) )
+           {
+               const std::string msg = std::string("Could not write into logfile ")+logfilename_+" !";
+               const std::string func = "summaryinfo::~summaryinfo()";
+               config_exception::print_error(func, msg);
+               delete logfile;
+               iswritten_=true;
+               return;
+           }
+           delete logfile;
+           iswritten_=true;
+       }
+       else
+       {
+           if ( !(logfile_) )
+           {
+               const std::string msg = "Could not open logfile!";
+               const std::string func = "summaryinfo::~summaryinfo()";
+               config_exception::print_error(func, msg);
+               iswritten_=true;
+               return;
+           }
+           if ( !(logfile_<<"Started at: "<<getstarttime_hr()<<"Running time: "<<getusertime()<<" seconds\n"<<"CPU time: "<<getcputime()<<" seconds\n"<<"Memory usage: "<<maxmemory_<<" MegaBytes\n"<<"Ended at: "<<gettime_hr()<<std::endl) )
+           {
+               const std::string msg = "Could not write into logfile!";
+               const std::string func = "summaryinfo::~summaryinfo()";
+               config_exception::print_error(func, msg);
+               iswritten_=true;
+               return;
+           }
+           iswritten_=true;
+       }
+    }
+}
+
+
+const summaryinfo& summaryinfo::operator=(const summaryinfo &other)
+{
+    return *this;
+}
+
+
+void summaryinfo::logfilename(const std::string &filename)
+{
+    logfilename_=filename;
+    acquiremaxmem();
+}
+
+
+std::string summaryinfo::logfilename()
+{
+    acquiremaxmem();
+    return logfilename_;
+}
+
+
+void summaryinfo::write_logfile(const bool flag)
+{
+    writelogfile_=flag;
+    acquiremaxmem();
+}
+
+
+bool summaryinfo::write_logfile()
+{
+    acquiremaxmem();
+    return writelogfile_;
+}
+
+
+unsigned int summaryinfo::getstarttime()
+{
+    return starttime_;
+}
+
+
+std::string summaryinfo::getstarttime_hr()
+{
+    return starttime_hr_;
+}
+
+
+void summaryinfo::acquiremaxmem()
+{
+    unsigned int memory=getmem();
+    if ( maxmemory_<memory ) maxmemory_=memory;
+}
+
+
+unsigned int summaryinfo::showmaxmem()
+{
+    acquiremaxmem();
+    return maxmemory_;
+}
+
+
+#ifndef __NO_AUTO_LOGGING
+const summaryinfo __summaryinfo;
+#endif
+
diff --git a/inputs/example/kgminkowski/analysis/src/intq.cc b/inputs/example/kgminkowski/analysis/src/intq.cc
new file mode 100644 (file)
index 0000000..0db2a54
--- /dev/null
@@ -0,0 +1,221 @@
+////////////////////////////////////////////////////////////////////////
+// intq.cc  Calculates intagrated (conserved) quantities for          //
+//          kgminkowski.                                              //
+////////////////////////////////////////////////////////////////////////
+
+
+#include "config.h"
+#include <gridripper/math.h>
+#include <gstream/gstream.h>
+#include <gridripper/multipole/ScalarFieldMP.h>
+#include <gridripper/io/BDataReader.h>
+#include <gridripper/phys/example/kgminkowski/KGMinkowski.h>
+#include <gridripper/math/O5Interpolation.h>
+
+
+using namespace std;
+using namespace gridripper::io;
+using namespace gridripper;
+using namespace gridripper::multipole;
+
+
+// O5 integration routine:
+GReal_t o5integrator(const tvector<GReal_t>& xarr, const tvector<GReal_t>& yarr)
+{
+    if ( xarr.size()<6 ) return 0.0;
+    gridripper::math::O5Interpolation o5int(xarr, yarr);
+    GReal_t sum=0.0;
+    for ( unsigned int i=0 ; i<xarr.size()-1 ; ++i )
+    {
+       sum+=o5int.integrateSection(i, 0);
+    }
+    return sum;
+}
+
+
+// The main program:
+int main (int argc, char *argv[])
+{
+
+    // Get the name of the configfile from the argument:
+    if ( argc!=2 ) { cerr<<"Usage:\tintq <configfile>\n"; exit(1); }
+    string configfile=argv[1];
+
+    // Set gstream debug level:
+    gstream_debug(1);
+
+    // Variable to store configuration file content:
+    config conf, confin;
+
+    // Read the content of configuration file:
+    conf.append(configfile);
+
+    // Get the input file name:
+    string inputFile=getconfig(conf, "inputFile");
+
+    // Append the input file parameters to configuration:
+    igpstream inputstream;
+    inputstream.open(inputFile, "grep \"^parameter\\..*=.*\" %f | tr '=' ' ' | tr '\t' ' ' | tr -s ' ' | cut -f -2 -d ' ' | cut -f 2- -d '.'");
+    confin.append(inputstream);
+    inputstream.close();
+    inputstream.clear();
+
+    // Get the bdata file names:
+    bool isongrid=false;
+    if ( inputFile.length()>=6 )
+    {
+       if ( inputFile.substr(0, 6)=="/grid/" ) isongrid=true;
+    }
+    vector<string> bdataFile;
+    string::size_type lastPerIndex=inputFile.rfind("/");
+    string bdataDir=(lastPerIndex!=string::npos ? inputFile.substr(0, lastPerIndex+1) : ".");
+    string bdataMask=(lastPerIndex!=string::npos ? inputFile.substr(lastPerIndex+1, inputFile.length()-(lastPerIndex+1)) : inputFile);
+    if ( bdataMask.length()>=3 )
+    {
+       if ( bdataMask.substr(bdataMask.length()-3, 3)==".in" ) bdataMask=bdataMask.substr(0, bdataMask.length()-3);
+    }
+    std::ipstream bdataList;
+    bdataList.open(string((isongrid ? "lfc-ls " : "ls "))+bdataDir+" | grep \""+bdataMask+".*\\.bdata\"");
+    string namebuff;
+    while ( getline(bdataList, namebuff) ) bdataFile.push_back(bdataDir+namebuff);
+    bdataList.close();
+    bdataList.clear();
+
+    // Get output file name for integrated quantities:
+    string intqFile=getconfig(conf, "intqFile");
+    // Get output file name for density quantities:
+    string densityFile=getconfig(conf, "densityFile");
+
+    // Get log file name:
+    string logFile=getconfig(conf, "logFile");
+    summaryinfo::logfilename(logFile);
+
+    // Get EPSILON:
+    GReal_t EPSILON=getconfig(conf, "EPSILON");
+
+    // Simaulation parameters:
+    GReal_t minr=0.0;
+    GReal_t minrIntq=minr;
+    GReal_t maxr=getconfig(confin, "maxr");
+    GReal_t maxrIntq=maxr;
+    int maxOrder=getconfig(confin, "maxOrder");
+    bool isTimeReversed=getconfig(confin, "isTimeReversed");
+
+    // Loop over time series:
+    ScalarFieldMP<GComplex_t> physf(maxOrder);
+    ScalarFieldMP<GComplex_t> physft(maxOrder);
+    ScalarFieldMP<GComplex_t> physfr(maxOrder);
+    tvector<GReal_t> tSequence;
+    tvector<GReal_t> EfluxInner, EfluxOuter;
+    tvector<GReal_t> LfluxInner, LfluxOuter;
+    ofstream intq;
+    intq.open(intqFile.c_str());
+    intq << "# " << "t\t" 
+          << "Etot\t" << "E("<<minrIntq<<","<<maxrIntq<<")\t" 
+          << "Eflow("<<minrIntq<<")\t" 
+          << "Eflow("<<maxrIntq<<")\t" 
+          << "Ltot\t" << "L("<<minrIntq<<","<<maxrIntq<<")\t" 
+          << "Lflow("<<minrIntq<<")\t" 
+          << "Lflow("<<maxrIntq<<")\t" 
+          <<endl;
+    intq.flush();
+    ofstream density;
+    density.open(densityFile.c_str());
+    density << "# " << "t\t" <<"r\t"<< "Refin\t" << "Edens\t" << "Ldens" <<endl;
+    for ( int i=0 ; i<(int)bdataFile.size() ; ++i )
+    {
+       BDataReader bdataReader(bdataFile[i]);
+       while ( bdataReader.read() )
+       {
+           GReal_t t=bdataReader.getTime();
+           int nc=bdataReader.getNumComponents();
+           if ( nc%2!=0 || (nc/2)%3!=0 || (maxOrder+1)*(maxOrder+1)*2*3!=nc ) { cerr<<"Inconsistent number of components!\n"; exit(1); }
+           GReal_t xMin=bdataReader.getXMin();
+           const tvalarray<int> &kXArray=bdataReader.getKxArray();
+           int maxi=(int)kXArray.size()-1;
+           GReal_t dx=bdataReader.getDeltaX();
+           int xMultiplier=bdataReader.getXMultiplier();
+           const GReal_t *theData=(const GReal_t*)bdataReader.getDoubleData();
+           tvector<GReal_t> totalR, partialR;
+           tvector<GReal_t> totalE, partialE;
+           tvector<GReal_t> totalL, partialL;
+           for ( int j=0 ; j<=maxi ; ++j )
+           {
+               GReal_t r=xMin+kXArray[j]*dx;
+               GReal_t refinementFactor=(GReal_t)xMultiplier/(GReal_t)(j<maxi ? kXArray[j+1]-kXArray[j] : kXArray[j]-kXArray[j-1]);
+               GReal_t *fieldAtX=(GReal_t*)(&theData[j*nc]);
+               int nc2=(maxOrder+1)*(maxOrder+1)*2;
+               ScalarFieldMP<GComplex_t> f(&fieldAtX[nc2*0], maxOrder);
+               ScalarFieldMP<GComplex_t> ft(&fieldAtX[nc2*1], maxOrder);
+               ScalarFieldMP<GComplex_t> fr(&fieldAtX[nc2*2], maxOrder);
+               physf=f/r;
+               physft=ft/r*(isTimeReversed ? -1.0 : 1.0);
+               physfr=fr-f/(r*r);
+               totalR.push_back(r);
+               if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialR.push_back(r);
+               // Calculate energy density:
+               GReal_t Edensity=(r<minr+EPSILON ? 0.0 :
+                                       real(0.5*(r*r)*L2S2(physft, physft)+0.5*(r*r)*L2S2(physfr, physfr)-0.5*L2S2(physf, LaplaceS2(physf)))
+                                );
+               totalE.push_back(Edensity);
+               if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialE.push_back(Edensity);
+               // Calculate angular momentum density:
+               GReal_t Ldensity=(r<minr+EPSILON ? 0.0 :
+                                       real((r*r)*L2S2(physft, partial_phi(physf)))
+                                );
+               totalL.push_back(Ldensity);
+               if ( minrIntq-EPSILON<r && r<maxrIntq+EPSILON ) partialL.push_back(Ldensity);
+               density << t << "\t" << r << "\t" << refinementFactor << "\t" << Edensity << "\t" << Ldensity <<endl;
+               if ( ::fabs(r-minrIntq)<EPSILON || ::fabs(r-maxrIntq)<EPSILON )
+               {
+                   if ( tSequence.size()==0 ) tSequence.push_back(t);
+                   else if ( ::fabs(tSequence[tSequence.size()-1]-t)>EPSILON ) tSequence.push_back(t);
+                   // Calculate energy flow:
+//                 GReal_t Eflow=(r<minr+EPSILON ? 0.0 :
+//                                     -Delta*real(L2S2(physfr, physft)
+//                                     -GI*q*e*r*L2S2(physf, divideByNeumann(scaledSigma, physfr, normBoundNeumann(scaledSigma(0, 0), H2S2SigmaSigma/norm(scale)), neumannTolerance))/scale)
+//                               );
+                   GReal_t Eflow=0.0;
+                   // Calculate angular momentum flow:
+//                 GReal_t Lflow=(r<minr+EPSILON ? 0.0 :
+//                                     -real(Delta*L2S2(physfr, partial_phi(physf))+GI*q*a*e*r*L2S2(physf, divideByNeumann(scaledSigma, mulSinThetaSqr(physfr), normBoundNeumann(scaledSigma(0, 0), H2S2SigmaSigma/norm(scale)), neumannTolerance)/scale))
+//                               );
+                   GReal_t Lflow=0.0;
+                   if ( ::fabs(r-minrIntq)<EPSILON )
+                   {
+                       EfluxInner.push_back(Eflow);
+                       LfluxInner.push_back(Lflow);
+                   }
+                   else
+                   {
+                       EfluxOuter.push_back(Eflow);
+                       LfluxOuter.push_back(Lflow);
+                   }
+               }
+           }
+           density << endl;
+           density.flush();
+           intq << t << "\t" 
+                  << o5integrator(totalR, totalE) << "\t" << o5integrator(partialR, partialE) << "\t" 
+                  << o5integrator(tSequence, EfluxInner) << "\t"
+                  << o5integrator(tSequence, EfluxOuter) << "\t"
+                  << o5integrator(totalR, totalL) << "\t" << o5integrator(partialR, partialL) << "\t" 
+                  << o5integrator(tSequence, LfluxInner) << "\t"
+                  << o5integrator(tSequence, LfluxOuter) << "\t"
+                  << endl;
+           intq.flush();
+       }
+    }
+    density.close();
+    density.clear();
+    intq.close();
+    intq.clear();
+
+    return 0;
+
+}
+
+
+/////////////
+// intq.cc //
+/////////////
diff --git a/inputs/example/kgminkowski/kgminkowski.in b/inputs/example/kgminkowski/kgminkowski.in
new file mode 100644 (file)
index 0000000..eab3a3c
--- /dev/null
@@ -0,0 +1,59 @@
+#################### Configuration file for kgminkowski ################
+
+### PDE and initial conditions
+pde=example.kgminkowski.KGMinkowski
+initCond=example.kgminkowski.Hunch
+#initCond=BData /usr/people/alaszlo/gridripperprod/inputs/example/kgminkowski/kgminkowski.009.bdata t0=0.0
+parameters=t0,amplitudeF,centerF,widthF,rampF,lF:%d,mF:%d,amplitudeFt,centerFt,widthFt,rampFt,lFt:%d,mFt:%d,omega,maxOrder:%d,maxr,isTimeReversed:%d
+parameter.t0=0.0
+parameter.amplitudeF=1.0
+parameter.centerF=31.9832359852999
+parameter.widthF=34.3219506001479
+parameter.rampF=34.3219506001479
+parameter.lF=2
+parameter.mF=2
+parameter.amplitudeFt=0.0
+parameter.centerFt=31.9832359852999
+parameter.widthFt=34.3219506001479
+parameter.rampFt=34.3219506001479
+parameter.lFt=2
+parameter.mFt=2
+parameter.omega=0.433804363739061
+parameter.maxOrder=12
+parameter.maxr=64.0
+parameter.isTimeReversed=0
+
+### Numerical method and AMR parameters
+gridInt=RK4 d=O421
+dtdx=0.125
+sigma=0.01
+maxi=512
+maxLevel=4
+refinementRatio=2
+amError=NormalizedL2ComponentError offset=constraint 0
+errorTolerance=1e-6
+errorCheckFreq=8
+regridFreq=16
+bufferZoneSize=2
+
+### Output file
+fileNo=0
+maxBDataSize=1G
+dtwrite=0.25
+tmax=t0+192.0
+maxRuntime=864000
+blankLines=2
+dataFormat=binary
+energyWritten=false
+writeIntq=false
+writeMarker=false
+
+### GUI only
+delay=0
+displayedComponents=0
+scale=5000
+memorySteps=2
+memoryMax=10
+
+### Debug
+#debug=true
index 947e0f7..f2654dd 100644 (file)
@@ -6,6 +6,8 @@ gridripperapp \- This GridRipperApp PDE solver application library.
 .SH SYNOPSIS
 To use as library, include the relevant library header files, like:
 
+.B #include <gridripper/io/BDataReader.h>
+
 .B #include <gridripper/phys/Hunch.h>
 
 Compile with the following flags:
@@ -16,43 +18,43 @@ Link with the following flags:
 
 .B `pkg-config gridripperapp --libs` `pkg-config gridripperapp --libs`
 
-(Please do apply the indicated repeating in the above linking option: 
+(Please do apply the indicated repeating in the above linking option:
 it is not a typo. It is needed to resolve circular dependency within.)
 
 To use the PDE solver application, use the following command:
 
 .B gridripper your_inputfile
 
-The produced binary output files are stored under name scheme *.bdata, 
-and can be converted to an easy-to-visualize ASCII format using the 
+The produced binary output files are stored under name scheme *.bdata,
+and can be converted to an easy-to-visualize ASCII format using the
 command line tool:
 
 .B gettslice
 
-In order to get help on the usage of this conversion tool, just start 
+In order to get help on the usage of this conversion tool, just start
 it on the command line without arguments.
 
 .SH DESCRIPTION
-The GridRipperApp is an application library for the GridRipper PDE 
-solver core library, written in C++ language. 
-The PDE solver is capable of solving 1+d dimensional problems with 
+The GridRipperApp is an application library for the GridRipper PDE
+solver core library, written in C++ language.
+The PDE solver is capable of solving 1+d dimensional problems with
 finite difference AMR + spectral methods.
 
-This library contains implementations of the partial differential 
-equations (PDE-s) passed to the GridRipper numerical core library. 
-Users may extend this collection library of PDE-s by defining new 
+This library contains implementations of the partial differential
+equations (PDE-s) passed to the GridRipper numerical core library.
+Users may extend this collection library of PDE-s by defining new
 ones.
 
-When this application library is compiled, it gets automatically linked 
-to its core library, generating the binary 
+When this application library is compiled, it gets automatically linked
+to its core library, generating the binary
 .B gridripper
 which is the actual PDE solver program.
 
 The tool
 .B gettslice
-is shipped with the gridripper+gridripperapp package and can be used to 
-pick certain time slices of the output binary files, called .bdata 
-files. This tool converts the time slices to an easy-to-visualize ASCII 
+is shipped with the gridripper+gridripperapp package and can be used to
+pick certain time slices of the output binary files, called .bdata
+files. This tool converts the time slices to an easy-to-visualize ASCII
 format.
 
 .SH REQUIREMENTS
@@ -61,25 +63,42 @@ GridRipperApp requires the GridRipper core library to be installed.
 .SH CONFIGURATION
 The PDE solver program
 .B gridripper
-must be started with an input file as its argument. The input file 
-contains the settings for the applied numerical methods and for the 
-actual PDE plus its initial conditions. This is described in detail 
-on the GridRipper homepage (see below). One can find example input files 
+must be started with an input file as its argument. The input file
+contains the settings for the applied numerical methods and for the
+actual PDE plus its initial conditions. This is described in detail
+on the GridRipper homepage (see below). One can find example input files
 under the directory reported by the following command:
 
+.B pkg-config gridripperapp --variable=inputsdir
+
+There is a particular educative working example input file under
+
+.B gridripperapp/inputs/example/kgminkowski
+
+with corresponding header and source code files under
+
+.B gridripperapp/include/gridripper/phys/example/kgminkowski
+
+and
+
+.B gridripperapp/src/phys/example/kgminkowski
+
+access path. This example solves the wave equation in 3+1 dimensional
+Minkowski spacetime for fields expanded with respect to spherical
+harmonics.
+
 The binary-to-ASCII converter program
 .B gettslice
-must be started with proper command line arguments. For a help on 
-these, just start it on the command line without arguments.
-
-.B pkg-config gridripperapp --variable=inputsdir
+which is shipped with gridridripper+gridripperapp, must be started
+with proper command line arguments. For a help on these, just start
+it on the command line without arguments.
 
 .SH ENVIRONMENT
 Make sure that the following environmental variables are set correctly:
 
 .B PATH, LD_LIBRARY_PATH, PKG_CONFIG_PATH, MANPATH
 
-Typical setting in your ~/.bashrc file (for bash shell) in case of an 
+Typical setting in your ~/.bashrc file (for bash shell) in case of an
 installation under /usr/local directory:
 
 .B export PATH=/usr/local/bin:$PATH
@@ -90,7 +109,7 @@ installation under /usr/local directory:
 
 .B export MANPATH=/usr/local/share/man:$MANPATH
 
-Typical setting in your ~/.tcshrc file (for tcsh shell) in case of an 
+Typical setting in your ~/.tcshrc file (for tcsh shell) in case of an
 installation under /usr/local directory:
 
 .B setenv PATH /usr/local/bin:$PATH
@@ -101,7 +120,7 @@ installation under /usr/local directory:
 
 .B setenv MANPATH /usr/local/share/man:$MANPATH
 
-Note that on some flavors of 64 bit operating systems the paths of 
+Note that on some flavors of 64 bit operating systems the paths of
 the form
 .B .../lib/...
 rather takes the form of
diff --git a/src/phys/example/kgminkowski/Factory.cxx b/src/phys/example/kgminkowski/Factory.cxx
new file mode 100644 (file)
index 0000000..c83f6b6
--- /dev/null
@@ -0,0 +1,83 @@
+#include <gridripper/phys/example/kgminkowski/Factory.h>
+#include <gridripper/phys/example/kgminkowski/Hunch.h>
+#include <gridripper/phys/example/kgminkowski/KGMinkowski.h>
+#include <sstream>
+#include <gridripper/Parameters.h>
+
+#include <gridripper/factory_ns.h>
+
+
+//////////////////// implementation of Factories ///////////////////////
+
+ParameterPreprocessor*
+gridripper::phys::example::kgminkowski::createPP(const string& name,
+           const string& arg, const Parameters& parameters)
+       throw(IllegalArgumentException&)
+{
+    ostringstream oss;
+    oss << "cannot create unknown ParameterPreprocessor \"" << name << "\"";
+    throw IllegalArgumentException(oss.str(),
+           "gridripper::phys::example::kgminkowski::createPP");
+}
+
+InitCond* gridripper::phys::example::kgminkowski::createInitCond(
+               const string& name,
+               string& arg, const Parameters* parameters, const PDE& pde)
+       throw(InitCond::Exception&, IllegalArgumentException&)
+{
+    const string location="gridripper::phys::example::kgminkowski::createInitCond(const string&, const string&, const Parameters*, const PDE&)";
+    try
+    {
+       if ( name=="example.kgminkowski.Hunch" )
+       {
+           return new Hunch(arg, parameters, pde);
+       }
+    }
+    catch(IllegalArgumentException& ex)
+    {
+       ex.addLocation(location);
+        throw ex;
+    }
+    ostringstream oss;
+    oss << "cannot create unknown initial condition \"" << name << "\"";
+    throw InitCond::Exception(oss.str(), location);
+}
+
+
+PDE* gridripper::phys::example::kgminkowski::createPDE(const string& name,
+                                                 const Parameters* parameters)
+       throw(IllegalArgumentException&)
+{
+    const string location = "gridripper::phys::example::kgminkowski::createPDE";
+    try
+    {
+       if ( name=="example.kgminkowski.KGMinkowski" )
+       {
+           if ( parameters==NULL ) throw(IllegalArgumentException("No parameters specified.", location));
+           int maxorder=parameters->getInt("maxOrder");
+           if ( maxorder<2 ) throw(IllegalArgumentException("Necessary: maxOrder>=2 !!!", location));
+           int arraysize=(maxorder+1)*(maxorder+1)*sizeof(GComplex_t)/sizeof(GReal_t);
+           return new KGMinkowski(parameters, maxorder, arraysize);
+       }
+    }
+    catch(IllegalArgumentException& ex)
+    {
+       ex.addLocation(location);
+       throw ex;
+    }
+    ostringstream oss;
+    oss << "cannot create unknown PDE \"" << name << "\"";
+    throw IllegalArgumentException(oss.str(), location);
+}
+
+
+//ODE* gridripper::phys::gr::fixmp::kerrhiggs::createODE(const string& name, 
+//                                               const Parameters* parameters)
+//     throw(IllegalArgumentException&)
+//{
+//    ostringstream oss;
+//    oss << "cannot create unknown ODE \"" << name << "\"";
+//    throw IllegalArgumentException(oss.str(),
+//         "gridripper::phys::gr::fixmp::kerrhiggs::createODE(const string&, const Parameters*");
+//}
+
diff --git a/src/phys/example/kgminkowski/Hunch.cxx b/src/phys/example/kgminkowski/Hunch.cxx
new file mode 100644 (file)
index 0000000..2cf95e7
--- /dev/null
@@ -0,0 +1,125 @@
+#include <gridripper/phys/example/kgminkowski/Hunch.h>
+#include <gridripper/phys/example/kgminkowski/KGMinkowski.h>
+#include <gridripper/amr1d/PDE.h>
+#include <gridripper/Parameters.h>
+#include <gridripper/math.h>
+
+using namespace gridripper;
+using namespace std;
+
+
+//////////////////// class Hunch ///////////////////////////////////////
+
+
+using namespace gridripper::phys::example::kgminkowski;
+
+
+Hunch::Hunch(string& args, const Parameters* p, const PDE& pde)
+            throw(InitCond::Exception&, IllegalArgumentException&):
+                FuncInitCond(args, p)
+{
+    const string location = "gridripper::phys::example::kgminkowski::Hunch::Hunch(const string&, const Parameters*, const PDE&)";
+    if(p != NULL)
+    {
+       try
+       {
+           t0=p->getReal("t0");
+           minr=0.0;
+           maxr=p->getReal("maxr");
+           maxOrder=p->getInt("maxOrder");
+           arraySize=(maxOrder+1)*(maxOrder+1)*sizeof(GComplex_t)/sizeof(GReal_t);
+           amplitudeF=p->getReal("amplitudeF");
+           centerF=p->getReal("centerF");
+           widthF=p->getReal("widthF");
+           rampF=p->getReal("rampF");
+           lF=p->getInt("lF");
+           mF=p->getInt("mF");
+           amplitudeFt=p->getReal("amplitudeFt");
+           centerFt=p->getReal("centerFt");
+           widthFt=p->getReal("widthFt");
+           rampFt=p->getReal("rampFt");
+           lFt=p->getInt("lFt");
+           mFt=p->getInt("mFt");
+           omega=p->getReal("omega");
+       }
+       catch(IllegalArgumentException& ex)
+       {
+           ex.addLocation(location);
+           throw ex;
+       }
+    }
+    else
+    {
+       throw IllegalArgumentException("No parameters specified!", location);
+    }
+    message="Hunch";
+}
+
+
+Hunch::~Hunch()
+{
+}
+
+
+GReal_t Hunch::getTime() const
+{
+    return t0;
+}
+
+
+string Hunch::getMessage() const
+{
+    return message;
+}
+
+
+// Compactly supported wave packet over real numbers.
+static GReal_t packet(const GReal_t x, const GReal_t center, const GReal_t width, const GReal_t ramp)
+{
+    return (center-0.5*width<x && x<center+0.5*width ? ::exp(-::fabs(ramp/(x-center+0.5*width))-::fabs(ramp/(x-center-0.5*width))+4.0*::fabs(ramp/width)) : 0.0);
+}
+
+
+// Derivative of compactly supported wave packet.
+static GReal_t packetprime(const GReal_t x, const GReal_t center, const GReal_t width, const GReal_t ramp)
+{
+    return (center-0.5*width<x && x<center+0.5*width ? (ramp/square(x-center+0.5*width)-ramp/square(x-center-0.5*width))*::exp(-::fabs(ramp/(x-center+0.5*width))-::fabs(ramp/(x-center-0.5*width))+4.0*::fabs(ramp/width)) : 0.0);
+}
+
+
+void Hunch::init(PDE* pde, GReal_t r, FieldWrapper& fw)
+    throw(InitCond::Exception&)
+{
+    // Convert to multipole expansion.
+    const GReal_t *vtmp=fw.data();
+    ScalarFieldMP<GComplex_t> f((GReal_t*)&vtmp[KGMinkowski::ind_f*arraySize], maxOrder);
+    ScalarFieldMP<GComplex_t> ft((GReal_t*)&vtmp[KGMinkowski::ind_ft*arraySize], maxOrder);
+    // Fill the spherical wave initial conditions for the physical field.
+    for ( int l=0 ; l<=maxOrder ; ++l )
+    {
+       for ( int m=-l ; m<=l ; ++m )
+       {
+           f(l, m)=(l==lF && m==mF ? amplitudeF*packet(r+t0, centerF, widthF, rampF)*exp(-GI*omega*(r+t0-centerF)) : 0.0);
+           ft(l, m)=((l==lFt && m==mFt ? amplitudeFt*packet(r+t0, centerFt, widthFt, rampFt)*exp(-GI*omega*(r+t0-centerFt)) : 0.0) + (l==lF && m==mF ? amplitudeF*(packetprime(r+t0, centerF, widthF, rampF)*exp(-GI*omega*(r+t0-centerF))-GI*omega*packet(r+t0, centerF, widthF, rampF)*exp(-GI*omega*(r+t0-centerF))) : 0.0));
+       }
+    }
+    // Transform physical field to field evolved on grid.
+    if ( minr+KGMinkowski::EPSILON<r && r<maxr-KGMinkowski::EPSILON )
+    {
+       if ( KGMinkowski::EPSILON<r && r<1.0/KGMinkowski::EPSILON )
+       {
+           f*=r;
+           ft*=r;
+       }
+       else
+       {
+           f.annulate();
+           ft.annulate();
+       }
+    }
+    else
+    {
+       f.annulate();
+       ft.annulate();
+    }
+}
diff --git a/src/phys/example/kgminkowski/KGMinkowski.cxx b/src/phys/example/kgminkowski/KGMinkowski.cxx
new file mode 100644 (file)
index 0000000..05ef983
--- /dev/null
@@ -0,0 +1,327 @@
+#include <gridripper/phys/example/kgminkowski/KGMinkowski.h>
+#include <gridripper/amr1d/Grid.h>
+#include <gridripper/amr1d/Grad.h>
+#include <gridripper/amr1d/Integrator.h>
+#include <gridripper/Parameters.h>
+#include <gsl/gsl_sf.h>
+#include <gridripper/math.h>
+#include <iostream>
+#include <time.h>
+
+
+using namespace gridripper;
+using namespace gridripper::multipole;
+using namespace std;
+
+
+//////////////////// class KGMinkowski /////////////////////////////////
+
+
+using namespace gridripper::phys::example::kgminkowski;
+
+// Error tolerance.
+const GReal_t KGMinkowski::EPSILON=1e-8;
+
+// Label arrays.
+const string KGMinkowski::COMP_NAMES[NUM_COMPS]={"f", "f_t", "f_r"};
+const string KGMinkowski::CONSTRAINT_NAMES[NUM_CONSTRAINTS]={"df/dr"};
+const string KGMinkowski::COORD_LABELS[2]={"t", "r"};
+
+
+// A helper template function to resize a tvector.
+template <typename T>
+static tvector<T>& resize(tvector<T> &v, const unsigned int newsize)
+{
+    v.resize(newsize);
+    return v;
+}
+
+
+tvector<string> KGMinkowski::componentNames(int maxOrder)
+{
+    tvector<string> compNames;
+    for ( int l=0 ; l<=maxOrder ; ++l )
+    {
+       for ( int m=-l ; m<=l ; ++m )
+       {
+           ostringstream oss;
+           oss<<"("<<l<<","<<m<<")";
+           for ( int comp=0 ; comp<NUM_COMPS ; ++comp )
+           {
+               compNames.push_back(COMP_NAMES[comp]+oss.str()+"re");
+               compNames.push_back(COMP_NAMES[comp]+oss.str()+"im");
+           }
+       }
+    }
+    return compNames;
+}
+
+tvector<string> KGMinkowski::constraintNames(int maxOrder)
+{
+    tvector<string> constrNames;
+    for ( int l=0 ; l<=maxOrder ; ++l )
+    {
+       for ( int m=-l ; m<=l ; ++m )
+       {
+           ostringstream oss;
+           oss<<"("<<l<<","<<m<<")";
+           for ( int cons=0 ; cons<NUM_CONSTRAINTS ; ++cons )
+           {
+               constrNames.push_back(CONSTRAINT_NAMES[cons]+oss.str()+"re");
+               constrNames.push_back(CONSTRAINT_NAMES[cons]+oss.str()+"im");
+           }
+       }
+    }
+    return constrNames;
+}
+
+
+// Initialize the PDE.
+KGMinkowski::KGMinkowski(const Parameters* p, const int maxorder, const int arraysize) throw(IllegalArgumentException&):
+    PDE(NUM_COMPS*arraysize, NUM_INDEP_COMPS*arraysize, 
+       componentNames(maxorder), constraintNames(maxOrder), COORD_LABELS),
+    maxOrder(maxorder), arraySize(arraysize),
+    I_f(ind_f*arraySize), I_ft(ind_ft*arraySize), I_fr(ind_f*arraySize),
+    df(NUM_COMPS*arraySize), fieldat(NUM_AUX_POINTS)
+{
+    const string location="gridripper::phys::example::kgminkowski::KGMinkowski::KGMinkowski(const Parameters*,.....)";
+    for ( int i=0 ; i<NUM_AUX_POINTS ; ++i ) fieldat[i].resize(NUM_COMPS*arraySize);
+    // Set the parameter values.
+    if ( p!=NULL )
+    {
+       try
+       {
+           t0=p->getReal("t0");
+           minr=0.0;
+           maxr=p->getReal("maxr");
+           if ( minr<0.0 || maxr<0.0 || maxr<=minr )
+               throw IllegalArgumentException("r domain should be non-negative for Minkowski!", location);
+           isTimeReversed=p->getInt("isTimeReversed");
+       }
+       catch(IllegalArgumentException& ex)
+       {
+           ex.addLocation(location);
+           throw ex;
+       }
+    }
+    else
+    {
+       throw IllegalArgumentException("No parameters specified!", location);
+    }
+    exclude=0x00000000;
+    if ( isTimeReversed ) { cout<<"Time is reversed.\n"; cout.flush(); }
+}
+
+
+KGMinkowski::~KGMinkowski()
+{
+}
+
+
+GReal_t KGMinkowski::getMinX() const
+{
+    return minr;
+}
+
+
+GReal_t KGMinkowski::getMaxX() const
+{
+    return maxr;
+}
+
+
+GReal_t KGMinkowski::getPhysicalMinX() const
+{
+    return minr;
+}
+
+
+GReal_t KGMinkowski::getPhysicalMaxX() const
+{
+    return maxr;
+}
+
+
+bool KGMinkowski::isXPeriodic() const
+{
+    return false;
+}
+
+
+// Evaluate the PDE.
+int KGMinkowski::eval(GReal_t* dF, GReal_t t, GReal_t r,
+                   const GReal_t* F, const GReal_t* Fr,
+                   const Grid& G, int i, Grad& d, GReal_t dt)
+{
+    // Interpret raw fields as multipole expanded fields.
+    ScalarFieldMP<GComplex_t> df((GReal_t*)&dF[I_f], maxOrder);
+    ScalarFieldMP<GComplex_t> dft((GReal_t*)&dF[I_ft], maxOrder);
+    ScalarFieldMP<GComplex_t> dfr((GReal_t*)&dF[I_fr], maxOrder);
+    ScalarFieldMP<GComplex_t> f((GReal_t*)&F[I_f], maxOrder);
+    ScalarFieldMP<GComplex_t> ft((GReal_t*)&F[I_ft], maxOrder);
+    ScalarFieldMP<GComplex_t> fr((GReal_t*)&F[I_fr], maxOrder);
+    ScalarFieldMP<GComplex_t> pf((GReal_t*)&Fr[I_f], maxOrder);
+    ScalarFieldMP<GComplex_t> pft((GReal_t*)&Fr[I_ft], maxOrder);
+    ScalarFieldMP<GComplex_t> pfr((GReal_t*)&Fr[I_fr], maxOrder);
+
+    // If we are at the origin, treat it specially.
+    if ( r<getMinX()+KGMinkowski::EPSILON )
+    {
+       // Field is constant zero at the origin in Minkowski case.
+       // Also the radial derivative is zero at the origin, 
+       // except for l==0 multipole component, which properly evolves.
+       df.annulate();
+       dft.annulate();
+       dfr.annulate();
+       dfr(0, 0)=pft(0, 0);
+    }
+    // If it is an ordinary point, treat it normally.
+    else
+    {
+       // The field equation.
+       df=(isTimeReversed ? -1.0 : 1.0)*ft;
+       dft=(isTimeReversed ? -1.0 : 1.0)*(pfr+(1.0/(r*r))*LaplaceS2(f));
+       dfr=(isTimeReversed ? -1.0 : 1.0)*pft;
+    }
+    return exclude;
+}
+
+
+// Evaluate the conserved constraints initially.
+void KGMinkowski::evalConstrainedComponents(Grid& G, int i, FieldWrapper& fw)
+{
+    GReal_t* v = fw.data();
+    Grad& d=G.getIntegrator()->getGrad();
+    d.d(G, this, i, df);
+    int offset=getNumIndependentComponents();
+    for ( int j=0 ; j<arraySize ; ++j ) v[offset+j]=df[I_f+j];
+}
+
+
+// Set nonevolving components.
+void KGMinkowski::setNonevolvingComponents(Grid& G, Grad& d)
+{
+    // Get the range of the grid and its step size.
+    GReal_t minrGrid=G.getLocationX(G.getLeftmostI());
+    GReal_t maxrGrid=G.getLocationX(G.getRightmostI());
+    GReal_t Dr=G.getDeltaX();
+    GReal_t r, ratio, difference;
+    // If the grid contains the origin, treat it specially.
+    r=0.0;
+    ratio=(r-minrGrid)/Dr;
+    difference=fabs(ratio-round(ratio));
+    if ( minrGrid-KGMinkowski::EPSILON<r && r<maxrGrid+KGMinkowski::EPSILON && difference<KGMinkowski::EPSILON )
+    {
+       // Get the index of the origin.
+       int i=(int)round(ratio);
+       // Get the field at the origin.
+       G.getField(i, fieldat[0]);
+       ScalarFieldMP<GComplex_t> f0((GReal_t*)&fieldat[0].data()[I_f], maxOrder);
+       ScalarFieldMP<GComplex_t> f0t((GReal_t*)&fieldat[0].data()[I_ft], maxOrder);
+       ScalarFieldMP<GComplex_t> f0r((GReal_t*)&fieldat[0].data()[I_fr], maxOrder);
+       // The field is exactly zero at the origin, and also its 
+       // radial derivative, except for l==0 multipole component, 
+       // which evolves properly.
+       f0.annulate();
+       f0t.annulate();
+       GComplex_t value=f0r(0, 0);
+       f0r.annulate();
+       f0r(0, 0)=value;
+       G.set(i, fieldat[0].data());
+    }
+    // If the grid contains the point next to origin, treat it specially.
+    r=0.0+Dr;
+    ratio=(r-minrGrid)/Dr;
+    difference=fabs(ratio-round(ratio));
+    if ( minrGrid-KGMinkowski::EPSILON<r && r<maxrGrid+KGMinkowski::EPSILON && difference<KGMinkowski::EPSILON )
+    {
+       // Get the index of the point next to origin.
+       int i=(int)round(ratio);
+       // Get the field at the point next to the origin.
+       G.getField(i, fieldat[1]);
+       ScalarFieldMP<GComplex_t> f1((GReal_t*)&fieldat[1].data()[I_f], maxOrder);
+       ScalarFieldMP<GComplex_t> f1t((GReal_t*)&fieldat[1].data()[I_ft], maxOrder);
+       ScalarFieldMP<GComplex_t> f1r((GReal_t*)&fieldat[1].data()[I_fr], maxOrder);
+       // Get the field at the point next to next to the origin.
+       G.getField(i+1, fieldat[2]);
+       ScalarFieldMP<GComplex_t> f2((GReal_t*)&fieldat[2].data()[I_f], maxOrder);
+       ScalarFieldMP<GComplex_t> f2t((GReal_t*)&fieldat[2].data()[I_ft], maxOrder);
+       ScalarFieldMP<GComplex_t> f2r((GReal_t*)&fieldat[2].data()[I_fr], maxOrder);
+       // Get the field at the point next to next to next to the origin.
+       G.getField(i+2, fieldat[3]);
+       ScalarFieldMP<GComplex_t> f3((GReal_t*)&fieldat[3].data()[I_f], maxOrder);
+       ScalarFieldMP<GComplex_t> f3t((GReal_t*)&fieldat[3].data()[I_ft], maxOrder);
+       ScalarFieldMP<GComplex_t> f3r((GReal_t*)&fieldat[3].data()[I_fr], maxOrder);
+       // The field has to be 2 times continuously differentiable at 
+       // the origin. This implies the following things.
+       for ( int l=0 ; l<=maxOrder ; ++l )
+       {
+           if ( l==0 ) continue;
+           for ( int m=-l ; m<=l ; ++m )
+           {
+               f1(l, m)=(l%2==0 ? 1.0/8.0 : 1.0/16.0)*f2(l, m); // <-This was necessary and sufficient for unigrid stability.
+               f1t(l, m)=(l%2==0 ? 1.0/8.0 : 1.0/16.0)*f2t(l, m);
+               f1r(l, m)=(1.0/Dr)*(l%2==0 ? (21.0/32.0)*f2(l, m)-(1.0/12.0)*f3(l, m) : (43.0/64.0)*f2(l, m)-(1.0/12.0)*f3(l, m));
+           }
+       }
+       G.set(i, fieldat[1].data());
+    }
+}
+
+
+FieldWrapper* KGMinkowski::createField(int level) const
+{
+    return new FieldComponents<GReal_t>(getNumComponents());
+}
+
+
+bool KGMinkowski::canBeExtended(int where) const
+{
+    return (where==LEFT);
+}
+
+
+// Extend field to the left to the nonphysical region.
+void KGMinkowski::mirrorLeft(GReal_t r, FieldWrapper& f) const
+{
+    // Convert fields to multipole expanded fields.
+    ScalarFieldMP<GComplex_t> v_f((GReal_t*)&f.data()[I_f], maxOrder);
+    ScalarFieldMP<GComplex_t> v_ft((GReal_t*)&f.data()[I_ft], maxOrder);
+    ScalarFieldMP<GComplex_t> v_fr((GReal_t*)&f.data()[I_fr], maxOrder);
+    v_f   = -reflect(v_f);
+    v_ft  = -reflect(v_ft);
+    v_fr = reflect(v_fr);
+}
+
+
+// Extend field to the right to the nonphysical region.
+void KGMinkowski::mirrorRight(GReal_t r, FieldWrapper& f) const
+{
+    cerr<<"KGMinkowski::mirrorRight should not be called.\n"; exit(1);
+}
+
+
+bool KGMinkowski::hasExtendedRefinedFieldData(const Grid& G, int i) const
+{
+    // Left bound:
+    if ( i<0 && G.containsLeftBound() )
+    {
+       int r=G.getRefinementRatio();
+       GReal_t dx=G.getDeltaX()/r;
+       GReal_t x=G.getLocationX(0)+i*dx;
+       GReal_t x_mirr=2.0*getMinX()-x;
+       GReal_t i_mirr=(x_mirr-G.getLocationX(0))/dx;
+       return G.hasRefinedData((int)(round(i_mirr)));
+    }
+    // Right bound:
+    else if ( i>G.getRefinedMaxI() && G.containsRightBound() )
+    {
+       return false;
+    }
+    // Otherwise:
+    else
+    {
+       return G.hasRefinedData(i);
+    }
+}
+