Changed table to this->table to prevent compile errors on solaris
[charm.git] / src / libs / ck-libs / fem / fem_mesh.h
1 /*
2 Charm++ Finite Element Framework:
3 C++ implementation file: Mesh representation and manipulation
4
5 This file lists the classes used to represent and manipulate 
6 Finite Element Meshes inside the FEM framework.
7
8 Orion Sky Lawlor, olawlor@acm.org, 1/3/2003
9 */
10 #ifndef __CHARM_FEM_MESH_H
11 #define __CHARM_FEM_MESH_H
12
13 #include "charm-api.h"
14
15 #ifndef __CHARM_IDXL_COMM_H
16 #  include "idxl_comm.h" /*for IDXL_Side */
17 #endif
18 #ifndef __CHARM_IDXL_LAYOUT_H
19 #  include "idxl_layout.h" /*for IDXL_Layout */
20 #endif
21
22
23 // Map the IDXL names to the old FEM names (FIXME: change all references, too)
24 typedef IDXL_Side FEM_Comm;
25 typedef IDXL_List FEM_Comm_List;
26 typedef IDXL_Rec FEM_Comm_Rec;
27 class IDXL_Chunk;
28
29 /// We want the FEM_Comm/IDXL_Side's to be accessible to *both* 
30 ///  FEM routines (via these data structures) and IDXL routines
31 ///  (via an idxl->addStatic registration).  Hence this class, which
32 ///  manages IDXL's view of FEM's data structures.
33 class FEM_Comm_Holder {
34         IDXL comm; //Our communication lists.
35         bool registered; //We're registered with IDXL
36         IDXL_Comm_t idx; //Our IDXL index, or -1 if we're unregistered
37         void registerIdx(IDXL_Chunk *c);
38 public:
39         FEM_Comm_Holder(FEM_Comm *sendComm, FEM_Comm *recvComm);
40         void pup(PUP::er &p);
41         ~FEM_Comm_Holder(void);
42         
43         /// Return our IDXL_Comm_t, registering with chunk if needed
44         inline IDXL_Comm_t getIndex(IDXL_Chunk *c) {
45                 if (idx==-1) registerIdx(c);
46                 return idx;
47         }
48 };
49
50
51 /** \addtogroup fem_mesh FEM Framework Library Mesh Representation */
52 /*\@{*/
53
54
55 /// This datatype is how the framework stores symmetries internally.
56 ///   Each bit of this type describes a different symmetry.
57 ///   There must be enough bits to accomidate several simulatanious symmetries.
58 typedef unsigned char FEM_Symmetries_t;
59
60
61 /// Describes one kind of symmetry condition
62 class FEM_Sym_Desc : public PUP::able {
63 public:
64         virtual ~FEM_Sym_Desc();
65
66         /// Apply this symmetry to this location vector
67         virtual CkVector3d applyLoc(const CkVector3d &loc) const =0;
68         
69         /// Apply this symmetry to this relative (vel or acc) vector
70         virtual CkVector3d applyVec(const CkVector3d &vec) const =0;
71         
72         /// Allows Desc's to be pup'd via | operator:
73         friend inline void operator|(PUP::er &p,FEM_Sym_Desc &a) {a.pup(p);}
74         friend inline void operator|(PUP::er &p,FEM_Sym_Desc* &a) {
75                 PUP::able *pa=a;  p(&pa);  a=(FEM_Sym_Desc *)pa;
76         }
77 };
78
79 /// Describes a linear-periodic (space shift) symmetry:
80 class FEM_Sym_Linear : public FEM_Sym_Desc {
81         CkVector3d shift; //Offset to add to locations
82 public:
83         FEM_Sym_Linear(const CkVector3d &shift_) :shift(shift_) {}
84         FEM_Sym_Linear(CkMigrateMessage *m) {}
85         
86         /// Apply this symmetry to this location vector
87         CkVector3d applyLoc(const CkVector3d &loc) const {return loc+shift;}
88         
89         /// Apply this symmetry to this relative (vel or acc) vector
90         virtual CkVector3d applyVec(const CkVector3d &vec) const {return vec;}
91         
92         virtual void pup(PUP::er &p);
93         PUPable_decl(FEM_Sym_Linear);
94 };
95
96 /**
97  * Describes all the different kinds of symmetries that apply to
98  * this mesh.
99  */
100 class FEM_Sym_List {
101         //This lists the different kinds of symmetry
102         CkPupAblePtrVec<FEM_Sym_Desc> sym; 
103         
104         FEM_Sym_List(const FEM_Sym_List &src); //NOT DEFINED: copy constructor
105 public:
106         FEM_Sym_List();
107         void operator=(const FEM_Sym_List &src); //Assignment operator
108         ~FEM_Sym_List();
109         
110         /// Add a new kind of symmetry to this list, returning
111         ///  the way objects with that symmetry should be marked.
112         FEM_Symmetries_t add(FEM_Sym_Desc *desc);
113         
114         /// Apply all the listed symmetries to this location
115         void applyLoc(CkVector3d *loc,FEM_Symmetries_t sym) const;
116         
117         /// Apply all the listed symmetries to this relative vector
118         void applyVec(CkVector3d *vec,FEM_Symmetries_t sym) const;
119         
120         void pup(PUP::er &p);
121 };
122
123
124 /**
125  * This is a simple 2D table.  The operations are mostly row-centric.
126  */
127 template <class T>
128 class BasicTable2d : public CkNoncopyable {
129 protected:
130         int rows; //Number of entries in table
131         int cols; //Size of each entry in table
132         T *table; //Data in table [rows * cols]
133 public:
134         BasicTable2d(T *src,int cols_,int rows_) 
135                 :rows(rows_), cols(cols_), table(src) {}
136         
137         /// "size" of the table is the number of rows:
138         inline int size(void) const {return rows;}
139         /// Width of the table is the number of columns:
140         inline int width(void) const {return cols;}
141         
142         T *getData(void) {return table;}
143         const T *getData(void) const {return table;}
144         
145 /// Element-by-element operations:
146         T operator() (int r,int c) const {return table[c+r*cols];}
147         T &operator() (int r,int c) {return table[c+r*cols];}
148         
149 /// Row-by-row operations
150         //Get a pointer to a row of the table:
151         inline T *getRow(int r) {return &table[r*cols];}
152         inline const T *getRow(int r) const {return &table[r*cols];}
153         inline void getRow(int r,T *dest,T idxBase=0) const {
154                 const T *src=getRow(r);
155                 for (int c=0;c<cols;c++) dest[c]=src[c]+idxBase;
156         }
157         inline T *operator[](int r) {return getRow(r);}
158         inline const T *operator[](int r) const {return getRow(r);}
159         inline void setRow(int r,const T *src,T idxBase=0) {
160                 T *dest=getRow(r);
161                 for (int c=0;c<cols;c++) dest[c]=src[c]-idxBase;
162         }
163         inline void setRow(int r,T value) {
164                 T *dest=getRow(r);
165                 for (int c=0;c<cols;c++) dest[c]=value;
166         }
167         
168 /// These affect the entire table:
169         void set(const T *src,T idxBase=0) {
170                 for (int r=0;r<rows;r++) 
171                 for (int c=0;c<cols;c++)
172                         table[c+r*cols]=src[c+r*cols]-idxBase;
173         }
174         void setTranspose(const T *srcT,int idxBase=0) {
175                 for (int r=0;r<rows;r++) 
176                 for (int c=0;c<cols;c++)
177                         table[c+r*cols]=srcT[r+c*rows]-idxBase;
178         }
179         void get(T *dest,T idxBase=0) const {
180                 for (int r=0;r<rows;r++) 
181                 for (int c=0;c<cols;c++)
182                         dest[c+r*cols]=table[c+r*cols]+idxBase;
183         }
184         void getTranspose(T *destT,int idxBase=0) const {
185                 for (int r=0;r<rows;r++) 
186                 for (int c=0;c<cols;c++)
187                         destT[r+c*rows]=table[c+r*cols]+idxBase;
188         }
189         void set(T value) {
190                 for (int r=0;r<rows;r++) setRow(r,value);
191         }
192 };
193
194 /**
195  * A heap-allocatable, resizable BasicTable2d.
196  * To be stored here, T must not require a copy constructor.
197  */
198 template <class T>
199 class AllocTable2d : public BasicTable2d<T> {
200         int max; //Maximum number of rows that can be used without reallocation
201         T fill; //Value to fill uninitialized regions with
202         T *allocTable; // the table that I allocated
203 public:
204         AllocTable2d(int cols_=0,int rows_=0,T fill_=0) 
205                 :BasicTable2d<T>(NULL,cols_,rows_), max(0), fill(fill_),allocTable(NULL)
206         {
207                 if (this->rows>0) allocate(this->rows);
208         }
209         ~AllocTable2d() {if(allocTable != NULL){delete[] allocTable;}}
210         /// Make room for this many rows
211         void allocate(int rows_) { 
212                 allocate(this->width(),rows_);
213         }
214         /// Make room for this many cols & rows
215         void allocate(int cols_,int rows_,int max_=0) { 
216                 if (cols_==this->cols && rows_<max) {
217                         //We have room--just update the size:
218                         this->rows=rows_;
219                         return;
220                 }
221                 if (max_==0) { //They gave no suggested size-- pick one:
222                         if (rows_==this->rows+1) //Growing slowly: grab a little extra
223                                 max_=10+rows_+(rows_>>2); //  125% plus 10
224                         else // for a big change, just go with the minimum needed: 
225                                 max_=rows_;
226                 }
227
228                 int oldRows=this->rows;
229                 this->cols=cols_;
230                 this->rows=rows_;
231                 this->max=max_;
232                 this->table=new T[max*this->cols];
233                 //Preserve old table entries (FIXME: assumes old cols is unchanged)
234                 int copyRows=0;
235                 if (allocTable!=NULL) { 
236                         copyRows=oldRows;
237                         if (copyRows>max) copyRows=max;
238                         memcpy(this->table,allocTable,sizeof(T)*this->cols*copyRows);
239                         delete[] allocTable;
240                 }else{
241                         for (int r=copyRows;r<max;r++)
242                                 setRow(r,fill);
243                 }
244                 allocTable = this->table;
245         }
246         
247         /// Pup routine and operator|:
248         void pup(PUP::er &p) {
249                 p|this->rows; p|this->cols;
250                 if (this->table==NULL) allocate(this->rows);
251                 p(this->table,this->rows*this->cols); //T better be a basic type, or this won't compile!
252         }
253
254         void pupSingle(PUP::er &p, int pupindx) {
255           p|this->table[pupindx];
256         }
257
258         friend void operator|(PUP::er &p,AllocTable2d<T> &t) {t.pup(p);}
259
260         /// Add a row to the table (by analogy with std::vector):
261         T *push_back(void) {
262                 if (this->rows>=max) 
263                 { //Not already enough room for the new row:
264                         int newMax=max+(max/4)+16; //Grow 25% longer
265                         allocate(this->cols,this->rows,newMax);
266                 }
267                 this->rows++;
268                 return getRow(this->rows-1);
269         }
270
271         /** to support replacement of attribute data by user 
272                         supplied data
273                         error checks have been performed at FEM_ATTRIB
274         */
275         void register_data(T *user,int len,int max){
276                 if(allocTable != NULL){
277                         delete [] allocTable;
278                         allocTable = NULL;
279                 }       
280                 this->table = user;
281                 this->rows = len;
282                 max = max;
283         }
284 };
285
286
287 class FEM_Entity; // Forward declaration
288 class FEM_Mesh; 
289
290
291 /// Return the human-readable version of this entity code, like "FEM_NODE".
292 ///  storage, which must be at least 80 bytes long, is used for
293 ///  non-static names, like the user tag "FEM_ELEM+2".
294 CDECL const char *FEM_Get_entity_name(int entity,char *storage);
295
296 /// Return the human-readable version of this attribute code, like "FEM_CONN".
297 ///  storage, which must be at least 80 bytes long, is used for
298 ///  non-static names, like the user tag "FEM_DATA+7".
299 CDECL const char *FEM_Get_attr_name(int attr,char *storage);
300
301
302
303 /**
304  * Describes an FEM entity's "attribute"--a user-visible, user-settable
305  * 2D table.  Common FEM_Attributes include: user data associated with nodes,
306  *  the element-to-node connectivity array, etc.
307  */
308 class FEM_Attribute {
309         FEM_Entity *e; //Owning entity (to get length, etc.)
310         FEM_Attribute *ghost; // Ghost attribute, which has the same width and datatype as us (or 0)
311         int attr; // My attribute code (e.g., FEM_DATA+7, FEM_CONN, etc.)
312         
313         int width; //Number of columns in our table of data (initially unknown)
314         int datatype; //Datatype of entries (initially unknown)
315         bool allocated; //True if subclass allocate has been called.
316         
317         //Abort with a nice error message saying: 
318         // Our <field> was previously set to <cur>; it cannot now be <next>
319         void bad(const char *field,bool forRead,int cur,int next,const char *caller) const;
320         
321 protected:
322         /**
323          * Allocate storage for at least length width-item records of type datatype.
324          * This routine is called after all three parameters are set,
325          * as a convenience for subclasses. 
326          */
327         virtual void allocate(int length,int width,int datatype) =0;
328 public:
329         FEM_Attribute(FEM_Entity *owner_,int myAttr_);
330         virtual void pup(PUP::er &p);
331         virtual void pupSingle(PUP::er &p, int pupindx);
332         virtual ~FEM_Attribute();
333         
334         /// Install this attribute as our ghost:
335         inline void setGhost(FEM_Attribute *ghost_) { ghost=ghost_;}
336         
337         /// Return true if we're a ghost
338         inline bool isGhost(void) const { return ghost!=NULL; }
339         
340         /// Return our attribute code
341         inline int getAttr(void) const {return attr;}
342
343         inline FEM_Entity *getEntity(void) {return e;}
344         
345         /// Return the number of rows in our table of data (0 if unknown).
346         /// This value is obtained directly from our owning Entity.
347         int getLength(void) const;
348         int getRealLength(void) const;
349
350         int getMax();
351         
352         /// Return the number of columns in our table of data (0 if unknown)
353         inline int getWidth(void) const { return width<0?0:width; }
354         inline int getRealWidth(void) const { return width; }
355         
356         /// Return our FEM_* datatype (-1 if unknown)
357         inline int getDatatype(void) const { return datatype; }
358         
359         /**
360          * Set our length (number of rows, or records) to this value.
361          * Default implementation calls setLength on our owning entity
362          * and tries to call allocate().
363          */
364         void setLength(int l,const char *caller="");
365         
366         /**
367          * Set our width (number of values per row) to this value.
368          * The default implementation sets width and tries to call allocate().
369          */
370         void setWidth(int w,const char *caller="");
371         
372         /**
373          * Set our datatype (e.g., FEM_INT, FEM_DOUBLE) to this value.
374          * The default implementation sets width and tries to call allocate().
375          */
376         void setDatatype(int dt,const char *caller="");
377         
378         /**
379          * Copy our width and datatype from this attribute.
380          * The default implementation calls setWidth and setDatatype.
381          * which should be enough for virtually any attribute.
382          */
383         virtual void copyShape(const FEM_Attribute &src);
384         
385         /// Check if all three of length, width, and datatype are set,
386         ///  but we're not yet allocated.  If so, call allocate; else ignore.
387         void tryAllocate(void);
388         
389         /// Our parent's length has changed: reallocate our storage
390         inline void reallocate(void) { allocated=false; tryAllocate(); }
391         
392         /// Return true if we've already had storage allocated.
393         inline bool isAllocated(void) const {return allocated;}
394         
395         /**
396          * Set our data to these (typically user-supplied, unchecked) values.
397          * Subclasses normally override this method as:
398          *    virtual void set( ...) {
399          *       super::set(...);
400          *       copy data from src.
401          *    }
402          */
403         virtual void set(const void *src, int firstItem,int length, 
404                 const IDXL_Layout &layout, const char *caller);
405         
406         /**
407          * Extract this quantity of user data.  Length and layout are
408          * parameter checked by the default implementation.
409          * Subclasses normally override this method as:
410          *    virtual void get( ...) {
411          *       super::get(...);
412          *       copy data to dest.
413          *    }
414          */
415         virtual void get(void *dest, int firstItem,int length,
416                 const IDXL_Layout &layout, const char *caller) const;
417         
418         /// Copy everything associated with src[srcEntity] into our dstEntity.
419         virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity) =0;
420
421         /** Register this user data for this attributre 
422                         Length, layout etc are checked by the default implementaion
423         */
424
425         virtual void register_data(void *dest, int length,int max,
426                 const IDXL_Layout &layout, const char *caller);
427         
428 };
429 PUPmarshall(FEM_Attribute);
430
431
432 /**
433  * Describes a single table of user data associated with an entity.
434  * Since the data can be of any type, it is stored as chars.
435  */
436 class FEM_DataAttribute : public FEM_Attribute {
437         typedef FEM_Attribute super;
438         AllocTable2d<unsigned char> *char_data; //Non-NULL for getDatatype==FEM_BYTE
439         AllocTable2d<int> *int_data; //Non-NULL for getDatatype==FEM_INT
440         AllocTable2d<float> *float_data; //Non-NULL for getDatatype==FEM_FLOAT
441         AllocTable2d<double> *double_data; //Non-NULL for getDatatype==FEM_DOUBLE
442 protected:
443         virtual void allocate(int length,int width,int datatype);
444 public:
445         FEM_DataAttribute(FEM_Entity *owner,int myAttr);
446         virtual void pup(PUP::er &p);
447         virtual void pupSingle(PUP::er &p, int pupindx);
448         ~FEM_DataAttribute();
449         
450         AllocTable2d<unsigned char> &getChar(void) {return *char_data;}
451         const AllocTable2d<unsigned char> &getChar(void) const {return *char_data;}
452         
453         AllocTable2d<int> &getInt(void) {return *int_data;}
454         const AllocTable2d<int> &getInt(void) const {return *int_data;}
455         
456         AllocTable2d<double> &getDouble(void) {return *double_data;}
457
458         
459         virtual void set(const void *src, int firstItem,int length, 
460                 const IDXL_Layout &layout, const char *caller);
461         
462         virtual void get(void *dest, int firstItem,int length, 
463                 const IDXL_Layout &layout, const char *caller) const;
464         
465         virtual void register_data(void *dest, int length,int max,
466                 const IDXL_Layout &layout, const char *caller);
467
468         
469         /// Copy src[srcEntity] into our dstEntity.
470         virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
471
472         /*used during refining to extrapolate the values of a node */
473         void interpolate(int A,int B,int D,double frac);
474         void interpolate(int *iNodes,int rNode,int k);
475 };
476 PUPmarshall(FEM_DataAttribute);
477
478 /**
479  * This table maps an entity to a set of integer indices.
480  * The canonical example of this is the element-node connectivity array.
481  */
482 class FEM_IndexAttribute : public FEM_Attribute {
483 public:
484         /// Checks incoming indices for validity.
485         class Checker {
486         public:
487                 virtual ~Checker();
488                 
489                 /**
490                  * Check this (newly set) row of our table for validity.
491                  * You're expected to abort or throw or exit if something is wrong.
492                  */
493                 virtual void check(int row,const BasicTable2d<int> &table,
494                         const char *caller) const =0;
495         };
496 private:
497         typedef FEM_Attribute super;
498         AllocTable2d<int> idx;
499         Checker *checker; //Checks indices (or NULL). This attribute will destroy this object
500 protected:
501         virtual void allocate(int length,int width,int datatype);
502 public:
503         FEM_IndexAttribute(FEM_Entity *owner,int myAttr, Checker *checker_=NULL);
504         virtual void pup(PUP::er &p);
505         virtual void pupSingle(PUP::er &p, int pupindx);
506         ~FEM_IndexAttribute();
507         
508         AllocTable2d<int> &get(void) {return idx;}
509         const AllocTable2d<int> &get(void) const {return idx;}
510         
511         virtual void set(const void *src, int firstItem,int length, 
512                 const IDXL_Layout &layout, const char *caller);
513         
514         virtual void get(void *dest, int firstItem,int length,
515                 const IDXL_Layout &layout, const char *caller) const;
516
517         virtual void register_data(void *dest, int length, int max,
518                 const IDXL_Layout &layout, const char *caller);
519         
520         /// Copy src[srcEntity] into our dstEntity.
521         virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
522 };
523 PUPmarshall(FEM_IndexAttribute);
524
525 /*
526         This table maps an entity to a list of integer indices 
527         of unknown length.
528         The node to element adjacency array is an example, where a node
529         is mapped to a list of element indices of unknown length.
530 */
531 class FEM_VarIndexAttribute : public FEM_Attribute{
532 public:
533  class ID{
534         //type is negative for ghost elements
535         //id refers to the index in the entity list
536         public:
537                 int type;
538                 int id;
539                 ID(){
540                         type=-1;
541                         id = -1;
542                 };
543                 ID(int _type,int _id){
544                   if(_id < 0) {
545                         type = -(_type+1);
546                         id = FEM_To_ghost_index(_id);
547                   }
548                   else {
549                         type = _type;
550                         id = _id;
551                   }
552                 };
553                 bool operator ==(const ID &rhs)const {
554                         return (type == rhs.type) && (id == rhs.id);
555                 }
556                 virtual void pup(PUP::er &p){
557                         p | type;
558                         p | id;
559                 };
560
561                 static ID createNodeID(int type,int node){
562                   ID temp(type, node);
563                   return temp;
564                 }
565                 int getSignedId() {
566                   if(type<0){
567                         return FEM_From_ghost_index(id);
568                   }
569                   else return id;
570                 }
571  };
572 private:
573         typedef FEM_Attribute super;
574         CkVec<CkVec<ID> > idx;
575         int oldlength;
576 protected:
577         virtual void allocate(int _length,int _width,int _datatype){
578           if(_length > oldlength){
579             setWidth(1,"allocate"); //there is 1 vector per entity 
580             oldlength = _length*2;
581             idx.reserve(oldlength);
582             for(int i=idx.size();i<oldlength;i++){
583               CkVec<ID> tempVec;
584               idx.insert(i,tempVec);
585             }
586           }
587         };
588 public:
589         FEM_VarIndexAttribute(FEM_Entity *owner,int myAttr);
590         ~FEM_VarIndexAttribute(){};
591         virtual void pup(PUP::er &p);
592         virtual void pupSingle(PUP::er &p, int pupindx);
593         CkVec<CkVec<ID> > &get(){return idx;};
594         const CkVec<CkVec<ID> > &get() const {return idx;};
595         
596         virtual void set(const void *src,int firstItem,int length,
597                 const IDXL_Layout &layout,const char *caller);
598
599         virtual void get(void *dest, int firstItem,int length,
600                 const IDXL_Layout &layout, const char *caller) const;
601         
602         virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
603
604         int findInRow(int row,const ID &data);
605         
606         void print();
607 };
608
609
610 class l2g_t;
611 /**
612  * Describes an entire class of "entities"--nodes, elements, or sparse 
613  *  data records. Basically consists of a length and a set of
614  *  FEM_Attributes. 
615  */
616 class FEM_Entity {
617         typedef CkVec<FEM_Symmetries_t> sym_t;
618         int length; // Number of entities of our type currently in the mesh
619         int max;    // Maximum number of entities of our type in the mesh that will be allowed
620         FEM_Mesh_alloc_fn resize; // should be a function pointer to the actual resize function later
621         void *args; // arguments to the resize function
622         
623         /**
624          * This is our main list of attributes-- everything about each of 
625          * our entities is in this list of attributes.  This list is searched
626          * by our "lookup" method and maintained by our subclasses "create" 
627          * method and calls to our "add" method.
628          * 
629          * It's a little funky having the superclass keep pointers to subclass
630          * objects (like element connectivity), but very nice to be able to
631          * easily loop over everything associated with an entity.
632          */
633         CkVec<FEM_Attribute *> attributes;
634         
635         /**
636          * Coordinates of each entity, from FEM_COORD.
637          * Datatype is always FEM_DOUBLE, width is always 2 or 3.
638          *  If NULL, coordinates are unknown.
639          */
640         FEM_DataAttribute *coord; 
641         void allocateCoord(void);
642         
643         /**
644          * Symmetries of each entity, from FEM_SYMMETRIES.  This bitvector per
645          * entity indicates which symmetry conditions the entity belongs to.
646          * Datatype is always FEM_BYTE (same as FEM_Symmetries_t), width 
647          * is always 1.  If NULL, all the symmetries are 0.
648          */
649         FEM_DataAttribute *sym; 
650         void allocateSym(void);
651         
652         /**
653          * Global numbers of each entity, from FEM_GLOBALNO.
654          * If NULL, the global numbers are unknown.
655          */
656         FEM_IndexAttribute *globalno;
657         void allocateGlobalno(void);
658         
659         /*
660                 used to allocate the integer array for storing the boundary
661                 values associated with an entity. 
662         */
663         void allocateBoundary();
664
665         /// Mesh sizing attribute for elements
666         /** Specifies a double edge length target for the mesh at each 
667             element; used in adaptivity algorithms */
668         FEM_DataAttribute *meshSizing;
669         void allocateMeshSizing(void);
670
671         /*
672             used to allocate the char array for storing whether each entity is valid
673                 When a node/element is deleted the flag in the valid table is set to 0.
674                 Additionally, we keep track of the first and last occurence in the array of
675                 invalid indices. This should make searching for slots to reuse quicker.
676         */
677         FEM_DataAttribute *valid;
678         unsigned int first_invalid, last_invalid;
679         
680 protected:
681         /**
682          * lookup of this attribute code has failed: check if it needs
683          * to be demand-created.  Subclasses should override this method
684          * to recognize a request for, and add their own attributes;
685          * otherwise call the default implementation.
686          *
687          * Every call to create must result in either a call to add()
688          * or a call to the superclass; but not both.
689          *
690          * Any entity with optional fields, that are created on demand, will
691          * have to override this method.  Entities with fixed fields that are
692          * known beforehand should just call add() from their constructor.
693          */
694         virtual void create(int attr,const char *caller);
695         
696         /// Add this attribute to this kind of Entity.
697         /// This superclass is responsible for eventually deleting the attribute.
698         /// This class also attaches the ghost attribute, so be sure to 
699         ///   call add before manipulating the attribute.
700         void add(FEM_Attribute *attribute);
701  public:
702
703         FEM_Entity *ghost; // Our ghost entity type, or NULL if we're the ghost
704                 
705         FEM_Comm ghostSend; //Non-ghosts we send out (only set for real entities)
706         FEM_Comm ghostRecv; //Ghosts we recv into (only set for ghost entities)
707
708         FEM_Entity(FEM_Entity *ghost_); //Default constructor
709         void pup(PUP::er &p);
710         virtual ~FEM_Entity();
711         
712         /// Return true if we're a ghost
713         bool isGhost(void) const {return ghost==NULL;}
714         
715         /// Switch from this, a real entity, to the ghosts:
716         FEM_Entity *getGhost(void) {return ghost;}
717         const FEM_Entity *getGhost(void) const {return ghost;}
718         
719         /// Return the number of entities of this type
720         inline int size(void) const {return length==-1?0:length;}
721         inline int realsize(void) const {return length;}
722
723         // return the maximum size 
724         inline int getMax() { if(max > 0) return max; else return length;}
725         
726         /// Return the human-readable name of this entity type, like "node"
727         virtual const char *getName(void) const =0;
728         
729         /// Copy all our attributes' widths and data types from this entity.
730         void copyShape(const FEM_Entity &src);
731         
732         /** 
733          *  The user is setting this many entities.  This reallocates
734          * all existing attributes to make room for the new entities.
735          */
736         void setLength(int newlen);
737
738         /** Support for registration API 
739          *  Set the current length and maximum length for this entity. 
740          *  If the current length exceeds the maximum length a resize
741          *  method is called .
742         */
743         void setMaxLength(int newLen,int newMaxLen,void *args,FEM_Mesh_alloc_fn fn);
744         
745         /// Copy everything associated with src[srcEntity] into our dstEntity.
746         /// dstEntity must have already been allocated, e.g., with setLength.
747         void copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity);
748
749         /// Add room for one more entity, with initial values from src[srcEntity],
750         /// and return the new entity's index.
751         int push_back(const FEM_Entity &src,int srcEntity);
752         
753         /**
754          * Find this attribute (from an FEM_ATTR code) of this entity, or 
755          * create the entity (using the create method below) or abort if it's
756          * not found.
757          */
758         FEM_Attribute *lookup(int attr,const char *caller);
759
760         
761         /**
762          * Get a list of the attribute numbers for this entity.
763          */
764         int getAttrs(int *attrs) const {
765                 int len=attributes.size();
766                 for (int i=0;i<len;i++) 
767                         attrs[i]=attributes[i]->getAttr();
768                 return len;
769         }
770         
771         /**
772          * Allocate or Modify the FEM_IS_VALID attribute data
773          */
774         void allocateValid();
775         void set_valid(unsigned int idx, bool isNode);
776         void set_invalid(unsigned int idx, bool isNode);
777         int is_valid(unsigned int idx);
778         unsigned int count_valid();
779         unsigned int get_next_invalid(FEM_Mesh *m, bool isNode, bool isGhost);
780         virtual bool hasConn(unsigned int idx)=0;
781         /**
782          * Set the coordinates for a single item
783          */
784         void set_coord(int idx, double x, double y);
785         void set_coord(int idx, double x, double y, double z);
786
787         /**expose the attribute vector for refining 
788                 . breaks modularity but more efficient
789         */
790         CkVec<FEM_Attribute *>* getAttrVec(){
791                 return &attributes;
792         }
793         
794         // Stupidest possible coordinate access
795         inline FEM_DataAttribute *getCoord(void) {return coord;}
796         inline const FEM_DataAttribute *getCoord(void) const {return coord;}
797         
798         //Symmetry array access:
799         const FEM_Symmetries_t *getSymmetries(void) const {
800                 if (sym==NULL) return NULL;
801                 else return (const FEM_Symmetries_t *)sym->getChar()[0];
802         }
803         FEM_Symmetries_t getSymmetries(int r) const { 
804                 if (sym==NULL) return FEM_Symmetries_t(0);
805                 else return sym->getChar()(r,0);
806         }
807         void setSymmetries(int r,FEM_Symmetries_t s);
808         
809         //Global numbering array access
810         bool hasGlobalno(void) const {return globalno!=0;}
811         int getGlobalno(int r) const {
812                 if (globalno==0) return -1; // Unknown global number
813                 return globalno->get()(r,0);
814         }
815         void setGlobalno(int r,int g);
816         void setAscendingGlobalno(void);
817         void setAscendingGlobalno(int base);
818         void copyOldGlobalno(const FEM_Entity &e);
819
820         // Mesh sizing array access
821         bool hasMeshSizing(void) const {return meshSizing!=0;}
822         double getMeshSizing(int r); 
823         void setMeshSizing(int r,double s);
824         void setMeshSizing(double *sf);
825         
826         //Ghost comm. list access
827         FEM_Comm &setGhostSend(void) { return ghostSend; }
828         const FEM_Comm &getGhostSend(void) const { return ghostSend; }
829         FEM_Comm &setGhostRecv(void) { 
830                 if (ghost==NULL) return ghostRecv;
831                 else return ghost->ghostRecv; 
832         }
833         const FEM_Comm &getGhostRecv(void) const { return ghost->ghostRecv; }
834         FEM_Comm_Holder ghostIDXL; //IDXL interface
835         
836         void addVarIndexAttribute(int code){
837                 FEM_VarIndexAttribute *varAttribute = new FEM_VarIndexAttribute(this,code);
838                 add(varAttribute);
839         }
840         
841         void print(const char *type,const IDXL_Print_Map &map);
842 };
843 PUPmarshall(FEM_Entity);
844
845 // Now that we have FEM_Entity, we can define attribute lenth, as entity length
846 inline int FEM_Attribute::getLength(void) const { return e->size(); }
847 inline int FEM_Attribute::getRealLength(void) const { return e->realsize(); }
848 inline int FEM_Attribute::getMax(){ return e->getMax();}
849
850 /**
851  * Describes a set of FEM Nodes--the FEM_NODE entity type. 
852  * Unlike elements, nodes have no connectivity; but they do have
853  * special shared-nodes communications and a "primary" flag.
854  */
855
856 class FEM_Elem;
857 class FEM_Node : public FEM_Entity {
858         typedef FEM_Entity super;
859
860         /**
861          * primary flag, from FEM_NODE_PRIMARY, indicates this chunk "owns" this node.
862          * Datatype is always FEM_BYTE (we need an FEM_BIT!), width is always 1,
863          * since there's only one such flag per node.
864          */
865         FEM_DataAttribute *primary; 
866         void allocatePrimary(void);
867         
868         void allocateElemAdjacency();
869         void allocateNodeAdjacency();
870
871         FEM_VarIndexAttribute *elemAdjacency; // stores the node to element adjacency vector 
872         FEM_VarIndexAttribute *nodeAdjacency; // stores the node to node adjacency vector 
873         typedef FEM_VarIndexAttribute::ID var_id;
874 protected:
875         virtual void create(int attr,const char *caller);
876 public:
877         FEM_Comm shared; //Shared nodes
878         FEM_Comm_Holder sharedIDXL; //IDXL interface to shared nodes
879         
880         FEM_Node(FEM_Node *ghost_);
881         void pup(PUP::er &p);
882         ~FEM_Node();
883         
884         virtual const char *getName(void) const;
885         
886         inline bool getPrimary(int nodeNo) const {
887                 if (primary==NULL) return true; //Everything must be primary
888                 else return primary->getChar()(nodeNo,0);
889         }
890         inline void setPrimary(int nodeNo,bool isPrimary) {
891                 if (primary==NULL) allocatePrimary();
892                 primary->getChar()(nodeNo,0)=isPrimary;
893         }
894         void fillElemAdjacencyTable(int type,const FEM_Elem &elem);
895         void setElemAdjacency(int type,const FEM_Elem &elem);
896         void fillNodeAdjacency(const FEM_Elem &elem);
897         void setNodeAdjacency(const FEM_Elem &elem);
898         void fillNodeAdjacencyForElement(int node,int nodesPerElem,const int *conn,FEM_VarIndexAttribute *adjacencyAttr);
899         bool hasConn(unsigned int idx);
900         void print(const char *type,const IDXL_Print_Map &map);
901 };
902 PUPmarshall(FEM_Node);
903
904 /**
905  * Describes one kind of FEM elements--the FEM_ELEM entity type.
906  * Elements are typically the central user-visible object in a FEM
907  * computation.
908  */
909 class FEM_Elem:public FEM_Entity {
910         typedef FEM_Entity super;
911 protected:
912         typedef AllocTable2d<int> conn_t;
913
914         int tuplesPerElem;
915
916         // The following are attributes that will commonly be used:
917         FEM_IndexAttribute *conn;                // FEM_CONN attribute: element-to-node mapping 
918         FEM_IndexAttribute *elemAdjacency;       // FEM_ELEM_ELEM_ADJACENCY attribute
919         FEM_IndexAttribute *elemAdjacencyTypes;  // FEM_ELEM_ELEM_ADJ_TYPES attribute
920
921 public:
922
923         FEM_Elem(const FEM_Mesh &mesh_, FEM_Elem *ghost_);
924         void pup(PUP::er &p);
925         ~FEM_Elem();
926         
927         virtual const char *getName(void) const;
928         
929         /// Directly access our connectivity table:
930         inline conn_t &setConn(void) {return conn->get();}
931         inline const conn_t &getConn(void) const {return conn->get();}
932         
933         void print(const char *type,const IDXL_Print_Map &map);
934
935         void create(int attr,const char *caller);
936
937         void allocateElemAdjacency();
938
939         FEM_IndexAttribute *getElemAdjacency(){return elemAdjacency;}
940
941         // Backward compatability routines:
942         int getConn(int elem,int nodeNo) const {return conn->get()(elem,nodeNo);}
943         int getNodesPer(void) const {return conn->get().width();}
944         int *connFor(int i) {return conn->get().getRow(i);}
945         const int *connFor(int i) const {return conn->get().getRow(i);}
946         void connIs(int i,const int *src) {conn->get().setRow(i,src);}
947         bool hasConn(unsigned int idx);
948 };
949 PUPmarshall(FEM_Elem);
950
951
952
953 /**
954  * FEM_Sparse describes a set of records of sparse data that are all the
955  * same size and all associated with the same number of nodes.
956  * Sparse data is associated with some subset of the nodes in the mesh,
957  * and gets copied to every chunk that has all those nodes.  The canonical
958  * use of sparse data is to describe boundary conditions.
959  */
960 class FEM_Sparse : public FEM_Elem {
961         typedef FEM_Elem super;
962         typedef AllocTable2d<int> elem_t;
963         
964         /**
965          * elem, from FEM_SPARSE_ELEM, is an optional (that is, possibly NULL) 
966          * array which changes the partitioning of sparse entities: if non-NULL,
967          * sparse entity t lives on the same chunk as FEM_ELEM+elem[2*t]
968          * local number elem[2*t+1].
969          *
970          * This attribute's width is always 2.
971          */
972         FEM_IndexAttribute *elem; //FEM_SPARSE_ELEM attribute
973         void allocateElem(void);
974         const FEM_Mesh &mesh; //Reference to enclosing mesh, for error checking
975 protected:
976         virtual void create(int attr,const char *caller);
977 public:
978         FEM_Sparse(const FEM_Mesh &mesh_, FEM_Sparse *ghost_);
979         void pup(PUP::er &p);
980         virtual ~FEM_Sparse();
981         
982         virtual const char *getName(void) const;
983         
984         /// Return true if we have an element partitioning table
985         bool hasElements(void) const {return elem!=NULL;}
986         
987         /// Directly access our element partitioning table (e.g., for re-numbering)
988         inline elem_t &setElem(void) {return elem->get();}
989         inline const elem_t &getElem(void) const {return elem->get();}
990 };
991 PUPmarshall(FEM_Sparse);
992
993 /** Describes a user function to pup a piece of mesh data 
994 */
995 class FEM_Userdata_pupfn {
996         FEM_Userdata_fn fn;
997         void *data;
998 public:
999         FEM_Userdata_pupfn(FEM_Userdata_fn fn_,void *data_)
1000                 :fn(fn_), data(data_) {}
1001         /// Call user's pup routine using this PUP::er
1002         void pup(PUP::er &p) {
1003                 (fn)((pup_er)&p,data);
1004         }
1005 };
1006
1007 /** Describes one piece of generic unassociated mesh data.
1008 */
1009 class FEM_Userdata_item {
1010         CkVec<char> data; ///< Serialized data from user's pup routine
1011 public:
1012         int tag; //User-assigned identifier
1013         FEM_Userdata_item(int tag_=-1) {tag=tag_;}
1014         
1015         /// Return true if we have stored data.
1016         bool hasStored(void) const {return data.size()!=0;}
1017         
1018         /// Store this userdata inside us:
1019         void store(FEM_Userdata_pupfn &f) {
1020                 data.resize(PUP::size(f));
1021                 PUP::toMemBuf(f,&data[0],data.size());
1022         }
1023         /// Extract this userdata from our stored copy:
1024         void restore(FEM_Userdata_pupfn &f) {
1025                 PUP::fromMemBuf(f,&data[0],data.size());
1026         }
1027         
1028         /// Save our stored data to this PUP::er
1029         void pup(PUP::er &p) {
1030                 p|tag;
1031                 p|data;
1032         }
1033 };
1034
1035 /** Describes all the unassociated data in a mesh. 
1036 */
1037 class FEM_Userdata_list {
1038         CkVec<FEM_Userdata_item> list;
1039 public:
1040         FEM_Userdata_item &find(int tag) {
1041                 for (int i=0;i<list.size();i++)
1042                         if (list[i].tag==tag)
1043                                 return list[i];
1044                 // It's not in the list-- add it.
1045                 list.push_back(FEM_Userdata_item(tag));
1046                 return list[list.size()-1];
1047         }
1048         int size(){
1049                 return list.size();
1050         }
1051         void pup(PUP::er &p) {p|list;}
1052 };
1053
1054
1055 void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType);
1056 void FEM_Is_NULL(const char *caller,const char *entityType,int type);
1057
1058 /**
1059  * This class describes several different types of a certain kind 
1060  * of entity.  For example, there might be a FEM_Entity_Types<FEM_Elem>
1061  * that lists the different kinds of element.
1062  *
1063  * This class exists to provide a nice "demand-creation" semantics,
1064  * where the user assigns array indices (the e in FEM_ELEM+e),
1065  * so we don't know that we're setting the first copy when we set it.
1066  *
1067  * It's not clear this class has any right to exist--it should either be
1068  * folded into FEM_Mesh or generalized into a "userNumberedVec" or some such.
1069  */
1070 template <class T>
1071 class FEM_Entity_Types {
1072         CkVec<T *> types; // Our main storage for different entity types
1073         const FEM_Mesh &mesh;
1074         const char *name; //FEM_SPARSE or FEM_ELEM, or some such.
1075
1076 public:
1077         FEM_Entity_Types(const FEM_Mesh &mesh_,const char *name_) 
1078                 :mesh(mesh_), name(name_) {}
1079         void pup(PUP::er &p) { 
1080                 // Can't just use p|types, because T has a funky constructor:
1081                 int n=types.size();
1082                 p|n;
1083                 for (int i=0;i<n;i++) {
1084                         int isNULL=0;
1085                         if (!p.isUnpacking()) isNULL=(types[i]==NULL);
1086                         p|isNULL;
1087                         if (!isNULL) set(i,"pup").pup(p);
1088                 }
1089         }
1090         ~FEM_Entity_Types() {
1091                 for (int i=0;i<types.size();i++)
1092                         if (types[i]) delete types[i];
1093         }
1094         
1095         /// Return the number of different entity types
1096         inline int size(void) const {return types.size();}
1097         
1098         /// Return a read-only copy of this type, or else abort if type isn't set.
1099         const T &get(int type,const char *caller="") const {
1100                 FEM_Index_Check(caller,name,type,types.size());
1101                 const T *ret=types[type];
1102                 if (ret==NULL) FEM_Is_NULL(caller,name,type);
1103                 return *ret;
1104         }
1105         
1106         /// Return true if we have a type t, and false otherwise
1107         bool has(int type) const { 
1108                 if (type>=types.size()) return false;
1109                 return types[type]!=NULL; 
1110         }
1111         
1112         /// Return a writable copy of this type, calling new T(mesh) if it's not there
1113         T &set(int type,const char *caller="") {
1114                 if (type<0) FEM_Index_Check(caller,name,type,types.size());
1115                 while (types.size()<=type) types.push_back(NULL); //Make room for new type
1116                 if (types[type]==NULL) { //Have to allocate a new T:
1117                         T *ghost=new T(mesh,NULL);
1118                         types[type]=new T(mesh,ghost);
1119                 }
1120                 return *types[type];
1121         }
1122         
1123         /// Read-only and write-only operator[]'s:
1124         inline T &operator[] (int type) { return set(type); }
1125         inline const T &operator[] (int type) const { return get(type); }
1126
1127 };
1128
1129 // Map fortran element (>=1) or node (0) marker to C version (>=1, -1)
1130 inline int zeroToMinusOne(int i) {
1131         if (i==0) return -1;
1132         return i;
1133 }
1134
1135
1136
1137
1138 /**
1139  * This class describes all the nodes and elements in
1140  * a finite-element mesh or submesh.
1141  */
1142 class FEM_ElemAdj_Layer;
1143 class femMeshModify;
1144 class FEM_Mesh : public CkNoncopyable {
1145   /// The symmetries in the mesh
1146   FEM_Sym_List symList;
1147   bool m_isSetting;
1148   femMeshModify *fmMM;
1149   
1150   void checkElemType(int elType,const char *caller) const;
1151   void checkSparseType(int uniqueID,const char *caller) const; 
1152   
1153   FEM_ElemAdj_Layer* lastElemAdjLayer;
1154
1155  public:
1156   void setFemMeshModify(femMeshModify *m);
1157   
1158   FEM_Mesh();
1159   void pup(PUP::er &p); //For migration
1160   ~FEM_Mesh();
1161   
1162   /// The nodes in this mesh:
1163   FEM_Node node; 
1164   
1165   /// The different element types in this mesh:
1166   FEM_Entity_Types<FEM_Elem> elem;
1167   
1168   /// The different sparse types in this mesh:
1169   FEM_Entity_Types<FEM_Sparse> sparse;
1170   
1171   /// The unassociated user data for this mesh:
1172   FEM_Userdata_list udata;
1173   
1174   /// The symmetries that apply to this mesh:
1175   void setSymList(const FEM_Sym_List &src) {symList=src;}
1176   const FEM_Sym_List &getSymList(void) const {return symList;}
1177   
1178   /// Set up the "shape" of our fields-- the number of element types,
1179   /// the datatypes for user data, etc--based on this mesh.
1180   void copyShape(const FEM_Mesh &src);
1181
1182   // Get the fem mesh modification object associated with this mesh or partition  
1183   femMeshModify *getfmMM();
1184
1185   //Return this type of element, given an element type
1186   FEM_Entity &setCount(int elTypeOrMinusOne) {
1187     if (elTypeOrMinusOne==-1) return node;
1188     else return elem[chkET(elTypeOrMinusOne)];
1189   }
1190   const FEM_Entity &getCount(int elTypeOrMinusOne) const {
1191     if (elTypeOrMinusOne==-1) return node;
1192     else return elem[chkET(elTypeOrMinusOne)];
1193   }
1194   FEM_Elem &setElem(int elType) {return elem[chkET(elType)];}
1195   const FEM_Elem &getElem(int elType) const {return elem[chkET(elType)];}
1196   int chkET(int elType) const; //Check this element type-- abort if it's bad
1197   
1198   /// Look up this FEM_Entity type in this mesh, or abort if it's not valid.
1199   FEM_Entity *lookup(int entity,const char *caller);
1200   const FEM_Entity *lookup(int entity,const char *caller) const;
1201   
1202   /// Set/get direction control:
1203   inline bool isSetting(void) const {return m_isSetting;}
1204   void becomeSetting(void) {m_isSetting=true;}
1205   void becomeGetting(void) {m_isSetting=false;}
1206   
1207   int nElems() const //Return total number of elements (of all types)
1208     {return nElems(elem.size());}
1209   /// Return the total number of elements before type t
1210   int nElems(int t) const;
1211   /// Return the "global" number of element elNo of type elType.
1212   int getGlobalElem(int elType,int elNo) const;
1213   /// Set our global numbers as 0...n-1 for nodes, elements, and sparse
1214   void setAscendingGlobalno(void);
1215   ///   The global numbers for elements runs across different types
1216   void setAbsoluteGlobalno();
1217   
1218   void copyOldGlobalno(const FEM_Mesh &m);
1219   void print(int idxBase);//Write a human-readable description to CkPrintf
1220   /// Extract a list of our entities:
1221   int getEntities(int *entites);
1222   
1223   
1224   
1225   /********** New methods ***********/
1226   /*
1227     This method creates the mapping from a node to all the elements that are 
1228     incident on it . It assumes the presence of one layer of ghost nodes that
1229     share a node.
1230   */
1231   void createNodeElemAdj();
1232   void createNodeNodeAdj();
1233   void createElemElemAdj();
1234   
1235   FEM_ElemAdj_Layer *getElemAdjLayer(void);
1236   
1237   // Terry's adjacency accessors & modifiers
1238
1239   //  ------- Element-to-element: preserve initial ordering relative to nodes
1240   /// Place all of element e's adjacent elements in neighbors; assumes
1241   /// neighbors allocated to correct size
1242   void e2e_getAll(int e, int *neighborss, int etype=0);
1243   /// Given id of element e, return the id of the idx-th adjacent element
1244   int e2e_getNbr(int e, short idx, int etype=0);
1245   /// Given id of element e and id of another element nbr, return i such that
1246   /// nbr is the i-th element adjacent to e
1247   int e2e_getIndex(int e, int nbr, int etype=0);
1248   /// Set the element adjacencies of element e to neighbors; assumes neighbors 
1249   /// has the correct size
1250   void e2e_setAll(int e, int *neighbors, int etype=0);
1251   /// Set the idx-th element adjacent to e to be newElem
1252   void e2e_setIndex(int e, short idx, int newElem, int etype=0);
1253   /// Find element oldNbr in e's adjacent elements and replace with newNbr
1254   void e2e_replace(int e, int oldNbr, int newNbr, int etype=0);
1255   /// Remove all neighboring elements in adjacency
1256   void e2e_removeAll(int e, int etype=0);
1257
1258
1259   //  ------- Element-to-node: preserve initial ordering
1260   /// Place all of element e's adjacent nodes in adjnodes; assumes
1261   /// adjnodes allocated to correct size
1262   void e2n_getAll(int e, int *adjnodes, int etype=0);
1263   /// Given id of element e, return the id of the idx-th adjacent node
1264   int e2n_getNode(int e, short idx, int etype=0);
1265   /// Given id of element e and id of a node n, return i such that
1266   /// n is the i-th node adjacent to e
1267   short e2n_getIndex(int e, int n, int etype=0);
1268   /// Set the node adjacencies of element e to adjnodes; assumes adjnodes 
1269   /// has the correct size
1270   void e2n_setAll(int e, int *adjnodes, int etype=0);
1271   /// Set the idx-th node adjacent to e to be newNode
1272   void e2n_setIndex(int e, short idx, int newNode, int etype=0);
1273   /// Find node oldNode in e's adjacent ndoes and replace with newNode
1274   void e2n_replace(int e, int oldNode, int newNode, int etype=0);
1275   /// Replace all entries with -1
1276   void e2n_removeAll(int e, int etype=0);
1277
1278
1279   //  ------- Node-to-node
1280   /// Place all of node n's adjacent nodes in adjnodes and the resulting 
1281   /// length of adjnodes in sz; assumes adjnodes is not allocated, but sz is
1282   void n2n_getAll(int n, int **adjnodes, int *sz);
1283   /// Adds newNode to node n's node adjacency list
1284   void n2n_add(int n, int newNode);
1285   /// Removes oldNode from n's node adjacency list
1286   void n2n_remove(int n, int oldNode);
1287   /// Finds oldNode in n's node adjacency list, and replaces it with newNode
1288   void n2n_replace(int n, int oldNode, int newNode);
1289   /// Remove all nodes from n's node adjacency list
1290   void n2n_removeAll(int n);
1291   /// Is queryNode in node n's adjacency vector?
1292   int n2n_exists(int n, int queryNode);
1293
1294   //  ------- Node-to-element
1295   /// Place all of node n's adjacent elements in adjelements and the resulting 
1296   /// length of adjelements in sz; assumes adjelements is not allocated, 
1297   /// but sz is
1298   void n2e_getAll(int n, int **adjelements, int *sz);
1299   /// Adds newElem to node n's element adjacency list
1300   void n2e_add(int n, int newElem);
1301   /// Removes oldElem from n's element adjacency list
1302   void n2e_remove(int n, int oldElem);
1303   /// Finds oldElem in n's element adjacency list, and replaces it with newElem
1304   void n2e_replace(int n, int oldElem, int newElem);
1305   /// Remove all elements from n's element adjacency list
1306   void n2e_removeAll(int n);
1307
1308   /// Get an element on edge (n1, n2) where n1, n2 are chunk-local
1309   /// node numberings and result is chunk-local element; return -1 in case 
1310   /// of failure
1311   int getElementOnEdge(int n1, int n2);
1312
1313   /// Get two elements adjacent to both n1 and n2
1314   void get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) ;
1315
1316
1317
1318
1319
1320 }; 
1321 PUPmarshall(FEM_Mesh);
1322 FEM_Mesh *FEM_Mesh_lookup(int fem_mesh,const char *caller);
1323 FEM_Entity *FEM_Entity_lookup(int fem_mesh,int entity,const char *caller);
1324 FEM_Attribute *FEM_Attribute_lookup(int fem_mesh,int entity,int attr,const char *caller);
1325
1326 void FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,     
1327         void *data, int firstItem,int length, const IDXL_Layout &layout);
1328
1329 //registration internal api
1330 void FEM_Register_array_layout(int fem_mesh,int entity,int attr,void *data,int firstItem,const IDXL_Layout &layout);
1331 void FEM_Register_entity_impl(int fem_mesh,int entity,void *args,int len,int max,FEM_Mesh_alloc_fn fn);
1332 /// Reassemble split chunks into a single mesh
1333 FEM_Mesh *FEM_Mesh_assemble(int nchunks,FEM_Mesh **chunks);
1334
1335 FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead);
1336 FEM_Mesh *FEM_readMesh(const char *prefix,int chunkNo,int nChunks);
1337 void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks);
1338
1339 /*\@}*/
1340
1341
1342 #endif