This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@3 1c20744c-e24a-42ec-9533-f5004cb800e5

12 files changed, 2511 insertions, 0 deletions

diff --git a/COPYRIGHT b/COPYRIGHT
new file mode 100644
index 0000000..57b2896
--- /dev/null
+++ b/COPYRIGHT
@@ -0,0 +1,352 @@
+
+
+This thorn is (c) Copyright the authors listed in the sources files. This
+thorn is distributed under the GNU GPL with the specific exception that
+the GNU GPL does not migrate to code linking to or using infrastructure
+from this thorn.
+
+- The Cactus Team <cactus@cactuscode.org>
+
+**************************************************************************
+
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                          675 Mass Ave, Cambridge, MA 02139, USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	Appendix: How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) 19yy  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) 19yy name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
+
+
diff --git a/README b/README
new file mode 100644
index 0000000..b2b21d1
--- /dev/null
+++ b/README
@@ -0,0 +1,18 @@
+Cactus Code Thorn PUGHInterp
+Authors    : Paul Walker and Thomas Radke, some optimization ideas from Erik Schnetter
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+   This thorn provides interpolation of arrays at arbitrary points
+   on a regular, uniform cartesian grid.
+
+
+2. Additional Information
+
+   More detailed information on thorn PUGHInterp is given in the
+   doc/documentation.tex file.
+
+   For information on how to invoke interpolation operators via the
+   general flesh interpolation API please refer to the flesh documentation.
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..b59a7a3
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,66 @@
+%version   $Header$
+\documentclass{article}
+\begin{document}
+
+\title{PUGHInterp}
+\author{Paul Walker, Thomas Radke, Erik Schnetter}
+\date{1997-2001}
+\maketitle
+
+\abstract{Thorn PUGHInterp provides interpolation of arrays
+          at arbitrary points.}
+
+\section{Purpose}
+Thorn PUGHInterp implements interpolation operators working on regular,
+uniform cartesian grids. They can be applied to interpolate both CCTK grid
+variables (grid functions and arrays distributed over all processors), and
+processor-local arrays at arbitrary points.\\
+The operators for distributed/local arrays are invoked via the flesh
+interpolation API routines {\tt CCTK\_InterpGV()} and {\tt CCTK\_InterpLocal()}
+resp. They are registered with the flesh under the name {\it ''uniform
+cartesian''} prepended by the interpolation order (eg. {\it ''second-order
+uniform cartesian''}). Currently there is first, second, and third-order
+interpolation implemented.\\
+
+\section{Implementation Notes}
+The interpolation operators registered for different orders are mapped
+via wrappers (in {\tt Startup.c}) onto a single routine (in {\tt Operator.c})
+just passing the order as an additional argument.\\
+The routine for distributed arrays will then map all points to interpolate
+at to the processors which own those points, and communicate the points'
+coordinates and corresponding input arrays ({\tt MPI\_Alltoall()} is used
+for this).\\
+Then the interpolation takes place in parallel on every processor, calling
+a core interpolation routine (located in {\tt Interpolate.c}). This one
+takes a list of input arrays and points and interpolates these to a
+list of output arrays (one output value per interpolation point).\\
+Again, for distributed arrays, the interpolation results for remote points
+are sent back to the requesting processors.\\[2ex]
+%
+Current limitations of the core interpolation routine's implementation are:
+%
+\begin{itemize}
+  \item arrays up to three ({\bf MAXDIM}) dimensions only can be handled
+  \item interpolation orders up to three ({\bf MAXORDER}) only are supported
+  \item coordinates must be given as {\bf CCTK\_REAL} types
+  \item input and output array types must be the same
+        (no type casting of interpolation results supported)
+\end{itemize}
+%
+Despite of these limitations, the code it was programmed almost generic
+in that it can easily be extended to support higher-dimensional arrays
+or more interpolation orders.\\
+Please see the NOTES in this source file for details.
+
+\section{Comments}
+%
+For more information on how to invoke interpolation operators please refer
+to the flesh documentation.
+%
+% Automatically created from the ccl files 
+% Do not worry for now.
+\include{interface}
+\include{param}
+\include{schedule}
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..565e2d8
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,4 @@
+# Interface definition for thorn PUGHInterp
+# $Header$
+
+implements: PUGHInterp
diff --git a/par/interp.par b/par/interp.par
new file mode 100644
index 0000000..f1ff599
--- /dev/null
+++ b/par/interp.par
@@ -0,0 +1,6 @@
+# interp.par
+#
+
+ActiveThorns = "PUGH PUGHInterp CartGrid3D"
+
+PUGHInterp::demo  = "yes"
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..9f1c756
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,9 @@
+# Parameter definitions for thorn PUGHInterp
+# $Header$
+
+#############################################################################
+### import CartGrid3D parameters
+#############################################################################
+shares: grid
+
+USES KEYWORD domain
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..5e5a7bc
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn PUGHInterp
+# $Header$
+
+schedule PUGHInterp_Startup at STARTUP after Driver_Startup
+{
+  LANG:C
+} "PUGHInterp startup routine"
diff --git a/src/Interpolate.c b/src/Interpolate.c
new file mode 100644
index 0000000..392cb18
--- /dev/null
+++ b/src/Interpolate.c
@@ -0,0 +1,427 @@
+ /*@@
+   @file      Interpolate.c
+   @date      Wed 17 Jan 2001
+   @author    Thomas Radke
+   @desc
+              Interpolation of arrays to arbitrary points
+
+              This interpolator is based on the Cactus 3.x Fortran version
+              written by Paul Walker. It also contains some nice optimization
+              features from Erik Schnetter.
+   @enddesc
+   @history
+   @date      Wed 17 Jan 2001
+   @author    Thomas Radke
+   @hdesc     Translation from Fortran to C
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Id$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c)
+
+
+/* the highest order of interpolation we support so far */
+#define MAXORDER  3
+
+/* the highest dimension for variables we can deal with (so far) */
+#define MAXDIM    3
+
+/* macro to do the interpolation of in_array[n] to out_array[n] at point[] */
+#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array,    \
+                    order, point, dims, n, coeff)                             \
+    {                                                                         \
+      int ii, jj, kk;                                                         \
+      cctk_type *fi;                                                          \
+      /* for some reason (probably loop unrolling & pipelining) the compiler  \
+         produces faster code if arrays are used instead of scalars */        \
+      cctk_subtype interp_result, fj[1], fk[1];                               \
+                                                                              \
+                                                                              \
+      interp_result = 0;                                                      \
+                                                                              \
+      /* NOTE: support >3D arrays by adding more loops */                     \
+      for (kk = 0; kk <= order; kk++)                                         \
+      {                                                                       \
+        fk[0] = 0;                                                            \
+        for (jj = 0; jj <= order; jj++)                                       \
+        {                                                                     \
+          /* NOTE: for >3D arrays adapt the index calculation here */         \
+          fi = (cctk_type *) in_array +                                       \
+               point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk));      \
+          fj[0] = 0;                                                          \
+          for (ii = 0; ii <= order; ii++)                                     \
+          {                                                                   \
+            fj[0] += fi[ii]subpart * coeff[0][ii];                            \
+          }                                                                   \
+          fk[0] += fj[0] * coeff[1][jj];                                      \
+        }                                                                     \
+        interp_result += fk[0] * coeff[2][kk];                                \
+      }                                                                       \
+                                                                              \
+      /* assign the result */                                                 \
+      ((cctk_type *) out_array)[n]subpart = interp_result;                    \
+    }
+
+/* this is needed by some preprocessors to pass into INTERPOLATE
+   as a dummy macro */
+#define NOTHING
+
+
+/*@@
+   @routine    PUGHInterp_Interpolate
+   @date       Wed 17 Jan 2001
+   @author     Thomas Radke
+   @desc
+               This routine interpolates a set of input arrays
+               to a set of output arrays (one-to-one) at arbitrary points\                     which are given by their coordinates and the underlying
+               regular, uniform grid.
+
+               Current limitations of this implementation are:
+
+                - arrays up to three (MAXDIM) dimensions only can be handled
+                - interpolation orders up to three (MAXORDER) only are supported
+                - coordinates must be given as CCTK_REAL types
+                - input and output array types must be the same
+                  (no type casting of interpolation results supported)
+
+               Despite of these limitations, the code was programmed almost
+               generic in that it can easily be extended to support higher-
+               dimensional arrays or more interpolation orders.
+               Please see the NOTES in this source file for details.
+   @enddesc
+
+   @var        num_points
+   @vdesc      number of points to interpolate at
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      dimensionality of the input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_arrays
+   @vdesc      number of input/output arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dims
+   @vdesc      dimensions of the input arrays
+   @vtype      int[ num_dims ]
+   @vio        in
+   @endvar
+   @var        coord
+   @vdesc      coordinates to interpolate at
+   @vtype      CCTK_REAL coord[ num_dims * num_points ]
+   @vio        in
+   @endvar
+   @var        origin
+   @vdesc      origin of the underlying grid
+   @vtype      CCTK_REAL origin[ num_dims ]
+   @vio        in
+   @endvar
+   @var        delta
+   @vdesc      deltas of the underlying grid
+   @vtype      CCTK_REAL delta[ num_dims ]
+   @vio        in
+   @endvar
+   @var        in_types
+   @vdesc      CCTK variable types of input arrays
+   @vtype      int in_types[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        in_arrays
+   @vdesc      list of input arrays
+   @vtype      void *in_arrays[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        out_types
+   @vdesc      CCTK variable types of output arrays
+   @vtype      int out_types[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        out_arrays
+   @vdesc      list of output arrays
+   @vtype      void *out_arrays[ num_arrays ]
+   @vio        out
+   @endvar
+
+   @returntype int
+   @returndesc
+               0  - successful interpolation
+               -1 - in case of any errors
+   @endreturndesc
+@@*/
+int PUGHInterp_Interpolate (int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            int dims[ /* num_dims */ ],
+                            CCTK_REAL coord[ /* num_dims * num_points */ ],
+                            CCTK_REAL origin[ /* num_dims */ ],
+                            CCTK_REAL delta[ /* num_dims */ ],
+                            int in_types[ /* num_arrays */ ],
+                            void *in_arrays[ /* num_arrays */ ],
+                            int out_types[ /* num_arrays */ ],
+                            void *out_arrays[ /* num_arrays */ ])
+{
+  int retval;
+  int i, a, n, out_of_bounds, shift;
+  int max_dims[MAXDIM], point[MAXDIM];
+  CCTK_REAL below[MAXDIM];
+  CCTK_REAL offset[MAXDIM];
+  CCTK_REAL coeff[MAXDIM][MAXORDER + 1];
+
+
+  /* verify parameters and check against our restrictions */
+  retval = -1;
+  if (num_dims < 1)
+  {
+    CCTK_WARN (1, "Number of dimensions must be positive");
+  }
+  else if (num_dims > MAXDIM)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Interpolation of %d-dimensional arrays not implemented",
+                num_dims);
+  }
+  else if (order < 1)
+  {
+    CCTK_WARN (1, "Inperpolation order must be positive");
+  }
+  else if (order > MAXORDER)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Interpolation order %d not implemented", order);
+  }
+  else if (num_points < 0)
+  {
+    CCTK_WARN (1, "Negative number of points given");
+  }
+  else
+  {
+    retval = 0;
+  }
+
+  /* immediately return in case of errors */
+  if (retval)
+  {
+    return (retval);
+  }
+
+  /* also immediately return if there's nothing to do */
+  if (num_points == 0)
+  {
+    return (retval);
+  }
+
+  /* duplicate the dims[] vector into one with MAXDIM-1 elements
+     (with the remaining elements zeroed out)
+     so that we can use nested loops over MAXDIM dimensions later on */
+  memset (max_dims, 0, sizeof (max_dims));
+  memcpy (max_dims, dims, (num_dims - 1) * sizeof (int));
+
+  /* zero out the coefficients and set the elements with index 'order' to one
+     so that we can use nested loops over MAXDIM dimensions later on */
+  memset (coeff, 0, sizeof (coeff));
+  for (i = num_dims; i < MAXDIM; i++)
+  {
+    coeff[i][0] = 1;
+  }
+
+  /* zero out the iterator */
+  memset (point, 0, sizeof (point));
+
+  /* loop over all points to interpolate at */
+  for (n = 0; n < num_points; n++)
+  {
+    /* reset the out-of-bounds flag */
+    out_of_bounds = 0;
+
+    /* loop over all dimensions */
+    for (i = 0; i < num_dims; i++)
+    {
+      /* nearest grid point for stencil */
+      point[i] = (coord[num_dims*n + i] - origin[i]) / delta[i];
+
+      /* if beyond upper bound shift the grid point to the left */
+      shift = point[i] + order - (dims[i] - 1);
+      if (shift > 0)
+      {
+        point[i] -= shift;
+      }
+
+      /* physical coordinate of that grid point */
+      below[i] = origin[i] + point[i] * delta[i];
+
+      /* offset from that grid point, in fractions of grid points */
+      offset[i] = (coord[num_dims*n + i] - below[i]) / delta[i];
+
+      /* test bounds */
+      out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i];
+    }
+
+    /* check bounds */
+    if (out_of_bounds)
+    {
+      char *msg;
+
+
+      /* put all information into a single message string for output */
+      msg = (char *) malloc (100 + num_dims*(10 + 4*20));
+      sprintf (msg, "Interpolation stencil out of bounds at grid point [%d",
+               point[0]);
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %d", msg, point[i]);
+      }
+      sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg,
+               (double) below[0], (double) (below[0] + offset[0]));
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (double) below[i], (double) (below[i] + offset[i]));
+      }
+      sprintf (msg, "%s]\ngrid is min/max        [%f / %f", msg,
+               (double) origin[0], (double) (origin[0] + dims[0]*delta[0]));
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (double) origin[i], (double) (origin[i] + dims[i]*delta[i]));
+      }
+      sprintf (msg, "%s]", msg);
+      CCTK_WARN (1, msg);
+      free (msg);
+
+      continue;
+    }
+
+    /* compute the interpolation coefficients according to the order
+       Thanks to Erik for formulating these so nicely. */
+    /* NOTE: support higher interpolation orders by adding the
+             appropriate coefficients in another else branch */
+    if (order == 1)
+    {
+      /* first order (linear) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = 1 - offset[i];
+        coeff[i][1] =     offset[i];
+      }
+    }
+    else if (order == 2)
+    {
+      /* second order (quadratic) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = (1-offset[i]) * (2-offset[i]) / (  2  *   1 );
+        coeff[i][1] = ( -offset[i]) * (2-offset[i]) / (  1  * (-1));
+        coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2));
+      }
+    }
+    else if (order == 3)
+    {
+      /* third order (cubic) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) /
+                      (  3  *   2  *   1 );
+        coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) /
+                      (  2  *   1  * (-1));
+        coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) /
+                      (  1  * (-1) * (-2));
+        coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) /
+                      ((-1) * (-2) * (-3));
+      }
+    }
+    else
+    {
+      CCTK_WARN (0, "Implementation error");
+    }
+
+    /* now loop over all arrays to interpolate at the current point */
+    for (a = 0; a < num_arrays; a++)
+    {
+      /* sorry, for now input and output arrays must be of same type */
+      if (in_types[a] != out_types[a])
+      {
+        CCTK_WARN (1, "Type casting of interpolation results not implemented");
+        continue;
+      }
+
+      /* now do the interpolation according to the array type
+         we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */
+      if (in_types[a] == CCTK_VARIABLE_REAL)
+      {
+        INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX)
+      {
+        INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#ifdef CCTK_REAL4
+      else if (in_types[a] == CCTK_VARIABLE_REAL4)
+      {
+        INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX8)
+      {
+        INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (in_types[a] == CCTK_VARIABLE_REAL8)
+      {
+        INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX16)
+      {
+        INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (in_types[a] == CCTK_VARIABLE_REAL16)
+      {
+        INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX32)
+      {
+        INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+      else
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Unsupported variable type %d", (int) in_types[a]);
+      }
+    } /* end of loop over all arrays */
+
+  } /* end of loop over all points to interpolate at */
+
+  /* we're done */
+  return (retval);
+}
diff --git a/src/Operator.c b/src/Operator.c
new file mode 100644
index 0000000..76a71f5
--- /dev/null
+++ b/src/Operator.c
@@ -0,0 +1,1293 @@
+ /*@@
+   @file      Operator.c
+   @date      Tue Apr 15 18:22:45 1997
+   @author    Paul Walker
+   @desc 
+              Definition of interpolation operators for regular uniform grids.
+   @enddesc 
+
+   @history
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @hdesc     conversion to Cactus 4.0 (copied from pughGetPoints.c)
+   @date      Wed 31 Jan 2001
+   @author    Thomas Radke
+   @hdesc     translation of fortran interpolators into C
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Operator_c)
+
+/* uncomment this to get some debugging output */
+/* #define PUGHINTERP_DEBUG 1 */
+
+/* macro do sort interpolation results from a single communication buffer
+   into their appropriate output arrays */
+#define SORT_TYPED_ARRAY(cctk_type, source, dest, points, index_array, offset)\
+          {                                                                   \
+            int i;                                                            \
+            cctk_type *src, *dst;                                             \
+                                                                              \
+                                                                              \
+            src = (cctk_type *) source;                                       \
+            dst = (cctk_type *) dest;                                         \
+            for (i = 0; i < points; i++)                                      \
+            {                                                                 \
+              dst[index_array[i + offset]] = *src++;                          \
+            }                                                                 \
+            source = (char *) src;                                            \
+          }
+
+
+
+#ifdef CCTK_MPI
+/* internal structure describing a handle for a single CCTK data type */
+typedef struct
+{
+  int vtypesize;              /* variable type's size in bytes */
+  MPI_Datatype mpi_type;      /* corresponding MPI datatype */
+  int num_arrays;             /* number of in/out arrays */
+  int array;                  /* iterator for num_arrays */
+  void *send_buffer;          /* communication send buffer for this type */
+  void *recv_buffer;          /* communication receive buffer for this type */
+  char *buffer;               /* work pointer for send_buffer */
+} type_desc_t;
+#endif
+
+
+/* prototypes of routines defined in this source file */
+static int PUGHInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      int interp_coord_array_types[]);
+#ifdef CCTK_MPI
+static int GetLocalCoords (cGH *GH,
+                           int num_points,
+                           const char *coord_system,
+                           pGExtras *extras,
+                           CCTK_REAL *coords[],
+                           int *num_local_points,
+                           CCTK_REAL **local_coords);
+#endif
+
+
+/*@@
+  @routine    PUGHInterp_InterpGV
+  @date       Sun Jul 04 1999
+  @author     Thomas Radke
+  @desc
+              The interpolation operator registered with the CCTK
+              under the name "regular uniform cartesian".
+
+              Interpolates a list of CCTK variables (domain-decomposed
+              grid functions or arrays) to a list of output arrays
+              (one-to-one) at a given number of interpolation points
+              (indicated by their coordinates). The points are located
+              on a coordinate system which is assumed to be a uniform
+              cartesian.
+  @enddesc
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH *
+  @vio        in
+  @endvar
+  @var        order
+  @vdesc      interpolation order
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coord_system
+  @vdesc      name of coordinate system to use for interpolation
+  @vtype      const char *
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_in_array_indices
+  @vdesc      number of input arrays (given by their indices)
+              to interpolate from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_out_arrays
+  @vdesc      number of output arrays to interpolate to
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        interp_coord_arrays
+  @vdesc      coordinates of points to interpolate at
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        interp_coord_array_types
+  @vdesc      CCTK data type of coordinate arrays
+  @vtype      int [num_dims]
+  @vio        in
+  @endvar
+  @var        in_arrays
+  @vdesc      list of input arrays to interpolate on
+  @vtype      void *[num_in_arrays]
+  @vio        in
+  @endvar
+  @var        in_array_types
+  @vdesc      CCTK data types of input arrays
+  @vtype      int [num_in_arrays]
+  @vio        in
+  @endvar
+  @var        out_arrays
+  @vdesc      list of output arrays to interpolate to
+  @vtype      void *[num_out_arrays]
+  @vio        out
+  @endvar
+  @var        out_array_types
+  @vdesc      CCTK data types of output arrays
+  @vtype      int [num_out_arrays]
+  @vio        in
+  @endvar
+
+  @returntype int
+  @returndesc
+              0  - successful interpolation
+              -1 - in case of any errors
+  @endreturndesc
+@@*/
+int PUGHInterp_InterpGV (cGH *GH,
+                         int order,
+                         const char *coord_system,
+                         int num_points,
+                         int num_in_array_indices,
+                         int num_out_arrays,
+                         void *interp_coord_arrays[],
+                         int interp_coord_array_types[],
+                         int in_array_indices[],
+                         void *out_arrays[],
+                         int out_array_types[])
+{
+  CCTK_REAL **data, *interp_local_coords;
+  CCTK_REAL *origin, *delta;
+  int nprocs;
+  int timelevel;
+  int dim, num_dims, *dims;
+  int array;
+  int point;
+  int retval;
+  void **in_arrays;
+  pGH *pughGH;
+  pGExtras *extras;
+  cGroupDynamicData group_data;
+#ifdef CCTK_MPI
+  int num_local_points;
+  int offset;
+  int proc;
+  int type;
+  void **local_out_arrays;
+  pughInterpGH *myGH;
+  int max_type;
+  type_desc_t *type_desc;
+#endif
+
+
+  /* get dimensionality of the coordinate system */
+  num_dims = CCTK_CoordSystemDim (coord_system);
+  if (num_dims <= 0)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Cannot get dimensions of coordinate system '%s'",
+                coord_system);
+    return (-1);
+  }
+
+  /* check other arguments */
+  retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+                                      num_in_array_indices, num_out_arrays,
+                                      interp_coord_array_types);
+  if (retval <= 0)
+  {
+    return (retval);
+  }
+
+  /* get dimensions, origin, and delta of the processor-local grid */
+  /* NOTE: getting the dimensions should be a flesh routine as well
+           for now we get the dimensions of every coordinate and take the
+           i'th element - this is inconsistent !! */
+  dims = (int *) malloc (2 * num_dims * sizeof (int));
+  origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+  delta  = origin + num_dims;
+  in_arrays = malloc (num_in_array_indices * sizeof (void *));
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    CCTK_GrouplshVI (GH, num_dims, dims + num_dims,
+                     CCTK_CoordIndex (dim + 1, NULL, coord_system));
+    dims[dim] = dims[dim + num_dims];
+
+    CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL,
+                          coord_system);
+    delta[dim] = (delta[dim] - origin[dim]) / dims[dim];
+  }
+
+  /* get extension handle for PUGH */
+  pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");
+  extras = NULL;
+
+  /* check that the input arrays dimensions match the coordinate system
+     (but their dimensionality can be less) */
+  for (array = 0; array < num_in_array_indices; array++)
+  {
+    if (CCTK_GroupDynamicData (GH,
+                               CCTK_GroupIndexFromVarI(in_array_indices[array]),
+                               &group_data) < 0)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Invalid input array index %d",
+                  in_array_indices[array]);
+      retval = -1;
+      continue;
+    }
+
+    if (group_data.dim > num_dims)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Input array variable with index %d has more dimensions "
+                  "than coordinate system '%s'",
+                  in_array_indices[array], coord_system);
+      retval = -1;
+      continue;
+    }
+
+    if (memcmp (group_data.lsh, dims, group_data.dim * sizeof (int)))
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Dimensions of input array variable with index %d "
+                  "doesn't match with coordinate system '%s'",
+                  in_array_indices[array], coord_system);
+      retval = -1;
+    }
+
+    /* the extras pointer will refer to the first input array with the
+       highest dimension (need this information later on) */
+    if (! extras || extras->dim < group_data.dim)
+    {
+      extras = ((pGA *) pughGH->variables[in_array_indices[array]][0])->extras;
+    }
+
+    /* get the data pointer to the input array (use current timelevel) */
+    timelevel = CCTK_NumTimeLevelsFromVarI (in_array_indices[array]) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    in_arrays[array] = CCTK_VarDataPtrI(GH, timelevel, in_array_indices[array]);
+  }
+  if (retval < 0)
+  {
+    free (in_arrays);
+    free (origin);
+    free (dims);
+    return (retval);
+  }
+
+  /* single-processor case is easy: no communication or buffering, just direct
+     interpolation of interp_coord_arrays from in_arrays into out_arrays */
+  nprocs = CCTK_nProcs (GH);
+  if (nprocs == 1)
+  {
+    /* sort the individual interpolation coordinate arrays into a single one */
+    interp_local_coords = (CCTK_REAL *) malloc (num_dims * num_points *
+                                                sizeof (CCTK_REAL));
+    data = (CCTK_REAL **) interp_coord_arrays;
+    for (point = 0; point < num_points; point++)
+    {
+      for (dim = 0; dim < num_dims; dim++)
+      {
+        *interp_local_coords++ = data[dim][point];
+      }
+    }
+    interp_local_coords -= num_dims * num_points;
+
+    /* call the interpolator function */
+    retval = PUGHInterp_Interpolate (order,
+                                     num_points, num_dims, num_out_arrays,
+                                     dims, interp_local_coords,
+                                     origin, delta,
+                                     out_array_types, in_arrays,
+                                     out_array_types, out_arrays);
+
+    /* free allocated resources */
+    free (interp_local_coords);
+    free (in_arrays);
+    free (origin);
+    free (dims);
+
+    return (retval);
+  }
+
+#ifdef CCTK_MPI
+  /*** Here follows the multi-processor case:
+       All processors locate their points to interpolate at
+       and exchange the coordinates so that every processor gets
+       those points which it can process locally.
+       After interpolation the results have to be send back to the
+       requesting processors.
+       For both communications MPI_Alltoallv() is used.
+
+       In order to minimize the total number of MPI_Alltoallv() calls
+       (which are quite expensive) we collect the interpolation results
+       for all output arrays of the same CCTK data type into a single
+       communication buffer. That means, after communication the data
+       needs to be resorted from the buffer into the output arrays.
+   ***/
+
+  /* first of all, set up a structure with information of the
+     CCTK data types we have to deal with */
+
+  /* get the maximum value of the output array CCTK data types
+     NOTE: we assume that CCTK data types are defined as consecutive
+           positive constants starting from zero */
+  max_type = 0;
+  for (array = 0; array < num_out_arrays; array++)
+  {
+    if (max_type < out_array_types[array])
+    {
+      max_type = out_array_types[array];
+    }
+  }
+
+  /* now allocate an array of structures for all potential types */
+  type_desc = (type_desc_t *) calloc (max_type + 1, sizeof (type_desc_t));
+
+  /* count the number of arrays of same type
+     (the num_arrays element was already initialized to zero by calloc() */
+  for (array = 0; array < num_out_arrays; array++)
+  {
+    type_desc[out_array_types[array]].num_arrays++;
+  }
+
+  /* fill in the type description information */
+  retval = 0;
+  for (type = 0; type <= max_type; type++)
+  {
+    if (type_desc[type].num_arrays > 0)
+    {
+      /* get the variable type size in bytes */
+      type_desc[type].vtypesize = CCTK_VarTypeSize (type);
+      if (type_desc[type].vtypesize <= 0)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Invalid variable type %d passed, "
+                    "arrays of such type will be skipped during interpolation",
+                    type);
+        type_desc[type].num_arrays = 0;
+        continue;
+      }
+
+      /* get the MPI data type to use for communicating such a CCTK data type */
+      type_desc[type].mpi_type = PUGH_MPIDataType (pughGH, type);
+      if (type_desc[type].mpi_type == MPI_DATATYPE_NULL)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "No MPI data type defined for variable type %d, "
+                    "arrays of such type will be skipped during interpolation",
+                    type);
+        type_desc[type].num_arrays = 0;
+        continue;
+      }
+
+      retval++;
+    }
+  }
+
+  /* check that there's at least one array with a valid CCTK data type */
+  if (retval <= 0)
+  {
+    free (in_arrays);
+    free (dims);
+    free (origin);
+    free (type_desc);
+    return (-1);
+  }
+
+  /* map the requested points to interpolate at onto the processors
+     they belong to and gather the coordinates of all the points
+     that this processor owns
+     the number of processor-local points is returned in num_local_points,
+     their coordinates in interp_local_coords */
+  retval = GetLocalCoords (GH, num_points, coord_system, extras,
+                           (CCTK_REAL **) interp_coord_arrays,
+                           &num_local_points, &interp_local_coords);
+  if (retval)
+  {
+    free (in_arrays);
+    free (origin);
+    free (dims);
+    free (type_desc);
+    return (retval);
+  }
+
+  /* allocate contiguous communication buffers for each CCTK data type
+     holding the local interpolation results from all input arrays
+     of that type
+     If there are no points to process on this processor
+     set the buffer pointer to an invalid but non-NULL value
+     otherwise we might get trouble with NULL pointers in MPI_Alltoallv () */
+  for (type = 0; type <= max_type; type++)
+  {
+    if (type_desc[type].num_arrays > 0 && num_local_points > 0)
+    {
+      type_desc[type].send_buffer = malloc (num_local_points *
+                                            type_desc[type].num_arrays *
+                                            type_desc[type].vtypesize);
+      type_desc[type].buffer = (char *) type_desc[type].send_buffer;
+      type_desc[type].array = 0;
+    }
+    else
+    {
+      /* dereferencing such an address should code crash on most systems */
+      type_desc[type].send_buffer = (void *) 1;
+    }
+  }
+
+  /* get extension handle for interp */
+  myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+
+  /* allocate new out_arrays array for local interpolation results
+     from this processor */
+  local_out_arrays = calloc (num_out_arrays, sizeof (void *));
+
+  /* now, in a loop over all processors, do the interpolation
+     and put the results in the communication buffer at the proper offset */
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (type = 0; type <= max_type; type++)
+    {
+      if (type_desc[type].num_arrays > 0)
+      {
+        for (array = 0; array < num_out_arrays; array++)
+        {
+          if (out_array_types[array] == type)
+          {
+            local_out_arrays[array] = type_desc[type].buffer;
+            type_desc[type].buffer += myGH->num_points_to[proc] *
+                                      type_desc[type].vtypesize;
+          }
+        }
+      }
+    }
+
+    /* call the interpolation operator to process all points of all
+       output arrays for this processor */
+    PUGHInterp_Interpolate (order,
+                            myGH->num_points_to[proc], num_dims, num_out_arrays,
+                            dims, interp_local_coords, origin, delta,
+                            out_array_types, in_arrays,
+                            out_array_types, local_out_arrays);
+
+    /* have to add offset for this processor to coordinates array */
+    interp_local_coords += myGH->num_points_to[proc] * num_dims;
+
+  } /* end of loop over all processors */
+
+  /* don't need these anymore */
+  if (num_local_points > 0)
+  {
+    interp_local_coords -= num_local_points * num_dims;
+    free (interp_local_coords);
+  }
+  free (local_out_arrays);
+  free (in_arrays);
+  free (origin);
+  free (dims);
+
+  /* now send the interpolation results back to the processors they were
+     requested from, also receive my own results that were computed
+     by other processors
+     Since all the locally computed results are in a single contiguous buffer
+     we need to call MPI_Alltoall() only once for each CCTK data type. */
+  for (type = 0; type <= max_type; type++)
+  {
+    /* skip unused types */
+    if (type_desc[type].num_arrays <= 0)
+    {
+      continue;
+    }
+
+    /* set up the communication (this is type-independent) */
+    myGH->send_count[0] = type_desc[type].num_arrays * myGH->num_points_to[0];
+    myGH->recv_count[0] = type_desc[type].num_arrays * myGH->num_points_from[0];
+    myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+    for (proc = 1; proc < pughGH->nprocs; proc++)
+    {
+      myGH->send_count[proc] = type_desc[type].num_arrays *
+                               myGH->num_points_to[proc];
+      myGH->recv_count[proc] = type_desc[type].num_arrays *
+                               myGH->num_points_from[proc];
+      myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+                               myGH->send_count[proc-1];
+      myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+                               myGH->recv_count[proc-1];
+    }
+
+    /* allocate buffer for receiving my own requested points */
+    /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */
+    if (num_points > 0)
+    {
+      type_desc[type].recv_buffer = malloc (num_points *
+                                            type_desc[type].num_arrays *
+                                            type_desc[type].vtypesize);
+    }
+    else
+    {
+      /* access to such a fake address should crash the code on most systems */
+      type_desc[type].recv_buffer = (void *) 1;
+    }
+
+    /* now exchange the data for this CCTK data type */
+    CACTUS_MPI_ERROR (MPI_Alltoallv (type_desc[type].send_buffer,
+                                     myGH->send_count, myGH->send_displ,
+                                     type_desc[type].mpi_type,
+                                     type_desc[type].recv_buffer,
+                                     myGH->recv_count, myGH->recv_displ,
+                                     type_desc[type].mpi_type,
+                                     pughGH->PUGH_COMM_WORLD));
+
+    /* now that the data is sent we don't need the buffer anymore */
+    if (num_local_points > 0)
+    {
+      free (type_desc[type].send_buffer);
+    }
+
+    /* no sort neccessary if there are no points */
+    if (num_points <= 0)
+    {
+      continue;
+    }
+
+    /* Go back from processor-sorted data to input-ordered data.
+       The creation of the indices array above makes this not so bad. */ 
+    offset = 0;
+    type_desc[type].buffer = (char *) type_desc[type].recv_buffer;
+    for (proc = 0; proc < nprocs; proc++)
+    {
+      for (array = 0; array < num_out_arrays; array++)
+      {
+        if (out_array_types[array] != type)
+        {
+          continue;
+        }
+
+        /* now do the sorting according to the CCTK data type */
+        if (out_array_types[array] == CCTK_VARIABLE_CHAR)
+        {
+          SORT_TYPED_ARRAY (CCTK_CHAR, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_INT)
+        {
+          SORT_TYPED_ARRAY (CCTK_INT, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#ifdef CCTK_REAL4
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL4)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL4, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX8)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX8, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL8)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL8, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX16)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX16, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL16)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL16, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX32)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX32, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+        else
+        {
+          CCTK_WARN (0, "Implementation error");
+        }
+
+      } /* end of loop over all output arrays */
+
+      /* advance the offset into the communication receive buffer */
+      offset += myGH->num_points_from[proc];
+
+    } /* end of loop over all processors */
+
+    /* this communication receive buffer isn't needed anymore */
+    free (type_desc[type].recv_buffer);
+
+  } /* end of loop over all types */
+
+  /* free remaining resources allocated within this run */
+  if (myGH->whichproc)
+  {
+    free (myGH->whichproc);
+    myGH->whichproc = NULL;
+  }
+  free (type_desc);
+#endif /* MPI */
+
+  return (0);
+}
+
+
+/*@@
+  @routine    PUGHInterp_InterpLocal
+  @date       Sun Jul 04 1999
+  @author     Thomas Radke
+  @desc
+              The interpolation operator registered with the CCTK
+              under the name "uniform cartesian".
+
+              Interpolates a list of non-distributed, processor-local
+              input arrays to a list of output arrays (one-to-one)
+              at a given number of interpolation points (indicated by
+              their coordinates). The points are located on a coordinate
+              system which is assumed to be a uniform cartesian.
+  @enddesc
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH *
+  @vio        in
+  @endvar
+  @var        order
+  @vdesc      interpolation order
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_dims
+  @vdesc      dimensionality of the underlying grid
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_in_arrays
+  @vdesc      number of input arrays to interpolate from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_out_arrays
+  @vdesc      number of output arrays to interpolate to
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coord_dims
+  @vdesc      dimensions of the underlying grid
+  @vtype      int [num_dims]
+  @vio        in
+  @endvar
+  @var        coord_arrays
+  @vdesc      list of grid coordinate arrays
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        coord_array_types
+  @vdesc      CCTK data type of coordinate arrays
+  @vtype      int [num_dims]
+  @vio        int
+  @endvar
+  @var        interp_coord_arrays
+  @vdesc      coordinates of points to interpolate at
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        interp_coord_array_types
+  @vdesc      CCTK data type of coordinate arrays to interpolate at
+  @vtype      int [num_dims]
+  @vio        out
+  @endvar
+  @var        in_arrays
+  @vdesc      list of input arrays to interpolate on
+  @vtype      void *[num_in_arrays]
+  @vio        in
+  @endvar
+  @var        in_array_types
+  @vdesc      CCTK data types of input arrays
+  @vtype      int [num_in_arrays]
+  @vio        in
+  @endvar
+  @var        out_arrays
+  @vdesc      list of output arrays to interpolate to
+  @vtype      void *[num_out_arrays]
+  @vio        out
+  @endvar
+  @var        out_array_types
+  @vdesc      CCTK data types of output arrays
+  @vtype      int [num_out_arrays]
+  @vio        in
+  @endvar
+
+  @returntype int
+  @returndesc
+              0  - successful interpolation
+              -1 - in case of any errors
+  @endreturndesc
+@@*/
+int PUGHInterp_InterpLocal (cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            int coord_dims[],
+                            void *coord_arrays[],
+                            int coord_array_types[],
+                            void *interp_coord_arrays[],
+                            int interp_coord_array_types[],
+                            void *in_arrays[],
+                            int in_array_types[],
+                            void *out_arrays[],
+                            int out_array_types[])
+{
+  int dim, point, retval;
+  CCTK_REAL *coords, *origin, *delta, **data;
+
+
+  /* check arguments */
+  retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+                                      num_in_arrays, num_out_arrays,
+                                      interp_coord_array_types);
+  if (retval <= 0)
+  {
+    return (retval);
+  }
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    if (coord_array_types[dim] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Coordinates should be of type CCTK_REAL");
+      return (-1);
+    }
+    if (interp_coord_array_types[dim] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Interpolation coordinates should be of type CCTK_REAL");
+      return (-1);
+    }
+  }
+
+  /* get the grid spacings - this assumes a cartesian grid */
+  origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+  delta  = origin + num_dims;
+  data   = (CCTK_REAL **) coord_arrays;
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    origin[dim] = data[dim][0];
+    delta[dim]  = data[dim][1] - data[dim][0];
+  }
+
+  /* sort the individual interpolation coordinate arrays into a single one */
+  coords = (CCTK_REAL *) malloc (num_dims * num_points * sizeof (CCTK_REAL));
+  data   = (CCTK_REAL **) interp_coord_arrays;
+  for (point = 0; point < num_points; point++)
+  {
+    for (dim = 0; dim < num_dims; dim++)
+    {
+      *coords++ = data[dim][point];
+    }
+  }
+  coords -= num_dims * num_points;
+
+  /* call the interpolator function */
+  retval = PUGHInterp_Interpolate (order,
+                                   num_points, num_dims, num_out_arrays,
+                                   coord_dims, coords, origin, delta,
+                                   in_array_types, in_arrays,
+                                   out_array_types, out_arrays);
+
+  /* free allocated resources */
+  free (coords);
+  free (origin);
+
+  return (retval);
+}
+
+
+/**************************************************************************/
+/*                            local routines                              */
+/**************************************************************************/
+
+#ifdef CCTK_MPI
+
+/*@@
+  @routine   GetLocalCoords
+  @date      Sun Jul 04 1999
+  @author    Thomas Radke
+  @desc 
+  Collect the coordinates of all points to be processed locally
+  into an array coords[num_dims][num_local_points].
+  <B>
+  This means for the single-processor case to sort
+
+   inCoords1[num_points], inCoords2[npoints], ..., inCoords<num_dims>[npoints]
+    into coords[num_dims][num_points]
+
+  where num_points == num_local_points.
+  <B>
+  In the multiprocessor case all processors map their points' coordinates
+  to the processor that owns this point and exchange this information
+  via MPI_Alltoall ().
+  num_local_points is then the number of all processors' points to be
+  interpolated locally on this processor.
+  <B>
+  This routine returns the number of points to be processed locally and
+  a pointer to the allocated array of their coordinates.
+  @enddesc 
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        isGlobal
+  @vdesc      flag indicating that coordinates are global (and need to be
+              collected over all processors)
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_dims
+  @vdesc      number of coordinate dimensions for each point
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coords
+  @vdesc      coordinates of each point to be interpolated on this processor
+  @vtype      CCTK_REAL array of size num_dims
+  @vio        in
+  @endvar
+  @var        num_local_points
+  @vdesc      number of points to be processed by this processor
+  @vtype      int *
+  @vio        out
+  @endvar
+  @var        local_coords
+  @vdesc      coordinates of each point to be processed by this processor
+  @vtype      pointer to CCTK_REAL array of dims[num_dims][num_local_points]
+  @vio        out
+  @endvar
+
+@@*/
+static int GetLocalCoords (cGH *GH,
+                           int num_points,
+                           const char *coord_system,
+                           pGExtras *extras,
+                           CCTK_REAL *coords[],
+                           int *num_local_points,
+                           CCTK_REAL **local_coords)
+{
+  DECLARE_CCTK_PARAMETERS
+  int dim, point;
+  int tmp, proc, nprocs;
+  pGH *pughGH;
+  pughInterpGH *myGH;
+  CCTK_REAL *range_min, *range_max;
+  CCTK_REAL *origin, *delta;
+  CCTK_REAL *proc_coords;
+#define FUDGE 0.00001
+
+
+  /* Do a fliperooni in octant mode */
+  /* This flips the sign around the origin only */
+  /* if the left side (cx0) is bigger than the position to interpolate at */
+  /* BoxInBox needs this for interpolating onto grids */ 
+  if (CCTK_Equals (domain, "octant"))
+  {
+    for (point = 0; point < num_points; point++)
+    {
+      tmp = 0;
+      for (dim = 0; dim < extras->dim; dim++)
+      {
+        if (coords[dim][point] <= GH->cctk_origin_space[dim])
+        {
+          tmp = 1;                 /* mark flipping */
+          break;
+        }
+      }
+/*** FIXME *** really change input arrays here ??? */
+      if (tmp)
+      {
+        for (dim = 0; dim < extras->dim; dim++)
+        {
+          coords[dim][point] = fabs (coords[dim][point]);
+        }
+      }
+    }
+  }
+
+  /* get GH extension handles for PUGHInterp and PUGH */
+  myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+  pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");
+
+  /* This holds the proccessor for *each* of num_points points */
+  if (num_points > 0)
+  {
+    myGH->whichproc = (int *) malloc (2 * num_points * sizeof (int));
+  }
+  else
+  {
+    myGH->whichproc = NULL;
+  }
+  /* indices[] is used to make the sorting easier
+     when receiving the output data */
+  myGH->indices = myGH->whichproc + num_points;
+
+  /* initialize whichproc with invalid processor number -1 */
+  for (point = 0; point < num_points; point++)
+  {
+    myGH->whichproc[point] = -1;
+  }
+  
+  /* initialize num_points_from to 0 for counting it up in the following loop */
+  nprocs = CCTK_nProcs (GH);
+  memset (myGH->num_points_from, 0, nprocs * sizeof (CCTK_INT));
+
+  /* allocate the ranges for my local coordinates */
+  range_min = (CCTK_REAL *) malloc (4 * extras->dim * sizeof (CCTK_REAL));
+  range_max = range_min + 1*extras->dim;
+  origin    = range_min + 2*extras->dim;
+  delta     = range_min + 3*extras->dim;
+
+  /* get the global origin and delta of the coordinate system */
+  for (dim = 0; dim < extras->dim; dim++)
+  {
+    CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system);
+    delta[dim] = (delta[dim] - origin[dim]) / extras->nsize[dim];
+  }
+
+  /* locate the points to interpolate at */
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (dim = 0; dim < extras->dim; dim++)
+    {
+      /* compute the coordinate ranges */
+      range_min[dim] = origin[dim] + (extras->lb[proc][dim] - FUDGE)*delta[dim];
+      range_max[dim] = origin[dim] + (extras->ub[proc][dim] + FUDGE)*delta[dim];
+    }
+
+    /* and now which point will be processed by what processor */
+    for (point = 0; point < num_points; point++)
+    {
+      /* skip points which have already been located */
+      if (myGH->whichproc[point] >= 0)
+      {
+        continue;
+      }
+
+      /* check whether the point belongs to this processor
+         (must be within min/max in all dimensions) */
+      tmp = 0;
+      for (dim = 0; dim < extras->dim; dim++)
+      {
+        if (coords[dim][point] >= range_min[dim] &&
+            coords[dim][point] <= range_max[dim])
+        {
+          tmp++;
+        }
+      }
+      if (tmp == extras->dim)
+      {
+        myGH->whichproc[point] = proc;
+        myGH->num_points_from[proc]++;
+      }
+    }
+  }
+  /* don't need this anymore */
+  free (range_min);
+
+  /* make sure that all points could be mapped onto a processor */
+  tmp = 0;
+  for (point = 0; point < num_points; point++)
+  {
+    if (myGH->whichproc[point] < 0)
+    {
+      int i;
+      char *msg = (char *) malloc (80 + extras->dim*20);
+
+
+      sprintf (msg, "Unable to locate point %d [%f",
+               point, coords[0][point]);
+      for (i = 1; i < extras->dim; i++)
+      {
+        sprintf (msg, "%s %f", msg, coords[i][point]);
+      }
+      sprintf (msg, "%s]", msg);
+      CCTK_WARN (1, msg);
+      free (msg);
+      tmp = 1;                       /* mark as error */
+    }
+  }
+  if (tmp)
+  {
+    if (myGH->whichproc)
+    {
+      free (myGH->whichproc);
+      myGH->whichproc = NULL;
+    }
+    return (-1);
+  }
+
+  /* Now we want to resolve the num_points_from[]. Currently this is
+     the form of ( in 2 proc mode )
+     P1:  Num from P1  NFP2
+     P2:  NFP1         NFP2
+
+     and this needs to become
+     P1:  Num to P1    NTP2
+     P2:  NTP1         NTP1
+
+     Since NTP1 = NFP2 (and this works in more proc mode too)
+     this is an all-to-all communication.
+     */
+  CACTUS_MPI_ERROR (MPI_Alltoall (myGH->num_points_from, 1, PUGH_MPI_INT,
+                                  myGH->num_points_to, 1, PUGH_MPI_INT,
+                                  pughGH->PUGH_COMM_WORLD));
+
+#ifdef PUGHINTERP_DEBUG
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    printf ("processor %d <-> %d  From: %d  To: %d\n",
+            CCTK_MyProc (GH), proc, myGH->num_points_from[proc],
+            myGH->num_points_to[proc]);
+  }
+#endif
+
+  /* Great. Now we know how many to expect from each processor,
+     and how many to send to each processor. So first we have
+     to send the locations to the processors which hold our data.
+     This means I send coords[dim][point] to whichproc[point].
+     I have num_points_from[proc] to send to each processor.
+     */
+  
+  /* This is backwards in the broadcast location; the number of points
+     we are getting is how many everyone else is sending to us,
+     eg, num_points_to, not how many we get back from everyone else,
+     eg, num_points_from. The number we are sending, of course, is
+     all of our locations, eg, num_points */
+  *num_local_points = 0;
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    *num_local_points += myGH->num_points_to[proc];
+  }
+
+#ifdef PUGHINTERP_DEBUG
+  printf ("processor %d gets %d points in total\n",
+           CCTK_MyProc (GH), *num_local_points);
+#endif
+
+  /* allocate the local coordinates array (sorted in processor order)
+     and the resulting coordinates array that I have to process */
+  proc_coords = (CCTK_REAL *) malloc (extras->dim * num_points *
+                                      sizeof (CCTK_REAL));
+  *local_coords = (CCTK_REAL *) malloc (extras->dim * *num_local_points *
+                                        sizeof (CCTK_REAL));
+
+  /* now sort my own coordinates as tupels of [extras->dim] */
+  tmp = 0;
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (point = 0; point < num_points; point++)
+    {
+      if (myGH->whichproc[point] == proc)
+      {
+        for (dim = 0; dim < extras->dim; dim++)
+        {
+          *proc_coords++ = coords[dim][point];
+        }
+        myGH->indices[tmp++] = point;
+      }
+    }
+  }
+  proc_coords -= tmp * extras->dim;
+
+  /* So load up the send and recv stuff */
+  /* Send extras->dim elements per data point */
+  myGH->send_count[0] = extras->dim * myGH->num_points_from[0];
+  myGH->recv_count[0] = extras->dim * myGH->num_points_to[0];
+  myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+  for (proc = 1; proc < nprocs; proc++)
+  {
+    myGH->send_count[proc] = extras->dim * myGH->num_points_from[proc];
+    myGH->recv_count[proc] = extras->dim * myGH->num_points_to[proc];
+    myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+                              myGH->send_count[proc-1];
+    myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+                              myGH->recv_count[proc-1];
+  }
+ 
+  /* Great, and now exchange the coordinates and collect the ones
+     that I have to process in *local_coords[] */
+  CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->send_count,
+                                   myGH->send_displ, PUGH_MPI_REAL,
+                                   *local_coords, myGH->recv_count,
+                                   myGH->recv_displ, PUGH_MPI_REAL,
+                                   pughGH->PUGH_COMM_WORLD));
+
+  /* don't need this anymore */
+  free (proc_coords);
+
+  return (0);
+}
+#endif /* CCTK_MPI */
+
+
+ /*@@
+   @routine    PUGHInterp_CheckArguments
+   @date       Thu 25 Jan 2001
+   @author     Thomas Radke
+   @desc
+               Checks the interpolation arguments passed in via
+               the flesh's general interpolation calling interface
+
+               This routine also verifies that the parameters meet
+               the limitations of PUGHInterp's interpolation operators.
+   @enddesc
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      dimensionality of the underlying grid
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_points
+   @vdesc      number of points to interpolate at
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_in_arrays
+   @vdesc      number of passed input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_out_arrays
+   @vdesc      number of passed input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        interp_coord_array_types
+   @vdesc      types of passed coordinates to interpolate at
+   @vtype      int [num_dims]
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               +1 for success
+                0 for success but nothing to do
+               -1 for failure (wrong parameters passed or limitations not met)
+   @endreturndesc
+@@*/
+static int PUGHInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      int interp_coord_array_types[])
+{
+  int i;
+
+
+  /* check for invalid arguments */
+  if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0)
+  {
+    return (-1);
+  }
+
+  /* check if there's anything to do at all */
+  /* NOTE: num_points can be 0 in a collective call */
+  if (num_dims == 0 || (CCTK_nProcs (GH) == 1 && num_points == 0) ||
+      num_in_arrays == 0 || num_out_arrays == 0)
+  {
+    return (0);
+  }
+
+  /* for now we can only deal with coordinates of type CCTK_REAL */
+  for (i = 0; i < num_dims; i++)
+  {
+    if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL");
+      return (-1);
+    }
+  }
+
+  /* PUGHInterp's interpolation operators compute one output array
+     per input array */
+  if (num_in_arrays != num_out_arrays)
+  {
+    CCTK_WARN (1, "Number of input arrays must match number of output arrays");
+    return (-1);
+  }
+
+  return (1);
+}
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..729e712
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,253 @@
+ /*@@
+   @file      Startup.c
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Startup routines for PUGHInterp
+   @enddesc
+   @history
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Interp.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Startup_c)
+
+
+/* prototypes of routines defined in this source file */
+void PUGHInterp_Startup (void);
+#ifdef CCTK_MPI
+static void *PUGHInterp_SetupGH (tFleshConfig *config,
+                                 int convergence_level,
+                                 cGH *GH);
+#endif
+
+
+#ifdef CCTK_MPI
+ /*@@
+   @routine   PUGHInterp_SetupGH
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Allocates the GH extension structure and - for the MPI case -
+              the count/displacement buffers used in MPI_Alltoall()
+              and MPI_Alltoallv()
+   @enddesc
+   @calls     CCTK_nProcs
+
+   @returntype void *
+   @returndesc
+               pointer to the allocated GH extension structure
+   @endreturndesc
+@@*/
+static void *PUGHInterp_SetupGH (tFleshConfig *config,
+                                 int convergence_level,
+                                 cGH *GH)
+{
+  int nprocs;
+  pughInterpGH *myGH;
+
+
+  /* suppress compiler warnings about unused variables */
+  config = config;
+  convergence_level = convergence_level;
+
+  myGH = (pughInterpGH *) malloc (sizeof (pughInterpGH));
+
+  /* allocate once for all fields of same type */
+  nprocs = CCTK_nProcs (GH);
+  myGH->send_count = (int *) malloc (4 * nprocs * sizeof (int));
+  myGH->send_displ = myGH->send_count + nprocs;
+  myGH->recv_count = myGH->send_displ + nprocs;
+  myGH->recv_displ = myGH->recv_count + nprocs;
+
+  myGH->num_points_from = (CCTK_INT *) malloc (2 * nprocs * sizeof (CCTK_INT));
+  myGH->num_points_to   = myGH->num_points_from + nprocs;
+
+  return (myGH);
+}
+#endif /* CCTK_MPI */
+
+
+ /*@@
+   @routine   PUGHInterp_InterpGV_NthOrder
+   @date      Wed 14 Feb 2001
+   @author    Thomas Radke
+   @desc
+              Wrappers for the different interpolation operators 
+              registered for first/second/third order interpolation.
+              These wrappers just call the common interpolation routine
+              passing all arguments plus the interpolation order.
+   @enddesc
+
+   @returntype int
+   @returndesc
+               the return code of the common interpolation routine
+   @endreturndesc
+@@*/
+static int PUGHInterp_InterpGV_1stOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 1, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpGV_2ndOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 2, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpGV_3rdOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 3, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpLocal_1stOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                                   int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 1, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+static int PUGHInterp_InterpLocal_2ndOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                            int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 2, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+static int PUGHInterp_InterpLocal_3rdOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                            int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 3, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+ /*@@
+   @routine   PUGHInterp_Startup
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              The startup registration routine for PUGHInterp.
+              Registers the PUGHInterp GH extensions setup routine
+              and the interpolation operators with the flesh.
+   @enddesc
+   @calls     CCTK_RegisterGHExtensionSetupGH
+              CCTK_InterpRegisterOperatorGV
+              CCTK_InterpRegisterOperatorLocal
+@@*/
+void PUGHInterp_Startup (void)
+{
+#ifdef CCTK_MPI
+  CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("PUGHInterp"),
+                                   PUGHInterp_SetupGH);
+#endif
+
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_1stOrder,
+                                 "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_1stOrder,
+                                    "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_2ndOrder,
+                                 "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_2ndOrder,
+                                    "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_3rdOrder,
+                                 "third-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_3rdOrder,
+                                    "third-order uniform cartesian");
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..b8d41c1
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,5 @@
+# Main make.code.defn file for thorn PUGHInterp
+# $Header$
+
+# Source files in this directory
+SRCS = Startup.c Operator.c Interpolate.c
diff --git a/src/pughInterpGH.h b/src/pughInterpGH.h
new file mode 100644
index 0000000..e6eef76
--- /dev/null
+++ b/src/pughInterpGH.h
@@ -0,0 +1,71 @@
+ /*@@
+   @file      pughInterpGH.h
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Definition of GH extensions of thorn PUGHInterp and
+              declaration of function prototypes
+   @enddesc
+   @history
+   @endhistory
+@@*/
+
+
+typedef struct
+{
+  int *send_count;             /* number of elements to be sent */
+  int *send_displ;             /* offset of each element to be sent */
+  int *recv_count;             /* number of elements to be received */
+  int *recv_displ;             /* offset of each element to be received */
+
+  CCTK_INT *num_points_to;     /* number of coordinate points to be sent */
+  CCTK_INT *num_points_from;   /* number of coordinate points to be received */
+
+  int *whichproc;              /* which processor does point i belong to */
+  int *indices;                /* indices to sort from processor-sorted into
+                                  point-sorted order */
+} pughInterpGH;
+
+
+/* prototypes of interpolation operators to be registered */
+int PUGHInterp_InterpGV (cGH *GH,
+                         int order,
+                         const char *coord_system,
+                         int num_points,
+                         int num_in_array_indices,
+                         int num_out_arrays,
+                         void *interp_coord_arrays[],
+                         int interp_coord_array_types[],
+                         int in_array_indices[],
+                         void *out_arrays[],
+                         int out_array_types[]);
+
+int PUGHInterp_InterpLocal (cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            int coord_dims[],
+                            void *coord_arrays[],
+                            int coord_array_types[],
+                            void *interp_coord_arrays[],
+                            int interp_coord_array_types[],
+                            void *in_arrays[],
+                            int in_array_types[],
+                            void *out_arrays[],
+                            int out_array_types[]);
+
+/* prototypes of routines used internally */
+int PUGHInterp_Interpolate (int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            int dims[ /* num_dims */ ],
+                            CCTK_REAL coord[ /* num_dims * num_points */ ],
+                            CCTK_REAL origin[ /* num_dims */ ],
+                            CCTK_REAL delta[ /* num_dims */ ],
+                            int in_types[ /* num_arrays */ ],
+                            void *in_arrays[ /* num_arrays */ ],
+                            int out_types[ /* num_arrays */ ],
+                            void *out_arrays[ /* num_arrays */ ]);