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