Implement vector groups correctly

darcs-hash:20040323114027-07bb3-1ac44b7b144e8936307942bc999d54ecb579719f.gz

5 files changed, 85 insertions, 28 deletions

diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 4c30b44cf..ec6c87858 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.46 2004/03/23 12:14:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.47 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 #include <limits.h>
@@ -28,17 +28,34 @@ using namespace std;
 
 // Constructors
 template<class T, int D>
-data<T,D>::data (const int varindex_, const operator_type transport_operator_)
+data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+                 const int vectorlength, const int vectorindex,
+                 data* const vectorleader)
   : gdata<D>(varindex_, transport_operator_),
-    _storage(0)
-{ }
+    _storage(0),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
+{
+  assert (vectorlength>=1);
+  assert (vectorindex>=0 && vectorindex<vectorlength);
+  assert ((vectorindex==0 && !vectorleader)
+          || (vectorindex!=0 && vectorleader));
+}
 
 template<class T, int D>
 data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+                 const int vectorlength, const int vectorindex,
+                 data* const vectorleader,
                  const ibbox& extent_, const int proc_)
   : gdata<D>(varindex_, transport_operator_),
-    _storage(0)
+    _storage(0),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
 {
+  assert (vectorlength>=1);
+  assert (vectorindex>=0 && vectorindex<vectorlength);
+  assert ((vectorindex==0 && !vectorleader)
+          || (vectorindex!=0 && vectorleader));
   allocate(extent_, proc_);
 }
 
@@ -80,7 +97,13 @@ void data<T,D>::allocate (const ibbox& extent_, const int proc_,
   if (rank==this->_proc) {
     this->_owns_storage = !mem;
     if (this->_owns_storage) {
-      this->_storage = new T[this->_size];
+      if (this->vectorindex == 0) {
+        assert (! this->vectorleader);
+        this->_storage = new T[this->vectorlength * this->_size];
+      } else {
+        assert (this->vectorleader);
+        this->_storage = this->vectorleader->vectordata (this->vectorindex);
+      }
     } else {
       this->_storage = (T*)mem;
     }
@@ -91,19 +114,33 @@ void data<T,D>::allocate (const ibbox& extent_, const int proc_,
 
 template<class T, int D>
 void data<T,D>::free () {
-  if (this->_storage && this->_owns_storage) delete [] _storage;
+  if (this->_storage && this->_owns_storage && this->vectorindex==0) {
+    delete [] _storage;
+  }
   _storage = 0;
   this->_has_storage = false;
 }
 
 template<class T, int D>
 void data<T,D>::transfer_from (gdata<D>* gsrc) {
+  assert (this->vectorlength==1);
   data* src = (data*)gsrc;
+  assert (src->vectorlength==1);
   assert (!_storage);
   *this = *src;
   *src = data(this->varindex, this->transport_operator);
 }
 
+template<class T, int D>
+T* data<T,D>::vectordata (const int vectorindex) const
+{
+  assert (this->vectorindex==0);
+  assert (! this->vectorleader);
+  assert (vectorindex>=0 && vectorindex<this->vectorlength);
+  assert (this->_storage && this->_owns_storage);
+  return this->_storage + vectorindex * this->_size;
+}
+
 
 
 // Processor management
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index a6658e23e..85e8a28fe 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.17 2004/03/23 12:14:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.18 2004/03/23 12:40:27 schnetter Exp $
 
 #ifndef DATA_HH
 #define DATA_HH
@@ -30,13 +30,21 @@ class data: public gdata<D> {
 
   // Fields
   T* _storage;			// the data (if located on this processor)
+   
+  int vectorlength;
+  int vectorindex;
+  data* vectorleader;
   
 public:
   
   // Constructors
   data (const int varindex = -1,
-        const operator_type transport_operator = op_error);
+        const operator_type transport_operator = op_error,
+        const int vectorlength = 1, const int vectorindex = 0,
+        data* const vectorleader = NULL);
   data (const int varindex, const operator_type transport_operator,
+        const int vectorlength, const int vectorindex,
+        data* const vectorleader,
         const ibbox& extent, const int proc);
 
   // Destructors
@@ -52,6 +60,10 @@ public:
   virtual void free ();
   virtual void transfer_from (gdata<D>* gsrc);
 
+private:
+  T* vectordata (const int vectorindex) const;
+public:
+
   // Processor management
   virtual void change_processor (comm_state<D>& state,
                                  const int newproc, void* const mem=0);
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 2fcec0d2f..b2b5198d4 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.15 2004/01/25 14:57:30 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.16 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 
@@ -13,15 +13,16 @@ using namespace std;
 
 
 // Constructors
-// VGF
 template<class T,int D>
 gf<T,D>::gf (const int varindex, const operator_type transport_operator,
              th<D>& t, dh<D>& d,
-	     const int tmin, const int tmax, const int prolongation_order_time)
+	     const int tmin, const int tmax, const int prolongation_order_time,
+             const int vectorlength, const int vectorindex,
+             gf* const vectorleader)
   : ggf<D>(varindex, transport_operator,
-           t, d, tmin, tmax, prolongation_order_time)
+           t, d, tmin, tmax, prolongation_order_time,
+           vectorlength, vectorindex, vectorleader)
 {
-  // VGF
   this->recompose (0, true);
 }
 
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index a1da9d3d6..877308389 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.34 2004/03/23 09:26:49 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.35 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -23,12 +23,16 @@ template<int D>
 ggf<D>::ggf (const int varindex, const operator_type transport_operator,
              th<D>& t, dh<D>& d,
              const int tmin, const int tmax,
-             const int prolongation_order_time)
+             const int prolongation_order_time,
+             const int vectorlength, const int vectorindex,
+             ggf* const vectorleader)
   : varindex(varindex), transport_operator(transport_operator), t(t),
     tmin(tmin), tmax(tmax),
     prolongation_order_time(prolongation_order_time),
     h(d.h), d(d),
-    storage(tmax-tmin+1)
+    storage(tmax-tmin+1),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
 {
   assert (&t.h == &d.h);
 
@@ -50,7 +54,6 @@ bool ggf<D>::operator== (const ggf<D>& f) const {
 
 
 // Modifiers
-// VGF
 template<int D>
 void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
   
@@ -68,7 +71,7 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
       for (int c=0; c<h.components(rl); ++c) {
       	storage[tl-tmin][rl][c].resize(h.mglevels(rl,c));
 	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-          storage[tl-tmin][rl][c][ml] = 0;
+          storage[tl-tmin][rl][c][ml] = NULL;
 	} // for ml
       } // for c
     } // for rl
@@ -80,10 +83,9 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
       for (int c=0; c<h.components(rl); ++c) {
       	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
 	  
-	  storage[tl-tmin][rl][c][ml] = typed_data();
+	  storage[tl-tmin][rl][c][ml] = typed_data(tl,rl,c,ml);
 	  
       	  // Allocate storage
-          // VGF
       	  storage[tl-tmin][rl][c][ml]->allocate
       	    (d.boxes[rl][c][ml].exterior, h.proc(rl,c));
 	  
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 7b7c5c159..64998c8a2 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.19 2004/03/23 09:26:49 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.20 2004/03/23 12:40:27 schnetter Exp $
 
 #ifndef GGF_HH
 #define GGF_HH
@@ -47,9 +47,9 @@ class ggf {
   typedef vector<mdata> cdata;  // ... for each component
   typedef vector<cdata> rdata;  // ... for each refinement level
   typedef vector<rdata> fdata;  // ... for each time level
-
+  
 public:				// should be readonly
-
+  
   // Fields
   int varindex;                 // Cactus variable index
   operator_type transport_operator;
@@ -63,14 +63,20 @@ public:				// should be readonly
 
 protected:
   fdata storage;		// storage
-
+  
+  int vectorlength;             // vector length
+  int vectorindex;              // index of *this
+  ggf* vectorleader;            // first vector element
+  
 public:
 
   // Constructors
   ggf (const int varindex, const operator_type transport_operator,
        th<D>& t, dh<D>& d,
        const int tmin, const int tmax,
-       const int prolongation_order_time);
+       const int prolongation_order_time,
+       const int vectorlength, const int vectorindex,
+       ggf* const vectorleader);
 
   // Destructors
   virtual ~ggf ();
@@ -81,7 +87,6 @@ public:
 
 
   // Modifiers
-  // VGF
   void recompose (const int initialise_from, const bool do_prolongate);
 
   // Cycle the time levels by rotating the data sets
@@ -96,7 +101,7 @@ public:
   
 protected:
   
-  virtual gdata<D>* typed_data() = 0;
+  virtual gdata<D>* typed_data (int tl, int rl, int c, int ml) = 0;