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