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