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