MSA & ParFUM: Pass templating of MSA type args through
[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     template<typename U>
2054     static inline void accumulate(T& a, const U& b) { a += b; }
2055     // identity for initializing at start of accumulate
2056     static inline T getIdentity() { return T(); }
2057     static inline bool pupEveryElement(){ return PUP_EVERY_ELEMENT; }
2058   };
2059
2060 extern double elemlistaccTime;
2061
2062 template <class T>
2063 class ElemList{
2064  public:
2065   CkVec<T> *vec;
2066   ElemList(){
2067     vec = new CkVec<T>();
2068   }
2069   ~ElemList(){
2070     delete vec;
2071   }
2072   ElemList(const ElemList &rhs){
2073     vec = new CkVec<T>();
2074     *this=rhs;
2075   }
2076   inline ElemList& operator=(const ElemList& rhs){
2077     //           vec = new CkVec<T>();
2078     *vec = *(rhs.vec);
2079                 return *this;
2080   }
2081   inline ElemList& operator+=(const ElemList& rhs){
2082     /*
2083       add the new unique elements to the List
2084     */
2085     double _start = CkWallTimer();
2086     for(int i=0;i<rhs.vec->length();i++){
2087       vec->push_back((*(rhs.vec))[i]);
2088     }
2089     //          uniquify();
2090     elemlistaccTime += (CkWallTimer() - _start);
2091     return *this;
2092   }
2093   ElemList(const T &val){
2094     vec =new CkVec<T>();
2095     vec->push_back(val);
2096   };
2097   inline virtual void pup(PUP::er &p){
2098     if(p.isUnpacking()){
2099       vec = new CkVec<T>();
2100     }
2101     pupCkVec(p,*vec);
2102   }
2103 };
2104 template <class T>
2105 class UniqElemList: public ElemList<T>{
2106 public:
2107   UniqElemList(const T &val):ElemList<T>(val){};
2108   UniqElemList():ElemList<T>(){};
2109   inline void uniquify(){
2110     CkVec<T> *lvec = this->vec;
2111     lvec->quickSort(8);
2112     if(lvec->length() != 0){
2113       int count=0;
2114       for(int i=1;i<lvec->length();i++){
2115         if((*lvec)[count] == (*lvec)[i]){
2116         }else{
2117           count++;
2118           if(i != count){
2119             (*lvec)[count] = (*lvec)[i];
2120           }
2121         }
2122       }
2123       lvec->resize(count+1);
2124     }
2125   }
2126 };
2127
2128 class NodeElem {
2129  public:
2130   //global number of this node
2131   int global;
2132   /*no of chunks that share this node
2133     owned by 1 element - numShared 0
2134     owned by 2 elements - numShared 2
2135   */
2136   int numShared;
2137   /*somewhat horrible semantics
2138     -if numShared == 0 shared is NULL
2139     - else stores all the chunks that share this node
2140   */
2141   int shared[MAX_SHARED_PER_NODE];
2142   //    int *shared;
2143   NodeElem(){
2144     global = -1;
2145     numShared = 0;
2146     //          shared = NULL;
2147   }
2148   NodeElem(int g_,int num_){
2149     global = g_; numShared= num_;
2150     /*          if(numShared != 0){
2151                 shared = new int[numShared];
2152                 }else{
2153                 shared = NULL;
2154                 }*/
2155     if(numShared > MAX_SHARED_PER_NODE){
2156       CkAbort("In Parallel partition the max number of chunks per node exceeds limit\n");
2157     }
2158   }
2159   NodeElem(int g_){
2160     global = g_;
2161     numShared = 0;
2162     //          shared = NULL;
2163   }
2164   NodeElem(const NodeElem &rhs){
2165     //          shared=NULL;
2166     *this = rhs;
2167   }
2168   inline NodeElem& operator=(const NodeElem &rhs){
2169     global = rhs.global;
2170     numShared = rhs.numShared;
2171     /*          if(shared != NULL){
2172                 delete [] shared;
2173                 }
2174                 shared = new int[numShared];*/
2175     memcpy(&shared[0],&((rhs.shared)[0]),numShared*sizeof(int));
2176     return *this;
2177   }
2178
2179   inline        bool operator == (const NodeElem &rhs){
2180     if(global == rhs.global){
2181       return true;
2182     }else{
2183       return false;
2184     }
2185   }
2186   inline        bool operator >=(const NodeElem &rhs){
2187     return (global >= rhs.global);
2188   }
2189
2190   inline bool operator <=(const NodeElem &rhs){
2191     return (global <= rhs.global);
2192   }
2193
2194   virtual void pup(PUP::er &p){
2195     p | global;
2196     p | numShared;
2197     /*          if(p.isUnpacking()){
2198                 if(numShared != 0){
2199                 shared = new int[numShared];
2200                 }else{
2201                 shared = NULL;
2202                 }
2203                 }*/
2204     p(shared,numShared);
2205   }
2206   ~NodeElem(){
2207     /*          if(shared != NULL){
2208                 delete [] shared;
2209                 }*/
2210   }
2211 };
2212
2213 /**
2214   This class is an MSA Entity. It is used for 2 purposes
2215   1 It is used for storing the mesh while creating the mesh
2216   for each chunk
2217   2     It is used for storing the ghost elements and nodes for
2218   a chunk.
2219 */
2220 class MeshElem{
2221  public:
2222   FEM_Mesh *m;
2223   CkVec<int> gedgechunk; // Chunk number of
2224   MeshElem(){
2225     m = new FEM_Mesh;
2226   }
2227   MeshElem(int dummy){
2228     m = new FEM_Mesh;
2229   }
2230   ~MeshElem(){
2231     delete m;
2232   }
2233   MeshElem(const MeshElem &rhs){
2234     m = NULL;
2235     *this = rhs;
2236   }
2237   inline MeshElem& operator=(const MeshElem &rhs){
2238     if(m != NULL){
2239       delete m;
2240     }
2241     m = new FEM_Mesh;
2242     m->copyShape(*(rhs.m));
2243     (*this) += rhs;
2244
2245     return *this;
2246   }
2247   inline MeshElem& operator+=(const MeshElem &rhs){
2248     int oldel = m->nElems();
2249     m->copyShape(*(rhs.m));
2250     for(int i=0;i<rhs.m->node.size();i++){
2251       m->node.push_back((rhs.m)->node,i);
2252     }
2253     if((rhs.m)->elem.size()>0){
2254       for(int t=0;t<(rhs.m)->elem.size();t++){
2255         if((rhs.m)->elem.has(t)){
2256           for(int e=0;e<(rhs.m)->elem.get(t).size();e++){
2257             m->elem[t].push_back((rhs.m)->elem.get(t),e);
2258           }
2259         }
2260       }
2261     }
2262
2263     return *this;
2264   }
2265   virtual void pup(PUP::er &p){
2266     if(p.isUnpacking()){
2267       m = new FEM_Mesh;
2268     }
2269     m->pup(p);
2270   }
2271 };
2272
2273
2274 class Hashnode{
2275 public:
2276         class tupledata{
2277                 public:
2278                 enum {MAX_TUPLE = 8};
2279                         int nodes[MAX_TUPLE];
2280                         tupledata(int _nodes[MAX_TUPLE]){
2281                                 memcpy(nodes,_nodes,sizeof(int)*MAX_TUPLE);
2282                         }
2283                         tupledata(tupledata &rhs){
2284                                 memcpy(nodes,rhs.nodes,sizeof(int)*MAX_TUPLE);
2285                         }
2286                         tupledata(const tupledata &rhs){
2287                                 memcpy(nodes,rhs.nodes,sizeof(int)*MAX_TUPLE);
2288                         }
2289                         tupledata(){};
2290                         //dont store the returned string
2291                         char *toString(int numnodes,char *str){
2292                                 str[0]='\0';
2293                                 for(int i=0;i<numnodes;i++){
2294                                         sprintf(&str[strlen(str)],"%d ",nodes[i]);
2295                                 }
2296                                 return str;
2297                         }
2298                         inline int &operator[](int i){
2299                                 return nodes[i];
2300                         }
2301                         inline const int &operator[](int i) const {
2302                                 return nodes[i];
2303                         }
2304                         virtual void pup(PUP::er &p){
2305                                 p(nodes,MAX_TUPLE);
2306                         }
2307         };
2308         int numnodes; //number of nodes in this tuple
2309         //TODO: replace *nodes with the above tupledata class
2310         tupledata nodes;        //the nodes in the tuple
2311         int chunk;              //the chunk number to which this element belongs
2312         int elementNo;          //local number of that element
2313         Hashnode(){
2314                 numnodes=0;
2315         };
2316         Hashnode(int _num,int _chunk,int _elNo,int _nodes[tupledata::MAX_TUPLE]): nodes(_nodes){
2317                 numnodes = _num;
2318                 chunk = _chunk;
2319                 elementNo = _elNo;
2320         }
2321         Hashnode(const Hashnode &rhs){
2322                 *this = rhs;
2323         }
2324         inline Hashnode &operator=(const Hashnode &rhs){
2325                 numnodes = rhs.numnodes;
2326                 for(int i=0;i<numnodes;i++){
2327                         nodes[i] = rhs.nodes[i];
2328                 }
2329                 chunk = rhs.chunk;
2330                 elementNo = rhs.elementNo;
2331                 return *this;
2332         }
2333         inline bool operator==(const Hashnode &rhs){
2334                 if(numnodes != rhs.numnodes){
2335                         return false;
2336                 }
2337                 for(int i=0;i<numnodes;i++){
2338                         if(nodes[i] != rhs.nodes[i]){
2339                                 return false;
2340                         }
2341                 }
2342                 if(chunk != rhs.chunk){
2343                         return false;
2344                 }
2345                 if(elementNo != rhs.elementNo){
2346                         return false;
2347                 }
2348                 return true;
2349         }
2350         inline bool operator>=(const Hashnode &rhs){
2351                 if(numnodes < rhs.numnodes){
2352                         return false;
2353                 };
2354                 if(numnodes > rhs.numnodes){
2355                         return true;
2356                 }
2357
2358     for(int i=0;i<numnodes;i++){
2359       if(nodes[i] < rhs.nodes[i]){
2360         return false;
2361       }
2362       if(nodes[i] > rhs.nodes[i]){
2363         return true;
2364       }
2365     }
2366     if(chunk < rhs.chunk){
2367       return false;
2368     }
2369     if(chunk > rhs.chunk){
2370       return true;
2371     }
2372     if(elementNo < rhs.elementNo){
2373       return false;
2374     }
2375     if(elementNo > rhs.elementNo){
2376       return true;
2377     }
2378     return true;
2379   }
2380
2381   inline bool operator<=(const Hashnode &rhs){
2382     if(numnodes < rhs.numnodes){
2383       return true;
2384     };
2385     if(numnodes > rhs.numnodes){
2386       return false;
2387     }
2388
2389     for(int i=0;i<numnodes;i++){
2390       if(nodes[i] < rhs.nodes[i]){
2391         return true;
2392       }
2393       if(nodes[i] > rhs.nodes[i]){
2394         return false;
2395       }
2396     }
2397     if(chunk < rhs.chunk){
2398       return true;
2399     }
2400     if(chunk > rhs.chunk){
2401       return false;
2402     }
2403     if(elementNo < rhs.elementNo){
2404       return true;
2405     }
2406     if(elementNo > rhs.elementNo){
2407       return false;
2408     }
2409     return true;
2410   }
2411
2412   inline bool equals(tupledata &tuple){
2413     for(int i=0;i<numnodes;i++){
2414       if(tuple.nodes[i] != nodes[i]){
2415         return false;
2416       }
2417     }
2418     return true;
2419   }
2420   virtual void pup(PUP::er &p){
2421     p | numnodes;
2422     p | nodes;
2423     p | chunk;
2424     p | elementNo;
2425   }
2426 };
2427
2428 typedef MSA2D<int, DefaultEntry<int>, MSA_DEFAULT_ENTRIES_PER_PAGE, MSA_ROW_MAJOR> MSA2DRM;
2429
2430 typedef MSA1D<int, DefaultEntry<int>, MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DINT;
2431
2432 typedef UniqElemList<int> IntList;
2433 typedef MSA1D<IntList, DefaultListEntry<IntList,true>,MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DINTLIST;
2434
2435 typedef UniqElemList<NodeElem> NodeList;
2436 typedef MSA1D<NodeList, DefaultListEntry<NodeList,true>,MSA_DEFAULT_ENTRIES_PER_PAGE> MSA1DNODELIST;
2437
2438 typedef MSA1D<MeshElem,DefaultEntry<MeshElem,true>,1> MSA1DFEMMESH;
2439
2440
2441
2442 struct conndata{
2443   int nelem;
2444   int nnode;
2445   MSA1DINT arr1;
2446   MSA1DINT arr2;
2447
2448   void pup(PUP::er &p){
2449     p|nelem;
2450     p|nnode;
2451     arr1.pup(p);
2452     arr2.pup(p);
2453   }
2454 };
2455
2456 /**
2457   Structure to store connectivity data after the
2458   global element partition has been returned by parmetis
2459 */
2460 struct partconndata{
2461   int nelem;
2462   int startindex;
2463   int *eptr,*eind;
2464   int *part;
2465   ~partconndata(){
2466     delete [] eptr;
2467     delete [] eind;
2468     delete [] part;
2469   };
2470 };
2471
2472 /**
2473   structure for storing the ghost layers
2474 */
2475 struct ghostdata{
2476   int numLayers;
2477   FEM_Ghost_Layer **layers;
2478   void pup(PUP::er &p){
2479     p | numLayers;
2480     if(p.isUnpacking()){
2481       layers = new FEM_Ghost_Layer *[numLayers];
2482       for(int i=0;i<numLayers;i++){
2483         layers[i] = new FEM_Ghost_Layer;
2484       }
2485
2486     }
2487     for(int i=0;i<numLayers;i++){
2488       layers[i]->pup(p);
2489     }
2490   }
2491   ~ghostdata(){
2492     //printf("destructor on ghostdata called \n");
2493     for(int i=0;i<numLayers;i++){
2494       delete layers[i];
2495     }
2496     delete [] layers;
2497   };
2498 };
2499
2500
2501
2502
2503 int FEM_master_parallel_part(int ,int ,FEM_Comm_t);
2504 int FEM_slave_parallel_part(int ,int ,FEM_Comm_t);
2505 struct partconndata* FEM_call_parmetis(int nelem, MSA1DINT::Read &rPtr, MSA1DINT::Read &rInd, FEM_Comm_t comm_context);
2506 void FEM_write_nodepart(MSA1DINTLIST::Accum &nodepart, struct partconndata *data, MPI_Comm comm_context);
2507 void FEM_write_part2node(MSA1DINTLIST::Read &nodepart, MSA1DNODELIST::Accum &part2node, struct partconndata *data, MPI_Comm comm_context);
2508 void FEM_write_part2elem(MSA1DINTLIST::Accum &part2elem, struct partconndata *data, MPI_Comm comm_context);
2509 FEM_Mesh * FEM_break_mesh(FEM_Mesh *m,int numElements,int numChunks);
2510 void sendBrokenMeshes(FEM_Mesh *mesh_array,FEM_Comm_t comm_context);
2511 void FEM_write_part2mesh(MSA1DFEMMESH::Accum &part2mesh, struct partconndata *partdata, struct conndata *data, MSA1DINTLIST::Read &nodepart, int numChunks, int myChunk, FEM_Mesh *mypiece);
2512 void addIDXLists(FEM_Mesh *m,NodeList &lnodes,int myChunk);
2513 struct ghostdata *gatherGhosts();
2514 void makeGhosts(FEM_Mesh *m,MPI_Comm comm,int masterRank,int numLayers,FEM_Ghost_Layer **layers);
2515 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);
2516 bool sharedWith(int lnode,int chunk,FEM_Mesh *m);
2517
2518 // End Parallel Partitioner
2519
2520
2521
2522 /*
2523   Some definitions of structures/classes used in map.C and possibly other FEM source files.
2524   Moved to this file by Isaac Dooley 4/5/05
2525
2526   FIXME: Clean up the organization of this file, and possibly move other structures from Map.C here
2527   Also all the stuff in this file might just belong in impl.h. I just moved it here so it could
2528   be included in fem.C for my element->element build adjacency table function.
2529 */
2530
2531 static FEM_Symmetries_t noSymmetries=(FEM_Symmetries_t)0;
2532
2533 // defined in map.C
2534 CkHashCode CkHashFunction_ints(const void *keyData,size_t keyLen);
2535 int CkHashCompare_ints(const void *k1,const void *k2,size_t keyLen);
2536 extern "C" int ck_fem_map_compare_int(const void *a, const void *b);
2537
2538 //A linked list of elements surrounding a tuple.
2539 //ElemLists are allocated one at a time, and hence a bit cleaner than chunklists:
2540 class elemList {
2541  public:
2542   int chunk;
2543   int tupleNo;//tuple number on this element for the tuple that this list sorrounds
2544   int localNo;//Local number of this element on this chunk (negative for a ghost)
2545   int type; //Kind of element
2546   FEM_Symmetries_t sym; //Symmetries this element was reached via
2547   elemList *next;
2548
2549   elemList(int chunk_,int localNo_,int type_,FEM_Symmetries_t sym_)
2550     :chunk(chunk_),localNo(localNo_),type(type_), sym(sym_)
2551     { next=NULL; }
2552   elemList(int chunk_,int localNo_,int type_,FEM_Symmetries_t sym_,int tupleNo_)
2553     :chunk(chunk_),localNo(localNo_),type(type_), sym(sym_) , tupleNo(tupleNo_)
2554     { next=NULL; }
2555
2556   ~elemList() {if (next) delete next;}
2557   void setNext(elemList *n) {next=n;}
2558 };
2559
2560
2561 //Maps node tuples to element lists
2562 class tupleTable : public CkHashtable {
2563   int tupleLen; //Nodes in a tuple
2564   CkHashtableIterator *it;
2565   static int roundUp(int val,int to) {
2566     return ((val+to-1)/to)*to;
2567   }
2568   static CkHashtableLayout makeLayout(int tupleLen) {
2569     int ks=tupleLen*sizeof(int);
2570     int oo=roundUp(ks+sizeof(char),sizeof(void *));
2571     int os=sizeof(elemList *);
2572     return CkHashtableLayout(ks,ks,oo,os,oo+os);
2573   }
2574
2575   //Make a canonical version of this tuple, so different
2576   // orderings of the same nodes don't end up in different lists.
2577   //I canonicalize by sorting:
2578   void canonicalize(const int *tuple,int *can)
2579     {
2580       switch(tupleLen) {
2581       case 1: //Short lists are easy to sort:
2582         can[0]=tuple[0]; break;
2583       case 2:
2584         if (tuple[0]<tuple[1])
2585           {can[0]=tuple[0]; can[1]=tuple[1];}
2586         else
2587           {can[0]=tuple[1]; can[1]=tuple[0];}
2588         break;
2589       default: //Should use std::sort here:
2590         memcpy(can,tuple,tupleLen*sizeof(int));
2591         qsort(can,tupleLen,sizeof(int),ck_fem_map_compare_int);
2592       };
2593     }
2594  public:
2595   enum {MAX_TUPLE=8};
2596
2597   tupleTable(int tupleLen_)
2598     :CkHashtable(makeLayout(tupleLen_),
2599                  137,0.5,
2600                  CkHashFunction_ints,
2601                  CkHashCompare_ints)
2602     {
2603       tupleLen=tupleLen_;
2604       if (tupleLen>MAX_TUPLE) CkAbort("Cannot have that many shared nodes!\n");
2605       it=NULL;
2606     }
2607   ~tupleTable() {
2608     beginLookup();
2609     elemList *doomed;
2610     while (NULL!=(doomed=(elemList *)lookupNext()))
2611       delete doomed;
2612   }
2613   //Lookup the elemList associated with this tuple, or return NULL
2614   elemList **lookupTuple(const int *tuple) {
2615     int can[MAX_TUPLE];
2616     canonicalize(tuple,can);
2617     return (elemList **)get(can);
2618   }
2619
2620   //Register this (new'd) element with this tuple
2621   void addTuple(const int *tuple,elemList *nu)
2622     {
2623       int can[MAX_TUPLE];
2624       canonicalize(tuple,can);
2625       //First try for an existing list:
2626       elemList **dest=(elemList **)get(can);
2627       if (dest!=NULL)
2628         { //A list already exists here-- link it into the new list
2629           nu->setNext(*dest);
2630         } else {//No pre-existing list-- initialize a new one.
2631           dest=(elemList **)put(can);
2632         }
2633       *dest=nu;
2634     }
2635   //Return all registered elemLists:
2636   void beginLookup(void) {
2637     it=iterator();
2638   }
2639   elemList *lookupNext(void) {
2640     void *ret=it->next();
2641     if (ret==NULL) {
2642       delete it;
2643       return NULL;
2644     }
2645     return *(elemList **)ret;
2646   }
2647 };
2648
2649
2650
2651
2652 /*This object maps a single entity to a list of the chunks
2653   that have a copy of it.  For the vast majority
2654   of entities, the list will contain exactly one element.
2655 */
2656 class chunkList : public CkNoncopyable {
2657  public:
2658   int chunk;//Chunk number; or -1 if the list is empty
2659   int localNo;//Local number of this entity on this chunk (if negative, is a ghost)
2660   FEM_Symmetries_t sym; //Symmetries this entity was reached via
2661   int layerNo; // -1 if real; if a ghost, our ghost layer number
2662   chunkList *next;
2663   chunkList() {chunk=-1;next=NULL;}
2664   chunkList(int chunk_,int localNo_,FEM_Symmetries_t sym_,int ln_=-1) {
2665     chunk=chunk_;
2666     localNo=localNo_;
2667     sym=sym_;
2668     layerNo=ln_;
2669     next=NULL;
2670   }
2671   ~chunkList() {delete next;}
2672   void set(int c,int l,FEM_Symmetries_t s,int ln) {
2673     chunk=c; localNo=l; sym=s; layerNo=ln;
2674   }
2675   //Is this chunk in the list?  If so, return false.
2676   // If not, add it and return true.
2677   bool addchunk(int c,int l,FEM_Symmetries_t s,int ln) {
2678     //Add chunk c to the list with local index l,
2679     // if it's not there already
2680     if (chunk==c && sym==s) return false;//Already in the list
2681     if (chunk==-1) {set(c,l,s,ln);return true;}
2682     if (next==NULL) {next=new chunkList(c,l,s,ln);return true;}
2683     else return next->addchunk(c,l,s,ln);
2684   }
2685   //Return this node's local number on chunk c (or -1 if none)
2686   int localOnChunk(int c,FEM_Symmetries_t s) const {
2687     const chunkList *l=onChunk(c,s);
2688     if (l==NULL) return -1;
2689     else return l->localNo;
2690   }
2691   const chunkList *onChunk(int c,FEM_Symmetries_t s) const {
2692     if (chunk==c && sym==s) return this;
2693     else if (next==NULL) return NULL;
2694     else return next->onChunk(c,s);
2695   }
2696   int isEmpty(void) const //Return 1 if this is an empty list
2697     {return (chunk==-1);}
2698   int isShared(void) const //Return 1 if this is a shared entity
2699     {return next!=NULL;}
2700   int isGhost(void) const  //Return 1 if this is a ghost entity
2701     {return FEM_Is_ghost_index(localNo); }
2702   int length(void) const {
2703     if (next==NULL) return isEmpty()?0:1;
2704     else return 1+next->length();
2705   }
2706   chunkList &operator[](int i) {
2707     if (i==0) return *this;
2708     else return (*next)[i-1];
2709   }
2710 };
2711
2712
2713 // end of map.h header file
2714
2715
2716 #include "ParFUM_locking.h"
2717
2718 /* File: util.h
2719  * Authors: Nilesh Choudhury
2720  *
2721  */
2722 /// A utility class with helper functions for adaptivity
2723 class FEM_MUtil {
2724   ///chunk index
2725   int idx;
2726   ///cross-pointer to the femMeshModify object for this chunk
2727   femMeshModify *mmod;
2728   ///a data structure to remember mappings from oldnode to newnode yet to be updated
2729   class tuple {
2730   public:
2731     ///old index of a node
2732     int oldIdx;
2733     ///the corresponding new index
2734     int newIdx;
2735     ///default constructor
2736     tuple() {
2737       oldIdx = -1;
2738       newIdx = -1;
2739     }
2740     ///constructor
2741     tuple(int o, int n) {
2742       oldIdx = o;
2743       newIdx = n;
2744     }
2745   };
2746   /// vector of pairs of 'to-be-updated' mappings
2747   CkVec<tuple> outStandingMappings;
2748
2749  public:
2750   ///default constructor
2751   FEM_MUtil() {}
2752   ///constructor
2753   FEM_MUtil(int i, femMeshModify *m);
2754   ///constructor
2755   FEM_MUtil(femMeshModify *m);
2756   ///destructor
2757   ~FEM_MUtil();
2758   ///Pup for this object
2759   void pup(PUP::er &p) {
2760     p|idx;
2761   }
2762   ///Return the chunk index this util object is attached to
2763   int getIdx() { return idx; }
2764   ///Get the shared ghost recv index in IDXL list for this element on one chunk
2765   int getRemoteIdx(FEM_Mesh *m, int elementid, int elemtype);
2766   ///Is this node shared among chunks
2767   bool isShared(int index);
2768   ///Returns the IDXL entry index of this 'entity' on IDXL list specified by type
2769   int exists_in_IDXL(FEM_Mesh *m, int localIdx, int chk, int type, int elemType=0);
2770   ///Read this shared entry 'sharedIdx' in the IDXL list specified by 'type'
2771   int lookup_in_IDXL(FEM_Mesh *m, int sharedIdx, int fromChk, int type, int elemType=0);
2772
2773   ///Get all chunks that own this element/node
2774   /** The entType signifies what type of entity to lock. node=0, elem=1;
2775       entNo signifies the local index of the entity
2776       numChunks is the number of chunks that need to be locked to lock that entity
2777       chunks identifies the chunks that need to be locked */
2778   void getChunkNos(int entType, int entNo, int *numChunks, IDXL_Share ***chunks, int elemType=0);
2779   ///build a table of mappings of the which chunk knows about which node (given by conn)
2780   void buildChunkToNodeTable(int *nodetype, int sharedcount, int ghostcount, int localcount, int *conn, int connSize, CkVec<int> ***allShared, int *numSharedChunks, CkVec<int> **allChunks, int ***sharedConn);
2781   ///The set of chunks to which this node is sent as a ghost
2782   chunkListMsg *getChunksSharingGhostNodeRemote(FEM_Mesh *m, int chk, int sharedIdx);
2783
2784   ///Add this node to all the shared idxl lists on all chunks that share this node
2785   void splitEntityAll(FEM_Mesh *m, int localIdx, int nBetween, int *between);
2786   ///adds this new shared node to chunks only which share this edge(n=2) or face(n=3)
2787   void splitEntitySharing(FEM_Mesh *m, int localIdx, int nBetween, int *between, int numChunks, int *chunks);
2788   ///Add the shared node to IDXL list and populate its data, remote helper of above function
2789   void splitEntityRemote(FEM_Mesh *m, int chk, int localIdx, int nBetween, int *between);
2790
2791   ///Delete this node on all remote chunks (present as shared or ghost)
2792   void removeNodeAll(FEM_Mesh *m, int localIdx);
2793   ///Remove this node from shared IDXL list and delete node (remote helper of above function)
2794   void removeNodeRemote(FEM_Mesh *m, int chk, int sharedIdx);
2795   ///Remove this ghost node on this chunk (called from remote chunk)
2796   void removeGhostNodeRemote(FEM_Mesh *m, int fromChk, int sharedIdx);
2797
2798   ///Add this element on this chunk (called from a remote chunk)
2799   int addElemRemote(FEM_Mesh *m, int chk, int elemtype, int connSize, int *conn, int numGhostIndex, int *ghostIndices);
2800   ///Add the element with this conn (indices, typeOfIndices) on this chunk (called from remote chunk)
2801   void addGhostElementRemote(FEM_Mesh *m, int chk, int elemType, int *indices, int *typeOfIndex, int connSize);
2802   ///Remove this element on this chunk (called from remote chunk)
2803   void removeElemRemote(FEM_Mesh *m, int chk, int elementid, int elemtype, int permanent);
2804   ///Remove this ghost element on this chunk, also delete some ghost nodes and elements
2805   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);
2806
2807   ///While acquiring a node (replace the oldIdx (ghost) with the newIdx)
2808   int Replace_node_local(FEM_Mesh *m, int oldIdx, int newIdx);
2809   ///Add this node to the shared Idxl list (called from remote chunk)
2810   void addToSharedList(FEM_Mesh *m, int fromChk, int sharedIdx);
2811   ///Acquire the element specified by this ghost index
2812   int eatIntoElement(int localIdx);
2813
2814   ///Does this chunk know about this node index (on any idxl list with 'chk')
2815   bool knowsAbtNode(int chk, int nodeId);
2816   ///get rid of unnecessary node sends to some chunks for a node
2817   void UpdateGhostSend(int nodeId, int *chunkl, int numchunkl);
2818   ///find the set of chunks where this node should be sent as a ghost
2819   void findGhostSend(int nodeId, int *&chunkl, int &numchunkl);
2820
2821   ///copies the elem data from elemid to newEl
2822   void copyElemData(int etype, int elemid, int newEl);
2823   ///Pack the data from this element/node and return it
2824   void packEntData(char **data, int *size, int *cnt, int localIdx, bool isnode, int elemType);
2825   ///update the attributes of 'newIndex' with this data
2826   void updateAttrs(char *data, int size, int newIndex, bool isnode, int elemType);
2827
2828   ///Return the owner of this lock (specified by nodeid)
2829   int getLockOwner(int nodeId);
2830
2831   ///Lock the idxl list with 'chk' (blocking call) (might be remote, if chk is smaller)
2832   void idxllock(FEM_Mesh *m, int chk, int type);
2833   ///Unlock the idxl list with 'chk' (might be remote if chk is smaller)
2834   void idxlunlock(FEM_Mesh *m, int chk, int type);
2835   ///Lock the idxl list with chk on this chunk
2836   void idxllockLocal(FEM_Mesh *m, int toChk, int type);
2837   ///Unlock the idxl list with chk on this chunk
2838   void idxlunlockLocal(FEM_Mesh *m, int toChk, int type);
2839
2840   ///Print the node-to-node adjacency for this node
2841   void FEM_Print_n2n(FEM_Mesh *m, int nodeid);
2842   ///Print the node-to-element adjacency for this node
2843   void FEM_Print_n2e(FEM_Mesh *m, int nodeid);
2844   ///Print the element-to-node adjacency for this element
2845   void FEM_Print_e2n(FEM_Mesh *m, int eid);
2846   ///Print the element-to-element adjacency for this element
2847   void FEM_Print_e2e(FEM_Mesh *m, int eid);
2848   ///Print the coords and boundary for this node
2849   void FEM_Print_coords(FEM_Mesh *m, int nodeid);
2850
2851   ///A bunch of tests to verify that connectivity and geometric structure of mesh is appropriate
2852   void StructureTest(FEM_Mesh *m);
2853   ///Test if the area of all elements is more than zero
2854   int AreaTest(FEM_Mesh *m);
2855   ///Test if the Idxl lists are consistent on all chunks
2856   int IdxlListTest(FEM_Mesh *m);
2857   ///Remote component of the above test (helper)
2858   void verifyIdxlListRemote(FEM_Mesh *m, int fromChk, int fsize, int type);
2859   ///Test that there are no remaining acquired locks on this mesh
2860   int residualLockTest(FEM_Mesh *m);
2861 };
2862
2863
2864
2865 // End util.h
2866 /* File: ParFUM_adapt_lock.h
2867  * Authors: Nilesh Choudhury, Terry Wilmarth
2868  *
2869  */
2870
2871
2872 #include "ParFUM_Mesh_Modify.h"
2873
2874 #include "ParFUM_Adapt_Algs.h"
2875
2876 #include "ParFUM_Adapt.h"
2877
2878 #include "ParFUM_Interpolate.h"
2879
2880 /*
2881   Orion's Standard Library
2882   Orion Sky Lawlor, 2/22/2000
2883   NAME:         vector2d.h
2884
2885   DESCRIPTION:  C++ 2-Dimentional vector library (no templates)
2886
2887   This file provides various utility routines for easily
2888   manipulating 2-D vectors-- included are arithmetic,
2889   dot product, magnitude and normalization terms.
2890   All routines are provided right in the header file (for inlining).
2891
2892   Converted from vector3d.h.
2893
2894 */
2895
2896 #ifndef __OSL_VECTOR_2D_H
2897 #define __OSL_VECTOR_2D_H
2898
2899 #include <math.h>
2900
2901 typedef double Real;
2902
2903 //vector2d is a cartesian vector in 2-space-- an x and y.
2904 class vector2d {
2905  public:
2906   Real x,y;
2907   vector2d(void) {}//Default consructor
2908   //Simple 1-value constructor
2909   explicit vector2d(const Real init) {x=y=init;}
2910   //Simple 1-value constructor
2911   explicit vector2d(int init) {x=y=init;}
2912   //2-value constructor
2913   vector2d(const Real Nx,const Real Ny) {x=Nx;y=Ny;}
2914   //Copy constructor
2915   vector2d(const vector2d &copy) {x=copy.x;y=copy.y;}
2916
2917   //Cast-to-Real * operators (treat vector as array)
2918   operator Real *() {return &x;}
2919   operator const Real *() const {return &x;}
2920
2921   /*Arithmetic operations: these are carefully restricted to just those
2922     that make unambiguous sense (to me... now...  ;-)
2923     Counterexamples: vector*vector makes no sense (use .dot()) because
2924     Real/vector is meaningless (and we'd want a*b/b==a for b!=0),
2925     ditto for vector&vector (dot?), vector|vector (projection?),
2926     vector^vector (cross?),Real+vector, vector+=real, etc.
2927   */
2928   vector2d &operator=(const vector2d &b) {x=b.x;y=b.y;return *this;}
2929   int operator==(const vector2d &b) const {return (x==b.x)&&(y==b.y);}
2930   int operator!=(const vector2d &b) const {return (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 vector2d &b) const {return vector2d(x-b.x,y-b.y);}
2933   vector2d operator*(const Real scale) const
2934     {return vector2d(x*scale,y*scale);}
2935   friend vector2d operator*(const Real scale,const vector2d &v)
2936     {return vector2d(v.x*scale,v.y*scale);}
2937   vector2d operator/(const Real &div) const
2938     {Real scale=1.0/div;return vector2d(x*scale,y*scale);}
2939   vector2d operator-(void) const {return vector2d(-x,-y);}
2940   void operator+=(const vector2d &b) {x+=b.x;y+=b.y;}
2941   void operator-=(const vector2d &b) {x-=b.x;y-=b.y;}
2942   void operator*=(const Real scale) {x*=scale;y*=scale;}
2943   void operator/=(const Real div) {Real scale=1.0/div;x*=scale;y*=scale;}
2944
2945   //Vector-specific operations
2946   //Return the square of the magnitude of this vector
2947   Real magSqr(void) const {return x*x+y*y;}
2948   //Return the magnitude (length) of this vector
2949   Real mag(void) const {return sqrt(magSqr());}
2950
2951   //Return the square of the distance to the vector b
2952   Real distSqr(const vector2d &b) const
2953     {return (x-b.x)*(x-b.x)+(y-b.y)*(y-b.y);}
2954   //Return the distance to the vector b
2955   Real dist(const vector2d &b) const {return sqrt(distSqr(b));}
2956
2957   //Return the dot product of this vector and b
2958   Real dot(const vector2d &b) const {return x*b.x+y*b.y;}
2959   //Return the cosine of the angle between this vector and b
2960   Real cosAng(const vector2d &b) const {return dot(b)/(mag()*b.mag());}
2961
2962   //Return the "direction" (unit vector) of this vector
2963   vector2d dir(void) const {return (*this)/mag();}
2964
2965   //Return the CCW perpendicular vector
2966   vector2d perp(void) const {return vector2d(-y,x);}
2967
2968   //Return this vector scaled by that
2969   vector2d &scale(const vector2d &b) {x*=b.x;y*=b.y;return *this;}
2970
2971   //Return the largest coordinate in this vector
2972   Real max(void) {return (x>y)?x:y;}
2973   //Make each of this vector's coordinates at least as big
2974   // as the given vector's coordinates.
2975   void enlarge(const vector2d &by)
2976     {if (by.x>x) x=by.x; if (by.y>y) y=by.y;}
2977 };
2978
2979 #endif //__OSL_VECTOR2D_H
2980
2981 /*
2982  * The former cktimer.h
2983  *
2984  */
2985
2986 #ifndef CMK_THRESHOLD_TIMER
2987 #define CMK_THRESHOLD_TIMER
2988
2989 /** Time a sequence of operations, printing out the
2990     names and times of any operations that exceed a threshold.
2991
2992     Use it with only the constructor and destructor like:
2993     void foo(void) {
2994     CkThresholdTimer t("foo");
2995     ...
2996     }
2997     (this times the whole execution of the routine,
2998     all the way to t's destructor is called on function return)
2999
3000     Or, you can start different sections like:
3001     void bar(void) {
3002     CkThresholdTimer t("first");
3003     ...
3004     t.start("second");
3005     ...
3006     t.start("third");
3007     ...
3008     }
3009
3010     This class *only* prints out the time if it exceeds
3011     a threshold-- by default, one millisecond.
3012 */
3013 class CkThresholdTimer {
3014   double threshold; // Print any times that exceed this (s).
3015   double lastStart; // Last activity started at this time (s).
3016   const char *lastWhat; // Last activity has this name.
3017
3018   void start_(const char *what) {
3019     lastStart=CmiWallTimer();
3020     lastWhat=what;
3021   }
3022   void done_(void) {
3023     double elapsed=CmiWallTimer()-lastStart;
3024     if (elapsed>threshold) {
3025       CmiPrintf("%s took %.2f s\n",lastWhat,elapsed);
3026     }
3027   }
3028  public:
3029   CkThresholdTimer(const char *what,double thresh=0.001)
3030     :threshold(thresh) { start_(what); }
3031   void start(const char *what) { done_(); start_(what); }
3032   ~CkThresholdTimer() {done_();}
3033 };
3034
3035 #endif
3036 // end cktimer or ParFUM_timer
3037
3038 #include "adapt_adj.h"
3039
3040 #endif
3041 // end ParFUM_internals.h
3042
3043 /*@}*/