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