c4b9cdf1cab4214f668ac3615eba73fe19331e04
[charm.git] / src / libs / ck-libs / ParFUM / ParFUM_internals.h
1
2 /**
3  \defgroup ParFUM  ParFUM Unstructured Mesh Framework
4 */
5 /*@{*/
6
7
8 /*** 
9    ParFUM_internals.h
10
11    This file should contain ALL required header code that is 
12    internal to ParFUM, but should not be visible to users. This
13    includes all the old fem_mesh.h type files.
14
15    Adaptivity Operations added by: Nilesh Choudhury
16 */
17
18 #ifndef __PARFUM_INTERNALS_H
19 #define __PARFUM_INTERNALS_H
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <assert.h>
24 #include <vector>
25
26 #include "charm-api.h" /* for Fortran name mangling FTN_NAME */
27 #include "ckvector3d.h"
28 #include "tcharm.h"
29 #include "charm++.h"
30 #include "converse.h" /* for CmiGetArg */
31 #include "cklists.h"
32 #include "mpi.h"
33 #include "pup_mpi.h"
34 #define checkMPI pup_checkMPI
35
36 #include "idxl_layout.h"
37 #include "idxl.h"
38 #include "idxl_comm.h"
39
40 #include "ParFUM.decl.h"
41 #include "msa/msa.h"
42 #include "cklists.h"
43 #include "pup.h"
44
45 #if CMK_STL_USE_DOT_H /* Pre-standard C++ */
46 #  include <iostream.h>
47 #else /* ISO C++ */
48 #  include <iostream>
49    using namespace std;
50 #endif
51
52 #include "ParFUM_Adapt.decl.h"
53
54 /* USE of this extern may be a BUG */
55 extern CProxy_femMeshModify meshMod;
56
57 #define MAX_CHUNK 1000000000
58
59 //this macro turns off prints in adaptivity added by nilesh
60 #define FEM_SILENT
61
62
63 CtvExtern(FEM_Adapt_Algs *, _adaptAlgs);
64
65 /*
66   Charm++ Finite Element Framework:
67   C++ implementation file: Mesh representation and manipulation
68
69   This file lists the classes used to represent and manipulate 
70   Finite Element Meshes inside the FEM framework.
71
72   Orion Sky Lawlor, olawlor@acm.org, 1/3/2003
73 */
74
75 // Map the IDXL names to the old FEM names (FIXME: change all references, too)
76 typedef IDXL_Side FEM_Comm;
77 typedef IDXL_List FEM_Comm_List;
78 typedef IDXL_Rec FEM_Comm_Rec;
79
80 /// We want the FEM_Comm/IDXL_Side's to be accessible to *both* 
81 ///  FEM routines (via these data structures) and IDXL routines
82 ///  (via an idxl->addStatic registration).  Hence this class, which
83 ///  manages IDXL's view of FEM's data structures.
84 class FEM_Comm_Holder {
85   IDXL comm; //Our communication lists.
86   bool registered; //We're registered with IDXL
87   IDXL_Comm_t idx; //Our IDXL index, or -1 if we're unregistered
88   void registerIdx(IDXL_Chunk *c);
89  public:
90   FEM_Comm_Holder(FEM_Comm *sendComm, FEM_Comm *recvComm);
91   void pup(PUP::er &p);
92   ~FEM_Comm_Holder(void);
93         
94   /// Return our IDXL_Comm_t, registering with chunk if needed
95   inline IDXL_Comm_t getIndex(IDXL_Chunk *c) {
96     if (idx==-1) registerIdx(c);
97     return idx;
98   }
99 };
100
101
102
103 /// This datatype is how the framework stores symmetries internally.
104 ///   Each bit of this type describes a different symmetry.
105 ///   There must be enough bits to accomidate several simulatanious symmetries.
106 typedef unsigned char FEM_Symmetries_t;
107
108
109 /// Describes one kind of symmetry condition
110 class FEM_Sym_Desc : public PUP::able {
111  public:
112   virtual ~FEM_Sym_Desc();
113
114   /// Apply this symmetry to this location vector
115   virtual CkVector3d applyLoc(const CkVector3d &loc) const =0;
116         
117   /// Apply this symmetry to this relative (vel or acc) vector
118   virtual CkVector3d applyVec(const CkVector3d &vec) const =0;
119         
120   /// Allows Desc's to be pup'd via | operator:
121   friend inline void operator|(PUP::er &p,FEM_Sym_Desc &a) {a.pup(p);}
122   friend inline void operator|(PUP::er &p,FEM_Sym_Desc* &a) {
123     PUP::able *pa=a;  p(&pa);  a=(FEM_Sym_Desc *)pa;
124   }
125 };
126
127 /// Describes a linear-periodic (space shift) symmetry:
128 class FEM_Sym_Linear : public FEM_Sym_Desc {
129   CkVector3d shift; //Offset to add to locations
130  public:
131   FEM_Sym_Linear(const CkVector3d &shift_) :shift(shift_) {}
132   FEM_Sym_Linear(CkMigrateMessage *m) {}
133         
134   /// Apply this symmetry to this location vector
135   CkVector3d applyLoc(const CkVector3d &loc) const {return loc+shift;}
136         
137   /// Apply this symmetry to this relative (vel or acc) vector
138   virtual CkVector3d applyVec(const CkVector3d &vec) const {return vec;}
139         
140   virtual void pup(PUP::er &p);
141   PUPable_decl(FEM_Sym_Linear);
142 };
143
144 /**
145  * Describes all the different kinds of symmetries that apply to
146  * this mesh.
147  */
148 class FEM_Sym_List {
149   //This lists the different kinds of symmetry
150   CkPupAblePtrVec<FEM_Sym_Desc> sym; 
151         
152   FEM_Sym_List(const FEM_Sym_List &src); //NOT DEFINED: copy constructor
153  public:
154   FEM_Sym_List();
155   void operator=(const FEM_Sym_List &src); //Assignment operator
156   ~FEM_Sym_List();
157         
158   /// Add a new kind of symmetry to this list, returning
159   ///  the way objects with that symmetry should be marked.
160   FEM_Symmetries_t add(FEM_Sym_Desc *desc);
161         
162   /// Apply all the listed symmetries to this location
163   void applyLoc(CkVector3d *loc,FEM_Symmetries_t sym) const;
164         
165   /// Apply all the listed symmetries to this relative vector
166   void applyVec(CkVector3d *vec,FEM_Symmetries_t sym) const;
167         
168   void pup(PUP::er &p);
169 };
170
171 /// This is a simple 2D table.  The operations are mostly row-centric.
172 template <class T>
173 class BasicTable2d : public CkNoncopyable {
174  protected:
175   int rows; //Number of entries in table
176   int cols; //Size of each entry in table
177   T *table; //Data in table [rows * cols]
178  public:
179   BasicTable2d(T *src,int cols_,int rows_) 
180     :rows(rows_), cols(cols_), table(src) {}
181         
182   /// "size" of the table is the number of rows:
183   inline int size(void) const {return rows;}
184   /// Width of the table is the number of columns:
185   inline int width(void) const {return cols;}
186         
187   T *getData(void) {return table;}
188   const T *getData(void) const {return table;}
189         
190   // Element-by-element operations:
191   T operator() (int r,int c) const {return table[c+r*cols];}
192   T &operator() (int r,int c) {return table[c+r*cols];}
193         
194   // Row-by-row operations
195   ///Get a pointer to a row of the table:
196   inline T *getRow(int r) {return &table[r*cols];}
197   ///Get a const pointer to a row of the table:
198   inline const T *getRow(int r) const {return &table[r*cols];}
199   inline void getRow(int r,T *dest,T idxBase=0) const {
200     const T *src=getRow(r);
201     for (int c=0;c<cols;c++) dest[c]=src[c]+idxBase;
202   }
203   inline T *operator[](int r) {return getRow(r);}
204   inline const T *operator[](int r) const {return getRow(r);}
205   inline void setRow(int r,const T *src,T idxBase=0) {
206     T *dest=getRow(r);
207     for (int c=0;c<cols;c++) dest[c]=src[c]-idxBase;
208   }
209   inline void setRow(int r,T value) {
210     T *dest=getRow(r);
211     for (int c=0;c<cols;c++) dest[c]=value;
212   }
213
214   // These affect the entire table:
215   void set(const T *src,T idxBase=0) {
216     for (int r=0;r<rows;r++) 
217       for (int c=0;c<cols;c++)
218         table[c+r*cols]=src[c+r*cols]-idxBase;
219   }
220   void setTranspose(const T *srcT,int idxBase=0) {
221     for (int r=0;r<rows;r++) 
222       for (int c=0;c<cols;c++)
223         table[c+r*cols]=srcT[r+c*rows]-idxBase;
224   }
225   void get(T *dest,T idxBase=0) const {
226     for (int r=0;r<rows;r++) 
227       for (int c=0;c<cols;c++)
228         dest[c+r*cols]=table[c+r*cols]+idxBase;
229   }
230   void getTranspose(T *destT,int idxBase=0) const {
231     for (int r=0;r<rows;r++) 
232       for (int c=0;c<cols;c++)
233         destT[r+c*rows]=table[c+r*cols]+idxBase;
234   }
235   void set(T value) {
236     for (int r=0;r<rows;r++) setRow(r,value);
237   }
238   void setRowLen(int rows_) {
239     rows = rows_;
240   }
241 };
242
243 /// A heap-allocatable, resizable BasicTable2d.
244 /** A heap-allocatable, resizable BasicTable2d.
245  * To be stored here, T must not require a copy constructor.
246  */
247 template <class T>
248 class AllocTable2d : public BasicTable2d<T> {
249   ///Maximum number of rows that can be used without reallocation
250   int max;
251   ///Value to fill uninitialized regions with
252   T fill; 
253   /// the table that I allocated
254   T *allocTable;
255  public:
256   ///default constructor
257   AllocTable2d(int cols_=0,int rows_=0,T fill_=0) 
258     :BasicTable2d<T>(NULL,cols_,rows_), max(0), fill(fill_)             
259     {
260                         allocTable = NULL;
261       if (this->rows>0) allocate(this->rows);
262     }
263   ///default destructor
264   ~AllocTable2d() {if(allocTable != NULL){delete[] allocTable;}}
265   /// Make room for this many rows
266   void allocate(int rows_) { 
267     allocate(this->width(),rows_);
268   }
269   /// Make room for this many cols & rows
270   void allocate(int cols_,int rows_,int max_=0) { 
271     if (cols_==this->cols && rows_<max) {
272       //We have room--just update the size:
273       this->rows=rows_;
274       return;
275     }
276     if (max_==0) { //They gave no suggested size-- pick one:
277       if (rows_==this->rows+1) //Growing slowly: grab a little extra
278         max_=10+rows_+(rows_>>2); //  125% plus 10
279       else // for a big change, just go with the minimum needed: 
280         max_=rows_;
281     }
282
283     int oldRows=this->rows;
284     this->cols=cols_;
285     this->rows=rows_;
286     this->max=max_;
287     this->table=new T[max*this->cols];
288     //Preserve old table entries (FIXME: assumes old cols is unchanged)
289     int copyRows=0;
290     if (allocTable!=NULL) { 
291       copyRows=oldRows;
292       if (copyRows>max) copyRows=max;
293       memcpy(this->table,allocTable,sizeof(T)*this->cols*copyRows);
294       delete[] allocTable;
295     }else{
296       for (int r=copyRows;r<max;r++)
297         setRow(r,fill);
298     }
299     allocTable = this->table;
300   }
301         
302   /// Pup routine and operator|:
303   void pup(PUP::er &p) {
304     p|this->rows; p|this->cols;
305     if (this->table==NULL) allocate(this->rows);
306     p(this->table,this->rows*this->cols); //T better be a basic type, or this won't compile!
307   }
308   /// Pup only one element which is at index 'pupindx' in the table
309   void pupSingle(PUP::er &p, int pupindx) {
310     int startIdx = this->cols*pupindx;
311     for(int i=0; i<this->cols; i++) {
312       p|this->table[startIdx+i];
313     }
314   }
315
316   friend void operator|(PUP::er &p,AllocTable2d<T> &t) {t.pup(p);}
317
318   /// Add a row to the table (by analogy with std::vector):
319   T *push_back(void) {
320     if (this->rows>=max) 
321       { //Not already enough room for the new row:
322         int newMax=max+(max/4)+16; //Grow 25% longer
323         allocate(this->cols,this->rows,newMax);
324       }
325     this->rows++;
326     return getRow(this->rows-1);
327   }
328
329   /** to support replacement of attribute data by user supplied data
330       error checks have been performed at FEM_ATTRIB
331   */
332   void register_data(T *user,int len,int max){
333     if(allocTable != NULL){
334       delete [] allocTable;
335       allocTable = NULL;
336     }   
337     this->table = user;
338     this->rows = len;
339     this->max = max;
340   }
341 };
342
343 /// Return the human-readable version of this entity code, like "FEM_NODE".
344 ///  storage, which must be at least 80 bytes long, is used for
345 ///  non-static names, like the user tag "FEM_ELEM+2".
346 CDECL const char *FEM_Get_entity_name(int entity,char *storage);
347
348 /// Return the human-readable version of this attribute code, like "FEM_CONN".
349 ///  storage, which must be at least 80 bytes long, is used for
350 ///  non-static names, like the user tag "FEM_DATA+7".
351 CDECL const char *FEM_Get_attr_name(int attr,char *storage);
352
353
354 ///A user-visible 2D table attached to a FEM_Entity
355 /** Describes an FEM entity's "attribute"--a user-visible, user-settable
356  * 2D table.  Common FEM_Attributes include: user data associated with nodes,
357  *  the element-to-node connectivity array, etc.
358  */
359 class FEM_Attribute {
360   ///Owning entity (to get length, etc.)
361   FEM_Entity *e;
362   /// Ghost attribute, which has the same width and datatype as us (or 0)
363   FEM_Attribute *ghost;
364   /// My attribute code (e.g., FEM_DATA+7, FEM_CONN, etc.)
365   int attr;
366   ///Number of columns in our table of data (initially unknown)
367   int width;
368   ///Datatype of entries (initially unknown)
369   int datatype;
370   ///True if subclass allocate has been called.
371   bool allocated;
372   
373   ///Abort with a nice error message saying: 
374   /** Our <field> was previously set to <cur>; it cannot now be <next>*/
375   void bad(const char *field,bool forRead,int cur,int next,const char *caller) const;
376   
377  protected:
378   /**
379    * Allocate storage for at least length width-item records of type datatype.
380    * This routine is called after all three parameters are set,
381    * as a convenience for subclasses. 
382    */
383   virtual void allocate(int length,int width,int datatype) =0;
384  public:
385   FEM_Attribute(FEM_Entity *owner_,int myAttr_);
386   virtual void pup(PUP::er &p);
387   virtual void pupSingle(PUP::er &p, int pupindx);
388   virtual ~FEM_Attribute();
389         
390   /// Install this attribute as our ghost:
391   inline void setGhost(FEM_Attribute *ghost_) { ghost=ghost_;}
392         
393   /// Return true if we're a ghost
394   inline bool isGhost(void) const { return ghost!=NULL; }
395         
396   /// Return our attribute code
397   inline int getAttr(void) const {return attr;}
398
399   inline FEM_Entity *getEntity(void) {return e;}
400         
401   /// Return the number of rows in our table of data (0 if unknown).
402   /// This value is obtained directly from our owning Entity.
403   int getLength(void) const;
404   int getRealLength(void) const;
405
406   int getMax();
407         
408   /// Return the number of columns in our table of data (0 if unknown)
409   inline int getWidth(void) const { return width<0?0:width; }
410   inline int getRealWidth(void) const { return width; }
411         
412   /// Return our FEM_* datatype (-1 if unknown)
413   inline int getDatatype(void) const { return datatype; }
414         
415   /**
416    * Set our length (number of rows, or records) to this value.
417    * Default implementation calls setLength on our owning entity
418    * and tries to call allocate().
419    */
420   void setLength(int l,const char *caller="");
421         
422   /**
423    * Set our width (number of values per row) to this value.
424    * The default implementation sets width and tries to call allocate().
425    */
426   void setWidth(int w,const char *caller="");
427         
428   /**
429    * Set our datatype (e.g., FEM_INT, FEM_DOUBLE) to this value.
430    * The default implementation sets width and tries to call allocate().
431    */
432   void setDatatype(int dt,const char *caller="");
433         
434   /**
435    * Copy our width and datatype from this attribute.
436    * The default implementation calls setWidth and setDatatype.
437    * which should be enough for virtually any attribute.
438    */
439   virtual void copyShape(const FEM_Attribute &src);
440         
441   /// Check if all three of length, width, and datatype are set,
442   ///  but we're not yet allocated.  If so, call allocate; else ignore.
443   void tryAllocate(void);
444         
445   /// Our parent's length has changed: reallocate our storage
446   inline void reallocate(void) { allocated=false; tryAllocate(); }
447         
448   /// Return true if we've already had storage allocated.
449   inline bool isAllocated(void) const {return allocated;}
450         
451   /**
452    * Set our data to these (typically user-supplied, unchecked) values.
453    * Subclasses normally override this method as:
454    *    virtual void set( ...) {
455    *       super::set(...);
456    *       copy data from src.
457    *    }
458    */
459   virtual void set(const void *src, int firstItem,int length, 
460                    const IDXL_Layout &layout, const char *caller);
461         
462   /**
463    * Extract this quantity of user data.  Length and layout are
464    * parameter checked by the default implementation.
465    * Subclasses normally override this method as:
466    *    virtual void get( ...) {
467    *       super::get(...);
468    *       copy data to dest.
469    *    }
470    */
471   virtual void get(void *dest, int firstItem,int length,
472                    const IDXL_Layout &layout, const char *caller) const;
473         
474   /// Copy everything associated with src[srcEntity] into our dstEntity.
475   virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity) =0;
476
477   /** Register this user data for this attributre 
478       Length, layout etc are checked by the default implementaion
479   */
480   virtual void register_data(void *dest, int length,int max,
481                              const IDXL_Layout &layout, const char *caller);
482         
483 };
484 PUPmarshall(FEM_Attribute);
485
486 ///A single table of user data associated with an entity.
487 /**
488  * Describes a single table of user data associated with an entity.
489  * Since the data can be of any type, it is stored as chars.
490  */
491 class FEM_DataAttribute : public FEM_Attribute {
492   typedef FEM_Attribute super;
493   ///Non-NULL for getDatatype==FEM_BYTE
494   AllocTable2d<unsigned char> *char_data;
495   ///Non-NULL for getDatatype==FEM_INT
496   AllocTable2d<int> *int_data;
497   ///Non-NULL for getDatatype==FEM_FLOAT
498   AllocTable2d<float> *float_data;
499   ///Non-NULL for getDatatype==FEM_DOUBLE
500   AllocTable2d<double> *double_data;
501  protected:
502   virtual void allocate(int length,int width,int datatype);
503  public:
504   FEM_DataAttribute(FEM_Entity *owner,int myAttr);
505   virtual void pup(PUP::er &p);
506   virtual void pupSingle(PUP::er &p, int pupindx);
507   ~FEM_DataAttribute();
508   
509   AllocTable2d<unsigned char> &getChar(void) {return *char_data;}
510   const AllocTable2d<unsigned char> &getChar(void) const {return *char_data;}
511   
512   AllocTable2d<int> &getInt(void) {return *int_data;}
513   const AllocTable2d<int> &getInt(void) const {return *int_data;}
514     
515   AllocTable2d<double> &getDouble(void) {return *double_data;}
516   const AllocTable2d<double> &getDouble(void) const {return *double_data;}
517
518   AllocTable2d<float> &getFloat(void) {return *float_data;}
519   const AllocTable2d<float> &getFloat(void) const {return *float_data;}
520
521   
522   virtual void set(const void *src, int firstItem,int length, 
523                    const IDXL_Layout &layout, const char *caller);
524   
525   virtual void get(void *dest, int firstItem,int length, 
526                    const IDXL_Layout &layout, const char *caller) const;
527   
528   virtual void register_data(void *dest, int length,int max,
529                              const IDXL_Layout &layout, const char *caller);
530   
531         
532   /// Copy src[srcEntity] into our dstEntity.
533   virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
534   
535   /// extrapolate the values of a node from 2 others
536   /** Extrapolate the value of D from A and B (using the fraction 'frac')
537    */
538   void interpolate(int A,int B,int D,double frac);
539   /// extrapolate the values of a node from some others
540   /** Same as above except the the original nodes are send in an array
541    */
542   void interpolate(int *iNodes,int rNode,int k);
543 };
544 PUPmarshall(FEM_DataAttribute);
545
546 ///The FEM_Attribute is of type integer indices
547 /**
548  * This table maps an entity to a set of integer indices.
549  * The canonical example of this is the element-node connectivity array.
550  */
551 class FEM_IndexAttribute : public FEM_Attribute {
552  public:
553   /// Checks incoming indices for validity.
554   class Checker {
555   public:
556     virtual ~Checker();
557     
558     /**
559      * Check this (newly set) row of our table for validity.
560      * You're expected to abort or throw or exit if something is wrong.
561      */
562     virtual void check(int row,const BasicTable2d<int> &table,
563                        const char *caller) const =0;
564   };
565  private:
566   typedef FEM_Attribute super;
567   AllocTable2d<int> idx;
568   ///Checks indices (or NULL). This attribute will destroy this object
569   Checker *checker;
570  protected:
571   virtual void allocate(int length,int width,int datatype);
572  public:
573   FEM_IndexAttribute(FEM_Entity *owner,int myAttr, Checker *checker_=NULL);
574   virtual void pup(PUP::er &p);
575   virtual void pupSingle(PUP::er &p, int pupindx);
576   ~FEM_IndexAttribute();
577   
578   AllocTable2d<int> &get(void) {return idx;}
579   const AllocTable2d<int> &get(void) const {return idx;}
580         
581   virtual void set(const void *src, int firstItem,int length, 
582                    const IDXL_Layout &layout, const char *caller);
583         
584   virtual void get(void *dest, int firstItem,int length,
585                    const IDXL_Layout &layout, const char *caller) const;
586
587   virtual void register_data(void *dest, int length, int max,
588                              const IDXL_Layout &layout, const char *caller);
589         
590   /// Copy src[srcEntity] into our dstEntity.
591   virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
592 };
593 PUPmarshall(FEM_IndexAttribute);
594
595 ///The FEM_Attribute is a variable set of integer indices
596 /**
597   This table maps an entity to a list of integer indices 
598   of unknown length.
599   The node to element adjacency array is an example, where a node
600   is mapped to a list of element indices of unknown length.
601 */
602 class FEM_VarIndexAttribute : public FEM_Attribute{
603  public:
604   class ID{
605   public:
606     ///type is negative for ghost elements
607     int type;
608     ///id refers to the index in the entity list
609     int id;
610
611     ///default constructor
612     ID(){
613       type=-1;
614       id = -1;
615     };
616     ///constructor - initializer
617     ID(int _type,int _id){
618       if(_id < 0) {
619         type = -(_type+1);
620         id = FEM_To_ghost_index(_id);
621       }
622       else {
623         type = _type;
624         id = _id;
625       }
626     };
627     bool operator ==(const ID &rhs)const {
628       return (type == rhs.type) && (id == rhs.id);
629     }
630     const ID& operator =(const ID &rhs) {
631       type = rhs.type;
632       id = rhs.id;
633       return *this;
634     }
635     virtual void pup(PUP::er &p){
636       p | type;
637       p | id;
638     };
639
640     static ID createNodeID(int type,int node){
641       ID temp(type, node);
642       return temp;
643     }
644     int getSignedId() {
645       if(type<0){
646         return FEM_From_ghost_index(id);
647       }
648       else return id;
649     }
650   };
651  private:
652   typedef FEM_Attribute super;
653   CkVec<CkVec<ID> > idx;
654   int oldlength;
655  protected:
656   virtual void allocate(int _length,int _width,int _datatype){
657     if(_length > oldlength){
658       setWidth(1,"allocate"); //there is 1 vector per entity 
659       oldlength = _length*2;
660       idx.reserve(oldlength);
661       for(int i=idx.size();i<oldlength;i++){
662         CkVec<ID> tempVec;
663         idx.insert(i,tempVec);
664       }
665     }
666   };
667  public:
668   FEM_VarIndexAttribute(FEM_Entity *owner,int myAttr);
669   ~FEM_VarIndexAttribute(){};
670   virtual void pup(PUP::er &p);
671   virtual void pupSingle(PUP::er &p, int pupindx);
672   CkVec<CkVec<ID> > &get(){return idx;};
673   const CkVec<CkVec<ID> > &get() const {return idx;};
674         
675   virtual void set(const void *src,int firstItem,int length,
676                    const IDXL_Layout &layout,const char *caller);
677
678   virtual void get(void *dest, int firstItem,int length,
679                    const IDXL_Layout &layout, const char *caller) const;
680         
681   virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
682
683   int findInRow(int row,const ID &data);
684         
685   void print();
686 };
687
688 ///FEM_Entity describes an entire entity (nodes/elements/sparse) which contains attributes
689 /**
690  * Describes an entire class of "entities"--nodes, elements, or sparse 
691  *  data records. Basically consists of a length and a set of
692  *  FEM_Attributes. 
693  */
694 class FEM_Entity {
695   ///symmetries
696   typedef CkVec<FEM_Symmetries_t> sym_t;
697   /// Number of entities of our type currently in the mesh
698   int length;
699   /// Maximum number of entities of our type in the mesh that will be allowed
700   int max;
701   /// should be a function pointer to the actual resize function later
702   FEM_Mesh_alloc_fn resize;
703   /// arguments to the resize function
704   void *args;
705         
706   /**
707    * This is our main list of attributes-- everything about each of 
708    * our entities is in this list of attributes.  This list is searched
709    * by our "lookup" method and maintained by our subclasses "create" 
710    * method and calls to our "add" method.
711    * 
712    * It's a little funky having the superclass keep pointers to subclass
713    * objects (like element connectivity), but very nice to be able to
714    * easily loop over everything associated with an entity.
715    */
716   CkVec<FEM_Attribute *> attributes;
717         
718   /**
719    * Coordinates of each entity, from FEM_COORD.
720    * Datatype is always FEM_DOUBLE, width is always 2 or 3.
721    *  If NULL, coordinates are unknown.
722    */
723   FEM_DataAttribute *coord; 
724   void allocateCoord(void);
725         
726   /**
727    * Symmetries of each entity, from FEM_SYMMETRIES.  This bitvector per
728    * entity indicates which symmetry conditions the entity belongs to.
729    * Datatype is always FEM_BYTE (same as FEM_Symmetries_t), width 
730    * is always 1.  If NULL, all the symmetries are 0.
731    */
732   FEM_DataAttribute *sym; 
733   void allocateSym(void);
734         
735   /**
736    * Global numbers of each entity, from FEM_GLOBALNO.
737    * If NULL, the global numbers are unknown.
738    */
739   FEM_IndexAttribute *globalno;
740   void allocateGlobalno(void);
741         
742   /**
743     used to allocate the integer array for storing the boundary
744     values associated with an entity. 
745   */
746   void allocateBoundary();
747
748   /// Mesh sizing attribute for elements
749   /** Specifies a double edge length target for the mesh at each 
750       element; used in adaptivity algorithms */
751   FEM_DataAttribute *meshSizing;
752   void allocateMeshSizing(void);
753
754   /**
755     used to allocate the char array for storing whether each entity is valid
756     When a node/element is deleted the flag in the valid table is set to 0.
757     Additionally, we keep track of the first and last occurence in the array of
758     invalid indices. This should make searching for slots to reuse quicker.
759   */
760   FEM_DataAttribute *valid;
761   int first_invalid, last_invalid;
762   /** This is a list of invalid entities.. its a last in first out queue,
763    * which can be implemented as simple array, where we add and remove from
764    * the last element, so that we always have a compact list, and we
765    * allocate only when we need to
766    */
767   int* invalidList;
768   ///length of invalid list 
769   int invalidListLen;
770   ///length of allocated invalid list
771   int invalidListAllLen;
772   
773  protected:
774   /**
775    * lookup of this attribute code has failed: check if it needs
776    * to be demand-created.  Subclasses should override this method
777    * to recognize a request for, and add their own attributes;
778    * otherwise call the default implementation.
779    *
780    * Every call to create must result in either a call to add()
781    * or a call to the superclass; but not both.
782    *
783    * Any entity with optional fields, that are created on demand, will
784    * have to override this method.  Entities with fixed fields that are
785    * known beforehand should just call add() from their constructor.
786    */
787   virtual void create(int attr,const char *caller);
788         
789   /// Add this attribute to this kind of Entity.
790   /// This superclass is responsible for eventually deleting the attribute.
791   /// This class also attaches the ghost attribute, so be sure to 
792   ///   call add before manipulating the attribute.
793   void add(FEM_Attribute *attribute);
794  public:
795
796   FEM_Entity *ghost; // Our ghost entity type, or NULL if we're the ghost
797                 
798   FEM_Comm ghostSend; //Non-ghosts we send out (only set for real entities)
799   FEM_Comm ghostRecv; //Ghosts we recv into (only set for ghost entities)
800   
801         FEM_Comm_Holder ghostIDXL; //IDXL interface
802
803   FEM_Entity(FEM_Entity *ghost_); //Default constructor
804   void pup(PUP::er &p);
805   virtual ~FEM_Entity();
806         
807   /// Return true if we're a ghost
808   bool isGhost(void) const {return ghost==NULL;}
809         
810   /// Switch from this, a real entity, to the ghosts:
811   FEM_Entity *getGhost(void) {return ghost;}
812   const FEM_Entity *getGhost(void) const {return ghost;}
813         
814         //should only be called on non-ghost elements that have ghosts.
815         //empty its corresponding ghost entity and clear the idxls
816         void clearGhost();
817         
818   /// Return the number of entities of this type
819   inline int size(void) const {return length==-1?0:length;}
820   inline int realsize(void) const {return length;}
821
822   // return the maximum size 
823   inline int getMax() { if(max > 0) return max; else return length;}
824         
825   /// Return the human-readable name of this entity type, like "node"
826   virtual const char *getName(void) const =0;
827         
828   /// Copy all our attributes' widths and data types from this entity.
829   void copyShape(const FEM_Entity &src);
830         
831   /** 
832    *  The user is setting this many entities.  This reallocates
833    * all existing attributes to make room for the new entities.
834    */
835   void setLength(int newlen, bool f=false);
836
837   /** Support for registration API 
838    *  Set the current length and maximum length for this entity. 
839    *  If the current length exceeds the maximum length a resize
840    *  method is called .
841    */
842   void setMaxLength(int newLen,int newMaxLen,void *args,FEM_Mesh_alloc_fn fn);
843         
844   /// Copy everything associated with src[srcEntity] into our dstEntity.
845   /// dstEntity must have already been allocated, e.g., with setLength.
846   void copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity);
847
848   /// Add room for one more entity, with initial values from src[srcEntity],
849   /// and return the new entity's index.
850   int push_back(const FEM_Entity &src,int srcEntity);
851         
852   /**
853    * Find this attribute (from an FEM_ATTR code) of this entity, or 
854    * create the entity (using the create method below) or abort if it's
855    * not found.
856    */
857   FEM_Attribute *lookup(int attr,const char *caller);
858
859         
860   /**
861    * Get a list of the attribute numbers for this entity.
862    */
863   int getAttrs(int *attrs) const {
864     int len=attributes.size();
865     for (int i=0;i<len;i++) 
866       attrs[i]=attributes[i]->getAttr();
867     return len;
868   }
869         
870   /**
871    * Allocate or Modify the FEM_IS_VALID attribute data
872    */
873   void allocateValid();
874   void set_valid(int idx, bool isNode);
875   void set_invalid(int idx, bool isNode);
876   int is_valid(int idx);         // will fail assertions for out of range indices
877   int is_valid_any_idx(int idx); // will not fail assertions for out of range indices
878   int count_valid();
879   int get_next_invalid(FEM_Mesh *m=NULL, bool isNode=false, bool isGhost=false);// the arguments are not really needed but Nilesh added them when he wrote a messy version and did not remove them when he fixed the implementation. Since changing the uses was too painful the default arguments were added.
880   virtual bool hasConn(int idx)=0;
881   /**
882    * Set the coordinates for a single item
883    */
884   void set_coord(int idx, double x, double y);
885   void set_coord(int idx, double x, double y, double z);
886
887   /** Expose the attribute vector for refining. breaks modularity but more efficient
888   */
889   CkVec<FEM_Attribute *>* getAttrVec(){
890     return &attributes;
891   }
892         
893   // Stupidest possible coordinate access
894   inline FEM_DataAttribute *getCoord(void) {return coord;}
895   inline const FEM_DataAttribute *getCoord(void) const {return coord;}
896         
897   //Symmetry array access:
898   const FEM_Symmetries_t *getSymmetries(void) const {
899     if (sym==NULL) return NULL;
900     else return (const FEM_Symmetries_t *)sym->getChar()[0];
901   }
902   FEM_Symmetries_t getSymmetries(int r) const { 
903     if (sym==NULL) return FEM_Symmetries_t(0);
904     else return sym->getChar()(r,0);
905   }
906   void setSymmetries(int r,FEM_Symmetries_t s);
907         
908   //Global numbering array access
909   bool hasGlobalno(void) const {return globalno!=0;}
910   int getGlobalno(int r) const {
911     if (globalno==0) return -1; // Unknown global number
912     return globalno->get()(r,0);
913   }
914   void setGlobalno(int r,int g);
915   void setAscendingGlobalno(void);
916   void setAscendingGlobalno(int base);
917   void copyOldGlobalno(const FEM_Entity &e);
918
919   // Mesh sizing array access
920   bool hasMeshSizing(void) const {return meshSizing!=0;}
921   double getMeshSizing(int r); 
922   void setMeshSizing(int r,double s);
923   void setMeshSizing(double *sf);
924         
925   //Ghost comm. list access
926   FEM_Comm &setGhostSend(void) { return ghostSend; }
927   const FEM_Comm &getGhostSend(void) const { return ghostSend; }
928   FEM_Comm &setGhostRecv(void) { 
929     if (ghost==NULL) return ghostRecv;
930     else return ghost->ghostRecv; 
931   }
932   const FEM_Comm &getGhostRecv(void) const { return ghost->ghostRecv; }
933         
934   void addVarIndexAttribute(int code){
935     FEM_VarIndexAttribute *varAttribute = new FEM_VarIndexAttribute(this,code);
936     add(varAttribute);
937   }
938         
939   void print(const char *type,const IDXL_Print_Map &map);
940 };
941 PUPmarshall(FEM_Entity);
942
943 // Now that we have FEM_Entity, we can define attribute lenth, as entity length
944 inline int FEM_Attribute::getLength(void) const { return e->size(); }
945 inline int FEM_Attribute::getRealLength(void) const { return e->realsize(); }
946 inline int FEM_Attribute::getMax(){ return e->getMax();}
947
948
949 ///FEM_Node is a type of FEM_Entity, which refers to nodes
950 /**
951  * Describes a set of FEM Nodes--the FEM_NODE entity type. 
952  * Unlike elements, nodes have no connectivity; but they do have
953  * special shared-nodes communications and a "primary" flag.
954  */
955 class FEM_Node : public FEM_Entity {
956   typedef FEM_Entity super;
957
958   /**
959    * primary flag, from FEM_NODE_PRIMARY, indicates this chunk "owns" this node.
960    * Datatype is always FEM_BYTE (we need an FEM_BIT!), width is always 1,
961    * since there's only one such flag per node.
962    */
963   FEM_DataAttribute *primary; 
964   void allocatePrimary(void);
965         
966   void allocateElemAdjacency();
967   void allocateNodeAdjacency();
968
969   FEM_VarIndexAttribute *elemAdjacency; ///< stores the node to element adjacency vector 
970   FEM_VarIndexAttribute *nodeAdjacency; ///< stores the node to node adjacency vector 
971   typedef FEM_VarIndexAttribute::ID var_id;
972  protected:
973   virtual void create(int attr,const char *caller);
974  public:
975   FEM_Comm shared; ///<Shared nodes
976   FEM_Comm_Holder sharedIDXL; ///<IDXL interface to shared nodes
977         
978   FEM_Node(FEM_Node *ghost_);
979   void pup(PUP::er &p);
980   ~FEM_Node();
981         
982   virtual const char *getName(void) const;
983         
984   inline bool getPrimary(int nodeNo) const {
985     if (primary==NULL) return true; //Everything must be primary
986     else return primary->getChar()(nodeNo,0);
987   }
988   inline void setPrimary(int nodeNo,bool isPrimary) {
989     if (primary==NULL) allocatePrimary();
990     primary->getChar()(nodeNo,0)=isPrimary;
991   }
992   void fillElemAdjacencyTable(int type,const FEM_Elem &elem);
993   void setElemAdjacency(int type,const FEM_Elem &elem);
994   void fillNodeAdjacency(const FEM_Elem &elem);
995   void setNodeAdjacency(const FEM_Elem &elem);
996   void fillNodeAdjacencyForElement(int node,int nodesPerElem,const int *conn,FEM_VarIndexAttribute *adjacencyAttr);
997   bool hasConn(int idx);
998   void print(const char *type,const IDXL_Print_Map &map);
999 };
1000 PUPmarshall(FEM_Node);
1001
1002
1003 ///FEM_Elem is a type of FEM_Entity, which refers to elems
1004 /**
1005  * Describes one kind of FEM elements--the FEM_ELEM entity type.
1006  * Elements are typically the central user-visible object in a FEM
1007  * computation.
1008  */
1009 class FEM_Elem:public FEM_Entity {
1010   typedef FEM_Entity super;
1011  protected:
1012   typedef AllocTable2d<int> conn_t;
1013
1014   int tuplesPerElem;
1015
1016   // The following are attributes that will commonly be used:
1017   FEM_IndexAttribute *conn;                ///< FEM_CONN attribute: element-to-node mapping 
1018   FEM_IndexAttribute *elemAdjacency;       ///< FEM_ELEM_ELEM_ADJACENCY attribute
1019   FEM_IndexAttribute *elemAdjacencyTypes;  ///< FEM_ELEM_ELEM_ADJ_TYPES attribute
1020   
1021  public:
1022   
1023   FEM_Elem(const FEM_Mesh &mesh_, FEM_Elem *ghost_);
1024   void pup(PUP::er &p);
1025   ~FEM_Elem();
1026   
1027   virtual const char *getName(void) const;
1028   
1029   // Directly access our connectivity table:
1030   inline conn_t &setConn(void) {return conn->get();}
1031   inline const conn_t &getConn(void) const {return conn->get();}
1032   
1033   void print(const char *type,const IDXL_Print_Map &map);
1034   
1035   void create(int attr,const char *caller);
1036   
1037   void allocateElemAdjacency();
1038   
1039   FEM_IndexAttribute *getElemAdjacency(){return elemAdjacency;}
1040   
1041   // Backward compatability routines:
1042   int getConn(int elem,int nodeNo) const {return conn->get()(elem,nodeNo);}
1043   int getNodesPer(void) const {return conn->get().width();}
1044   int *connFor(int i) {return conn->get().getRow(i);}
1045   const int *connFor(int i) const {return conn->get().getRow(i);}
1046   void connIs(int i,const int *src) {conn->get().setRow(i,src);}
1047   bool hasConn(int idx);
1048 };
1049 PUPmarshall(FEM_Elem);
1050
1051
1052
1053 ///FEM_Sparse is a type of FEM_Entity, which refers to edges
1054 /**
1055  * FEM_Sparse describes a set of records of sparse data that are all the
1056  * same size and all associated with the same number of nodes.
1057  * Sparse data is associated with some subset of the nodes in the mesh,
1058  * and gets copied to every chunk that has all those nodes.  The canonical
1059  * use of sparse data is to describe boundary conditions.
1060  */
1061 class FEM_Sparse : public FEM_Elem {
1062   typedef FEM_Elem super;
1063   typedef AllocTable2d<int> elem_t;
1064         
1065   /**
1066    * elem, from FEM_SPARSE_ELEM, is an optional (that is, possibly NULL) 
1067    * array which changes the partitioning of sparse entities: if non-NULL,
1068    * sparse entity t lives on the same chunk as FEM_ELEM+elem[2*t]
1069    * local number elem[2*t+1].
1070    *
1071    * This attribute's width is always 2.
1072    */
1073   FEM_IndexAttribute *elem; //FEM_SPARSE_ELEM attribute
1074   void allocateElem(void);
1075   const FEM_Mesh &mesh; //Reference to enclosing mesh, for error checking
1076  protected:
1077   virtual void create(int attr,const char *caller);
1078  public:
1079   FEM_Sparse(const FEM_Mesh &mesh_, FEM_Sparse *ghost_);
1080   void pup(PUP::er &p);
1081   virtual ~FEM_Sparse();
1082         
1083   virtual const char *getName(void) const;
1084         
1085   /// Return true if we have an element partitioning table
1086   bool hasElements(void) const {return elem!=NULL;}
1087         
1088   /// Directly access our element partitioning table (e.g., for re-numbering)
1089   inline elem_t &setElem(void) {return elem->get();}
1090   inline const elem_t &getElem(void) const {return elem->get();}
1091 };
1092 PUPmarshall(FEM_Sparse);
1093
1094 /** Describes a user function to pup a piece of mesh data 
1095  */
1096 class FEM_Userdata_pupfn {
1097   FEM_Userdata_fn fn;
1098   void *data;
1099  public:
1100   FEM_Userdata_pupfn(FEM_Userdata_fn fn_,void *data_)
1101     :fn(fn_), data(data_) {}
1102   /// Call user's pup routine using this PUP::er
1103   void pup(PUP::er &p) {
1104     (fn)((pup_er)&p,data);
1105   }
1106 };
1107
1108 /// Describes one piece of generic unassociated mesh data.
1109 class FEM_Userdata_item {
1110   CkVec<char> data; ///< Serialized data from user's pup routine
1111  public:
1112   int tag; //User-assigned identifier
1113   FEM_Userdata_item(int tag_=-1) {tag=tag_;}
1114         
1115   /// Return true if we have stored data.
1116   bool hasStored(void) const {return data.size()!=0;}
1117         
1118   /// Store this userdata inside us:
1119   void store(FEM_Userdata_pupfn &f) {
1120     data.resize(PUP::size(f));
1121     PUP::toMemBuf(f,&data[0],data.size());
1122   }
1123   /// Extract this userdata from our stored copy:
1124   void restore(FEM_Userdata_pupfn &f) {
1125     PUP::fromMemBuf(f,&data[0],data.size());
1126   }
1127         
1128   /// Save our stored data to this PUP::er
1129   void pup(PUP::er &p) {
1130     p|tag;
1131     p|data;
1132   }
1133 };
1134
1135 /// Describes all the unassociated data in a mesh. 
1136 class FEM_Userdata_list {
1137   CkVec<FEM_Userdata_item> list;
1138  public:
1139   FEM_Userdata_item &find(int tag) {
1140     for (int i=0;i<list.size();i++)
1141       if (list[i].tag==tag)
1142         return list[i];
1143     // It's not in the list-- add it.
1144     list.push_back(FEM_Userdata_item(tag));
1145     return list[list.size()-1];
1146   }
1147   int size(){
1148     return list.size();
1149   }
1150   void pup(PUP::er &p) {p|list;}
1151 };
1152
1153
1154 void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType);
1155 void FEM_Is_NULL(const char *caller,const char *entityType,int type);
1156
1157 /**
1158  * This class describes several different types of a certain kind 
1159  * of entity.  For example, there might be a FEM_Entity_Types<FEM_Elem>
1160  * that lists the different kinds of element.
1161  *
1162  * This class exists to provide a nice "demand-creation" semantics,
1163  * where the user assigns array indices (the e in FEM_ELEM+e),
1164  * so we don't know that we're setting the first copy when we set it.
1165  *
1166  * It's not clear this class has any right to exist--it should either be
1167  * folded into FEM_Mesh or generalized into a "userNumberedVec" or some such.
1168  */
1169 template <class T>
1170 class FEM_Entity_Types {
1171   CkVec<T *> types; // Our main storage for different entity types
1172   const FEM_Mesh &mesh;
1173   const char *name; //FEM_SPARSE or FEM_ELEM, or some such.
1174
1175  public:
1176   FEM_Entity_Types(const FEM_Mesh &mesh_,const char *name_) 
1177     :mesh(mesh_), name(name_) {}
1178   void pup(PUP::er &p) { 
1179     // Can't just use p|types, because T has a funky constructor:
1180     int n=types.size();
1181     p|n;
1182     for (int i=0;i<n;i++) {
1183       int isNULL=0;
1184       if (!p.isUnpacking()) isNULL=(types[i]==NULL);
1185       p|isNULL;
1186       if (!isNULL) set(i,"pup").pup(p);
1187     }
1188   }
1189   ~FEM_Entity_Types() {
1190     for (int i=0;i<types.size();i++)
1191       if (types[i]) delete types[i];
1192   }
1193         
1194   /// Return the number of different entity types
1195   inline int size(void) const {return types.size();}
1196         
1197   /// Return a read-only copy of this type, or else abort if type isn't set.
1198   const T &get(int type,const char *caller="") const {
1199     FEM_Index_Check(caller,name,type,types.size());
1200     const T *ret=types[type];
1201     if (ret==NULL) FEM_Is_NULL(caller,name,type);
1202     return *ret;
1203   }
1204         
1205   /// Return true if we have a type t, and false otherwise
1206   bool has(int type) const { 
1207     if (type>=types.size()) return false;
1208     return types[type]!=NULL; 
1209   }
1210         
1211   /// Return a writable copy of this type, calling new T(mesh) if it's not there
1212   T &set(int type,const char *caller="") {
1213     if (type<0) FEM_Index_Check(caller,name,type,types.size());
1214     while (types.size()<=type) types.push_back(NULL); //Make room for new type
1215     if (types[type]==NULL) { //Have to allocate a new T:
1216       T *ghost=new T(mesh,NULL);
1217       types[type]=new T(mesh,ghost);
1218     }
1219     return *types[type];
1220   }
1221         
1222   /// Read-only and write-only operator[]'s:
1223   inline T &operator[] (int type) { return set(type); }
1224   inline const T &operator[] (int type) const { return get(type); }
1225
1226 };
1227
1228 /// Map fortran element (>=1) or node (0) marker to C version (>=1, -1)
1229 inline int zeroToMinusOne(int i) {
1230   if (i==0) return -1;
1231   return i;
1232 }
1233
1234
1235 class ParFUMShadowArray;
1236
1237 ///A FEM_Mesh is a collection of entities.
1238 /**
1239  * This class describes all the nodes and elements in
1240  * a finite-element mesh or submesh.
1241  */
1242 class FEM_Mesh : public CkNoncopyable {
1243   /// The symmetries in the mesh
1244   FEM_Sym_List symList;
1245   bool m_isSetting;
1246   
1247   void checkElemType(int elType,const char *caller) const;
1248   void checkSparseType(int uniqueID,const char *caller) const; 
1249
1250  public:
1251   femMeshModify *fmMM;
1252   ParFUMShadowArray *parfumSA;
1253   bool lastLayerSet;
1254   FEM_ElemAdj_Layer* lastElemAdjLayer;
1255   void setFemMeshModify(femMeshModify *m);
1256   void setParfumSA(ParFUMShadowArray *m);
1257   
1258   FEM_Mesh();
1259   void pup(PUP::er &p); //For migration
1260   ~FEM_Mesh();
1261
1262   /// The nodes in this mesh:
1263   FEM_Node node; 
1264   
1265   /// The different element types in this mesh:
1266   FEM_Entity_Types<FEM_Elem> elem;
1267   
1268   /// The different sparse types in this mesh:
1269   FEM_Entity_Types<FEM_Sparse> sparse;
1270   
1271   /// The unassociated user data for this mesh:
1272   FEM_Userdata_list udata;
1273   
1274   /// The symmetries that apply to this mesh:
1275   void setSymList(const FEM_Sym_List &src) {symList=src;}
1276   const FEM_Sym_List &getSymList(void) const {return symList;}
1277   
1278   /// Set up the "shape" of our fields-- the number of element types,
1279   /// the datatypes for user data, etc--based on this mesh.
1280   void copyShape(const FEM_Mesh &src);
1281
1282   // Get the fem mesh modification object associated with this mesh or partition  
1283   femMeshModify *getfmMM();
1284
1285   //Return this type of element, given an element type
1286   FEM_Entity &setCount(int elTypeOrMinusOne) {
1287     if (elTypeOrMinusOne==-1) return node;
1288     else return elem[chkET(elTypeOrMinusOne)];
1289   }
1290   const FEM_Entity &getCount(int elTypeOrMinusOne) const {
1291     if (elTypeOrMinusOne==-1) return node;
1292     else return elem[chkET(elTypeOrMinusOne)];
1293   }
1294   FEM_Elem &setElem(int elType) {return elem[chkET(elType)];}
1295   const FEM_Elem &getElem(int elType) const {return elem[chkET(elType)];}
1296   int chkET(int elType) const; //Check this element type-- abort if it's bad
1297   
1298   /// Look up this FEM_Entity type in this mesh, or abort if it's not valid.
1299   FEM_Entity *lookup(int entity,const char *caller);
1300   const FEM_Entity *lookup(int entity,const char *caller) const;
1301   
1302   /// Set/get direction control:
1303   inline bool isSetting(void) const {return m_isSetting;}
1304   void becomeSetting(void) {m_isSetting=true;}
1305   void becomeGetting(void) {m_isSetting=false;}
1306   
1307   int nElems() const //Return total number of elements (of all types)
1308     {return nElems(elem.size());}
1309   /// Return the total number of elements before type t
1310   int nElems(int t) const;
1311   /// Return the "global" number of element elNo of type elType.
1312   int getGlobalElem(int elType,int elNo) const;
1313   /// Set our global numbers as 0...n-1 for nodes, elements, and sparse
1314   void setAscendingGlobalno(void);
1315   ///   The global numbers for elements runs across different types
1316   void setAbsoluteGlobalno();
1317   
1318   void copyOldGlobalno(const FEM_Mesh &m);
1319   void print(int idxBase);//Write a human-readable description to CkPrintf
1320   /// Extract a list of our entities:
1321   int getEntities(int *entites);
1322   
1323   /**
1324          * clearing the idxl and data for shared nodes and ghost nodes and shared elements
1325          * 
1326          * */
1327
1328         void clearSharedNodes();
1329         void clearGhostNodes();
1330         void clearGhostElems();
1331   
1332   /********** New methods ***********/
1333   /*
1334     This method creates the mapping from a node to all the elements that are 
1335     incident on it . It assumes the presence of one layer of ghost nodes that
1336     share a node.
1337   */
1338   void createNodeElemAdj();
1339   void createNodeNodeAdj();
1340   void createElemElemAdj();
1341   
1342   FEM_ElemAdj_Layer *getElemAdjLayer(void);
1343   
1344   // Terry's adjacency accessors & modifiers
1345
1346   //  ------- Element-to-element: preserve initial ordering relative to nodes
1347   /// Place all of element e's adjacent elements in neighbors; assumes
1348   /// neighbors allocated to correct size
1349   void e2e_getAll(int e, int *neighborss, int etype=0);
1350   /// Given id of element e, return the id of the idx-th adjacent element
1351   int e2e_getNbr(int e, short idx, int etype=0);
1352   /// Given id of element e and id of another element nbr, return i such that
1353   /// nbr is the i-th element adjacent to e
1354   int e2e_getIndex(int e, int nbr, int etype=0);
1355   /// Set the element adjacencies of element e to neighbors; assumes neighbors 
1356   /// has the correct size
1357   void e2e_setAll(int e, int *neighbors, int etype=0);
1358   /// Set the idx-th element adjacent to e to be newElem
1359   void e2e_setIndex(int e, short idx, int newElem, int etype=0);
1360   /// Find element oldNbr in e's adjacent elements and replace with newNbr
1361   void e2e_replace(int e, int oldNbr, int newNbr, int etype=0);
1362   /// Remove all neighboring elements in adjacency
1363   void e2e_removeAll(int e, int etype=0);
1364
1365
1366   //  ------- Element-to-node: preserve initial ordering
1367   /// Place all of element e's adjacent nodes in adjnodes; assumes
1368   /// adjnodes allocated to correct size
1369   void e2n_getAll(int e, int *adjnodes, int etype=0);
1370   /// Given id of element e, return the id of the idx-th adjacent node
1371   int e2n_getNode(int e, short idx, int etype=0);
1372   /// Given id of element e and id of a node n, return i such that
1373   /// n is the i-th node adjacent to e
1374   short e2n_getIndex(int e, int n, int etype=0);
1375   /// Set the node adjacencies of element e to adjnodes; assumes adjnodes 
1376   /// has the correct size
1377   void e2n_setAll(int e, int *adjnodes, int etype=0);
1378   /// Set the idx-th node adjacent to e to be newNode
1379   void e2n_setIndex(int e, short idx, int newNode, int etype=0);
1380   /// Find node oldNode in e's adjacent ndoes and replace with newNode
1381   void e2n_replace(int e, int oldNode, int newNode, int etype=0);
1382   /// Replace all entries with -1
1383   void e2n_removeAll(int e, int etype=0);
1384
1385
1386   //  ------- Node-to-node
1387   int n2n_getLength(int n);
1388   /// Place all of node n's adjacent nodes in adjnodes and the resulting 
1389   /// length of adjnodes in sz; assumes adjnodes is not allocated, but sz is
1390   void n2n_getAll(int n, int *&adjnodes, int &sz);
1391   /// Adds newNode to node n's node adjacency list
1392   void n2n_add(int n, int newNode);
1393   /// Removes oldNode from n's node adjacency list
1394   void n2n_remove(int n, int oldNode);
1395   /// Finds oldNode in n's node adjacency list, and replaces it with newNode
1396   void n2n_replace(int n, int oldNode, int newNode);
1397   /// Remove all nodes from n's node adjacency list
1398   void n2n_removeAll(int n);
1399   /// Is queryNode in node n's adjacency vector?
1400   int n2n_exists(int n, int queryNode);
1401
1402   //  ------- Node-to-element
1403   int n2e_getLength(int n);
1404   /// Place all of node n's adjacent elements in adjelements and the resulting 
1405   /// length of adjelements in sz; assumes adjelements is not allocated, 
1406   /// but sz is
1407   void n2e_getAll(int n, int *&adjelements, int &sz);
1408   /// Adds newElem to node n's element adjacency list
1409   void n2e_add(int n, int newElem);
1410   /// Removes oldElem from n's element adjacency list
1411   void n2e_remove(int n, int oldElem);
1412   /// Finds oldElem in n's element adjacency list, and replaces it with newElem
1413   void n2e_replace(int n, int oldElem, int newElem);
1414   /// Remove all elements from n's element adjacency list
1415   void n2e_removeAll(int n);
1416
1417   /// Get an element on edge (n1, n2) where n1, n2 are chunk-local
1418   /// node numberings and result is chunk-local element; return -1 in case 
1419   /// of failure
1420   int getElementOnEdge(int n1, int n2);
1421
1422   /// Get two elements adjacent to both n1 and n2
1423   void get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) ;
1424 }; 
1425 PUPmarshall(FEM_Mesh);
1426
1427 FEM_Mesh *FEM_Mesh_lookup(int fem_mesh,const char *caller);
1428 FEM_Entity *FEM_Entity_lookup(int fem_mesh,int entity,const char *caller);
1429 FEM_Attribute *FEM_Attribute_lookup(int fem_mesh,int entity,int attr,const char *caller);
1430
1431 void FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,     
1432                           void *data, int firstItem,int length, const IDXL_Layout &layout);
1433
1434 //registration internal api
1435 void FEM_Register_array_layout(int fem_mesh,int entity,int attr,void *data,int firstItem,const IDXL_Layout &layout);
1436 void FEM_Register_entity_impl(int fem_mesh,int entity,void *args,int len,int max,FEM_Mesh_alloc_fn fn);
1437 /// Reassemble split chunks into a single mesh
1438 FEM_Mesh *FEM_Mesh_assemble(int nchunks,FEM_Mesh **chunks);
1439
1440 FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead);
1441 FEM_Mesh *FEM_readMesh(const char *prefix,int chunkNo,int nChunks);
1442 void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks);
1443
1444
1445 /*Charm++ Finite Element Framework:
1446   C++ implementation file
1447
1448   This is the main internal implementation file for FEM.
1449   It includes all the other headers, and contains quite
1450   a lot of assorted leftovers and utility routines.
1451
1452   Orion Sky Lawlor, olawlor@acm.org, 9/28/00
1453 */
1454
1455 /* A stupid, stupid number: the maximum value of user-defined "elType" fields. 
1456    This should be dynamic, so any use of this should be considered a bug.
1457 */
1458 #define FEM_MAX_ELTYPE 20
1459
1460 // Verbose abort routine used by FEM framework:
1461 void FEM_Abort(const char *msg);
1462 void FEM_Abort(const char *caller,const char *sprintf_msg,int int0=0,int int1=0,int int2=0);
1463
1464
1465 ///This class describes a local-to-global index mapping, used in FEM_Print.
1466 /**The default is the identity mapping.*/
1467 class l2g_t {
1468  public:
1469   //Return the global number associated with this local element
1470   virtual int el(int t,int localNo) const {return localNo;}
1471   //Return the global number associated with this local node
1472   virtual int no(int localNo) const {return localNo;}
1473 };
1474
1475 /// Map (user-assigned) numbers to T's
1476 template <class T>
1477 class NumberedVec {
1478   CkPupPtrVec<T, CkPupAlwaysAllocatePtr<T> > vec;
1479         
1480  public:
1481   //Extend the vector to have up to this element
1482   void makeLonger(int toHaveElement)
1483     {
1484       int oldSize=vec.size(), newSize=toHaveElement+1;
1485       if (oldSize>=newSize) return; //Nothing to do
1486       vec.resize(newSize);
1487       for (int j=oldSize;j<newSize;j++)
1488         vec[j]=new T;
1489     }
1490   //Reinitialize element i:
1491   void reinit(int doomedEl) {
1492     vec[doomedEl].destroy();
1493     vec[doomedEl]=new T;
1494   }
1495         
1496   int size(void) const {return vec.size();}
1497         
1498   //Same old bracket operators, but return the actual object, not a pointer:
1499   T &operator[](int i) {
1500     if (i>=vec.size()) makeLonger(i);
1501     return *( vec[i] );
1502   }
1503   const T &operator[](int i) const {return *( vec[i] );}
1504         
1505   void pup(PUP::er &p) {
1506     vec.pup(p);
1507   }
1508   friend void operator|(PUP::er &p,NumberedVec<T> &v) {v.pup(p);}
1509 };
1510
1511
1512 ///Smart pointer-to-new[]'d array-of-T
1513 template <class T>
1514 class ArrayPtrT : public CkNoncopyable {
1515   T *sto;
1516  public:
1517   ArrayPtrT() {sto=NULL;}
1518   ArrayPtrT(int *src) {sto=src;}
1519   ~ArrayPtrT() {if (sto) delete[] sto;}
1520   void operator=(T *src) {
1521     if (sto) delete[] sto;
1522     sto=src;
1523   }
1524   operator T *(void) {return sto;}
1525   operator const T *(void) const {return sto;}
1526   T& operator[](int i) {return sto[i];}
1527   const T& operator[](int i) const {return sto[i];}
1528 };
1529 typedef ArrayPtrT<int> intArrayPtr;
1530
1531
1532
1533 /// Unmarshall into a heap-allocated copy
1534 template<class T>
1535 class marshallNewHeapCopy {
1536   T *cur;
1537  public:
1538   //Used on send side:
1539   marshallNewHeapCopy(T *readFrom) :cur(readFrom) {}
1540   marshallNewHeapCopy(const marshallNewHeapCopy &h) :cur(h.cur) {}
1541   marshallNewHeapCopy(void) { //Used on recv side:
1542     cur=new T;
1543   }
1544         
1545   void pup(PUP::er &p) {
1546     cur->pup(p);
1547   }
1548   operator T *() {return cur;}
1549   friend void operator|(PUP::er &p,marshallNewHeapCopy<T> &h) {h.pup(p);}
1550 };
1551 typedef marshallNewHeapCopy<FEM_Mesh> marshallMeshChunk;
1552
1553
1554 /// Keeps a list of dynamically-allocated T objects, indexed by a user-carried, persistent "int".
1555 template<class T>
1556 class FEM_T_List {
1557   CkPupPtrVec<T> list; // Vector of T's
1558  protected:
1559   int FIRST_DT; // User index of first T
1560   int size(void) const {return list.size();}
1561         
1562   /// If this isn't a valid, allocated index, abort.
1563   inline void check(int l,const char *caller) const {
1564     if (l<FIRST_DT || l>=FIRST_DT+list.size() || list[l-FIRST_DT]==NULL) 
1565       badIndex(l,caller);
1566   }
1567         
1568   void badIndex(int l,const char *caller) const {
1569     if (l<FIRST_DT || l>FIRST_DT+list.size()) bad(l,0,caller);
1570     else bad(l,1,caller);
1571   }
1572  public:
1573   FEM_T_List(int FIRST_DT_) :FIRST_DT(FIRST_DT_) {}
1574   virtual ~FEM_T_List() {}
1575   void pup(PUP::er &p) { p|list; }
1576         
1577   /// This routine is called when we're passed an invalid T index.
1578   virtual void bad(int l,int bad_code,const char *caller) const =0;
1579         
1580   /// Insert a new T (allocated with "new"), returning the user index:
1581   int put(T *t) {
1582     for (int i=0;i<list.size();i++) 
1583       if (list[i]==NULL) {
1584         list[i]=t;
1585         return FIRST_DT+i;
1586       }
1587     int ret=list.size();
1588     list.push_back(t);
1589     return FIRST_DT+ret;
1590   }
1591         
1592   /// Get this T given its user index.
1593   inline T *lookup(int l,const char *caller) const {
1594     check(l,caller);
1595     return list[l-FIRST_DT];
1596   }
1597         
1598   /// Free this T
1599   void destroy(int l,const char *caller) {
1600     check(l,caller);
1601     list[l-FIRST_DT].destroy();
1602   }
1603         
1604   /// Clear all stored T's:
1605   void empty(void) {
1606     for (int i=0;i<list.size();i++) list[i].destroy();
1607   }
1608 };
1609 class FEM_Mesh_list : public FEM_T_List<FEM_Mesh> {
1610   typedef FEM_T_List<FEM_Mesh> super;
1611  public:
1612   FEM_Mesh_list() :super(FEM_MESH_FIRST) { }
1613         
1614   virtual void bad(int l,int bad_code,const char *caller) const;
1615 };
1616
1617 #define CHK(p) do{if((p)==0)CkAbort("FEM>Memory Allocation failure.");}while(0)
1618
1619
1620 ///A collection of meshes on one processor
1621 /**
1622    FEM global data object.  Keeps track of the global
1623    list of meshes, and the default read and write meshes.
1624   
1625    This class was once an array element, back when the 
1626    FEM framework was built directly on Charm++.
1627   
1628    There's only one of this object per thread, and it's
1629    kept in a thread-private variable.
1630 */
1631 class FEM_chunk 
1632 {
1633  public:
1634   FEM_Mesh_list meshes; ///< Global list of meshes.
1635   int default_read; ///< Index of default read mesh.
1636   int default_write; ///< Index of default write mesh.
1637   
1638   /// Default communicator to use
1639   FEM_Comm_t defaultComm;
1640
1641   /// Global index (rank) in default communicator
1642   int thisIndex;
1643
1644 #ifdef CMK_OPTIMIZE /* Skip the check, for speed. */
1645   inline void check(const char *where) { }
1646 #else /* Do an extensive self-check */
1647   void check(const char *where);
1648 #endif
1649
1650  private:
1651   CkVec<int> listTmp;//List of local entities, for ghost list exchange
1652  
1653   void initFields(void);
1654
1655  public:
1656   FEM_chunk(FEM_Comm_t defaultComm_);
1657   FEM_chunk(CkMigrateMessage *msg);
1658   void pup(PUP::er &p);
1659   ~FEM_chunk();
1660   
1661   /// Return this thread's single static FEM_chunk instance:
1662   static FEM_chunk *get(const char *caller);
1663   
1664   inline FEM_Mesh *lookup(int fem_mesh,const char *caller) {
1665     return meshes.lookup(fem_mesh,caller);
1666   }
1667
1668   inline FEM_Mesh *getMesh(const char *caller) 
1669     {return meshes.lookup(default_read,caller);}
1670   inline FEM_Mesh *setMesh(const char *caller) 
1671     {return meshes.lookup(default_write,caller);}
1672
1673   void print(int fem_mesh,int idxBase);
1674   int getPrimary(int nodeNo) { return getMesh("getPrimary")->node.getPrimary(nodeNo); }
1675   const FEM_Comm &getComm(void) {return getMesh("getComm")->node.shared;}
1676
1677   // Basically everything below here should be moved to IDXL:
1678   void exchangeGhostLists(int elemType,int inLen,const int *inList,int idxbase);
1679   void recvList(int elemType,int fmChk,int nIdx,const int *idx);
1680   const CkVec<int> &getList(void) {return listTmp;}
1681   void emptyList(void) {listTmp.length()=0;}
1682   
1683   void reduce_field(int idxl_datatype, const void *nodes, void *outbuf, int op);
1684   void reduce(int idxl_datatype, const void *inbuf, void *outbuf, int op);
1685   void readField(int idxl_datatype, void *nodes, const char *fname);  
1686 };
1687
1688 /// Describes a single layer of ghost elements.
1689 class FEM_Ghost_Layer : public CkNoncopyable {
1690  public:
1691   int nodesPerTuple; //Number of shared nodes needed to connect elements
1692   bool addNodes; //Add ghost nodes to the chunks
1693   class elemGhostInfo {
1694   public:
1695     bool add; //Add this kind of ghost element to the chunks
1696     int tuplesPerElem; //# of tuples surrounding this element
1697     intArrayPtr elem2tuple; //The tuples around this element [nodesPerTuple * tuplesPerElem]
1698     elemGhostInfo(void) {add=false;tuplesPerElem=0;}
1699     ~elemGhostInfo(void) {}
1700     void pup(PUP::er &p) {//CkAbort("FEM> Shouldn't call elemGhostInfo::pup!\n");
1701     }
1702   };
1703   elemGhostInfo elem[FEM_MAX_ELTYPE];
1704   virtual void pup(PUP::er &p){
1705     p | nodesPerTuple;
1706     p | addNodes;
1707     for(int i=0;i<FEM_MAX_ELTYPE;i++){
1708       p | elem[i].add;
1709       p | elem[i].tuplesPerElem;
1710       if(elem[i].tuplesPerElem == 0){
1711         continue;
1712       }
1713       int *arr;
1714       if(p.isUnpacking()){
1715         arr = new int[nodesPerTuple*elem[i].tuplesPerElem];
1716       }else{
1717         arr = elem[i].elem2tuple;
1718       }
1719       p(arr,nodesPerTuple*elem[i].tuplesPerElem);
1720       if(p.isUnpacking()){
1721         elem[i].elem2tuple = arr;
1722       }
1723     }
1724   }
1725 };
1726
1727 /// Describes a set of required adjacent elements for this kind of element,
1728 ///  stored as an explicit adjacency list.
1729 class FEM_Ghost_Stencil {
1730   /// Our element type.
1731   int elType;
1732   /// Number of elements we describe.
1733   int n;
1734   /// If true, add ghost nodes as well as elements
1735   bool addNodes;
1736         
1737   /// Last adjacency entry (plus one), indexed by element.
1738   ///  That is, element i's data is at [ends[i-1],ends[i])
1739   intArrayPtr ends;
1740         
1741   /** Adjacency entries for each element.
1742       Stored as a series of pairs: elType, elNum.
1743       The first pair for element i starts at
1744       2*(ends[i-1]) 
1745       the last pair for element i starts at
1746       2*(ends[i]-1)
1747       This array then has, in total, 2*ends[n-1] elements.
1748   */
1749   intArrayPtr adj;
1750  public:
1751   /**
1752      Create a stencil with this number of elements, 
1753      and these adjecent elements.
1754   */
1755   FEM_Ghost_Stencil(int elType_, int n_,bool addNodes_,
1756                     const int *ends_,
1757                     const int *adj_,
1758                     int idxBase);
1759         
1760   /// Make sure this stencil makes sense for this mesh.
1761   void check(const FEM_Mesh &mesh) const;
1762         
1763   /// Return the type of element we describe
1764   inline int getType(void) const {return elType;}
1765         
1766   /**
1767      Return a pair consisting of the i'th element's
1768      j'th neighbor: the return value's first int is an element type,
1769      the second int is an element number.  Returns NULL if i doesn't
1770      have a j'th neighbor.
1771   */
1772   inline const int *getNeighbor(int i,int j) const {
1773     int start=0;
1774     if (i>0) start=ends[i-1];
1775     if (j>=ends[i]-start) return 0;
1776     return &adj[2*(start+j)];
1777   }
1778         
1779   inline bool wantNodes(void) const {return addNodes;}
1780 };
1781
1782 /// Describes a way to grow a set of ghosts.
1783 class FEM_Ghost_Region {
1784  public:        
1785   FEM_Ghost_Layer *layer;
1786   FEM_Ghost_Stencil *stencil;
1787         
1788   FEM_Ghost_Region() {layer=0; stencil=0;}
1789   FEM_Ghost_Region(FEM_Ghost_Layer *l) {layer=l; stencil=0;}
1790   FEM_Ghost_Region(FEM_Ghost_Stencil *s) {layer=0; stencil=s;}
1791 };
1792
1793
1794 //Accumulates all symmetries of the mesh before splitting:
1795 class FEM_Initial_Symmetries; /*Defined in symmetries.C*/
1796
1797 /// Describes all the data needed for partitioning a mesh.
1798 class FEM_Partition : public CkNoncopyable {
1799   /// Maps element number to (0-based) chunk number, allocated with new[]
1800   int *elem2chunk;
1801         
1802   /// Describes the different regions of ghost elements:
1803   CkVec<FEM_Ghost_Region> regions;
1804   FEM_Ghost_Layer *lastLayer;
1805         
1806   /// Describes the problem domain's spatial symmetries.
1807   FEM_Initial_Symmetries *sym;
1808  public:
1809   FEM_Partition();
1810   ~FEM_Partition();
1811         
1812   // Manipulate partitioning information
1813   void setPartition(const int *elem2chunk, int nElem, int idxBase);
1814   const int *getPartition(FEM_Mesh *src,int nChunks) const;
1815         
1816   // Manipulate ghost layers
1817   FEM_Ghost_Layer *addLayer(void) {
1818     lastLayer=new FEM_Ghost_Layer();
1819     regions.push_back(lastLayer);
1820     return lastLayer;
1821   }
1822   FEM_Ghost_Layer *curLayer(void) {
1823     if (lastLayer==0) CkAbort("Must call FEM_Add_ghost_layer before FEM_Add_ghost_elem\n");
1824     return lastLayer;
1825   }
1826         
1827   // Manipulate ghost stencils
1828   void addGhostStencil(FEM_Ghost_Stencil *s) {
1829     regions.push_back(s);
1830     lastLayer=0;
1831   }
1832   void markGhostStencilLayer(void) {
1833     regions.push_back(FEM_Ghost_Region());
1834   }
1835         
1836   // Read back ghost regions
1837   int getRegions(void) const {return regions.size();}
1838   const FEM_Ghost_Region &getRegion(int regNo) const {return regions[regNo];}
1839         
1840   // Manipulate spatial symmetries:
1841   void setSymmetries(int nNodes_,int *new_can,const int *sym_src);
1842   void addLinearPeriodic(int nFaces_,int nPer,
1843                          const int *facesA,const int *facesB,int idxBase,
1844                          int nNodes_,const CkVector3d *nodeLocs);
1845   const int *getCanon(void) const;
1846   const FEM_Symmetries_t *getSymmetries(void) const;
1847   const FEM_Sym_List &getSymList(void) const;
1848
1849
1850 };
1851 // Access the latest partition:
1852 FEM_Partition &FEM_curPartition(void);
1853
1854 //Declare this at the start of every API routine:
1855 #define FEMAPI(routineName) TCHARM_API_TRACE(routineName,"fem")
1856 //#define FEMAPI(routineName) printf("%s\n", routineName);
1857
1858
1859 /*Partition this mesh's elements into n chunks,
1860   writing each element's 0-based chunk number to elem2chunk.
1861 */
1862 void FEM_Mesh_partition(const FEM_Mesh *mesh,int nchunks,int *elem2chunk);
1863
1864 /*A way to stream out partitioned chunks of a mesh.
1865   By streaming, we can send the chunks as they are built,
1866   dramatically reducing the memory needed by the framework.
1867 */
1868 class FEM_Mesh_Output {
1869  public:
1870   virtual ~FEM_Mesh_Output() {} /*<- for whining compilers*/
1871   //Transfer ownership of this mesh chunk
1872   virtual void accept(int chunkNo,FEM_Mesh *msg) =0;
1873 };
1874
1875 /*After partitioning, create a sub-mesh for each chunk's elements,
1876   including communication lists between chunks.
1877 */
1878 void FEM_Mesh_split(FEM_Mesh *mesh,int nchunks,
1879                     const FEM_Partition &partition,FEM_Mesh_Output *out);
1880
1881
1882 //Make a new[]'d copy of this (len-entry) array, changing the index as spec'd
1883 int *CkCopyArray(const int *src,int len,int indexBase);
1884
1885
1886 // Isaac's stuff:
1887 // Describes Element Faces. For use with finding element->element adjacencies
1888 // based on FEM_Ghost_Layer
1889 class FEM_ElemAdj_Layer : public CkNoncopyable {
1890  public:
1891   int initialized;
1892   int nodesPerTuple; //Number of shared nodes for a pair of elements
1893   
1894   class elemAdjInfo {
1895   public:
1896     //  int recentElType; // should not be here, but if it is it should be pup'ed
1897     int tuplesPerElem; //# of tuples surrounding this element, i.e. number of faces on an element
1898     intArrayPtr elem2tuple; //The tuples around this element [nodesPerTuple * tuplesPerElem]
1899     elemAdjInfo(void) {/*add=false;*/tuplesPerElem=0;}
1900     ~elemAdjInfo(void) {}
1901     void pup(PUP::er &p) {//CkAbort("FEM> Shouldn't call elemGhostInfo::pup!\n");
1902     }
1903   };
1904   
1905   elemAdjInfo elem[FEM_MAX_ELTYPE];
1906   
1907   FEM_ElemAdj_Layer() {initialized=0;}
1908   
1909   virtual void pup(PUP::er &p){
1910     p | nodesPerTuple;
1911     p | initialized;
1912     for(int i=0;i<FEM_MAX_ELTYPE;i++){
1913       p | elem[i].tuplesPerElem;
1914       if(elem[i].tuplesPerElem == 0){
1915         continue;
1916       }
1917       int *arr;
1918       if(p.isUnpacking()){
1919         arr = new int[nodesPerTuple*elem[i].tuplesPerElem];
1920       }else{
1921         arr = elem[i].elem2tuple;
1922       }
1923       p(arr,nodesPerTuple*elem[i].tuplesPerElem);
1924       if(p.isUnpacking()){
1925         elem[i].elem2tuple = arr;
1926       }
1927     }
1928   }
1929 };
1930
1931
1932 /*\@}*/
1933
1934 // End impl.h
1935
1936
1937
1938 /*
1939   File containing the data structures, function declaration 
1940   and msa array declarations used during parallel partitioning
1941   of the mesh.
1942   Author Sayantan Chakravorty
1943   05/30/2004
1944 */
1945
1946
1947 //#define PARALLEL_DEBUG
1948
1949 #ifdef PARALLEL_DEBUG
1950 #define DEBUG(x) x
1951 #else
1952 #define DEBUG(x)
1953 #endif
1954
1955 #define MESH_CHUNK_TAG 3000
1956 #define MAX_SHARED_PER_NODE 20
1957
1958 template <class T, bool PUP_EVERY_ELEMENT=true >
1959   class DefaultListEntry {
1960     public:
1961     inline void accumulate(T& a, const T& b) { a += b; }
1962     // identity for initializing at start of accumulate
1963     inline T getIdentity() { return T(); }
1964     inline bool pupEveryElement(){ return PUP_EVERY_ELEMENT; }
1965   };
1966
1967 extern double elemlistaccTime;
1968
1969 template <class T>
1970 class ElemList{
1971  public:
1972   CkVec<T> *vec;
1973   ElemList(){
1974     vec = new CkVec<T>();
1975   }
1976   ~ElemList(){
1977     delete vec;
1978   }
1979   ElemList(const ElemList &rhs){
1980     vec = new CkVec<T>();
1981     *this=rhs;
1982   }
1983   inline ElemList& operator=(const ElemList& rhs){
1984     //           vec = new CkVec<T>();
1985     *vec = *(rhs.vec);
1986                 return *this;
1987   }
1988   inline ElemList& operator+=(const ElemList& rhs){
1989     /*
1990       add the new unique elements to the List
1991     */
1992     double _start = CkWallTimer();
1993     int len = vec->size();
1994     for(int i=0;i<rhs.vec->length();i++){
1995       vec->push_back((*(rhs.vec))[i]);
1996     }
1997     //          uniquify();
1998     elemlistaccTime += (CkWallTimer() - _start);
1999     return *this;
2000   }
2001   ElemList(const T &val){
2002     vec =new CkVec<T>();
2003     vec->push_back(val);
2004   };
2005   inline virtual void pup(PUP::er &p){
2006     if(p.isUnpacking()){
2007       vec = new CkVec<T>();
2008     }
2009     pupCkVec(p,*vec);
2010   }
2011 };
2012 template <class T>
2013 class UniqElemList: public ElemList<T>{
2014 public:
2015   UniqElemList(const T &val):ElemList<T>(val){};
2016   UniqElemList():ElemList<T>(){};
2017   inline void uniquify(){
2018     CkVec<T> *lvec = this->vec;
2019     lvec->quickSort(8);
2020     if(lvec->length() != 0){
2021       int count=0;
2022       for(int i=1;i<lvec->length();i++){
2023         if((*lvec)[count] == (*lvec)[i]){       
2024         }else{                                  
2025           count++;
2026           if(i != count){
2027             (*lvec)[count] = (*lvec)[i];
2028           }
2029         }
2030       }
2031       lvec->resize(count+1);
2032     }   
2033   }
2034 };
2035
2036 class NodeElem {
2037  public:
2038   //global number of this node
2039   int global;
2040   /*no of chunks that share this node
2041     owned by 1 element - numShared 0
2042     owned by 2 elements - numShared 2
2043   */
2044   int numShared;
2045   /*somewhat horrible semantics
2046     -if numShared == 0 shared is NULL
2047     - else stores all the chunks that share this node
2048   */
2049   int shared[MAX_SHARED_PER_NODE];
2050   //    int *shared;
2051   NodeElem(){
2052     global = -1;
2053     numShared = 0;
2054     //          shared = NULL;
2055   }
2056   NodeElem(int g_,int num_){
2057     global = g_; numShared= num_;
2058     /*          if(numShared != 0){
2059                 shared = new int[numShared];
2060                 }else{
2061                 shared = NULL;
2062                 }*/
2063     if(numShared > MAX_SHARED_PER_NODE){
2064       CkAbort("In Parallel partition the max number of chunks per node exceeds limit\n");
2065     }
2066   }
2067   NodeElem(int g_){
2068     global = g_;
2069     numShared = 0;
2070     //          shared = NULL;
2071   }
2072   NodeElem(const NodeElem &rhs){
2073     //          shared=NULL;
2074     *this = rhs;
2075   }
2076   inline NodeElem& operator=(const NodeElem &rhs){
2077     global = rhs.global;
2078     numShared = rhs.numShared;
2079     /*          if(shared != NULL){
2080                 delete [] shared;
2081                 }
2082                 shared = new int[numShared];*/
2083     memcpy(&shared[0],&((rhs.shared)[0]),numShared*sizeof(int));
2084     return *this;
2085   }
2086
2087   inline        bool operator == (const NodeElem &rhs){
2088     if(global == rhs.global){
2089       return true;
2090     }else{
2091       return false;
2092     }
2093   }
2094   inline        bool operator >=(const NodeElem &rhs){
2095     return (global >= rhs.global);
2096   }
2097
2098   inline bool operator <=(const NodeElem &rhs){
2099     return (global <= rhs.global);
2100   }
2101
2102   virtual void pup(PUP::er &p){
2103     p | global;
2104     p | numShared;
2105     /*          if(p.isUnpacking()){
2106                 if(numShared != 0){
2107                 shared = new int[numShared];
2108                 }else{
2109                 shared = NULL;
2110                 }
2111                 }*/
2112     p(shared,numShared);
2113   }
2114   ~NodeElem(){
2115     /*          if(shared != NULL){
2116                 delete [] shared;
2117                 }*/
2118   }
2119 };
2120
2121 /**
2122   This class is an MSA Entity. It is used for 2 purposes
2123   1 It is used for storing the mesh while creating the mesh
2124   for each chunk
2125   2     It is used for storing the ghost elements and nodes for
2126   a chunk.
2127 */
2128 class MeshElem{
2129  public: 
2130   FEM_Mesh *m;
2131   CkVec<int> gedgechunk; // Chunk number of 
2132   MeshElem(){
2133     m = new FEM_Mesh;
2134   }
2135   MeshElem(int dummy){
2136     m = new FEM_Mesh;
2137   }
2138   ~MeshElem(){
2139     delete m;
2140   }
2141   MeshElem(const MeshElem &rhs){
2142     m = NULL;
2143     *this = rhs;
2144   }
2145   inline MeshElem& operator=(const MeshElem &rhs){
2146     if(m != NULL){
2147       delete m;
2148     }   
2149     m = new FEM_Mesh;
2150     m->copyShape(*(rhs.m));
2151     (*this) += rhs;
2152                 
2153     return *this;
2154   }
2155   inline MeshElem& operator+=(const MeshElem &rhs){
2156     int oldel = m->nElems();
2157     m->copyShape(*(rhs.m));
2158     for(int i=0;i<rhs.m->node.size();i++){
2159       m->node.push_back((rhs.m)->node,i);
2160     }   
2161     if((rhs.m)->elem.size()>0){
2162       for(int t=0;t<(rhs.m)->elem.size();t++){
2163         if((rhs.m)->elem.has(t)){
2164           for(int e=0;e<(rhs.m)->elem.get(t).size();e++){
2165             m->elem[t].push_back((rhs.m)->elem.get(t),e);
2166           }     
2167         }       
2168       }
2169     }
2170
2171     return *this;
2172   }
2173   virtual void pup(PUP::er &p){
2174     if(p.isUnpacking()){
2175       m = new FEM_Mesh;
2176     }
2177     m->pup(p);
2178   }
2179 };
2180
2181
2182 class Hashnode{
2183 public:
2184         class tupledata{
2185                 public:
2186                 enum {MAX_TUPLE = 8};
2187                         int nodes[MAX_TUPLE];
2188                         tupledata(int _nodes[MAX_TUPLE]){
2189                                 memcpy(nodes,_nodes,sizeof(int)*MAX_TUPLE);
2190                         }
2191                         tupledata(tupledata &rhs){
2192                                 memcpy(nodes,rhs.nodes,sizeof(int)*MAX_TUPLE);
2193                         }
2194                         tupledata(const tupledata &rhs){
2195                                 memcpy(nodes,rhs.nodes,sizeof(int)*MAX_TUPLE);
2196                         }
2197                         tupledata(){};
2198                         //dont store the returned string
2199                         char *toString(int numnodes,char *str){
2200                                 str[0]='\0';
2201                                 for(int i=0;i<numnodes;i++){
2202                                         sprintf(&str[strlen(str)],"%d ",nodes[i]);
2203                                 }
2204                                 return str;
2205                         }
2206                         inline int &operator[](int i){
2207                                 return nodes[i];
2208                         }
2209                         inline const int &operator[](int i) const {
2210                                 return nodes[i];
2211                         }
2212                         virtual void pup(PUP::er &p){
2213                                 p(nodes,MAX_TUPLE);
2214                         }
2215         };
2216         int numnodes; //number of nodes in this tuple
2217         //TODO: replace *nodes with the above tupledata class
2218         tupledata nodes;        //the nodes in the tuple
2219         int chunk;              //the chunk number to which this element belongs 
2220         int elementNo;          //local number of that element
2221         Hashnode(){
2222                 numnodes=0;
2223         };
2224         Hashnode(int _num,int _chunk,int _elNo,int _nodes[tupledata::MAX_TUPLE]): nodes(_nodes){
2225                 numnodes = _num;
2226                 chunk = _chunk;
2227                 elementNo = _elNo;
2228         }
2229         Hashnode(const Hashnode &rhs){
2230                 *this = rhs;
2231         }
2232         inline Hashnode &operator=(const Hashnode &rhs){
2233                 numnodes = rhs.numnodes;
2234                 for(int i=0;i<numnodes;i++){
2235                         nodes[i] = rhs.nodes[i];
2236                 }
2237                 chunk = rhs.chunk;
2238                 elementNo = rhs.elementNo;
2239                 return *this;
2240         }
2241         inline bool operator==(const Hashnode &rhs){
2242                 if(numnodes != rhs.numnodes){
2243                         return false;
2244                 }
2245                 for(int i=0;i<numnodes;i++){
2246                         if(nodes[i] != rhs.nodes[i]){
2247                                 return false;
2248                         }
2249                 }
2250                 if(chunk != rhs.chunk){
2251                         return false;
2252                 }
2253                 if(elementNo != rhs.elementNo){
2254                         return false;
2255                 }
2256                 return true;
2257         }
2258         inline bool operator>=(const Hashnode &rhs){
2259                 if(numnodes < rhs.numnodes){
2260                         return false;
2261                 };
2262                 if(numnodes > rhs.numnodes){
2263                         return true;
2264                 }
2265
2266     for(int i=0;i<numnodes;i++){
2267       if(nodes[i] < rhs.nodes[i]){
2268         return false;
2269       }
2270       if(nodes[i] > rhs.nodes[i]){
2271         return true;
2272       }
2273     }
2274     if(chunk < rhs.chunk){
2275       return false;
2276     }
2277     if(chunk > rhs.chunk){
2278       return true;
2279     }
2280     if(elementNo < rhs.elementNo){
2281       return false;
2282     }
2283     if(elementNo > rhs.elementNo){
2284       return true;
2285     }
2286     return true;
2287   }
2288         
2289   inline bool operator<=(const Hashnode &rhs){
2290     if(numnodes < rhs.numnodes){
2291       return true;
2292     };
2293     if(numnodes > rhs.numnodes){
2294       return false;
2295     }
2296                 
2297     for(int i=0;i<numnodes;i++){
2298       if(nodes[i] < rhs.nodes[i]){
2299         return true;
2300       }
2301       if(nodes[i] > rhs.nodes[i]){
2302         return false;
2303       }
2304     }
2305     if(chunk < rhs.chunk){
2306       return true;
2307     }
2308     if(chunk > rhs.chunk){
2309       return false;
2310     }
2311     if(elementNo < rhs.elementNo){
2312       return true;
2313     }
2314     if(elementNo > rhs.elementNo){
2315       return false;
2316     }
2317     return true;
2318   }
2319         
2320   inline bool equals(tupledata &tuple){
2321     for(int i=0;i<numnodes;i++){
2322       if(tuple.nodes[i] != nodes[i]){
2323         return false;
2324       }
2325     }
2326     return true;
2327   }
2328   virtual void pup(PUP::er &p){
2329     p | numnodes;
2330     p | nodes;
2331     p | chunk;
2332     p | elementNo;
2333   }
2334 };
2335
2336 template <class T>
2337 ostream& operator << (ostream& os, const ElemList<T> & s){
2338 };
2339
2340 template <class T>
2341 CkOutStream& operator << (CkOutStream& os, const ElemList<T>& s) {
2342 };
2343
2344 typedef MSA2D<int, DefaultEntry<int>, MSA_DEFAULT_ENTRIES_PER_PAGE, MSA_ROW_MAJOR> MSA2DRM;
2345
2346 typedef MSA1D<int, DefaultEntry<int>, MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DINT;
2347
2348 typedef UniqElemList<int> IntList;
2349 typedef MSA1D<IntList, DefaultListEntry<IntList,true>,MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DINTLIST;
2350
2351 typedef UniqElemList<NodeElem> NodeList;
2352 typedef MSA1D<NodeList, DefaultListEntry<NodeList,true>,MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DNODELIST;
2353
2354 typedef MSA1D<MeshElem,DefaultEntry<MeshElem,true>,1> MSA1DFEMMESH;
2355
2356 typedef UniqElemList<Hashnode> Hashtuple;
2357 typedef MSA1D<Hashtuple,DefaultListEntry<Hashtuple,true>,MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DHASH;
2358
2359
2360
2361 struct conndata{
2362   int nelem;
2363   int nnode;
2364   MSA1DINT arr1;
2365   MSA1DINT arr2;
2366
2367   void pup(PUP::er &p){
2368     p|nelem;
2369     p|nnode;
2370     arr1.pup(p);
2371     arr2.pup(p);
2372   }
2373 };
2374
2375 /**
2376   Structure to store connectivity data after the 
2377   global element partition has been returned by parmetis
2378 */
2379 struct partconndata{
2380   int nelem;
2381   int startindex;
2382   int *eptr,*eind;
2383   int *part;
2384   ~partconndata(){
2385     delete [] eptr;
2386     delete [] eind;
2387     delete [] part;
2388   };
2389 };
2390
2391 /**
2392   structure for storing the ghost layers
2393 */
2394 struct ghostdata{
2395   int numLayers;
2396   FEM_Ghost_Layer **layers;
2397   void pup(PUP::er &p){
2398     p | numLayers;
2399     if(p.isUnpacking()){
2400       layers = new FEM_Ghost_Layer *[numLayers];
2401       for(int i=0;i<numLayers;i++){
2402         layers[i] = new FEM_Ghost_Layer;
2403       }
2404
2405     }
2406     for(int i=0;i<numLayers;i++){
2407       layers[i]->pup(p);
2408     }
2409   }
2410   ~ghostdata(){
2411     //printf("destructor on ghostdata called \n");
2412     for(int i=0;i<numLayers;i++){
2413       delete layers[i];
2414     }
2415     delete [] layers;
2416   };
2417 };
2418
2419
2420 class MsaHashtable{
2421  public:
2422   int numSlots;
2423   MSA1DHASH table;
2424   MsaHashtable(int _numSlots,int numWorkers):numSlots(_numSlots),table(_numSlots,numWorkers){
2425   }
2426   MsaHashtable(){};
2427
2428   virtual void pup(PUP::er &p){
2429     p | numSlots;
2430     p | table;
2431   }
2432   int addTuple(int *tuple,int nodesPerTuple,int chunk,int elementNo){
2433     // sort the tuples to get a canonical form
2434     // bubble sort should do just as well since the number
2435     // of nodes is less than 10.
2436     for(int i=0;i<nodesPerTuple-1;i++){
2437       for(int j=i+1;j<nodesPerTuple;j++){
2438         if(tuple[j] < tuple[i]){
2439           int t = tuple[j];
2440           tuple[j] = tuple[i];
2441           tuple[i] = t;
2442         }
2443       }
2444     }
2445     //find out the index
2446     long long sum = 0;
2447     for(int i=0;i<nodesPerTuple;i++){
2448       sum = sum *numSlots + tuple[i];
2449     }
2450     int index = (int )(sum %(long )numSlots);
2451     Hashnode entry(nodesPerTuple,chunk,elementNo,tuple);
2452         
2453     Hashtuple &list=table.accumulate(index);
2454     list.vec->push_back(entry);
2455     char str[100];
2456     DEBUG(printf("[%d] adding tuple %s element %d to index %d \n",chunk,entry.nodes.toString(nodesPerTuple,str),elementNo,index));
2457     return index;
2458   }
2459
2460   void print(){
2461     char str[100];
2462     for(int i=0;i<numSlots;i++){
2463       const Hashtuple &t = table.get(i);
2464       for(int j=0;j<t.vec->size();j++){
2465         Hashnode &tuple = (*t.vec)[j];
2466         printf("ghost element chunk %d element %d index %d tuple < %s>\n",tuple.chunk,tuple.elementNo,i,tuple.nodes.toString(tuple.numnodes,str));
2467       }
2468     }
2469   }
2470   void sync(){
2471     table.sync();
2472   }
2473   const Hashtuple &get(int i){
2474     return table.get(i);
2475   }
2476         
2477 };
2478
2479
2480 int FEM_master_parallel_part(int ,int ,FEM_Comm_t);
2481 int FEM_slave_parallel_part(int ,int ,FEM_Comm_t);
2482 struct partconndata* FEM_call_parmetis(struct conndata &data,FEM_Comm_t comm_context);
2483 void FEM_write_nodepart(MSA1DINTLIST    &nodepart,struct partconndata *data,MPI_Comm comm_context);
2484 void FEM_write_part2node(MSA1DINTLIST   &nodepart,MSA1DNODELIST &part2node,struct partconndata *data,MPI_Comm comm_context);
2485 void FEM_write_part2elem(MSA1DINTLIST &part2elem,struct partconndata *data,MPI_Comm comm_context);
2486 FEM_Mesh * FEM_break_mesh(FEM_Mesh *m,int numElements,int numChunks);
2487 void sendBrokenMeshes(FEM_Mesh *mesh_array,FEM_Comm_t comm_context);
2488 void    FEM_write_part2mesh(MSA1DFEMMESH &part2mesh,struct partconndata *partdata,struct conndata *data,MSA1DINTLIST &nodepart,int numChunks,int myChunk,FEM_Mesh *mypiece);
2489 void addIDXLists(FEM_Mesh *m,NodeList &lnodes,int myChunk);
2490 struct ghostdata *gatherGhosts();
2491 void makeGhosts(FEM_Mesh *m,MPI_Comm comm,int masterRank,int numLayers,FEM_Ghost_Layer **layers);
2492 void makeGhost(FEM_Mesh *m,MPI_Comm comm,int masterRank,int totalShared,FEM_Ghost_Layer *layer, CkHashtableT<CkHashtableAdaptorT<int>,char> &sharedNode,CkHashtableT<CkHashtableAdaptorT<int>,int> &global2local);
2493 bool sharedWith(int lnode,int chunk,FEM_Mesh *m);
2494
2495 // End Parallel Partitioner
2496
2497
2498
2499 /*
2500   Some definitions of structures/classes used in map.C and possibly other FEM source files.
2501   Moved to this file by Isaac Dooley 4/5/05
2502
2503   FIXME: Clean up the organization of this file, and possibly move other structures from Map.C here
2504   Also all the stuff in this file might just belong in impl.h. I just moved it here so it could
2505   be included in fem.C for my element->element build adjacency table function.
2506 */
2507
2508 static FEM_Symmetries_t noSymmetries=(FEM_Symmetries_t)0;
2509
2510 // defined in map.C
2511 CkHashCode CkHashFunction_ints(const void *keyData,size_t keyLen);
2512 int CkHashCompare_ints(const void *k1,const void *k2,size_t keyLen);
2513 extern "C" int ck_fem_map_compare_int(const void *a, const void *b);
2514
2515 //A linked list of elements surrounding a tuple.
2516 //ElemLists are allocated one at a time, and hence a bit cleaner than chunklists:
2517 class elemList {
2518  public:
2519   int chunk;
2520   int tupleNo;//tuple number on this element for the tuple that this list sorrounds 
2521   int localNo;//Local number of this element on this chunk (negative for a ghost)
2522   int type; //Kind of element
2523   FEM_Symmetries_t sym; //Symmetries this element was reached via
2524   elemList *next;
2525
2526   elemList(int chunk_,int localNo_,int type_,FEM_Symmetries_t sym_)
2527     :chunk(chunk_),localNo(localNo_),type(type_), sym(sym_) 
2528     { next=NULL; }
2529   elemList(int chunk_,int localNo_,int type_,FEM_Symmetries_t sym_,int tupleNo_)
2530     :chunk(chunk_),localNo(localNo_),type(type_), sym(sym_) , tupleNo(tupleNo_)
2531     { next=NULL; }
2532         
2533   ~elemList() {if (next) delete next;}
2534   void setNext(elemList *n) {next=n;}
2535 };
2536
2537
2538 //Maps node tuples to element lists
2539 class tupleTable : public CkHashtable {
2540   int tupleLen; //Nodes in a tuple
2541   CkHashtableIterator *it;
2542   static int roundUp(int val,int to) {
2543     return ((val+to-1)/to)*to;
2544   }
2545   static CkHashtableLayout makeLayout(int tupleLen) {
2546     int ks=tupleLen*sizeof(int);
2547     int oo=roundUp(ks+sizeof(char),sizeof(void *));
2548     int os=sizeof(elemList *);
2549     return CkHashtableLayout(ks,ks,oo,os,oo+os);
2550   }
2551   
2552   //Make a canonical version of this tuple, so different
2553   // orderings of the same nodes don't end up in different lists.
2554   //I canonicalize by sorting:
2555   void canonicalize(const int *tuple,int *can)
2556     {
2557       switch(tupleLen) {
2558       case 1: //Short lists are easy to sort:
2559         can[0]=tuple[0]; break;
2560       case 2:
2561         if (tuple[0]<tuple[1])
2562           {can[0]=tuple[0]; can[1]=tuple[1];}
2563         else
2564           {can[0]=tuple[1]; can[1]=tuple[0];}
2565         break;
2566       default: //Should use std::sort here:
2567         memcpy(can,tuple,tupleLen*sizeof(int));
2568         qsort(can,tupleLen,sizeof(int),ck_fem_map_compare_int);
2569       };
2570     }
2571  public:
2572   enum {MAX_TUPLE=8};
2573
2574   tupleTable(int tupleLen_)
2575     :CkHashtable(makeLayout(tupleLen_),
2576                  137,0.5,
2577                  CkHashFunction_ints,
2578                  CkHashCompare_ints)
2579     {
2580       tupleLen=tupleLen_;
2581       if (tupleLen>MAX_TUPLE) CkAbort("Cannot have that many shared nodes!\n");
2582       it=NULL;
2583     }
2584   ~tupleTable() {
2585     beginLookup();
2586     elemList *doomed;
2587     while (NULL!=(doomed=(elemList *)lookupNext()))
2588       delete doomed;
2589   }
2590   //Lookup the elemList associated with this tuple, or return NULL
2591   elemList **lookupTuple(const int *tuple) {
2592     int can[MAX_TUPLE];
2593     canonicalize(tuple,can);
2594     return (elemList **)get(can);
2595   }     
2596         
2597   //Register this (new'd) element with this tuple
2598   void addTuple(const int *tuple,elemList *nu)
2599     {
2600       int can[MAX_TUPLE];
2601       canonicalize(tuple,can);
2602       //First try for an existing list:
2603       elemList **dest=(elemList **)get(can);
2604       if (dest!=NULL) 
2605         { //A list already exists here-- link it into the new list
2606           nu->setNext(*dest);
2607         } else {//No pre-existing list-- initialize a new one.
2608           dest=(elemList **)put(can);
2609         }
2610       *dest=nu;
2611     }
2612   //Return all registered elemLists:
2613   void beginLookup(void) {
2614     it=iterator();
2615   }
2616   elemList *lookupNext(void) {
2617     void *ret=it->next();
2618     if (ret==NULL) {
2619       delete it; 
2620       return NULL;
2621     }
2622     return *(elemList **)ret;
2623   }
2624 };
2625
2626
2627
2628
2629 /*This object maps a single entity to a list of the chunks
2630   that have a copy of it.  For the vast majority
2631   of entities, the list will contain exactly one element.
2632 */
2633 class chunkList : public CkNoncopyable {
2634  public:
2635   int chunk;//Chunk number; or -1 if the list is empty
2636   int localNo;//Local number of this entity on this chunk (if negative, is a ghost)
2637   FEM_Symmetries_t sym; //Symmetries this entity was reached via
2638   int layerNo; // -1 if real; if a ghost, our ghost layer number
2639   chunkList *next;
2640   chunkList() {chunk=-1;next=NULL;}
2641   chunkList(int chunk_,int localNo_,FEM_Symmetries_t sym_,int ln_=-1) {
2642     chunk=chunk_;
2643     localNo=localNo_;
2644     sym=sym_;
2645     layerNo=ln_;
2646     next=NULL;
2647   }
2648   ~chunkList() {delete next;}
2649   void set(int c,int l,FEM_Symmetries_t s,int ln) {
2650     chunk=c; localNo=l; sym=s; layerNo=ln;
2651   }
2652   //Is this chunk in the list?  If so, return false.
2653   // If not, add it and return true.
2654   bool addchunk(int c,int l,FEM_Symmetries_t s,int ln) {
2655     //Add chunk c to the list with local index l,
2656     // if it's not there already
2657     if (chunk==c && sym==s) return false;//Already in the list
2658     if (chunk==-1) {set(c,l,s,ln);return true;}
2659     if (next==NULL) {next=new chunkList(c,l,s,ln);return true;}
2660     else return next->addchunk(c,l,s,ln);
2661   }
2662   //Return this node's local number on chunk c (or -1 if none)
2663   int localOnChunk(int c,FEM_Symmetries_t s) const {
2664     const chunkList *l=onChunk(c,s);
2665     if (l==NULL) return -1;
2666     else return l->localNo;
2667   }
2668   const chunkList *onChunk(int c,FEM_Symmetries_t s) const {
2669     if (chunk==c && sym==s) return this;
2670     else if (next==NULL) return NULL;
2671     else return next->onChunk(c,s);
2672   }
2673   int isEmpty(void) const //Return 1 if this is an empty list 
2674     {return (chunk==-1);}
2675   int isShared(void) const //Return 1 if this is a shared entity
2676     {return next!=NULL;}
2677   int isGhost(void) const  //Return 1 if this is a ghost entity
2678     {return FEM_Is_ghost_index(localNo); }
2679   int length(void) const {
2680     if (next==NULL) return isEmpty()?0:1;
2681     else return 1+next->length();
2682   }
2683   chunkList &operator[](int i) {
2684     if (i==0) return *this;
2685     else return (*next)[i-1];
2686   }
2687 };
2688
2689
2690 // end of map.h header file
2691
2692
2693 #include "ParFUM_locking.h"
2694
2695 /* File: util.h
2696  * Authors: Nilesh Choudhury
2697  * 
2698  */
2699 /// A utility class with helper functions for adaptivity
2700 class FEM_MUtil {
2701   ///chunk index
2702   int idx;
2703   ///cross-pointer to the femMeshModify object for this chunk
2704   femMeshModify *mmod;
2705   ///a data structure to remember mappings from oldnode to newnode yet to be updated
2706   class tuple {
2707   public:
2708     ///old index of a node
2709     int oldIdx;
2710     ///the corresponding new index
2711     int newIdx;
2712     ///default constructor
2713     tuple() {
2714       oldIdx = -1;
2715       newIdx = -1;
2716     }
2717     ///constructor
2718     tuple(int o, int n) {
2719       oldIdx = o;
2720       newIdx = n;
2721     }
2722   };
2723   /// vector of pairs of 'to-be-updated' mappings
2724   CkVec<tuple> outStandingMappings;
2725
2726  public:
2727   ///default constructor
2728   FEM_MUtil() {}
2729   ///constructor
2730   FEM_MUtil(int i, femMeshModify *m);
2731   ///constructor
2732   FEM_MUtil(femMeshModify *m);
2733   ///destructor
2734   ~FEM_MUtil();
2735   ///Pup for this object
2736   void pup(PUP::er &p) {
2737     p|idx;
2738   }
2739   ///Return the chunk index this util object is attached to
2740   int getIdx() { return idx; }
2741   ///Get the shared ghost recv index in IDXL list for this element on one chunk
2742   int getRemoteIdx(FEM_Mesh *m, int elementid, int elemtype);
2743   ///Is this node shared among chunks
2744   bool isShared(int index);
2745   ///Returns the IDXL entry index of this 'entity' on IDXL list specified by type
2746   int exists_in_IDXL(FEM_Mesh *m, int localIdx, int chk, int type, int elemType=0);
2747   ///Read this shared entry 'sharedIdx' in the IDXL list specified by 'type'
2748   int lookup_in_IDXL(FEM_Mesh *m, int sharedIdx, int fromChk, int type, int elemType=0);
2749
2750   ///Get all chunks that own this element/node
2751   /** The entType signifies what type of entity to lock. node=0, elem=1;
2752       entNo signifies the local index of the entity
2753       numChunks is the number of chunks that need to be locked to lock that entity
2754       chunks identifies the chunks that need to be locked */
2755   void getChunkNos(int entType, int entNo, int *numChunks, IDXL_Share ***chunks, int elemType=0);
2756   ///build a table of mappings of the which chunk knows about which node (given by conn)
2757   void buildChunkToNodeTable(int *nodetype, int sharedcount, int ghostcount, int localcount, int *conn, int connSize, CkVec<int> ***allShared, int *numSharedChunks, CkVec<int> **allChunks, int ***sharedConn);
2758   ///The set of chunks to which this node is sent as a ghost
2759   chunkListMsg *getChunksSharingGhostNodeRemote(FEM_Mesh *m, int chk, int sharedIdx);
2760
2761   ///Add this node to all the shared idxl lists on all chunks that share this node
2762   void splitEntityAll(FEM_Mesh *m, int localIdx, int nBetween, int *between);
2763   ///adds this new shared node to chunks only which share this edge(n=2) or face(n=3)
2764   void splitEntitySharing(FEM_Mesh *m, int localIdx, int nBetween, int *between, int numChunks, int *chunks);
2765   ///Add the shared node to IDXL list and populate its data, remote helper of above function
2766   void splitEntityRemote(FEM_Mesh *m, int chk, int localIdx, int nBetween, int *between);
2767
2768   ///Delete this node on all remote chunks (present as shared or ghost)
2769   void removeNodeAll(FEM_Mesh *m, int localIdx);
2770   ///Remove this node from shared IDXL list and delete node (remote helper of above function)
2771   void removeNodeRemote(FEM_Mesh *m, int chk, int sharedIdx);
2772   ///Remove this ghost node on this chunk (called from remote chunk)
2773   void removeGhostNodeRemote(FEM_Mesh *m, int fromChk, int sharedIdx);
2774
2775   ///Add this element on this chunk (called from a remote chunk)
2776   int addElemRemote(FEM_Mesh *m, int chk, int elemtype, int connSize, int *conn, int numGhostIndex, int *ghostIndices);
2777   ///Add the element with this conn (indices, typeOfIndices) on this chunk (called from remote chunk)
2778   void addGhostElementRemote(FEM_Mesh *m, int chk, int elemType, int *indices, int *typeOfIndex, int connSize);
2779   ///Remove this element on this chunk (called from remote chunk)
2780   void removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent);
2781   ///Remove this ghost element on this chunk, also delete some ghost nodes and elements
2782   void removeGhostElementRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int numGhostIndex, int *ghostIndices, int numGhostRNIndex, int *ghostRNIndices, int numGhostREIndex, int *ghostREIndices, int numSharedIndex, int *sharedIndices);
2783
2784   ///While acquiring a node (replace the oldIdx (ghost) with the newIdx)
2785   int Replace_node_local(FEM_Mesh *m, int oldIdx, int newIdx);
2786   ///Add this node to the shared Idxl list (called from remote chunk)
2787   void addToSharedList(FEM_Mesh *m, int fromChk, int sharedIdx);
2788   ///Acquire the element specified by this ghost index
2789   int eatIntoElement(int localIdx);
2790
2791   ///Does this chunk know about this node index (on any idxl list with 'chk')
2792   bool knowsAbtNode(int chk, int nodeId);
2793   ///get rid of unnecessary node sends to some chunks for a node
2794   void UpdateGhostSend(int nodeId, int *chunkl, int numchunkl);
2795   ///find the set of chunks where this node should be sent as a ghost
2796   void findGhostSend(int nodeId, int *&chunkl, int &numchunkl);
2797
2798   ///copies the elem data from elemid to newEl
2799   void copyElemData(int etype, int elemid, int newEl);
2800   ///Pack the data from this element/node and return it
2801   void packEntData(char **data, int *size, int *cnt, int localIdx, bool isnode, int elemType);
2802   ///update the attributes of 'newIndex' with this data
2803   void updateAttrs(char *data, int size, int newIndex, bool isnode, int elemType);
2804
2805   ///Return the owner of this lock (specified by nodeid)
2806   int getLockOwner(int nodeId);
2807
2808   ///Lock the idxl list with 'chk' (blocking call) (might be remote, if chk is smaller)
2809   void idxllock(FEM_Mesh *m, int chk, int type);
2810   ///Unlock the idxl list with 'chk' (might be remote if chk is smaller)
2811   void idxlunlock(FEM_Mesh *m, int chk, int type);
2812   ///Lock the idxl list with chk on this chunk 
2813   void idxllockLocal(FEM_Mesh *m, int toChk, int type);
2814   ///Unlock the idxl list with chk on this chunk 
2815   void idxlunlockLocal(FEM_Mesh *m, int toChk, int type);
2816
2817   ///Print the node-to-node adjacency for this node
2818   void FEM_Print_n2n(FEM_Mesh *m, int nodeid);
2819   ///Print the node-to-element adjacency for this node
2820   void FEM_Print_n2e(FEM_Mesh *m, int nodeid);
2821   ///Print the element-to-node adjacency for this element
2822   void FEM_Print_e2n(FEM_Mesh *m, int eid);
2823   ///Print the element-to-element adjacency for this element
2824   void FEM_Print_e2e(FEM_Mesh *m, int eid);
2825   ///Print the coords and boundary for this node
2826   void FEM_Print_coords(FEM_Mesh *m, int nodeid);
2827
2828   ///A bunch of tests to verify that connectivity and geometric structure of mesh is appropriate
2829   void StructureTest(FEM_Mesh *m);
2830   ///Test if the area of all elements is more than zero
2831   int AreaTest(FEM_Mesh *m);
2832   ///Test if the Idxl lists are consistent on all chunks
2833   int IdxlListTest(FEM_Mesh *m);
2834   ///Remote component of the above test (helper)
2835   void verifyIdxlListRemote(FEM_Mesh *m, int fromChk, int fsize, int type);
2836   ///Test that there are no remaining acquired locks on this mesh
2837   int residualLockTest(FEM_Mesh *m);
2838 };
2839
2840
2841
2842 // End util.h
2843 /* File: ParFUM_adapt_lock.h
2844  * Authors: Nilesh Choudhury, Terry Wilmarth
2845  *
2846  */
2847
2848
2849 #include "ParFUM_Mesh_Modify.h"
2850
2851 #include "ParFUM_Adapt_Algs.h"
2852
2853 #include "ParFUM_Adapt.h"
2854
2855 #include "ParFUM_Interpolate.h"
2856
2857 /*
2858   Orion's Standard Library
2859   Orion Sky Lawlor, 2/22/2000
2860   NAME:         vector2d.h
2861
2862   DESCRIPTION:  C++ 2-Dimentional vector library (no templates)
2863
2864   This file provides various utility routines for easily
2865   manipulating 2-D vectors-- included are arithmetic,
2866   dot product, magnitude and normalization terms. 
2867   All routines are provided right in the header file (for inlining).
2868
2869   Converted from vector3d.h.
2870
2871 */
2872
2873 #ifndef __OSL_VECTOR_2D_H
2874 #define __OSL_VECTOR_2D_H
2875
2876 #include <math.h>
2877
2878 typedef double Real;
2879
2880 //vector2d is a cartesian vector in 2-space-- an x and y.
2881 class vector2d {
2882  public:
2883   Real x,y;
2884   vector2d(void) {}//Default consructor
2885   //Simple 1-value constructor
2886   explicit vector2d(const Real init) {x=y=init;}
2887   //Simple 1-value constructor
2888   explicit vector2d(int init) {x=y=init;}
2889   //2-value constructor
2890   vector2d(const Real Nx,const Real Ny) {x=Nx;y=Ny;}
2891   //Copy constructor
2892   vector2d(const vector2d &copy) {x=copy.x;y=copy.y;}
2893         
2894   //Cast-to-Real * operators (treat vector as array)
2895   operator Real *() {return &x;}
2896   operator const Real *() const {return &x;}
2897         
2898   /*Arithmetic operations: these are carefully restricted to just those
2899     that make unambiguous sense (to me... now...  ;-)
2900     Counterexamples: vector*vector makes no sense (use .dot()) because
2901     Real/vector is meaningless (and we'd want a*b/b==a for b!=0), 
2902     ditto for vector&vector (dot?), vector|vector (projection?), 
2903     vector^vector (cross?),Real+vector, vector+=real, etc.
2904   */
2905   vector2d &operator=(const vector2d &b) {x=b.x;y=b.y;return *this;}
2906   int operator==(const vector2d &b) const {return (x==b.x)&&(y==b.y);}
2907   int operator!=(const vector2d &b) const {return (x!=b.x)||(y!=b.y);}
2908   vector2d operator+(const vector2d &b) const {return vector2d(x+b.x,y+b.y);}
2909   vector2d operator-(const vector2d &b) const {return vector2d(x-b.x,y-b.y);}
2910   vector2d operator*(const Real scale) const 
2911     {return vector2d(x*scale,y*scale);}
2912   friend vector2d operator*(const Real scale,const vector2d &v)
2913     {return vector2d(v.x*scale,v.y*scale);}
2914   vector2d operator/(const Real &div) const
2915     {Real scale=1.0/div;return vector2d(x*scale,y*scale);}
2916   vector2d operator-(void) const {return vector2d(-x,-y);}
2917   void operator+=(const vector2d &b) {x+=b.x;y+=b.y;}
2918   void operator-=(const vector2d &b) {x-=b.x;y-=b.y;}
2919   void operator*=(const Real scale) {x*=scale;y*=scale;}
2920   void operator/=(const Real div) {Real scale=1.0/div;x*=scale;y*=scale;}
2921
2922   //Vector-specific operations
2923   //Return the square of the magnitude of this vector
2924   Real magSqr(void) const {return x*x+y*y;}
2925   //Return the magnitude (length) of this vector
2926   Real mag(void) const {return sqrt(magSqr());}
2927         
2928   //Return the square of the distance to the vector b
2929   Real distSqr(const vector2d &b) const 
2930     {return (x-b.x)*(x-b.x)+(y-b.y)*(y-b.y);}
2931   //Return the distance to the vector b
2932   Real dist(const vector2d &b) const {return sqrt(distSqr(b));}
2933         
2934   //Return the dot product of this vector and b
2935   Real dot(const vector2d &b) const {return x*b.x+y*b.y;}
2936   //Return the cosine of the angle between this vector and b
2937   Real cosAng(const vector2d &b) const {return dot(b)/(mag()*b.mag());}
2938         
2939   //Return the "direction" (unit vector) of this vector
2940   vector2d dir(void) const {return (*this)/mag();}
2941
2942   //Return the CCW perpendicular vector
2943   vector2d perp(void) const {return vector2d(-y,x);}
2944
2945   //Return this vector scaled by that
2946   vector2d &scale(const vector2d &b) {x*=b.x;y*=b.y;return *this;}
2947         
2948   //Return the largest coordinate in this vector
2949   Real max(void) {return (x>y)?x:y;}
2950   //Make each of this vector's coordinates at least as big
2951   // as the given vector's coordinates.
2952   void enlarge(const vector2d &by)
2953     {if (by.x>x) x=by.x; if (by.y>y) y=by.y;}
2954 };
2955
2956 #endif //__OSL_VECTOR2D_H
2957
2958 /*
2959  * The former cktimer.h
2960  *
2961  */ 
2962
2963 #ifndef CMK_THRESHOLD_TIMER
2964 #define CMK_THRESHOLD_TIMER
2965
2966 /** Time a sequence of operations, printing out the
2967     names and times of any operations that exceed a threshold. 
2968
2969     Use it with only the constructor and destructor like:
2970     void foo(void) {
2971     CkThresholdTimer t("foo");
2972     ...
2973     }
2974     (this times the whole execution of the routine,
2975     all the way to t's destructor is called on function return)
2976
2977     Or, you can start different sections like:
2978     void bar(void) {
2979     CkThresholdTimer t("first");
2980     ...
2981     t.start("second");
2982     ...
2983     t.start("third");
2984     ...
2985     }
2986   
2987     This class *only* prints out the time if it exceeds
2988     a threshold-- by default, one millisecond.
2989 */
2990 class CkThresholdTimer {
2991   double threshold; // Print any times that exceed this (s).
2992   double lastStart; // Last activity started at this time (s).
2993   const char *lastWhat; // Last activity has this name.
2994         
2995   void start_(const char *what) {
2996     lastStart=CmiWallTimer();
2997     lastWhat=what;
2998   }
2999   void done_(void) {
3000     double elapsed=CmiWallTimer()-lastStart;
3001     if (elapsed>threshold) {
3002       CmiPrintf("%s took %.2f s\n",lastWhat,elapsed);
3003     }
3004   }
3005  public:
3006   CkThresholdTimer(const char *what,double thresh=0.001) 
3007     :threshold(thresh) { start_(what); }
3008   void start(const char *what) { done_(); start_(what); }
3009   ~CkThresholdTimer() {done_();}
3010 };
3011
3012 #endif
3013 // end cktimer or ParFUM_timer
3014
3015 #include "adapt_adj.h"
3016
3017 #endif
3018 // end ParFUM_internals.h
3019
3020 /*@}*/