Add normalize routine to normalise bboxsets. This isn't perfect yet,

Add normalize routine to normalise bboxsets.  This isn't perfect yet,
i.e. it still misses some cases.

Fix bug in calculating the regions that should be restricted to:
exclude all regions used as source for boundary prolongations, even if
these prolongations go to a different processor.

darcs-hash:20030318163025-07bb3-1cfed59a57c2db71d3dae5528e93a570ccb13101.gz

5 files changed, 102 insertions, 12 deletions

diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index b440bb59d..4b7b7c831 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -1,9 +1,10 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.12 2003/02/25 22:57:00 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.13 2003/03/18 17:30:25 schnetter Exp $
 
 #include <assert.h>
 
 #include <iostream>
 #include <set>
+#include <stack>
 
 #include "defs.hh"
 
@@ -58,8 +59,67 @@ bool bboxset<T,D>::invariant () const {
 // Normalisation
 template<class T, int D>
 void bboxset<T,D>::normalize () {
-  // TODO
   assert (invariant());
+  const int num_initial_boxes = bs.size();
+  int num_combined_boxes = 0;
+  stack<box> todo, done;
+  for (typename set<box>::const_iterator elt = bs.begin(); elt != bs.end(); ++elt) {
+    done.push (*elt);
+  }
+  // TODO: This will not catch all cases where bboxes can be combined.
+  for (int d=0; d<D; ++d) {
+    todo = done;
+    done = stack<box>();
+    while (! todo.empty()) {
+    restart:;
+      box item = todo.top();
+      todo.pop();
+      stack<box> work = done;
+      done = stack<box>();
+      while (! work.empty()) {
+        box comp = work.top();
+        work.pop();
+        {
+          assert (all(comp.stride() == item.stride()));
+          if (comp.upper()[d] + item.stride()[d] == item.lower()[d]) {
+            if (all((comp.lower() == item.lower()
+                     && comp.upper() == item.upper()).replace (d, true))) {
+              box newbox = box(comp.lower(), item.upper(), item.stride());
+              todo.push (newbox);
+              while (! work.empty()) {
+                done.push (work.top());
+                work.pop();
+              }
+              ++num_combined_boxes;
+              goto restart;
+            }
+          }
+          if (item.upper()[d] + item.stride()[d] == comp.lower()[d]) {
+            if (all((comp.lower() == item.lower()
+                     && comp.upper() == item.upper()).replace (d, true))) {
+              box newbox = box(item.lower(), comp.upper(), item.stride());
+              todo.push (newbox);
+              while (! work.empty()) {
+                done.push (work.top());
+                work.pop();
+              }
+              ++num_combined_boxes;
+              goto restart;
+            }
+          }
+        }
+        done.push (comp);
+      } // while work
+      done.push (item);
+    } // while todo
+  } // for d
+  bs.clear();
+  while (! done.empty()) {
+    bs.insert (done.top());
+    done.pop();
+  }
+  const int num_final_boxes = bs.size();
+  assert (num_initial_boxes - num_combined_boxes == num_final_boxes);
 }
 
 
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 3309fb31c..19c79ae7d 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.14 2003/02/27 22:23:19 shawley Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.15 2003/03/18 17:30:25 schnetter Exp $
 
 #include <assert.h>
 #include <ctype.h>
@@ -6,6 +6,7 @@
 #include <iostream>
 #include <list>
 #include <set>
+#include <stack>
 #include <vector>
 
 #include "defs.hh"
@@ -71,6 +72,18 @@ ostream& output (ostream& os, const set<T>& s) {
   return os;
 }
 
+// Stack output
+template<class T>
+ostream& output (ostream& os, const stack<T>& s) {
+  stack<T> s2 (s);
+  list<T> l;
+  while (! s2.empty()) {
+    l.insert (l.begin(), s2.top());
+    s2.pop();
+  }
+  return output (os, l);
+}
+
 // Vector output
 template<class T>
 ostream& output (ostream& os, const vector<T>& v) {
@@ -96,6 +109,7 @@ template istream& input (istream& os, vector<vector<vect<vect<bool,2>,3> > >& v)
 template ostream& output (ostream& os, const list<bbox<int,3> >& l);
 template ostream& output (ostream& os, const set<bbox<int,3> >& s);
 template ostream& output (ostream& os, const set<bboxset<int,3> >& s);
+template ostream& output (ostream& os, const stack<bbox<int,3> >& s);
 template ostream& output (ostream& os, const vector<int>& v);
 template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
 template ostream& output (ostream& os, const vector<list<bbox<int,3> > >& v);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 4e3cab8d8..b948d7a9a 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.10 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.11 2003/03/18 17:30:25 schnetter Exp $
 
 #ifndef DEFS_HH
 #define DEFS_HH
@@ -14,6 +14,7 @@
 #include <iostream>
 #include <list>
 #include <set>
+#include <stack>
 #include <vector>
 
 using namespace std;
@@ -121,6 +122,7 @@ inline istream& operator>> (istream& is, vector<T>& v) {
 // Container output
 template<class T> ostream& output (ostream& os, const list<T>& l);
 template<class T> ostream& output (ostream& os, const set<T>& s);
+template<class T> ostream& output (ostream& os, const stack<T>& s);
 template<class T> ostream& output (ostream& os, const vector<T>& v);
 
 template<class T>
@@ -134,6 +136,11 @@ inline ostream& operator<< (ostream& os, const set<T>& s) {
 }
 
 template<class T>
+inline ostream& operator<< (ostream& os, const stack<T>& s) {
+  return output(os,s);
+}
+
+template<class T>
 inline ostream& operator<< (ostream& os, const vector<T>& v) {
   return output(os,v);
 }
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 7b856565a..02bfd4b03 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.26 2003/02/25 22:57:00 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.27 2003/03/18 17:30:25 schnetter Exp $
 
 #include <assert.h>
 
@@ -260,13 +260,15 @@ void dh<D>::recompose () {
 	      // (the restriction must not fill points that are used
 	      // to prolongate the boundaries)
               ibset recvs = intrf.contracted_for(intr) & intr;
-              const iblist& sendlist
-                = boxes[rl][c][ml].send_ref_bnd_fine[cc];
-              for (typename iblist::const_iterator sli = sendlist.begin();
-                   sli != sendlist.end();
-                   ++sli) {
-                recvs -= *sli;
+              for (int ccc=0; ccc<h.components(rl); ++ccc) {
+                const iblist& sendlist
+                  = boxes[rl][ccc][ml].send_ref_bnd_fine[cc];
+                for (typename iblist::const_iterator sli = sendlist.begin();
+                     sli != sendlist.end(); ++sli) {
+                  recvs -= *sli;
+                }
               }
+//               recvs.normalize();
               assert (recvs.setsize() <= 1);
               if (recvs.setsize() == 1) {
                 const ibbox recv = *recvs.begin();
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index f41971ae5..814b1bf48 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.13 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.14 2003/03/18 17:30:26 schnetter Exp $
 
 #ifndef VECT_HH
 #define VECT_HH
@@ -221,6 +221,13 @@ public:
   }
   
   // Non-modifying operators
+  vect replace (const int d, const T x) {
+    assert (d>=0 && d<D);
+    vect r;
+    for (int dd=0; dd<D; ++dd) r[dd]=dd==d?x:(*this)[dd];
+    return r;
+  }
+  
   vect operator+ () const {
     vect r;
     for (int d=0; d<D; ++d) r[d]=+r[d];