Acheived backward compatability for fem_mesh files by
[charm.git] / src / libs / ck-libs / fem / fem_impl.h
1 /*Charm++ Finite Element Framework:
2 C++ implementation file
3
4 This is the under-the-hood implementation file for FEM.
5 Orion Sky Lawlor, olawlor@acm.org, 9/28/00
6 */
7 #ifndef _FEM_IMPL_H
8 #define _FEM_IMPL_H
9
10 #include <stdio.h>
11
12 #include "charm-api.h"
13 #include "tcharm.h"
14 #include "ckvector3d.h"
15 #include "fem.h"
16
17 // temporary Datatype representation
18 // will go away once MPI user-defined datatypes are ready
19 struct DType {
20   int base_type; //FEM_* datatype
21   int vec_len; //Number of items of this datatype
22   int init_offset; // offset of field in bytes from the beginning of data
23   int distance; // distance in bytes between successive field values
24   DType(void) {}
25   DType( const int b,  const int v=1,  const int i=0,  const int d=0)
26     : base_type(b), vec_len(v), init_offset(i) 
27   {
28     distance = (d ? d : length());
29   }
30   //Default copy constructor, assignment operator
31
32   //Return the total number of bytes required by this FEM_* data type
33   static int type_size(int dataType) {
34     switch(dataType) {
35       case FEM_BYTE : return 1; break;
36       case FEM_INT : return sizeof(int); break;
37       case FEM_REAL : return sizeof(float); break;
38       case FEM_DOUBLE : return sizeof(double); break;
39       default: CkAbort("Unrecognized data type field passed to FEM framework!\n");
40     }
41     return -1;
42   }
43   
44   //Return the total number of bytes required by the data stored in this DType
45   int length(const int nitems=1) const {
46     return type_size(base_type) * vec_len * nitems;
47   }
48 };
49
50 //This datatype is how the framework stores symmetries internally.
51 //  Each bit of this type describes a different symmetry.
52 //  There must be enough bits to accomidate several simulatanious symmetries.
53 typedef unsigned char FEM_Symmetries_t;
54
55 class FEMinit {
56  public:
57         int numElements;
58         CkArrayID threads;
59         int flags;
60         CkChareID coordinator;
61         FEMinit() {}
62         FEMinit(int ne_,const CkArrayID &t_,int f_,const CkChareID &c_)
63                 :numElements(ne_),threads(t_),flags(f_),coordinator(c_) {}
64         void pup(PUP::er &p) {
65                 p|numElements;
66                 p|threads;
67                 p|flags;
68                 p|coordinator;
69         }
70 };
71 PUPmarshall(FEMinit);
72
73 /*This class describes a local-to-global index mapping, used in FEM_Print.
74 The default is the identity mapping.*/
75 class l2g_t {
76 public:
77         //Return the global number associated with this local element
78         virtual int el(int localNo) const {return localNo;}
79         //Return the global number associated with this local node
80         virtual int no(int localNo) const {return localNo;}
81 };
82
83 /* Map (user-assigned) numbers to T's */
84 template <class T>
85 class NumberedVec {
86         CkPupPtrVec<T, CkPupAlwaysAllocatePtr<T> > vec;
87         
88 public:
89         //Extend the vector to have up to this element
90         void makeLonger(int toHaveElement)
91         {
92                 int oldSize=vec.size(), newSize=toHaveElement+1;
93                 if (oldSize>=newSize) return; //Nothing to do
94                 vec.setSize(newSize);
95                 vec.length()=newSize;
96                 for (int j=oldSize;j<newSize;j++)
97                         vec[j]=new T;
98         }
99         //Reinitialize element i:
100         void reinit(int doomedEl) {
101                 vec[doomedEl].destroy();
102                 vec[doomedEl]=new T;
103         }
104         
105         int size(void) const {return vec.size();}
106         
107         //Same old bracket operators, but return the actual object, not a pointer:
108         T &operator[](int i) {return *( vec[i] );}
109         const T &operator[](int i) const {return *( vec[i] );}
110         
111         void pup(PUP::er &p) {
112                 vec.pup(p);
113         }
114         friend void operator|(PUP::er &p,NumberedVec<T> &v) {v.pup(p);}
115 };
116
117
118 /**
119  * FEM_Comm_Share describes how one item is shared with one other chunk.
120  * It lists the chunk and the item's location in the communication list;
121  * this location in the communication list is useful because it may be
122  * the *only* way to specify a particular item to another chunk
123  * if we don't have global numbers to identify things with.
124  */
125 class FEM_Comm_Share {
126  public:
127   int chk; //Chunk we're shared with
128   int idx; //Our index in the local comm. list for that chunk
129   FEM_Comm_Share(int x=0) {chk=idx=-1;}
130   FEM_Comm_Share(int c,int i) :chk(c), idx(i) {}
131   void pup(PUP::er &p) {p(chk); p(idx);}
132 };
133 PUPmarshall(FEM_Comm_Share);
134
135 /**
136  * FEM_Comm_Rec lists all the chunks that share an item.
137  */
138 class FEM_Comm_Rec {
139         int item; //Index of item we describe
140         CkVec<FEM_Comm_Share> shares;
141 public:
142         FEM_Comm_Rec(int item_=-1);
143         ~FEM_Comm_Rec();
144
145         void pup(PUP::er &p);
146         
147         inline int getItem(void) const {return item;}
148         inline int getShared(void) const {return shares.size();}
149         inline int getChk(int shareNo) const {return shares[shareNo].chk;}
150         inline int getIdx(int shareNo) const {return shares[shareNo].idx;}
151         bool hasChk(int chk) const {
152                 for (int i=0;i<getShared();i++)
153                         if (getChk(i)==chk) return true;
154                 return false;
155         }
156         void add(int chk,int idx);
157 };
158
159 /**
160  * Map an item to its FEM_Comm_Rec.  We only bother listing items
161  * that are actually shared with other chunks; completely local items
162  * are not present.
163  */
164 class FEM_Comm_Map {
165         CkHashtableT<CkHashtableAdaptorT<int>,FEM_Comm_Rec *> map;
166 public:
167         FEM_Comm_Map();
168         void pup(PUP::er &p);
169         ~FEM_Comm_Map();
170
171         //Add a FEM_Comm_. entry for this item
172         void add(int item,int chk,int idx);
173         
174         //Look up this item's FEM_Comm_Rec.  Returns NULL if item is not shared.
175         const FEM_Comm_Rec *get(int item) const;
176 };
177
178 /**
179  * FEM_Comm_List lists the items we share with one other chunk. 
180  * This list is used to build outgoing and interpret incoming messages--
181  * for an outgoing message, the listed items are copied into the message;
182  * for an incoming message, the listed items are copied out of the message.
183  */
184 class FEM_Comm_List {
185         int pe; //Global number of other chunk  
186         CkPupBasicVec<int> shared; //Local indices of shared items
187 public:
188         FEM_Comm_List();
189         FEM_Comm_List(int otherPe);
190         ~FEM_Comm_List();
191         int getDest(void) const {return pe;}
192         int size(void) const {return shared.size();}
193         int operator[](int idx) const {return shared[idx]; }
194         const int *getVec(void) const {return &shared[0];}
195         int push_back(int localIdx) {
196                 int ret=shared.size();
197                 shared.push_back(localIdx);
198                 return ret;
199         }
200         void pup(PUP::er &p);
201 };
202 PUPmarshall(FEM_Comm_List)
203
204 /**
205  * FEM_Comm describes all the shared items of a given chunk.
206  * It provides both item->chunks that share it (map)
207  * and chunk->items shared with it (comm)
208  */
209 class FEM_Comm : public CkNoncopyable {
210         //Indexed by local item number
211         FEM_Comm_Map map;
212         //Indexed by (local) chunk number
213         CkPupPtrVec<FEM_Comm_List, CkPupAlwaysAllocatePtr<FEM_Comm_List> > comm; 
214         
215         //Return the Comm_List associated with processor pe
216         FEM_Comm_List *getListN(int pe) { 
217                 for (int i=0;i<comm.size();i++)
218                         if (comm[i]->getDest()==pe)
219                                 return comm[i];
220                 return NULL; 
221         }
222 public:
223         FEM_Comm(void);
224         ~FEM_Comm();
225         int totalShared() const;//Return total number of shared nodes
226         void pup(PUP::er &p); //For migration
227         
228         int size(void) const {return comm.size();}
229         //Return the i'th (local) chunk we communicate with
230         const FEM_Comm_List &getLocalList(int idx) const { return *comm[idx]; }
231         int findLocalList(int pe) const {
232                 for (int i=0;i<comm.size();i++) 
233                         if (comm[i]->getDest()==pe)
234                                 return i;
235                 return -1;
236         }
237         
238         //Return the Comm_List associated with processor pe
239         const FEM_Comm_List &getList(int pe) const { 
240                 const FEM_Comm_List *ret=((FEM_Comm *)this)->getListN(pe);
241                 if (ret==NULL) CkAbort("FEM> Communication lists corrupted (unexpected message)");
242                 return *ret; 
243         }
244         //Return the FEM_Comm_List for processor pe, adding if needed
245         FEM_Comm_List &addList(int pe) {
246                 FEM_Comm_List *ret=getListN(pe);
247                 if (ret==NULL) { //Have to add a new list:
248                         ret=new FEM_Comm_List(pe);
249                         comm.push_back(ret);
250                 }
251                 return *ret;
252         }
253         
254         //Look up an item's FEM_Comm_Rec
255         const FEM_Comm_Rec *getRec(int item) const {return map.get(item);}
256         
257         //This item is shared with the given (local) chunk
258         void addNode(int localNo,int sharedWithChk) {
259                 map.add(localNo,sharedWithChk,
260                         comm[sharedWithChk]->push_back(localNo));
261         }
262         
263         //Used in creating comm. lists:
264         //Add this local number to both lists:
265         void add(int myChunk,int myLocalNo,
266                  int hisChunk,int hisLocalNo,FEM_Comm &hisList);
267         
268         void print(const l2g_t &l2g);
269 };
270
271
272 /*
273 This is a simple 2D table.  The operations are mostly row-centric.
274 */
275 template <class T>
276 class BasicTable2d : public CkNoncopyable {
277 protected:
278         int rows; //Number of entries in table
279         int cols; //Size of each entry in table
280         T *table; //Data in table [rows * cols]
281 public:
282         BasicTable2d(T *src,int cols_,int rows_) 
283                 :rows(rows_), cols(cols_), table(src) {}
284         
285         //"size" of the table is the number of rows:
286         inline int size(void) const {return rows;}
287         //Width of the table is the number of columns:
288         inline int width(void) const {return cols;}
289         
290         T *getData(void) {return table;}
291         const T *getData(void) const {return table;}
292         
293 //Element-by-element operations:
294         T operator() (int r,int c) const {return table[c+r*cols];}
295         T &operator() (int r,int c) {return table[c+r*cols];}
296         
297 //Row-by-row operations
298         //Get a pointer to a row of the table:
299         inline T *getRow(int r) {return &table[r*cols];}
300         inline const T *getRow(int r) const {return &table[r*cols];}
301         inline T *operator[](int r) {return getRow(r);}
302         inline const T *operator[](int r) const {return getRow(r);}
303         inline void setRow(int r,const T *src,T idxBase=0) {
304                 T *dest=getRow(r);
305                 for (int c=0;c<cols;c++) dest[c]=src[c]-idxBase;
306         }
307         inline void setRow(int r,T value) {
308                 T *dest=getRow(r);
309                 for (int c=0;c<cols;c++) dest[c]=value;
310         }
311         
312 //These affect the entire table:
313         void set(const T *src,T idxBase=0) {
314                 for (int r=0;r<rows;r++) 
315                 for (int c=0;c<cols;c++)
316                         table[c+r*cols]=src[c+r*cols]-idxBase;
317         }
318         void setTranspose(const T *srcT,int idxBase=0) {
319                 for (int r=0;r<rows;r++) 
320                 for (int c=0;c<cols;c++)
321                         table[c+r*cols]=srcT[r+c*rows]-idxBase;
322         }
323         void get(T *dest,T idxBase=0) const {
324                 for (int r=0;r<rows;r++) 
325                 for (int c=0;c<cols;c++)
326                         dest[c+r*cols]=table[c+r*cols]+idxBase;
327         }
328         void getTranspose(T *destT,int idxBase=0) const {
329                 for (int r=0;r<rows;r++) 
330                 for (int c=0;c<cols;c++)
331                         destT[r+c*rows]=table[c+r*cols]+idxBase;
332         }
333         void set(T value) {
334                 for (int r=0;r<rows;r++) setRow(r,value);
335         }
336 };
337
338 //As above, but heap-allocatable and resizable.
339 // T must not require a copy constructor.
340 template <class T>
341 class AllocTable2d : public BasicTable2d<T> {
342         int max; //Maximum number of rows that can be used without reallocation
343 public:
344         AllocTable2d(int cols_=0,int rows_=0) 
345                 :BasicTable2d<T>(NULL,cols_,rows_), max(0)
346         {
347                 if (rows>0) allocate(rows);
348         }
349         ~AllocTable2d() {delete[] table;}
350         void allocate(int rows_) { //Make room for this many rows
351                 allocate(width(),rows_,rows_);
352         }
353         void allocate(int cols_,int rows_,int max_=0) { //Make room for this many cols & rows
354                 int oldRows=rows;
355                 T *oldTable=table;
356                 if (max_==0) max_=rows_;
357                 cols=cols_;
358                 rows=rows_;
359                 max=max_;
360                 table=new T[max*cols];
361                 if (oldTable!=NULL) { //Preserve old table entries:
362                         memcpy(table,oldTable,sizeof(T)*cols*oldRows);
363                         delete[] oldTable;
364                 }
365         }
366         
367         //Pup routine and operator|:
368         void pup(PUP::er &p) {
369                 p|rows; p|cols;
370                 if (table==NULL) allocate(rows);
371                 p(table,rows*cols); //T better be a basic type, or this won't compile!
372         }
373         friend void operator|(PUP::er &p,AllocTable2d<T> &t) {t.pup(p);}
374
375         //Add a row to the table (by analogy with std::vector):
376         T *push_back(void) {
377                 if (rows>=max) 
378                 { //Not already enough room for the new row:
379                         int newMax=max+(max/4)+16; //Grow 25% longer
380                         allocate(cols,rows,newMax);
381                 }
382                 rows++;
383                 return getRow(rows-1);
384         }
385 };
386
387 //Smart pointer-to-new[]'d array-of-T
388 template <class T>
389 class ArrayPtrT : public CkNoncopyable {
390         T *sto;
391 public:
392         ArrayPtrT() {sto=NULL;}
393         ArrayPtrT(int *src) {sto=src;}
394         ~ArrayPtrT() {if (sto) delete[] sto;}
395         void operator=(T *src) {
396                 if (sto) delete[] sto;
397                 sto=src;
398         }
399         operator T *(void) {return sto;}
400         operator const T *(void) const {return sto;}
401 };
402 typedef ArrayPtrT<int> intArrayPtr;
403
404
405 /* Describes an FEM item-- a set of nodes or elements */
406 class FEM_Item : public CkNoncopyable {
407 public:
408         typedef AllocTable2d<double> udata_t;
409         typedef CkPupBasicVec<FEM_Symmetries_t> sym_t;
410 protected:
411         int ghostStart; //Index of first ghost object
412         udata_t udata; //Uninterpreted item-associated user data (one row per item)
413         sym_t *sym; //Symmetries of each item (or NULL if all 0)
414 public:
415         FEM_Comm ghostSend; //Non-ghosts we send out
416         FEM_Comm ghostRecv; //Ghosts we recv into
417         
418         FEM_Item() //Default constructor
419           {ghostStart=0;sym=NULL;}
420         ~FEM_Item() {if (sym) delete sym;}
421         void pup(PUP::er &p);
422         void print(const char *type,const l2g_t &l2g);
423         
424         //Manipulate the user data array:
425         void allocate(int nItems,int dataPer);
426         udata_t &setUdata(void) {return udata;}
427         const udata_t &getUdata(void) const {return udata;}
428         inline int size(void) const {return udata.size();}
429         inline int getDataPer(void) const {return udata.width();}
430         void udataIs(int r,const double *src) {udata.setRow(r,src);}
431         const double *udataFor(int r) const {return udata.getRow(r);}
432         
433         //Get ghost info:
434         void startGhosts(int atIndex) { ghostStart=atIndex; }
435         int getGhostStart(void) const { return ghostStart; }
436         int isGhostIndex(int idx) const { return idx>=ghostStart; }
437         
438         //Symmetry array:
439         const FEM_Symmetries_t *getSymmetries(void) const {
440                 if (sym==NULL) return NULL;
441                 else return sym->getVec();
442         }
443         FEM_Symmetries_t getSymmetries(int r) const { 
444                 if (sym==NULL) return FEM_Symmetries_t(0);
445                 else return (*sym)[r]; 
446         }
447         void setSymmetries(int r,FEM_Symmetries_t s);
448 };
449 PUPmarshall(FEM_Item);
450
451 /* Describes one kind of FEM elements */
452 class FEM_Elem:public FEM_Item {
453 public:
454         typedef AllocTable2d<int> conn_t;
455 private:
456         conn_t conn; //Connectivity data (one row per element)
457 public:
458         FEM_Elem():FEM_Item() {}
459         
460         void pup(PUP::er &p);
461         void print(const char *type,const l2g_t &l2g);
462         
463         void allocate(int nItems,int dataPer,int nodesPer);
464         conn_t &setConn(void) {return conn;}
465         const conn_t &getConn(void) const {return conn;}
466         int getConn(int elem,int nodeNo) const {return conn(elem,nodeNo);}
467         int getNodesPer(void) const {return conn.width();}
468         int *connFor(int i) {return conn.getRow(i);}
469         const int *connFor(int i) const {return conn.getRow(i);}
470         void connIs(int i,const int *src) {conn.setRow(i,src);}
471 };
472 PUPmarshall(FEM_Elem);
473
474 /**
475  * FEM_Sparse describes a set of records of sparse data that are all the
476  * same size and all associated with the same number of nodes.
477  * Sparse data is associated with some subset of the nodes in the mesh,
478  * and gets copied to every chunk that has all those nodes.  The canonical
479  * use of sparse data is to describe boundary conditions.
480  */
481 class FEM_Sparse : public CkNoncopyable {
482         AllocTable2d<int> nodes; //Each row is the nodes surrounding a tuple
483         AllocTable2d<char> data; //Each row is the user data for a tuple
484         intArrayPtr elem; // *OPTIONAL* partitioning based on elements (2*size() ints)
485 public:
486         void allocate(int n_); //Allocate storage for data and nodes of n tuples
487         
488         FEM_Sparse() { } //Used during pup
489         FEM_Sparse(int nodesPer_,int bytesPer_) 
490                 :nodes(nodesPer_), data(bytesPer_) { }
491         
492         //Return the number of records:
493         inline int size(void) const {return data.size();}
494         //Return the size of each record:
495         inline int getNodesPer(void) const {return nodes.width();}
496         inline int getDataPer(void) const {return data.width();}
497         
498         //Examine/change a single record:
499         inline const int *getNodes(int i) const {return nodes.getRow(i);}
500         inline const char *getData(int i) const {return data.getRow(i);}
501         inline int *getNodes(int i) {return nodes.getRow(i);}
502         inline char *getData(int i) {return data.getRow(i);}
503         inline void setNodes(int i,const int *d,int idxBase) {nodes.setRow(i,d,idxBase);}
504         inline void setData(int i,const char *d) {data.setRow(i,d);}
505         
506         //Allocate and set the entire table:
507         void set(int records,const int *n,int idxBase,const char *d);
508         
509         //Get the entire table:
510         void get(int *n,int idxBase,char *d) const;
511         
512         //Set the optional element-partitioning array
513         void setElem(int *elem_) {elem=elem_;}
514         int *getElem(void) {return elem;}
515         const int *getElem(void) const {return elem;}
516         
517         void pup(PUP::er &p) {p|nodes; p|data;}
518 };
519 PUPmarshall(FEM_Sparse);
520
521
522 /*Describes one kind of symmetry condition
523 */
524 class FEM_Sym_Desc : public PUP::able {
525 public:
526         virtual ~FEM_Sym_Desc();
527
528         //Apply this symmetry to this location vector
529         virtual CkVector3d applyLoc(const CkVector3d &loc) const =0;
530         
531         //Apply this symmetry to this relative (vel or acc) vector
532         virtual CkVector3d applyVec(const CkVector3d &vec) const =0;
533         
534         //Make a new copy of this class:
535         virtual FEM_Sym_Desc *clone(void) const =0;
536         
537         //Allows Desc's to be pup'd via | operator:
538         friend inline void operator|(PUP::er &p,FEM_Sym_Desc &a) {a.pup(p);}
539         friend inline void operator|(PUP::er &p,FEM_Sym_Desc* &a) {
540                 PUP::able *pa=a;  p(&pa);  a=(FEM_Sym_Desc *)pa;
541         }
542 };
543
544 //Describes a linear-periodic (space shift) symmetry:
545 class FEM_Sym_Linear : public FEM_Sym_Desc {
546         CkVector3d shift; //Offset to add to locations
547 public:
548         FEM_Sym_Linear(const CkVector3d &shift_) :shift(shift_) {}
549         FEM_Sym_Linear(CkMigrateMessage *m) {}
550         
551         //Apply this symmetry to this location vector
552         CkVector3d applyLoc(const CkVector3d &loc) const {return loc+shift;}
553         
554         //Apply this symmetry to this relative (vel or acc) vector
555         virtual CkVector3d applyVec(const CkVector3d &vec) const {return vec;}
556         
557         virtual FEM_Sym_Desc *clone(void) const {
558                 return new FEM_Sym_Linear(shift);
559         }
560         
561         virtual void pup(PUP::er &p);
562         PUPable_decl(FEM_Sym_Linear);
563 };
564
565 /*
566 Describes all the different kinds of symmetries that apply to
567 this mesh.
568 */
569 class FEM_Sym_List {
570         //This lists the different kinds of symmetry
571         CkPupAblePtrVec<FEM_Sym_Desc> sym; 
572         
573         FEM_Sym_List(const FEM_Sym_List &src); //NOT DEFINED: copy constructor
574 public:
575         FEM_Sym_List();
576         void operator=(const FEM_Sym_List &src); //Assignment operator
577         ~FEM_Sym_List();
578         
579         //Add a new kind of symmetry to this list, returning
580         // the way objects with that symmetry should be marked.
581         FEM_Symmetries_t add(FEM_Sym_Desc *desc);
582         
583         //Apply all the listed symmetries to this location
584         void applyLoc(CkVector3d *loc,FEM_Symmetries_t sym) const;
585         
586         //Apply all the listed symmetries to this relative vector
587         void applyVec(CkVector3d *vec,FEM_Symmetries_t sym) const;
588         
589         void pup(PUP::er &p);
590 };
591
592 /*This class describes the nodes and elements in
593   a finite-element mesh or submesh*/
594 class FEM_Mesh : public CkNoncopyable {
595         CkPupPtrVec<FEM_Sparse, CkPupAlwaysAllocatePtr<FEM_Sparse> > sparse;
596         FEM_Sym_List symList;
597 public:
598         void setSymList(const FEM_Sym_List &src) {symList=src;}
599         const FEM_Sym_List &getSymList(void) const {return symList;}
600
601         int nSparse(void) const {return sparse.size();}
602         void setSparse(int uniqueID,FEM_Sparse *s);
603         FEM_Sparse &setSparse(int uniqueID);
604         const FEM_Sparse &getSparse(int uniqueID) const;
605         
606         FEM_Item node; //Describes the nodes in the mesh
607         NumberedVec<FEM_Elem> elem; //Describes the different types of elements in the mesh
608         
609         //Set up our fields based on this mesh:
610         void makeRoom(const FEM_Mesh &src) {
611                 elem.makeLonger(src.elem.size()-1);
612                 setSymList(src.getSymList());
613         }
614         
615         //Return this type of element, given an element type
616         FEM_Item &setCount(int elTypeOrMinusOne) {
617                 if (elTypeOrMinusOne==-1) return node;
618                 else return elem[chkET(elTypeOrMinusOne)];
619         }
620         const FEM_Item &getCount(int elTypeOrMinusOne) const {
621                 if (elTypeOrMinusOne==-1) return node;
622                 else return elem[chkET(elTypeOrMinusOne)];
623         }
624         FEM_Elem &setElem(int elType) {return elem[chkET(elType)];}
625         const FEM_Elem &getElem(int elType) const {return elem[chkET(elType)];}
626         int chkET(int elType) const; //Check this element type-- abort if it's bad
627         
628         FEM_Mesh();
629         ~FEM_Mesh();
630         
631         int nElems() const //Return total number of elements (of all types)
632           {return nElems(elem.size());}
633         int nElems(int t) const;//Return total number of elements before type t
634         int getGlobalElem(int elType,int elNo) const;
635         void pup(PUP::er &p); //For migration
636         void print(const l2g_t &l2g);//Write human-readable description to CkPrintf
637 };
638
639 //Describes a single chunk of the finite-element mesh
640 class MeshChunk : public CkNoncopyable {
641  public:
642         FEM_Mesh m; //The chunk mesh
643         FEM_Comm comm; //Shared nodes
644         int *elemNums; // Maps local elem#-> global elem#  [m.nElems()]
645         int *nodeNums; // Maps local node#-> global node#  [m.node.n]
646         int *isPrimary; // Indicates us as owner of node  [m.node.n]
647
648         MeshChunk(void);
649         ~MeshChunk();
650         //Allocates elemNums, nodeNums, and isPrimary
651         void allocate() {
652           elemNums=new int[m.nElems()];
653           nodeNums=new int[m.node.size()];
654           isPrimary=new int[m.node.size()];
655         }
656         void pup(PUP::er &p); //For send/recv
657
658         void read(FILE *fp);
659         void write(FILE *fp);
660 };
661
662 class CallMeshUpdated {
663         int val; //Value to pass to function below
664         FEM_Update_mesh_fn cfn; //if 0, skip meshUpdated call
665         FEM_Update_mesh_fortran_fn ffn; //if 0, skip f90 meshUpdated call
666 public:
667         CallMeshUpdated() 
668                 :val(0), cfn(0), ffn(0) {}
669         CallMeshUpdated(FEM_Update_mesh_fn cfn_,int val_) 
670                 :val(val_), cfn(cfn_), ffn(0) {}
671         CallMeshUpdated(FEM_Update_mesh_fortran_fn ffn_,int val_) 
672                 :val(val_), cfn(0), ffn(ffn_) {}
673         void call(void) 
674         { //Call the user's meshUpdated function:
675                 if (cfn) { cfn(val); }
676                 if (ffn) { ffn(&val); }
677         }
678 };
679 PUPmarshallBytes(CallMeshUpdated);
680
681 class UpdateMeshChunk : public MeshChunk {
682 public:
683         //These fields are (only) used during an updateMesh
684         int updateCount; //Mesh update serial number
685         int fromChunk; //Source chunk
686         CallMeshUpdated meshUpdated;
687         int doWhat; //If 0, do nothing; if 1, repartition; if 2, resume
688         
689         UpdateMeshChunk() {
690                 updateCount=fromChunk=doWhat=0;
691         }
692         void pup(PUP::er &p) {
693                 MeshChunk::pup(p);
694                 p.comment(" UpdateMesh data: ");        
695                 p(updateCount); p(fromChunk);
696                 p|meshUpdated;
697                 p(doWhat);
698         }
699 };
700
701 /* Unmarshall into a heap-allocated copy */
702 template<class T>
703 class marshallNewHeapCopy {
704         T *cur;
705 public:
706         marshallNewHeapCopy(T *readFrom) :cur(readFrom) {}
707         marshallNewHeapCopy(const marshallNewHeapCopy &h) :cur(h.cur) {}
708         marshallNewHeapCopy(void) {
709                 cur=new T;
710         }
711         
712         void pup(PUP::er &p) {
713                 cur->pup(p);
714         }
715         operator T *() {return cur;}
716         friend void operator|(PUP::er &p,marshallNewHeapCopy<T> &h) {h.pup(p);}
717 };
718 typedef marshallNewHeapCopy<UpdateMeshChunk> marshallUpdateMeshChunk;
719 typedef marshallNewHeapCopy<MeshChunk> marshallMeshChunk;
720
721
722 #include "fem.decl.h"
723
724 #define CHK(p) do{if((p)==0)CkAbort("FEM>Memory Allocation failure.");}while(0)
725
726 class FEM_DataMsg : public CMessage_FEM_DataMsg
727 {
728  public:
729   int from; //Source's chunk number
730   int dtype; //Field ID of data
731   int length; //Length in bytes of below array
732   int tag; //Sequence number
733   void *data;
734   double alignPad; //Makes sure this structure is double-aligned
735   
736   FEM_DataMsg(int t, int f, int d,int l) : 
737     from(f), dtype(d), length(l), tag(t) { data = (void*) (this+1); }
738   FEM_DataMsg(void) { data = (void*) (this+1); }
739   static void *pack(FEM_DataMsg *);
740   static FEM_DataMsg *unpack(void *);
741   static void *alloc(int, size_t, int*, int);
742 };
743
744 /* Maximum number of fields that can be registered */
745 #define MAXDT 20
746
747 class FEMchunk : public ArrayElement1D
748 {
749 public:
750 // updated_mesh keeps the still-being-assembled next mesh chunk.
751 // It is created and written by the FEM_Set routines called from driver.
752   UpdateMeshChunk *updated_mesh;
753   int updateCount; //Number of mesh updates
754
755   //The current finite-element mesh
756   MeshChunk *cur_mesh;
757   
758 private:
759   FEMinit init;
760
761   CProxy_FEMchunk thisproxy;
762   TCharm *thread;
763   
764   DType dtypes[MAXDT];
765   int ntypes;
766
767   CmmTable messages; // update messages to be processed
768   int updateSeqnum; // sequence number for last update operation
769
770   //Describes the current data we're waiting for:
771   typedef enum {INVALID_UPDATE,NODE_UPDATE,GHOST_UPDATE} updateType_t;
772   updateType_t updateType;
773
774   const FEM_Comm *updateComm; //Communicator we're blocked on
775   int nRecd; //Number of messages received for this update
776   void *updateBuf; //User data addr for current update
777
778   void beginUpdate(void *buf,int fid,
779                   const FEM_Comm *sendComm,const FEM_Comm *recvComm,updateType_t t);
780   void recvUpdate(FEM_DataMsg *);
781   void update_node(FEM_DataMsg *);
782   void update_ghost(FEM_DataMsg *);
783   void waitForUpdate(void);
784
785   void *reductionBuf; //Place to return reduction result
786
787   CkVec<int> listTmp;//List of local items 
788   int listCount; //Number of lists received
789   bool listSuspended;
790   bool finishListExchange(const FEM_Comm &l);
791
792   void initFields(void);
793   void setMesh(MeshChunk *msg=0);
794
795  public:
796
797   int tsize;
798   int doneCalled;
799
800   FEMchunk(const FEMinit &init);
801   FEMchunk(CkMigrateMessage *msg);
802   ~FEMchunk();
803
804   void ckJustMigrated(void);
805
806   void run(void);
807   void run(marshallMeshChunk &);
808   void reductionResult(FEM_DataMsg *);
809   void updateMesh(int doWhat);
810   void meshUpdated(marshallMeshChunk &);
811   void meshUpdatedComplete(void) {thread->resume();}
812
813   int new_DT(const DType &d) {
814     if(ntypes>=MAXDT) {
815       CkAbort("FEM_Create_Field> Too many registered datatypes!");
816     }
817     dtypes[ntypes] = d;
818     ntypes++;
819     return ntypes-1;
820   }
821   
822   void update(int fid, void *nodes);
823   void updateGhost(int fid, int elemType, void *nodes);
824   void recv(FEM_DataMsg *);
825   
826   void exchangeGhostLists(int elemType,int inLen,const int *inList,int idxbase);
827   void recvList(int elemType,int fmChk,int nIdx,const int *idx);
828   const CkVec<int> &getList(void) {return listTmp;}
829   void emptyList(void) {listTmp.length()=0;}
830   
831   void reduce_field(int fid, const void *nodes, void *outbuf, int op);
832   void reduce(int fid, const void *inbuf, void *outbuf, int op);
833   void readField(int fid, void *nodes, const char *fname);
834   void print(void);
835   FEM_Mesh &getMesh(void) { return cur_mesh->m; }
836   int *getNodeNums(void) { return cur_mesh->nodeNums; }
837   int *getElemNums(void) { return cur_mesh->elemNums; }
838   int getPrimary(int nodeNo) { return cur_mesh->isPrimary[nodeNo]; }
839
840   const FEM_Comm &getComm(void) const {return cur_mesh->comm;}
841
842   void pup(PUP::er &p);
843 };
844
845 /*Partition this mesh's elements into n chunks,
846  writing each element's 0-based chunk number to elem2chunk.
847 */
848 void fem_partition(const FEM_Mesh *mesh,int nchunks,int *elem2chunk);
849
850 //Describes a single layer of ghost elements
851 class ghostLayer : public CkNoncopyable {
852 public:
853         int nodesPerTuple; //Number of shared nodes needed to connect elements
854         bool addNodes; //Add ghost nodes to the chunks
855         class elemGhostInfo {
856         public:
857                 bool add; //Add this kind of ghost element to the chunks
858                 int tuplesPerElem; //# of tuples surrounding this element
859                 intArrayPtr elem2tuple; //The tuples around this element [nodesPerTuple * tuplesPerElem]
860                 elemGhostInfo(void) {add=false;tuplesPerElem=0;}
861                 ~elemGhostInfo(void) {}
862                 void pup(PUP::er &p) {CkAbort("FEM> Shouldn't call elemGhostInfo::pup!\n");}
863         };
864         NumberedVec<elemGhostInfo> elem;
865 };
866
867 //Accumulates all symmetries of the mesh before splitting:
868 class FEM_Initial_Symmetries; /*Defined in symmetries.C*/
869
870 //Describes all ghost elements
871 class FEM_Ghost : public CkNoncopyable {
872         CkVec<ghostLayer *> layers;
873         
874         FEM_Initial_Symmetries *sym;
875 public:
876         FEM_Ghost() {sym=NULL;}
877         ~FEM_Ghost() {for (int i=0;i<getLayers();i++) delete layers[i];}
878         
879         int getLayers(void) const {return layers.size();}
880         ghostLayer *addLayer(void) {
881                 ghostLayer *l=new ghostLayer();
882                 layers.push_back(l);
883                 return l;
884         }
885         const ghostLayer &getLayer(int layerNo) const {return *layers[layerNo];}
886         
887         void setSymmetries(int nNodes_,int *new_can,const int *sym_src);
888         void addLinearPeriodic(int nFaces_,int nPer,
889                 const int *facesA,const int *facesB,int idxBase,
890                 int nNodes_,const CkVector3d *nodeLocs);
891         const int *getCanon(void) const;
892         const FEM_Symmetries_t *getSymmetries(void) const;
893         const FEM_Sym_List &getSymList(void) const;
894 };
895
896 //Declare this at the start of every API routine:
897 #define FEMAPI(routineName) TCHARM_API_TRACE(routineName,"fem")
898
899 /*A way to stream out partitioned chunks of a mesh.
900   By streaming, we can send the chunks as they are built,
901   dramatically reducing the memory needed by the framework.
902 */
903 class MeshChunkOutput {
904  public:
905         virtual ~MeshChunkOutput() {} /*<- for whining compilers*/
906         //Transfer ownership of this mesh chunk
907         virtual void accept(int chunkNo,MeshChunk *msg) =0;
908 };
909
910 /*After partitioning, create a sub-mesh for each chunk's elements,
911 including communication lists between chunks.
912 */
913 void fem_split(FEM_Mesh *mesh,int nchunks,const int *elem2chunk,
914                const FEM_Ghost &ghosts,MeshChunkOutput *out);
915
916 /*The inverse of fem_split: reassemble split chunks into a single mesh*/
917 FEM_Mesh *fem_assemble(int nchunks,MeshChunk **msgs);
918
919
920 //Make a new[]'d copy of this (len-item) array, changing the index as spec'd
921 int *CkCopyArray(const int *src,int len,int indexBase);
922 const FEM_Mesh *FEM_Get_FEM_Mesh(void);
923 FEM_Ghost &FEM_Set_FEM_Ghost(void);
924
925 #endif
926
927