1e16ce83aa8b4dfc3caac745826c6eaea14bee4f
[charm.git] / src / libs / ck-libs / ParFUM / mesh.C
1 /*
2 Finite Element Method (FEM) Framework for Charm++
3 Parallel Programming Lab, Univ. of Illinois 2002
4 Orion Sky Lawlor, olawlor@acm.org, 12/20/2002
5
6 FEM Implementation file: mesh creation and user-data manipulation.
7 */
8
9 #include "ParFUM.h"
10 #include "ParFUM_internals.h"
11
12 extern int femVersion;
13
14 /*** IDXL Interface ****/
15 FEM_Comm_Holder::FEM_Comm_Holder(FEM_Comm *sendComm, FEM_Comm *recvComm)
16         :comm(sendComm,recvComm)
17 {
18         registered=false;
19         idx=-1; 
20 }
21 void FEM_Comm_Holder::registerIdx(IDXL_Chunk *c) {
22         assert(!registered);
23         IDXL_Chunk *owner=c;
24         if (idx!=-1) // had an old index: try to get it back
25                 idx=owner->addStatic(&comm,idx);
26         else //No old index:
27                 idx=owner->addStatic(&comm);
28 }
29 void FEM_Comm_Holder::pup(PUP::er &p) {
30         p|idx;
31         if (p.isUnpacking() && idx!=-1) 
32         { // Try to grab the same index we had on our old processor:
33                 registerIdx(IDXL_Chunk::get("FEM_Comm_Holder::pup"));
34         }
35 }
36
37 FEM_Comm_Holder::~FEM_Comm_Holder(void)
38 {
39         if (registered) 
40         { // Try to unregister from IDXL:
41                 const char *caller="FEM_Comm_Holder::~FEM_Comm_Holder";
42                 IDXL_Chunk *owner=IDXL_Chunk::getNULL();
43                 if (owner) owner->destroy(idx,caller); 
44         }
45 }
46
47 /******* FEM_Mesh API ******/
48
49 static void checkIsSet(int fem_mesh,bool wantSet,const char *caller) {
50         if (FEM_Mesh_is_set(fem_mesh)!=wantSet) {
51                 const char *msg=wantSet?"This mesh (%d) is not a setting mesh":
52                         "This mesh (%d) is not a getting mesh";
53                 FEM_Abort(caller,msg,fem_mesh);
54         }
55 }
56
57 /* Connectivity: Map calls to appropriate version of FEM_Mesh_data */
58 CDECL void 
59 FEM_Mesh_conn(int fem_mesh,int entity,
60         int *conn, int firstItem, int length, int width) 
61 {
62         FEM_Mesh_data(fem_mesh,entity,FEM_CONN, conn, firstItem,length, FEM_INDEX_0, width);
63 }
64 FDECL void 
65 FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(int *fem_mesh,int *entity,
66         int *conn, int *firstItem,int *length, int *width)
67 {
68         //Can't just call the C version of this routine, because we use 1-based indices:
69         FEM_Mesh_data(*fem_mesh,*entity,FEM_CONN, conn, *firstItem-1,*length, FEM_INDEX_1, *width);
70 }
71
72 CDECL void
73 FEM_Mesh_set_conn(int fem_mesh,int entity,
74         const int *conn, int firstItem,int length, int width)
75 {
76         checkIsSet(fem_mesh,true,"FEM_Mesh_set_conn");
77         FEM_Mesh_conn(fem_mesh,entity,(int *)conn,firstItem,length,width);
78 }
79 FDECL void
80 FTN_NAME(FEM_MESH_SET_CONN,fem_mesh_set_conn)(int *fem_mesh,int *entity,
81         const int *conn, int *firstItem,int *length, int *width)
82 {
83         checkIsSet(*fem_mesh,true,"fem_mesh_set_conn");
84         FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(fem_mesh,entity,(int *)conn,firstItem,length,width);
85 }
86
87 CDECL void
88 FEM_Mesh_get_conn(int fem_mesh,int entity,
89         int *conn, int firstItem,int length, int width)
90 {
91         checkIsSet(fem_mesh,false,"FEM_Mesh_get_conn");
92         FEM_Mesh_conn(fem_mesh,entity,conn,firstItem,length,width);
93 }
94 FDECL void
95 FTN_NAME(FEM_MESH_GET_CONN,fem_mesh_get_conn)(int *fem_mesh,int *entity,
96         int *conn, int *firstItem,int *length, int *width)
97 {
98         checkIsSet(*fem_mesh,false,"fem_mesh_get_conn");
99         FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(fem_mesh,entity,conn,firstItem,length,width);
100 }
101
102
103 /* Data: map to FEM_Mesh_offset */
104 CDECL void
105 FEM_Mesh_data(int fem_mesh,int entity,int attr,         
106         void *data, int firstItem,int length, int datatype,int width)
107 {
108         IDXL_Layout lo(datatype,width);
109         FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,lo);
110 }
111
112 FORTRAN_AS_C(FEM_MESH_DATA,FEM_Mesh_data,fem_mesh_data,
113         (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
114         (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
115 )
116
117
118 CDECL void
119 FEM_Mesh_set_data(int fem_mesh,int entity,int attr,     
120         const void *data, int firstItem,int length, int datatype,int width)
121 {
122         checkIsSet(fem_mesh,true,"FEM_Mesh_set_data");
123         FEM_Mesh_data(fem_mesh,entity,attr,(void *)data,firstItem,length,datatype,width);
124 }
125 FORTRAN_AS_C(FEM_MESH_SET_DATA,FEM_Mesh_set_data,fem_mesh_set_data,
126         (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
127         (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
128 )
129
130 CDECL void
131 FEM_Mesh_get_data(int fem_mesh,int entity,int attr,     
132         void *data, int firstItem,int length, int datatype,int width)
133 {
134         checkIsSet(fem_mesh,false,"FEM_Mesh_get_data");
135         FEM_Mesh_data(fem_mesh,entity,attr,data,firstItem,length,datatype,width);
136 }
137 FORTRAN_AS_C(FEM_MESH_GET_DATA,FEM_Mesh_get_data,fem_mesh_get_data,
138         (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
139         (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
140 )
141
142 CDECL void
143 FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,  
144         void *data, int firstItem,int length, IDXL_Layout_t layout)
145 {
146         const char *caller="FEM_Mesh_data_layout";
147         FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,
148                 IDXL_Layout_List::get().get(layout,caller));
149 }
150
151 FORTRAN_AS_C(FEM_MESH_DATA_LAYOUT,FEM_Mesh_data_layout,fem_mesh_data_layout,
152         (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *layout),
153         (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*layout)
154 )
155
156 CDECL void
157 FEM_Mesh_data_offset(int fem_mesh,int entity,int attr,  
158         void *data, int firstItem,int length,
159         int type,int width, int offsetBytes,int distanceBytes,int skewBytes)
160 {
161         const char *caller="FEM_Mesh_data_offset";
162         FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,
163                 IDXL_Layout(type,width,offsetBytes,distanceBytes,skewBytes));
164 }
165 FORTRAN_AS_C(FEM_MESH_DATA_OFFSET,FEM_Mesh_data_offset,fem_mesh_data_offset,
166         (int *fem_mesh,int *entity,int *attr,
167          void *data,int *firstItem,int *length,
168          int *type,int *width,int *offset,int *distance,int *skew),
169         (*fem_mesh,*entity,*attr,
170          data,*firstItem-1,*length,
171          *type,*width,*offset,*distance,*skew)
172 )
173
174
175 void FEM_Register_array(int fem_mesh,int entity,int attr,
176         void *data, int datatype,int width,int firstItem){
177         IDXL_Layout lo(datatype,width);
178 /*      if(attr == FEM_CONN){
179                 printf("CONN width %d \n",width);
180                 int len = FEM_Mesh_get_length(fem_mesh,entity);
181                 int *connd = (int *)data;
182                 for(int i=0;i<len;i++){
183                         printf("%d -> (%d %d %d) \n",i+1,connd[3*i],connd[3*i+1],connd[3*i+2]);
184                 }
185         }
186         printf("firstItem %d \n",firstItem);*/
187         FEM_Register_array_layout(fem_mesh,entity,attr,data,firstItem,lo);
188 }
189
190 void FEM_Register_array_layout(int fem_mesh,int entity,int attr,        
191         void *data, IDXL_Layout_t layout,int firstItem){
192         const char *caller="FEM_Register_array_layout";
193         FEM_Register_array_layout(fem_mesh,entity,attr,data,firstItem, 
194                 IDXL_Layout_List::get().get(layout,caller));
195
196 }
197
198 /*registration api */
199 CDECL void 
200 FEM_Register_array(int fem_mesh,int entity,int attr,
201         void *data, int datatype,int width)
202 {       
203         FEM_Register_array(fem_mesh,entity,attr,data,datatype,width,0);
204 }
205
206 CDECL void
207 FEM_Register_array_layout(int fem_mesh,int entity,int attr,     
208         void *data, IDXL_Layout_t layout){
209         FEM_Register_array_layout(fem_mesh,entity,attr,data,layout,0);
210 }
211
212
213 CDECL void 
214 FEM_Register_entity(int fem_mesh,int entity,void *data,
215                 int len,int max,FEM_Mesh_alloc_fn fn) {
216                 FEM_Register_entity_impl(fem_mesh,entity,data,len,max,fn);
217 }
218
219 /**TODO: add the fortran api for registration*/
220
221 FORTRAN_AS_C(FEM_REGISTER_ARRAY,FEM_Register_array,fem_register_array,
222         (int *fem_mesh,int *entity,int *attr,void *data,int *datatype,int *width),(*fem_mesh,*entity,*attr,data,*datatype,*width,0))
223
224
225 FORTRAN_AS_C(FEM_REGISTER_ARRAY_LAYOUT,FEM_Register_array_layout,fem_register_array_layout,
226         (int *fem_mesh,int *entity,int *attr,void *data,int *layout),(*fem_mesh,*entity,*attr,data,*layout,0))
227
228 FORTRAN_AS_C(FEM_REGISTER_ENTITY,FEM_Register_entity,fem_register_entity,
229         (int *fem_mesh,int *entity,void *data,int *len,int *max,FEM_Mesh_alloc_fn fn),(*fem_mesh,*entity,data,*len,*max,fn))
230
231
232 // User data API:
233 CDECL void 
234 FEM_Mesh_pup(int fem_mesh,int dataTag,FEM_Userdata_fn fn,void *data) {
235         const char *caller="FEM_Mesh_pup"; FEMAPI(caller);
236         FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
237         FEM_Userdata_item &i=m->udata.find(dataTag);
238         FEM_Userdata_pupfn f(fn,data);
239         if (m->isSetting()) i.store(f);
240         else /* m->isGetting() */ {
241                 if (!i.hasStored())
242                         FEM_Abort(caller,"Never stored any user data at tag %d",dataTag);
243                 i.restore(f);
244         }
245 }
246 FORTRAN_AS_C(FEM_MESH_PUP,FEM_Mesh_pup,fem_mesh_pup,
247         (int *m,int *t,FEM_Userdata_fn fn,void *data), (*m,*t,fn,data))
248
249 // Accessor API:
250 CDECL void 
251 FEM_Mesh_set_length(int fem_mesh,int entity,int newLength) {
252         const char *caller="FEM_Mesh_set_length"; FEMAPI(caller);
253         checkIsSet(fem_mesh,true,caller);
254         FEM_Entity_lookup(fem_mesh,entity,caller)->setLength(newLength);
255 }
256 FORTRAN_AS_C(FEM_MESH_SET_LENGTH,FEM_Mesh_set_length,fem_mesh_set_length,
257         (int *fem_mesh,int *entity,int *newLength),
258         (*fem_mesh,*entity,*newLength)
259 )
260
261
262 CDECL int 
263 FEM_Mesh_get_length(int fem_mesh,int entity) {
264         const char *caller="FEM_Mesh_get_length"; FEMAPI(caller);
265         int len=FEM_Entity_lookup(fem_mesh,entity,caller)->size();
266         return len;
267 }
268 FORTRAN_AS_C_RETURN(int,
269         FEM_MESH_GET_LENGTH,FEM_Mesh_get_length,fem_mesh_get_length,
270         (int *fem_mesh,int *entity),(*fem_mesh,*entity)
271 )
272
273
274 CDECL void 
275 FEM_Mesh_set_width(int fem_mesh,int entity,int attr,int newWidth) {
276         const char *caller="FEM_Mesh_set_width";
277         FEMAPI(caller);
278         checkIsSet(fem_mesh,true,caller);
279         FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->setWidth(newWidth,caller);
280 }
281 FORTRAN_AS_C(FEM_MESH_SET_WIDTH,FEM_Mesh_set_width,fem_mesh_set_width,
282         (int *fem_mesh,int *entity,int *attr,int *newWidth),
283         (*fem_mesh,*entity,*attr,*newWidth)
284 )
285
286 CDECL int 
287 FEM_Mesh_get_width(int fem_mesh,int entity,int attr) {
288         const char *caller="FEM_Mesh_get_width";
289         FEMAPI(caller);
290         return FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->getWidth();
291 }
292 FORTRAN_AS_C_RETURN(int,
293         FEM_MESH_GET_WIDTH,FEM_Mesh_get_width,fem_mesh_get_width,
294         (int *fem_mesh,int *entity,int *attr),(*fem_mesh,*entity,*attr)
295 )
296
297 CDECL int 
298 FEM_Mesh_get_datatype(int fem_mesh,int entity,int attr) {
299         const char *caller="FEM_Mesh_get_datatype";
300         FEMAPI(caller);
301         return FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->getDatatype();
302 }
303 FORTRAN_AS_C_RETURN(int,
304         FEM_MESH_GET_DATATYPE,FEM_Mesh_get_datatype,fem_mesh_get_datatype,
305         (int *fem_mesh,int *entity,int *attr),(*fem_mesh,*entity,*attr)
306 )
307
308 CDECL int 
309 FEM_Mesh_is_set(int fem_mesh) /* return 1 if this is a writing mesh */
310 {
311         return (FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_is_get")->isSetting())?1:0;
312 }
313 FORTRAN_AS_C_RETURN(int,
314         FEM_MESH_IS_SET,FEM_Mesh_is_set,fem_mesh_is_set,
315         (int *fem_mesh),(*fem_mesh)
316 )
317
318 CDECL int 
319 FEM_Mesh_is_get(int fem_mesh) /* return 1 if this is a readable mesh */
320 {
321         return (!FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_is_get")->isSetting())?1:0;
322 }
323 FORTRAN_AS_C_RETURN(int,
324         FEM_MESH_IS_GET,FEM_Mesh_is_get,fem_mesh_is_get,
325         (int *fem_mesh),(*fem_mesh)
326 )
327
328 CDECL void 
329 FEM_Mesh_become_get(int fem_mesh) /* Make this a readable mesh */
330 { FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_become_get")->becomeGetting(); }
331 FORTRAN_AS_C(FEM_MESH_BECOME_GET,FEM_Mesh_become_get,fem_mesh_become_get, (int *m),(*m))
332
333 CDECL void 
334 FEM_Mesh_become_set(int fem_mesh)
335 { FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_become_get")->becomeSetting(); }
336 FORTRAN_AS_C(FEM_MESH_BECOME_SET,FEM_Mesh_become_set,fem_mesh_become_set, (int *m),(*m))
337
338
339 CDECL IDXL_t 
340 FEM_Comm_shared(int fem_mesh,int entity) {
341         const char *caller="FEM_Comm_shared";
342         FEMAPI(caller); 
343         if (entity!=FEM_NODE) FEM_Abort(caller,"Only shared nodes supported");
344         return FEM_Mesh_lookup(fem_mesh,caller)->node.
345                 sharedIDXL.getIndex(IDXL_Chunk::get(caller));
346 }
347 FORTRAN_AS_C_RETURN(int,
348         FEM_COMM_SHARED,FEM_Comm_shared,fem_comm_shared,
349         (int *fem_mesh,int *entity),(*fem_mesh,*entity)
350 )
351
352 CDECL IDXL_t 
353 FEM_Comm_ghost(int fem_mesh,int entity) {
354         const char *caller="FEM_Comm_ghost";
355         FEMAPI(caller);
356         FEM_Entity *e=FEM_Entity_lookup(fem_mesh,entity,caller);
357         if (e->isGhost()) FEM_Abort(caller,"Can only call FEM_Comm_ghost on real entity type");
358         return e->ghostIDXL.getIndex(IDXL_Chunk::get(caller));
359 }
360 FORTRAN_AS_C_RETURN(int,
361         FEM_COMM_GHOST,FEM_Comm_ghost,fem_comm_ghost,
362         (int *fem_mesh,int *entity),(*fem_mesh,*entity)
363 )
364
365
366 // Internal API:
367 void FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,     
368         void *data, int firstItem,int length, const IDXL_Layout &layout) 
369 {
370         if (femVersion == 0 && length==0) return;
371         const char *caller="FEM_Mesh_data";
372         FEMAPI(caller);
373         FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
374         FEM_Attribute *a=m->lookup(entity,caller)->lookup(attr,caller);
375         
376         if (m->isSetting()) 
377                 a->set(data,firstItem,length,layout,caller);
378         else /* m->isGetting()*/
379                 a->get(data,firstItem,length,layout,caller);
380 }
381
382 /** the internal registration function */
383 void FEM_Register_array_layout(int fem_mesh,int entity,int attr,void *data,int firstItem,const IDXL_Layout &layout){
384         const char *caller="FEM_Register_array";
385         FEMAPI(caller);
386         FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
387         FEM_Entity *e = m->lookup(entity,caller);
388         int length = e->size();
389         //should actually be a call on the entity
390         int max = e->getMax();
391         FEM_Attribute *a = e->lookup(attr,caller);
392         if(a->getWidth() == 0){
393                 /** While registering an attribute for the first time
394                  * Its width should get set correctly
395                  * **/
396         //      a->setDatatype(layout.type);
397                 a->setWidth(layout.width);
398         }
399         
400         
401         if(m->isSetting()){
402           a->register_data(data,length,max,layout,caller);
403         }else{
404           a->get(data,firstItem,length,layout,caller);
405           //replace the attribute's data array with the user's data
406           a->register_data(data,max/*length*/,max,layout,caller);
407         }
408 }
409 void FEM_Register_entity_impl(int fem_mesh,int entity,void *args,int len,int max,FEM_Mesh_alloc_fn fn){
410         char *caller = "FEM_Register_entity";
411         FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
412 /*      if(!m->isSetting()){
413                 CmiAbort("Register entity called on mesh that can't be written into");
414         }
415 */
416         FEM_Entity *e = m->lookup(entity,caller);
417         e->setMaxLength(len,max,args,fn);
418 }
419
420 FEM_Entity *FEM_Entity_lookup(int fem_mesh,int entity,const char *caller) {
421         return FEM_Mesh_lookup(fem_mesh,caller)->lookup(entity,caller);
422 }
423 FEM_Attribute *FEM_Attribute_lookup(int fem_mesh,int entity,int attr,const char *caller) {
424         return FEM_Entity_lookup(fem_mesh,entity,caller)->lookup(attr,caller);
425 }
426
427 CDECL int FEM_Mesh_get_entities(int fem_mesh, int *entities) {
428         const char *caller="FEM_Mesh_get_entities";
429         FEMAPI(caller); 
430         return FEM_Mesh_lookup(fem_mesh,caller)->getEntities(entities);
431 }
432 FORTRAN_AS_C_RETURN(int,
433         FEM_MESH_GET_ENTITIES,FEM_Mesh_get_entities,fem_mesh_get_entities,
434         (int *mesh,int *ent), (*mesh,ent)
435 )
436
437 CDECL int FEM_Mesh_get_attributes(int fem_mesh,int entity,int *attributes) {
438         const char *caller="FEM_Mesh_get_attributes";
439         FEMAPI(caller);
440         return FEM_Entity_lookup(fem_mesh,entity,caller)->getAttrs(attributes);
441 }
442 FORTRAN_AS_C_RETURN(int,
443         FEM_MESH_GET_ATTRIBUTES,FEM_Mesh_get_attributes,fem_mesh_get_attributes,
444         (int *mesh,int *ent,int *attrs), (*mesh,*ent,attrs)
445 )
446
447 /************** FEM_Attribute ****************/
448
449 CDECL const char *FEM_Get_datatype_name(int datatype,char *storage) {
450         switch(datatype) {
451         case FEM_BYTE: return "FEM_BYTE";
452         case FEM_INT: return "FEM_INT";
453         case FEM_FLOAT: return "FEM_FLOAT";
454         case FEM_DOUBLE: return "FEM_DOUBLE";
455         case FEM_INDEX_0: return "FEM_INDEX_0";
456         case FEM_INDEX_1: return "FEM_INDEX_1";
457         };
458         sprintf(storage,"unknown datatype code (%d)",datatype);
459         return storage;
460 }
461
462 /// Return the human-readable version of this FEM_ATTR code.
463 ///  For example, FEM_attr2name(FEM_CONN)=="FEM_CONN".
464 CDECL const char *FEM_Get_attr_name(int attr,char *storage) 
465 {
466         if (attr<FEM_ATTRIB_TAG_MAX) 
467         { //It's a user tag:
468                 sprintf(storage,"FEM_DATA+%d",attr-FEM_DATA);
469                 return storage;
470         }
471         switch(attr) {
472         case FEM_CONN: return "FEM_CONN"; 
473         case FEM_SPARSE_ELEM: return "FEM_SPARSE_ELEM";
474         case FEM_COOR: return "FEM_COOR";
475         case FEM_GLOBALNO: return "FEM_GLOBALNO";
476         case FEM_PARTITION: return "FEM_PARTITION";
477         case FEM_SYMMETRIES: return "FEM_SYMMETRIES";
478         case FEM_NODE_PRIMARY: return "FEM_NODE_PRIMARY";
479         case FEM_NODE_ELEM_ADJACENCY: return "FEM_NODE_ELEM_ADJACENCY";
480         case FEM_NODE_NODE_ADJACENCY: return "FEM_NODE_NODE_ADJACENCY";
481         case FEM_ELEM_ELEM_ADJACENCY: return "FEM_ELEM_ELEM_ADJACENCY";
482         case FEM_ELEM_ELEM_ADJ_TYPES: return "FEM_ELEM_ELEM_ADJ_TYPES";
483         case FEM_IS_VALID_ATTR: return "FEM_IS_VALID_ATTR";
484         case FEM_MESH_SIZING: return "FEM_MESH_SIZING";
485
486         default: break;
487         };
488         sprintf(storage,"unknown attribute code (%d)",attr);
489         return storage;
490 }
491
492 //Abort with a nice error message saying: 
493 // Our <field> was previously set to <cur>; it cannot now be <operation> <next>
494 void FEM_Attribute::bad(const char *field,bool forRead,int cur,int next,const char *caller) const
495 {
496         char nameStorage[256];
497         const char *name=FEM_Get_attr_name(attr,nameStorage);
498         char errBuf[1024];
499         const char *cannotBe=NULL;
500         if (forRead) {
501                 if (cur==-1) {
502                         sprintf(errBuf,"The %s %s %s was never set-- it cannot now be read",
503                                 e->getName(),name,field);
504                 }
505                 else /* already had a value */
506                         cannotBe="read as";
507         }
508         else /* for write */ {
509                 cannotBe="set to";
510         }
511         if (cannotBe!=NULL) /* Use standard ... <something> cannot be <something>'d... error message */
512                 sprintf(errBuf,"The %s %s %s was previously set to %d; it cannot now be %s %d",
513                         e->getName(),name,field,cur,cannotBe,next);
514         
515         FEM_Abort(caller,errBuf);
516 }
517
518
519 FEM_Attribute::FEM_Attribute(FEM_Entity *e_,int attr_)
520                 :e(e_),ghost(0),attr(attr_),width(0),datatype(-1), allocated(false)
521 {
522         tryAllocate();
523         if (femVersion == 0) width=-1;
524 }
525 void FEM_Attribute::pup(PUP::er &p) {
526         // e, attr, and ghost are always set by the constructor
527         p|width;
528         if (p.isUnpacking() && femVersion > 0 && width<0)  width=0;
529         p|datatype;
530         if (p.isUnpacking()) tryAllocate();
531 }
532 void FEM_Attribute::pupSingle(PUP::er &p, int pupindx) {
533   return;
534 }
535 FEM_Attribute::~FEM_Attribute() {}
536
537 void FEM_Attribute::setLength(int next,const char *caller) {
538         int cur=getLength();
539         if (next==cur) return; //Already set--nothing to do 
540         if (cur>0) bad("length",false,cur,next, caller);
541         e->setLength(next);
542         tryAllocate();
543 }
544         
545 void FEM_Attribute::setWidth(int next,const char *caller) {
546         int cur=getWidth();
547         if (next==cur) return; //Already set--nothing to do 
548         if (cur>0) bad("width",false,cur,next, caller);
549         width=next;
550         tryAllocate();
551         if (ghost) ghost->setWidth(width,caller);
552 }
553
554 void FEM_Attribute::setDatatype(int next,const char *caller) {
555         int cur=getDatatype();
556         if (next==cur) return; //Already set--nothing to do 
557         if (cur!=-1) bad("datatype",false,cur,next, caller);
558         datatype=next;
559         tryAllocate();
560         if (ghost) ghost->setDatatype(datatype,caller);
561 }
562
563 void FEM_Attribute::copyShape(const FEM_Attribute &src) {
564         setWidth(src.getWidth());
565         if (src.getDatatype()!=-1)
566           setDatatype(src.getDatatype()); //Automatically calls tryAllocate
567 }
568 void FEM_Attribute::set(const void *src, int firstItem,int length, 
569                 const IDXL_Layout &layout, const char *caller) 
570 {
571         if (firstItem!=0) { /* If this isn't the start... */
572                 if (length!=1) /* And we're not setting one at a time */
573                         CmiAbort("FEM_Mesh_data: unexpected firstItem");
574         }
575
576         if (femVersion == 0 && getRealLength() == -1) setLength(length);
577         else if (getLength()==0) setLength(length);
578         else if (length!=1 && length!=getLength()) 
579                 bad("length",false,getLength(),length, caller);
580         
581         int width=layout.width;
582         if (femVersion==0 && getRealWidth()==-1) setWidth(width);
583         else if (getWidth()==0) setWidth(width);
584         else if (width!=getWidth()) 
585                 bad("width",false,getWidth(),width, caller);
586         
587         int datatype=layout.type;
588         if (getDatatype()==-1) setDatatype(datatype);
589         else if (datatype!=getDatatype()) 
590                 bad("datatype",false,getDatatype(),datatype, caller);
591         
592         /* Assert: our storage should be allocated now.
593            Our subclass will actually copy user data */
594 }
595
596 void FEM_Attribute::get(void *dest, int firstItem,int length, 
597                 const IDXL_Layout &layout, const char *caller)  const
598 {
599         if (length==0) return; //Nothing to get
600         if (length!=1 && length!=getLength()) 
601                 bad("length",true,getLength(),length, caller);
602         
603         int width=layout.width;
604         if (width!=getWidth()) 
605                 bad("width",true,getWidth(),width, caller);
606         
607         int datatype=layout.type;
608         if (datatype!=getDatatype()) 
609                 bad("datatype",true,getDatatype(),datatype, caller);
610         
611         /* our subclass will actually copy into user data */
612 }
613
614 /*check if the layout is the same as earlier */
615
616 void FEM_Attribute::register_data(void *user, int length,int max,
617         const IDXL_Layout &layout, const char *caller){
618         
619                 int width=layout.width;
620                 if (femVersion == 0 && getRealWidth()==-1) setWidth(width);
621                 else if (getWidth()==0){
622                         setWidth(width);
623                 }else{
624                         if (width!=getWidth()){
625                                 bad("width",false,getWidth(),width, caller);
626                         }
627                 }       
628         
629                 int datatype=layout.type;
630                 if (getDatatype()==-1){
631                         setDatatype(datatype);
632                 }else{
633                         if (datatype!=getDatatype()){ 
634                                 bad("datatype",false,getDatatype(),datatype, caller);
635                         }
636                 }       
637                 
638 }
639
640 //Check if all three of length, width, and datatype are set.
641 // If so, call allocate.
642 void FEM_Attribute::tryAllocate(void) {
643         int lenNull, widthNull;
644         if (femVersion == 0) {
645           // version 0 takes -1 as empty
646           lenNull = (getRealLength()==-1);
647           widthNull = (getRealWidth()==-1);
648         }
649         else {
650           lenNull = (getLength()==0);
651           widthNull = (getWidth()==0);
652         }
653         if ((!allocated) && !lenNull && !widthNull && getDatatype()!=-1) {
654           allocated=true;
655                 allocate(getMax(),getWidth(),getDatatype());
656         }
657 }
658
659 /*********************** DataAttribute *******************/
660 FEM_DataAttribute::FEM_DataAttribute(FEM_Entity *e,int myAttr)
661         :FEM_Attribute(e,myAttr), 
662          char_data(0),int_data(0),float_data(0),double_data(0)
663 {
664 }
665 void FEM_DataAttribute::pup(PUP::er &p) {
666         super::pup(p);
667         switch(getDatatype()) {
668         case -1: /* not allocated yet */ break;
669         case FEM_BYTE:   
670           if (char_data) {
671             /*if(p.isSizing()) {
672               char_data->setRowLen(getEntity()->size());
673               }*/
674             char_data->pup(p);
675           } 
676           break;
677         case FEM_INT:    
678           if (int_data) {
679             /*if(p.isSizing()) {
680               int_data->setRowLen(getEntity()->size());
681               }*/
682             int_data->pup(p);
683           }
684           break;
685         case FEM_FLOAT:  
686           if (float_data) {
687             /*if(p.isSizing()) {
688               float_data->setRowLen(getEntity()->size());
689               }*/
690             float_data->pup(p);
691           }
692           break;
693         case FEM_DOUBLE: 
694           if (double_data) {
695             /*if(p.isSizing()) {
696               double_data->setRowLen(getEntity()->size());
697               }*/
698             double_data->pup(p);
699           }
700           break;
701         default: CkAbort("Invalid datatype in FEM_DataAttribute::pup");
702         }
703 }
704 void FEM_DataAttribute::pupSingle(PUP::er &p, int pupindx) {
705         super::pupSingle(p,pupindx);
706         switch(getDatatype()) {
707         case -1: /* not allocated yet */ break;
708         case FEM_BYTE:   if (char_data) char_data->pupSingle(p,pupindx); break;
709         case FEM_INT:    if (int_data) int_data->pupSingle(p,pupindx); break;
710         case FEM_FLOAT:  if (float_data) float_data->pupSingle(p,pupindx); break;
711         case FEM_DOUBLE: if (double_data) double_data->pupSingle(p,pupindx); break;
712         default: CkAbort("Invalid datatype in FEM_DataAttribute::pupSingle");
713         }
714 }
715 FEM_DataAttribute::~FEM_DataAttribute() {
716         if (char_data) delete char_data;
717         if (int_data) delete int_data;
718         if (float_data) delete float_data;
719         if (double_data) delete double_data;
720         
721 }
722
723 /// Copy this data out of the user's (layout-formatted) array:
724 template <class T>
725 inline void setTableData(const void *user, int firstItem, int length, 
726         IDXL_LAYOUT_PARAM, AllocTable2d<T> *table) 
727 {
728         for (int r=0;r<length;r++) {
729                 register T *tableRow=table->getRow(firstItem+r);
730                 for (int c=0;c<width;c++)
731                         tableRow[c]=IDXL_LAYOUT_DEREF(T,user,r,c);
732         }
733 }
734
735 /// Copy this data into the user's (layout-formatted) array:
736 template <class T>
737 inline void getTableData(void *user, int firstItem, int length, 
738         IDXL_LAYOUT_PARAM, const AllocTable2d<T> *table) 
739 {
740         for (int r=0;r<length;r++) {
741                 register const T *tableRow=table->getRow(firstItem+r);
742                 for (int c=0;c<width;c++)
743                         IDXL_LAYOUT_DEREF(T,user,r,c)=tableRow[c];
744         }
745
746 }
747
748 void FEM_DataAttribute::set(const void *u, int f,int l, 
749                 const IDXL_Layout &layout, const char *caller)
750 {
751         super::set(u,f,l,layout,caller);
752         switch(getDatatype()) {
753         case FEM_BYTE:  setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),char_data); break;
754         case FEM_INT: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),int_data); break;
755         case FEM_FLOAT: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),float_data); break;
756         case FEM_DOUBLE: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),double_data); break;
757         }
758 }
759         
760 void FEM_DataAttribute::get(void *u, int f,int l,
761                 const IDXL_Layout &layout, const char *caller) const
762 {
763         super::get(u,f,l,layout,caller);
764         switch(getDatatype()) {
765         case FEM_BYTE:  getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),char_data); break;
766         case FEM_INT: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),int_data); break;
767         case FEM_FLOAT: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),float_data); break;
768         case FEM_DOUBLE: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),double_data); break;
769         }
770 }
771
772 void FEM_DataAttribute::register_data(void *u,int l,int max,
773             const IDXL_Layout &layout, const char *caller)
774 {
775         super::register_data(u,l,max,layout,caller);
776         switch(getDatatype()){
777                 case FEM_BYTE: char_data->register_data((unsigned char *)u,l,max); break;
778                 case FEM_INT:    int_data->register_data((int *)u,l,max); break;
779                 case FEM_FLOAT: float_data->register_data((float *)u,l,max);break;
780                 case FEM_DOUBLE: double_data->register_data((double *)u,l,max);break;
781         }
782 }
783
784 template<class T>
785 inline AllocTable2d<T> *allocTablePtr(AllocTable2d<T> *src,int len,int wid) {
786         if (src==NULL) src=new AllocTable2d<T>;
787         src->allocate(wid,len);
788         return src;
789 }
790 void FEM_DataAttribute::allocate(int l,int w,int datatype)
791 {
792         switch(datatype) {
793         case FEM_BYTE:  char_data=allocTablePtr(char_data,l,w); break;
794         case FEM_INT: int_data=allocTablePtr(int_data,l,w); break;
795         case FEM_FLOAT: float_data=allocTablePtr(float_data,l,w); break;
796         case FEM_DOUBLE: double_data=allocTablePtr(double_data,l,w); break;
797         default: CkAbort("Invalid datatype in FEM_DataAttribute::allocate");
798         };
799 }
800 void FEM_DataAttribute::copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity)
801 {
802         const FEM_DataAttribute *dsrc=(const FEM_DataAttribute *)&src;
803         switch(getDatatype()) {
804         case FEM_BYTE:  char_data->setRow(dstEntity,dsrc->char_data->getRow(srcEntity)); break;
805         case FEM_INT: 
806                         int_data->setRow(dstEntity,dsrc->int_data->getRow(srcEntity)); break;
807         case FEM_FLOAT: float_data->setRow(dstEntity,dsrc->float_data->getRow(srcEntity)); break;
808         case FEM_DOUBLE: double_data->setRow(dstEntity,dsrc->double_data->getRow(srcEntity)); break;
809         }
810 }
811
812 template<class T>
813 inline void interpolateAttrs(AllocTable2d<T> *data,int A,int B,int D,double frac,int width){
814   T *rowA = data->getRow(A);
815   T *rowB = data->getRow(B);
816   T *rowD = data->getRow(D);
817   for(int i=0;i<width;i++){
818     double val = (double )rowA[i];
819     val *= (frac);
820     val += (1-frac) *((double )rowB[i]);
821     rowD[i] = (T )val;
822   }
823 }
824
825 template<class T>
826 inline void interpolateAttrs(AllocTable2d<T> *data,int *iNodes,int rNode,int k,int width){
827   T *row[8];
828   for (int i=0; i<k; i++) {
829     row[i] = data->getRow(iNodes[i]);
830   }
831   T *rowR = data->getRow(rNode);
832   for(int i=0;i<width;i++){
833     double val = 0.0;
834     for (int j=0; j<k; j++) {
835       val += (double)row[j][i];
836     }
837     val = val/k;
838     rowR[i] = (T )val;
839   }
840 }
841
842 template<class T>
843 inline void minAttrs(AllocTable2d<T> *data,int A,int B,int D,double frac,int width){
844   T *rowA = data->getRow(A);
845   T *rowB = data->getRow(B);
846   T *rowD = data->getRow(D);
847   for(int i=0;i<width;i++){
848     if(rowA[i] < rowB[i]){
849       rowD[i] = rowA[i];
850     }else{
851       rowD[i] = rowB[i];
852     }
853   }
854 }
855
856 template<class T>
857 inline void minAttrs(AllocTable2d<T> *data,int *iNodes,int rNode,int k,int width){
858   T *row[8];
859   for (int i=0; i<k; i++) {
860     row[i] = data->getRow(iNodes[i]);
861   }
862   T *rowR = data->getRow(rNode);
863   for(int i=0;i<width;i++){
864     rowR[i] = row[0][i]; 
865     for (int j=1; j<k; j++) {
866       if (row[j][i] < rowR[i]) {
867         rowR[i] = row[j][i];
868       }
869     }
870   }
871 }
872
873 void FEM_DataAttribute::interpolate(int A,int B,int D,double frac){
874   switch(getDatatype()){
875   case FEM_BYTE:
876     minAttrs(char_data,A,B,D,frac,getWidth());          
877     break;
878   case FEM_INT:
879     minAttrs(int_data,A,B,D,frac,getWidth());           
880     break;
881   case FEM_FLOAT:
882     interpolateAttrs(float_data,A,B,D,frac,getWidth());         
883     break;
884   case FEM_DOUBLE:
885     interpolateAttrs(double_data,A,B,D,frac,getWidth());                
886     break;
887   }
888 }
889
890 void FEM_DataAttribute::interpolate(int *iNodes,int rNode,int k){
891   switch(getDatatype()){
892   case FEM_BYTE:
893     minAttrs(char_data,iNodes,rNode,k,getWidth());              
894     break;
895   case FEM_INT:
896     minAttrs(int_data,iNodes,rNode,k,getWidth());               
897     break;
898   case FEM_FLOAT:
899     interpolateAttrs(float_data,iNodes,rNode,k,getWidth());             
900     break;
901   case FEM_DOUBLE:
902     interpolateAttrs(double_data,iNodes,rNode,k,getWidth());            
903     break;
904   }
905 }
906
907 /*********************** FEM_IndexAttribute *******************/
908 FEM_IndexAttribute::Checker::~Checker() {}
909
910 FEM_IndexAttribute::FEM_IndexAttribute(FEM_Entity *e,int myAttr,FEM_IndexAttribute::Checker *checker_)
911         :FEM_Attribute(e,myAttr), idx(0,0,-1), checker(checker_)
912 {
913         setDatatype(FEM_INT);
914 }
915 void FEM_IndexAttribute::pup(PUP::er &p) {
916         super::pup(p);
917         /*if(p.isSizing()) {
918           idx.setRowLen(getEntity()->size());
919           }*/
920         p|idx;
921 }
922 void FEM_IndexAttribute::pupSingle(PUP::er &p, int pupindx) {
923         super::pupSingle(p,pupindx);
924         idx.pupSingle(p,pupindx);
925 }
926 FEM_IndexAttribute::~FEM_IndexAttribute() {
927         if (checker) delete checker;
928 }
929
930 void FEM_IndexAttribute::allocate(int length,int width,int datatype)
931 {
932         idx.allocate(width,length);
933 }
934
935 /**
936  * Convert a datatype, which must be FEM_INDEX_0 or FEM_INDEX_1, to 
937  * the first valid index (base index) of that type.  Otherwise 
938  * call FEM_Abort, because the datatype is wrong.
939  */
940 static int type2base(int base_type,const char *caller) {
941         if (base_type==FEM_INDEX_0) return 0;
942         if (base_type==FEM_INDEX_1) return 1;
943         FEM_Abort(caller,"You must use the datatype FEM_INDEX_0 or FEM_INDEX_1 with FEM_CONN, not %d",
944                 base_type);
945         return 0; //< for whining compilers
946 }
947
948 /// Copy this data out of the user's (layout-formatted, indexBase) array:
949 void setIndexTableData(const void *user, int firstItem, int length, 
950         IDXL_LAYOUT_PARAM, AllocTable2d<int> *table,int indexBase) 
951 {
952         for (int r=0;r<length;r++) {
953                 register int *tableRow=table->getRow(firstItem+r);
954                 for (int c=0;c<width;c++)
955                         tableRow[c]=IDXL_LAYOUT_DEREF(int,user,r,c)-indexBase;
956         }
957 }
958
959 /// Copy this data into the user's (layout-formatted, indexBase) array:
960 void getIndexTableData(void *user, int firstItem, int length, 
961         IDXL_LAYOUT_PARAM, const AllocTable2d<int> *table,int indexBase) 
962 {
963         for (int r=0;r<length;r++) {
964                 register const int *tableRow=table->getRow(firstItem+r);
965                 for (int c=0;c<width;c++)
966                         IDXL_LAYOUT_DEREF(int,user,r,c)=tableRow[c]+indexBase;
967         }
968 }
969
970 void FEM_IndexAttribute::set(const void *src, int firstItem,int length,
971                 const IDXL_Layout &layout,const char *caller)
972 {
973         IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
974         super::set(src,firstItem,length,lo,caller);
975         
976         int indexBase=type2base(layout.type,caller);
977         setIndexTableData(src,firstItem,length,IDXL_LAYOUT_CALL(layout),&idx,indexBase);
978         
979         if (checker) 
980                 for (int r=0;r<length;r++)
981                         checker->check(firstItem+r,idx,caller);
982 }
983
984 void FEM_IndexAttribute::get(void *dest, int firstItem,int length, 
985                 const IDXL_Layout &layout,const char *caller) const
986 {
987         IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
988         super::get(dest,firstItem,length,lo,caller);
989         
990         int indexBase=type2base(layout.type,caller);
991         getIndexTableData(dest,firstItem,length,IDXL_LAYOUT_CALL(layout),&idx,indexBase);
992 }
993
994 void FEM_IndexAttribute::register_data(void *user, int length,int max,
995                 const IDXL_Layout &layout, const char *caller){
996         IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
997         super::register_data(user,length,max,lo,caller);
998
999         idx.register_data((int *)user,length,max);
1000 }
1001
1002 void FEM_IndexAttribute::copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity)
1003 {
1004         const FEM_IndexAttribute *csrc=(const FEM_IndexAttribute *)&src;
1005         idx.setRow(dstEntity,csrc->idx.getRow(srcEntity));
1006 }
1007
1008 /*************FEM_VarIndexAttribute***************/
1009
1010 FEM_VarIndexAttribute::FEM_VarIndexAttribute(FEM_Entity *e,int myAttr)
1011         :FEM_Attribute(e,myAttr)
1012 {
1013   oldlength = 0;
1014         allocate(getMax(),getWidth(),getDatatype());
1015         setDatatype(FEM_INT);
1016 };
1017
1018 void FEM_VarIndexAttribute::pup(PUP::er &p){
1019         super::pup(p);
1020         p | idx;
1021         p|oldlength;
1022 };
1023
1024 void FEM_VarIndexAttribute::pupSingle(PUP::er &p, int pupindx){
1025         super::pupSingle(p,pupindx);
1026         p|idx[pupindx];
1027 };
1028
1029 void FEM_VarIndexAttribute::set(const void *src,int firstItem,int length,
1030                 const IDXL_Layout &layout,const char *caller){
1031                 printf("set not yet implemented for FEM_VarIndexAttribute \n");
1032 }
1033
1034 void FEM_VarIndexAttribute::get(void *dest, int firstItem,int length,
1035                 const IDXL_Layout &layout, const char *caller) const{
1036          printf("get not yet implemented for FEM_VarIndexAttribute \n");
1037                         
1038 }
1039
1040 void FEM_VarIndexAttribute::copyEntity(int dstEntity,const FEM_Attribute &_src,int srcEntity){
1041         FEM_VarIndexAttribute &src = (FEM_VarIndexAttribute &)_src;
1042         const CkVec<CkVec<ID> > &srcTable = src.get();
1043         CkVec<ID> temp = srcTable[srcEntity];
1044         idx.insert(dstEntity, temp);
1045 }
1046
1047 void FEM_VarIndexAttribute::print(){
1048         for(int i=0;i<idx.size();i++){
1049                 printf("%d -> ",i);
1050                 for(int j=0;j<idx[i].size();j++){
1051                         printf("(%d %d) ",((idx[i])[j]).type,((idx[i])[j]).id);
1052                 }
1053                 printf("\n");
1054         }
1055 };
1056
1057 int FEM_VarIndexAttribute::findInRow(int row,const ID &data){
1058         if(row >= idx.length()){
1059                 return -1;
1060         }
1061         CkVec<ID> &rowVec = idx[row];
1062         for(int i=0;i<rowVec.length();i++){
1063                 if(data == rowVec[i]){
1064                         return i;
1065                 }
1066         }
1067         return -1;
1068 }
1069
1070 /********************** Entity **************************/
1071
1072 /// Return the human-readable version of this entity code.
1073 CDECL const char *FEM_Get_entity_name(int entity,char *storage) 
1074 {
1075         char *dest=storage;
1076         if (entity<FEM_ENTITY_FIRST || entity>=FEM_ENTITY_LAST) {
1077                 sprintf(dest,"unknown entity code (%d)",entity);
1078         }
1079         else {
1080                 if (entity>FEM_ENTITY_FIRST+FEM_GHOST) {
1081                         sprintf(dest,"FEM_GHOST+");
1082                         dest+=strlen(dest); /* we want "FEM_GHOST+foo" */
1083                         entity-=FEM_GHOST;
1084                 }
1085                 if (entity==FEM_NODE)
1086                         sprintf(dest,"FEM_NODE");
1087                 else if (entity>=FEM_SPARSE)
1088                         sprintf(dest,"FEM_SPARSE+%d",entity-FEM_SPARSE);
1089                 else /* entity>=FEM_ELEM */
1090                         sprintf(dest,"FEM_ELEM+%d",entity-FEM_ELEM);
1091         }
1092         return storage;
1093 }
1094
1095 FEM_Entity::FEM_Entity(FEM_Entity *ghost_) //Default constructor
1096   :length(0), max(0),ghost(ghost_), coord(0), sym(0), globalno(0), valid(0), meshSizing(0),
1097    ghostIDXL(ghost?&ghostSend:NULL, ghost?&ghost->ghostRecv:NULL),resize(NULL),
1098    invalidList(NULL),first_invalid(0),last_invalid(0),invalidListLen(0),invalidListAllLen(0)
1099 {
1100         //No attributes initially
1101         if (femVersion == 0) {
1102                 length=-1;
1103                 max=-1;
1104         }
1105
1106 void FEM_Entity::pup(PUP::er &p) {
1107         p|length;
1108         p|max;
1109         if (p.isUnpacking() && femVersion > 0 && length<0)  length=0;
1110         p.comment(" Ghosts to send out: ");
1111         ghostSend.pup(p);
1112         p.comment(" Ghosts to recv: ");
1113         ghostRecv.pup(p);
1114         p.comment(" Ghost IDXL tag: ");
1115         ghostIDXL.pup(p);
1116         
1117         int nAttributes=attributes.size();
1118         p|nAttributes;
1119         for (int a=0;a<nAttributes;a++) 
1120         {
1121         /* Beautiful hack: the FEM_Attribute objects are normally hideously cross-linked
1122            with their owning classes.  Thus instead of trying to rebuild the FEM_Attributes
1123            from scratch here, we just use the existing "lookup" method to demand-create those
1124            that need it, or just *find* those that are already there.
1125            
1126            This is a much better fit for this situation than using the general PUP::able.
1127          */
1128           int attr=0; // The attribute type we're pupping
1129           FEM_Attribute *r=NULL;
1130           if (!p.isUnpacking()) { //Send side: we already know the source
1131             r=attributes[a];
1132             attr=r->getAttr();
1133           }
1134           p|attr;
1135           if (p.isUnpacking()) { //Recv side: create (or recycle) the destination
1136             r=lookup(attr,"FEM_Entity::pup");
1137           }
1138           
1139           { //Put the human-readable attribute name in the output file:
1140             char attrNameStorage[256];
1141             p.comment(FEM_Get_attr_name(attr,attrNameStorage));
1142           }
1143           r->pup(p);
1144         }
1145         p|first_invalid;
1146         p|last_invalid;
1147         p|invalidListLen;
1148         p|invalidListAllLen;
1149         if(p.isUnpacking()) {
1150           if(invalidList!=NULL) delete[] invalidList; //was just created in allocateValid
1151           if(invalidListAllLen>0) invalidList = new int[invalidListAllLen];
1152         }
1153         if(invalidListAllLen>0) {
1154           PUParray(p, (int*)invalidList, invalidListAllLen);
1155         }
1156         if(p.isDeleting()) {
1157           if(invalidList!=NULL) delete[] invalidList;
1158         }
1159         
1160         if (ghost!=NULL) {
1161                 p.comment(" ---- Ghost attributes ---- ");
1162                 ghost->pup(p);
1163         }
1164 }
1165
1166 FEM_Entity::~FEM_Entity() 
1167 {
1168         delete ghost;
1169         for (int a=0;a<attributes.size();a++)
1170                 delete attributes[a];
1171 }
1172
1173 /// Copy our attributes' widths and data types from this entity.
1174 void FEM_Entity::copyShape(const FEM_Entity &src) {
1175         for (int a=0;a<src.attributes.size();a++) 
1176         { // We need each of his attributes:
1177                 const FEM_Attribute *Asrc=src.attributes[a];
1178                 FEM_Attribute *Adst=lookup(Asrc->getAttr(),"FEM_Entity::copyShape");
1179                 Adst->copyShape(*Asrc);
1180         }
1181         if (ghost) ghost->copyShape(*src.ghost);
1182 }
1183
1184 void FEM_Entity::setLength(int newlen, bool f) 
1185 {
1186     if (!resize) {
1187         if (size() != newlen) {
1188             length = newlen;
1189             // Each of our attributes need to be expanded for our new length:
1190             for (int a=0; a<attributes.size(); a++) {
1191                 attributes[a]->reallocate();
1192             }
1193         }
1194     }
1195     else {
1196       int otherlen = length;
1197       if(!f) {
1198         otherlen = newlen;
1199       }
1200       length = newlen;
1201         if (length > max) {
1202             if (max > 4) {
1203                 max = max + (max >> 2);
1204             } else {
1205                 max = max + 10;
1206             }
1207             for (int a=0;a<attributes.size();a++){
1208                 int code = attributes[a]->getAttr();
1209                 if (!(code <= FEM_ATTRIB_TAG_MAX || 
1210                             code == FEM_CONN || code == FEM_COORD ||
1211                       code == FEM_BOUNDARY)) { //user can store bpundary & connectivity also
1212                     attributes[a]->reallocate();
1213                 }
1214             }   
1215             // call resize with args max n;
1216             //CkPrintf("Resize called \n");
1217             resize(args,&otherlen,&max); //resets length to otherlen
1218             //resize(args,&length,&max);
1219             length = newlen;
1220         }
1221     }
1222 }
1223
1224 void FEM_Entity::allocateValid(void) {
1225   if (!valid){
1226     valid=new FEM_DataAttribute(this,FEM_IS_VALID_ATTR);
1227     add(valid);
1228     valid->setWidth(1); //Only 1 flag per node
1229     valid->setLength(size());
1230     valid->setDatatype(FEM_BYTE);
1231     valid->reallocate();
1232     
1233     // Set all to valid initially
1234     for(int i=0;i<size();i++) {
1235       valid->getChar()(i,0)=1;
1236     }
1237     if(true) { //maintains a list of invalid elements. FASTER
1238       invalidListLen=0; invalidListAllLen=16;
1239       invalidList = new int[invalidListAllLen]; 
1240       //its ok to waste a few bytes for an entity, especially when it can 
1241       //bring up the performance of the average case
1242       for(int i=0; i<invalidListAllLen; i++) {
1243         invalidList[i] = -1; //means this position is empty
1244       }
1245     }
1246     else {
1247       first_invalid = last_invalid = 0;
1248     }
1249   }
1250 }
1251
1252 void FEM_Entity::set_valid(int idx, bool isNode){
1253   if(true) { //maintains a list of invalid elements. FASTER
1254     CkAssert(idx < size() && idx >=0);
1255     valid->getChar()(idx,0)=1;
1256     if(invalidListLen>0) {
1257       if(invalidList[invalidListLen-1]==idx) {
1258         invalidListLen--; //pull out the last entry from the invalid list
1259       }
1260     }
1261   }
1262   else {
1263     int size1 = size();
1264     if(size1==0) {
1265       valid->getChar()(idx,0)=1;
1266     }
1267     else {
1268       CkAssert(idx < size1 && idx >=0 && first_invalid<=last_invalid);
1269       valid->getChar()(idx,0)=1;
1270       if(idx == first_invalid) {
1271         // Move first_invalid to the next invalid entry 
1272         while((first_invalid<last_invalid) && is_valid(first_invalid)){
1273           first_invalid++;
1274         }
1275       }
1276       else if(idx == last_invalid) {
1277         // Move last_invalid to the previous invalid entry      
1278         while((first_invalid<last_invalid) && is_valid(last_invalid)) {
1279           last_invalid--;
1280         }
1281       }
1282     }
1283     // If we have no invalid elements left, then put both pointers to 0
1284     if( first_invalid == last_invalid && is_valid(first_invalid) )
1285       first_invalid = last_invalid = 0;
1286   }
1287 }
1288
1289 void FEM_Entity::set_invalid(int idx, bool isNode){
1290   if(true) { //maintains a list of invalid elements. FASTER
1291     CkAssert(idx < size() && idx >=0);
1292     valid->getChar()(idx,0)=0;
1293     if(invalidListLen==invalidListAllLen) {
1294       invalidListAllLen = invalidListAllLen<<1;
1295       int* tmp = invalidList;
1296       invalidList = new int[invalidListAllLen];
1297       memcpy(invalidList,tmp,invalidListLen*sizeof(int));
1298       delete[] tmp;
1299     }
1300     invalidList[invalidListLen] = idx;
1301     invalidListLen++;
1302   }
1303   else {
1304     CkAssert(idx < size() && idx >=0 && first_invalid<=last_invalid);
1305     valid->getChar()(idx,0)=0;
1306     
1307     // If there are currently no invalid entities
1308     if(first_invalid==0 && last_invalid==0 && is_valid(0)){
1309       first_invalid = last_invalid = idx;
1310       return;
1311     }
1312     if(idx < first_invalid){
1313       first_invalid = idx;
1314     }
1315     if(idx > last_invalid){
1316       last_invalid = idx;
1317     }
1318     // TODO:
1319     // We should probably have an algorithm for shrinking the entire attribute 
1320     // array if we invalidate the final element. In this case we should scan backwards
1321     // to find the largest indexed valid entity and resize down to it.
1322     // 
1323     // It may be necessary to modify the idxl lists if we do this type of shrinking.
1324     // Someone needs to confirm whether that is necessary. If not, then it should be 
1325     // simple to allow shrinking of the number of nodes or elements.
1326   }
1327 }
1328
1329 int FEM_Entity::is_valid(int idx){
1330     CkAssert(idx < size() && idx >=0);
1331     return valid->getChar()(idx,0);
1332 }
1333
1334 int FEM_Entity::is_valid_any_idx(int idx){
1335   if(idx == -1 || idx >= size())
1336         return false;
1337   else if(idx < 0)
1338         return ghost->is_valid_any_idx(FEM_To_ghost_index(idx));
1339   else 
1340     return valid->getChar()(idx,0);
1341 }
1342
1343 int FEM_Entity::count_valid(){
1344   int count=0;
1345   for(int i=0;i<size();i++)
1346         if(is_valid(i)) count++;
1347   return count;
1348 }
1349
1350 /// Get an entry(entity index) that corresponds to an invalid entity
1351 /// The invalid slot in the tables can then be reused when "creating" a new element or node
1352 /// We either return an empty slot, or resize the array and return a value at the end
1353 /// If someone has a better name for this function, please change it.
1354 int FEM_Entity::get_next_invalid(FEM_Mesh *m, bool isNode, bool isGhost){
1355   int retval=0;
1356 //  if(true) { //maintains a list of invalid elements. FASTER
1357     if(invalidListLen>0) {
1358       retval = invalidList[invalidListLen-1];
1359     }
1360     else {
1361       retval = size();
1362       setLength(retval+1,true);  
1363     }
1364 //  }
1365 /*  else {
1366     int size1 = size();
1367     if(size1==0) { //special case because for size=0, first_invalid=last_invalid=0
1368       retval= 0;
1369       setLength(1);
1370     }
1371     else {
1372       CkAssert(!is_valid(first_invalid) || first_invalid==0);
1373       // if we have an invalid entity to return
1374       bool flag1 = false;
1375       if(!is_valid(first_invalid)){
1376         retval = first_invalid;
1377         if(isNode && !isGhost) { //it is a node & the entity is not a ghost entity
1378           while(retval <= last_invalid) {
1379             if(!is_valid(retval)) {
1380               if(m->getfmMM()->fmLockN[retval].haslocks()) {
1381                 retval++;
1382               }
1383               else if(hasConn(retval)) { //has some connectivity
1384                 retval++;
1385               }
1386               else {
1387                 flag1 = true;
1388                 break;
1389               }
1390             }
1391             else retval++;
1392           }
1393         }
1394         else if(isNode) {
1395           while(retval <= last_invalid) {
1396             if(!is_valid(retval)) {
1397               if(hasConn(retval)) { //has some connectivity
1398                 retval++;
1399               }
1400               else {
1401                 flag1 = true;
1402                 break;
1403               }
1404             }
1405             else retval++;
1406           }
1407         }  
1408         else{
1409           // resize array and return new entity
1410           flag1 = true;
1411         }
1412       }
1413       if(!flag1) {
1414         retval = size1;
1415         setLength(retval+1);
1416       }
1417     }
1418   }*/
1419   set_valid(retval,isNode);
1420   return retval;
1421 }
1422
1423 void FEM_Entity::setMaxLength(int newLen,int newMaxLen,void *pargs,FEM_Mesh_alloc_fn fn){
1424   //CkPrintf("resize fn %p \n",fn);
1425         max = newMaxLen;
1426         resize = fn;
1427         args = pargs;
1428         setLength(newLen);
1429 };
1430
1431
1432 /// Copy src[srcEntity] into our dstEntity.
1433 void FEM_Entity::copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity) {
1434         FEM_Entity *dp=this; //Destination entity
1435         const FEM_Entity *sp=&src;
1436         for (int a=0;a<sp->attributes.size();a++) 
1437         { //We copy each of his attributes:
1438                 const FEM_Attribute *Asrc=sp->attributes[a];
1439                 FEM_Attribute *Adst=dp->lookup(Asrc->getAttr(),"FEM_Entity::copyEntity");
1440                 Adst->copyEntity(dstEntity,*Asrc,srcEntity);
1441         }
1442 }
1443
1444 /// Add room for one more entity, with initial values from src[srcEntity],
1445 /// and return the new entity's index.
1446 int FEM_Entity::push_back(const FEM_Entity &src,int srcEntity) {
1447         int dstEntity=size();
1448         setLength(dstEntity+1);
1449         copyEntity(dstEntity,src,srcEntity);
1450         return dstEntity;
1451 }
1452
1453 /// Add this attribute to this kind of Entity.
1454 /// This method is normally called by the default lookup method.
1455 void FEM_Entity::add(FEM_Attribute *attribute) {
1456         if (ghost!=NULL) 
1457         { //Look up (or create) the ghost attribute, too:
1458                 attribute->setGhost(ghost->lookup(attribute->getAttr(),"FEM_Entity::add"));
1459         }
1460         attributes.push_back(attribute);
1461 }
1462
1463 /**
1464  * Find this attribute (from an FEM_ATTR code) of this entity.
1465  * The default implementation searches the list of userdata attributes;
1466  * subclasses with other attributes should override this routine.
1467  */
1468 FEM_Attribute *FEM_Entity::lookup(int attr,const char *caller) {
1469         //Try to find an existing attribute (FIXME: keep attributes in a map, to speed this up)
1470         for (int a=0;a<attributes.size();a++) {
1471                 if (attributes[a]->getAttr()==attr)
1472                         return attributes[a];
1473         }
1474         
1475         //If we get here, no existing attribute fits the bill: create one
1476         create(attr,caller);
1477         
1478         // If create did its job, the next lookup should succeed:
1479         return lookup(attr,caller);
1480 }
1481
1482 /**
1483  * Create a new attribute from an FEM_ATTR code.
1484  * The default implementation handles FEM_DATA tags; entity-specific
1485  * attributes (like FEM_CONN) need to be overridden and created 
1486  * by subclasses.
1487  */
1488 void FEM_Entity::create(int attr,const char *caller) {
1489   if (attr<=FEM_ATTRIB_TAG_MAX) //It's a valid user data tag
1490         add(new FEM_DataAttribute(this,attr));
1491   else if (attr==FEM_COORD) 
1492         allocateCoord();
1493   else if (attr==FEM_SYMMETRIES) 
1494         allocateSym();
1495   else if (attr==FEM_GLOBALNO) 
1496         allocateGlobalno();
1497   else if (attr==FEM_IS_VALID_ATTR)
1498         allocateValid();
1499   else if (attr==FEM_MESH_SIZING) 
1500         allocateMeshSizing();
1501   else if(attr == FEM_CHUNK){
1502         FEM_IndexAttribute *chunkNo= new FEM_IndexAttribute(this,FEM_CHUNK,NULL);
1503         add(chunkNo);
1504         chunkNo->setWidth(1);
1505   }
1506   else if(attr == FEM_BOUNDARY){
1507         allocateBoundary();
1508   }else if(attr == FEM_ADAPT_ADJ){
1509                 FEM_DataAttribute *adaptAdj = new FEM_DataAttribute(this,attr);
1510                 adaptAdj->setDatatype(FEM_BYTE);
1511                 add(adaptAdj);
1512         }else if(attr == FEM_ADAPT_LOCK){
1513                 FEM_DataAttribute *adaptLock = new FEM_DataAttribute(this,attr);
1514                 adaptLock->setDatatype(FEM_INT);
1515                 adaptLock->setWidth(2);
1516                 add(adaptLock);
1517         }
1518   else {
1519         //It's an unrecognized tag: abort
1520         char attrNameStorage[256], msg[1024];
1521         sprintf(msg,"Could not locate the attribute %s for entity %s",
1522                         FEM_Get_attr_name(attr,attrNameStorage), getName());
1523         FEM_Abort(caller,msg);
1524   }
1525 }
1526
1527 void FEM_Entity::allocateCoord(void) {
1528         if (coord) CkAbort("FEM_Entity::allocateCoord called, but already allocated");
1529         coord=new FEM_DataAttribute(this,FEM_COORD);
1530         add(coord); // coord will be deleted by FEM_Entity now
1531         coord->setDatatype(FEM_DOUBLE);
1532 }
1533
1534 void FEM_Entity::allocateSym(void) {
1535         if (sym) CkAbort("FEM_Entity::allocateSym called, but already allocated");
1536         sym=new FEM_DataAttribute(this,FEM_SYMMETRIES);
1537         add(sym); // sym will be deleted via attributes list now
1538         sym->setWidth(1);
1539         sym->setDatatype(FEM_BYTE); //Same as FEM_Symmetries_t
1540 }
1541
1542 void FEM_Entity::setSymmetries(int r,FEM_Symmetries_t s)
1543 {
1544         if (!sym) {
1545                 if (s==0) return; //Don't bother allocating just for 0
1546                 allocateSym();
1547         }
1548         sym->getChar()(r,0)=s;
1549 }
1550
1551
1552 /*
1553  * Set the coordinates for a node or other entity.
1554  * Use the appropriate 2d or 3d version.
1555  * 
1556  * Note that first the function attempts to find coordinates in
1557  * the FEM_COORD attribute's array "*coord". If this fails, then
1558  * it will use the user's FEM_DATA field. Little error checking is
1559  * done, so the functions may crash if used inappropriately.
1560  */
1561 inline void FEM_Entity::set_coord(int idx, double x, double y){
1562   if(coord){
1563         coord->getDouble()(idx,0)=x;
1564         coord->getDouble()(idx,1)=y;
1565   }
1566   else {
1567         FEM_DataAttribute* attr =       (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
1568         attr->getDouble()(idx,0)=x;
1569         attr->getDouble()(idx,1)=y;
1570   }
1571 }
1572
1573 inline void FEM_Entity::set_coord(int idx, double x, double y, double z){
1574   if(coord){
1575         coord->getDouble()(idx,0)=x;
1576         coord->getDouble()(idx,1)=y;
1577         coord->getDouble()(idx,2)=z;
1578   }
1579   else {
1580         FEM_DataAttribute* attr =       (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
1581         attr->getDouble()(idx,0)=x;
1582         attr->getDouble()(idx,1)=y;
1583         attr->getDouble()(idx,2)=z;
1584   }
1585 }
1586
1587
1588 void FEM_Entity::allocateGlobalno(void) {
1589         if (globalno) CkAbort("FEM_Entity::allocateGlobalno called, but already allocated");
1590         globalno=new FEM_IndexAttribute(this,FEM_GLOBALNO,NULL);
1591         add(globalno); // globalno will be deleted via attributes list now
1592         globalno->setWidth(1);
1593 }
1594
1595 void FEM_Entity::allocateMeshSizing(void) {
1596   if (meshSizing) 
1597     CkAbort("FEM_Entity::allocateMeshSizing called, but already allocated");
1598   meshSizing=new FEM_DataAttribute(this,FEM_MESH_SIZING);
1599   add(meshSizing); // meshSizing will be deleted via attributes list now
1600   meshSizing->setWidth(1);
1601   meshSizing->setDatatype(FEM_DOUBLE);
1602 }
1603
1604 double FEM_Entity::getMeshSizing(int r) {
1605   if (!meshSizing) {
1606     allocateMeshSizing();
1607     return -1.0;
1608   }
1609   if (r >= 0)  return meshSizing->getDouble()(r,0);
1610   else  return ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0);
1611 }
1612
1613 void FEM_Entity::setMeshSizing(int r,double s)
1614 {
1615   if (!meshSizing) allocateMeshSizing();
1616   if (s <= 0.0) return;
1617   if (r >= 0)  meshSizing->getDouble()(r,0)=s;
1618   else ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0)=s;
1619 }
1620
1621 void FEM_Entity::setMeshSizing(double *sf)
1622 {
1623   if (!meshSizing) allocateMeshSizing();
1624   int len = size();
1625   for (int i=0; i<len; i++)
1626     meshSizing->getDouble()(i,0)=sf[i];
1627 }
1628
1629 void FEM_Entity::allocateBoundary(){
1630         FEM_DataAttribute *bound = new FEM_DataAttribute(this,FEM_BOUNDARY);
1631         add(bound);
1632         bound->setWidth(1);
1633         bound->setDatatype(FEM_INT);
1634 }
1635
1636 void FEM_Entity::setGlobalno(int r,int g) {
1637         if (!globalno) allocateGlobalno();
1638         globalno->get()(r,0)=g;
1639 }
1640 void FEM_Entity::setAscendingGlobalno(void) {
1641         if (!globalno) {
1642                 allocateGlobalno();
1643                 int len=size();
1644                 for (int i=0;i<len;i++) globalno->get()(i,0)=i;
1645         }
1646 }
1647 void FEM_Entity::setAscendingGlobalno(int base) {
1648         if (!globalno) {
1649                 allocateGlobalno();
1650                 int len=size();
1651                 for (int i=0;i<len;i++) globalno->get()(i,0)=i+base;
1652         }
1653 }
1654 void FEM_Entity::copyOldGlobalno(const FEM_Entity &e) {
1655         if ((!hasGlobalno()) && e.hasGlobalno() && size()>=e.size()) {
1656                 for (int i=0;i<size();i++) 
1657                         setGlobalno(i,e.getGlobalno(i));
1658         }
1659 }
1660
1661 /*
1662  * method to clear ghosts associated with this entity
1663  * */
1664
1665 void FEM_Entity::clearGhost(){
1666         CkAssert(ghost != NULL);
1667         ghostSend.clear();
1668         ghost->setLength(0);
1669         ghost->ghostRecv.clear();
1670 };
1671
1672 /********************** Node *****************/
1673 FEM_Node::FEM_Node(FEM_Node *ghost_) 
1674   :FEM_Entity(ghost_), primary(0), sharedIDXL(&shared,&shared),
1675    elemAdjacency(0),nodeAdjacency(0)
1676 {}
1677
1678 void FEM_Node::allocatePrimary(void) {
1679   if (primary) CkAbort("FEM_Node::allocatePrimary called, but already allocated");
1680   primary=new FEM_DataAttribute(this,FEM_NODE_PRIMARY);
1681   add(primary); // primary will be deleted by FEM_Entity now
1682   primary->setWidth(1); //Only 1 flag per node
1683   primary->setDatatype(FEM_BYTE);
1684 }
1685
1686 bool FEM_Node::hasConn(int idx) {
1687   if((elemAdjacency->get()[idx].length() > 0)||(nodeAdjacency->get()[idx].length() > 0))
1688     return true;
1689   else return false;
1690 }
1691
1692 void FEM_Node::pup(PUP::er &p) {
1693         p.comment(" ---------------- Nodes ------------------ ");       
1694         super::pup(p);
1695         p.comment(" ---- Shared nodes ----- "); 
1696         shared.pup(p);
1697         p.comment(" shared nodes IDXL ");
1698         sharedIDXL.pup(p);
1699 }
1700 FEM_Node::~FEM_Node() {
1701 }
1702
1703
1704 const char *FEM_Node::getName(void) const {return "FEM_NODE";}
1705
1706 void FEM_Node::create(int attr,const char *caller) {
1707   if (attr==FEM_NODE_PRIMARY)
1708         allocatePrimary();
1709   else if(attr == FEM_NODE_ELEM_ADJACENCY)
1710         allocateElemAdjacency();
1711   else if(attr == FEM_NODE_NODE_ADJACENCY)
1712         allocateNodeAdjacency();
1713   else
1714         super::create(attr,caller);
1715 }
1716
1717
1718 /********************** Elem *****************/
1719 /// This checker verifies that FEM_Elem::conn's entries are valid node indices.
1720 class FEM_Elem_Conn_Checker : public FEM_IndexAttribute::Checker {
1721         const FEM_Entity &sizeSrc;
1722         const FEM_Entity *sizeSrc2;
1723 public:
1724         FEM_Elem_Conn_Checker(const FEM_Entity &sizeSrc_,const FEM_Entity *sizeSrc2_) 
1725                 :sizeSrc(sizeSrc_), sizeSrc2(sizeSrc2_) {}
1726         
1727         void check(int row,const BasicTable2d<int> &table,const char *caller) const {
1728                 const int *idx=table.getRow(row);
1729                 int n=table.width();
1730                 int max=sizeSrc.size();
1731                 if (sizeSrc2) max+=sizeSrc2->size();
1732                 for (int i=0;i<n;i++) 
1733                         if ((idx[i]<0) || (idx[i]>=max))
1734                         { /* This index is out of bounds: */
1735                                 if (idx[i]<0)
1736                                         FEM_Abort(caller,"Connectivity entry %d's value, %d, is negative",row,idx[i]);
1737                                 else /* (idx[i]>=max) */
1738                                         FEM_Abort(caller,
1739                                                 "Connectivity entry %d's value, %d, should be less than the number of nodes, %d",
1740                                                 row,idx[i],max);
1741                         }
1742         }
1743 };
1744
1745 FEM_Elem::FEM_Elem(const FEM_Mesh &mesh, FEM_Elem *ghost_) 
1746   :FEM_Entity(ghost_), elemAdjacency(0), elemAdjacencyTypes(0)
1747 {
1748         FEM_IndexAttribute::Checker *c;
1749         if (isGhost()) // Ghost elements can point to both real as well as ghost nodes
1750                 c=new FEM_Elem_Conn_Checker(mesh.node, mesh.node.getGhost());
1751         else /* is real */ //Real elements only point to real nodes
1752                 c=new FEM_Elem_Conn_Checker(mesh.node, NULL);
1753         conn=new FEM_IndexAttribute(this,FEM_CONN,c);
1754         add(conn); // conn will be deleted by FEM_Entity now
1755 }
1756 void FEM_Elem::pup(PUP::er &p) {
1757         p.comment(" ------------- Element data ---------- ");
1758         FEM_Entity::pup(p);
1759 }
1760 FEM_Elem::~FEM_Elem() {
1761 }
1762
1763
1764 void FEM_Elem::create(int attr,const char *caller) {
1765   // We need to catch both FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES, 
1766   // since if either one falls through to 
1767   // the super::create(), it will not know what to do, and will fail
1768   //
1769   // Note: allocateElemAdjacency() will create both attribute fields since they
1770   //       should always be used together.
1771
1772   if(attr == FEM_CONN) ;
1773     //who allocates the conn?
1774   else if(attr == FEM_ELEM_ELEM_ADJACENCY)
1775     allocateElemAdjacency();
1776   else if(attr == FEM_ELEM_ELEM_ADJ_TYPES)
1777     allocateElemAdjacency();
1778   else
1779     super::create(attr,caller);
1780 }
1781
1782
1783 const char *FEM_Elem::getName(void) const {
1784         return "FEM_ELEM";
1785 }
1786
1787 bool FEM_Elem::hasConn(int idx) {
1788   return false;
1789 }
1790
1791 /********************* Sparse ******************/
1792 /**
1793  * This checker makes sure FEM_Sparse::elem's two element indices
1794  * (element type, element index) are valid.
1795  */
1796 class FEM_Sparse_Elem_Checker : public FEM_IndexAttribute::Checker {
1797         const FEM_Mesh &mesh;
1798 public:
1799         FEM_Sparse_Elem_Checker(const FEM_Mesh &mesh_) :mesh(mesh_) {}
1800         
1801         void check(int row,const BasicTable2d<int> &table,const char *caller) const {
1802                 //assert: table.getWidth==2
1803                 const int *elem=table.getRow(row);
1804                 int maxT=mesh.elem.size();
1805                 if ((elem[0]<0) || (elem[1]<0))
1806                         FEM_Abort(caller,"Sparse element entry %d's values, %d and %d, are negative",
1807                                 row,elem[0],elem[1]);
1808                 int t=elem[0];
1809                 if (t>=maxT)
1810                         FEM_Abort(caller,"Sparse element entry %d's element type, %d, is too big",
1811                                 row,elem[0]);
1812                 if (elem[1]>=mesh.elem[t].size())
1813                         FEM_Abort(caller,"Sparse element entry %d's element index, %d, is too big",
1814                                 row,elem[1]);
1815         }
1816 };
1817
1818 FEM_Sparse::FEM_Sparse(const FEM_Mesh &mesh_,FEM_Sparse *ghost_) 
1819         :FEM_Elem(mesh_,ghost_), elem(0), mesh(mesh_)
1820 {
1821 }
1822 void FEM_Sparse::allocateElem(void) {
1823         if (elem) CkAbort("FEM_Sparse::allocateElem called, but already allocated");
1824         FEM_IndexAttribute::Checker *checker=new FEM_Sparse_Elem_Checker(mesh);
1825         elem=new FEM_IndexAttribute(this,FEM_SPARSE_ELEM,checker);
1826         add(elem); //FEM_Entity will delete elem now
1827         elem->setWidth(2); //SPARSE_ELEM consists of pairs: element type, element number
1828 }
1829 void FEM_Sparse::pup(PUP::er &p) {
1830         p.comment(" ------------- Sparse Element ---------- ");
1831         super::pup(p);
1832 }
1833 FEM_Sparse::~FEM_Sparse() {
1834 }
1835
1836 const char *FEM_Sparse::getName(void) const { return "FEM_SPARSE"; }
1837
1838 void FEM_Sparse::create(int attr,const char *caller) {
1839         if (attr==FEM_SPARSE_ELEM)
1840                 allocateElem();
1841         else /*super*/ FEM_Entity::create(attr,caller);
1842 }
1843
1844
1845 /******************* Mesh *********************/
1846 FEM_Mesh::FEM_Mesh() 
1847         :node(new FEM_Node(NULL)),
1848          elem(*this,"FEM_ELEM"),
1849          sparse(*this,"FEM_SPARSE"),
1850      lastElemAdjLayer(NULL),
1851          fmMM(NULL)
1852 {
1853         m_isSetting=true; //Meshes start out setting
1854         lastLayerSet = false;
1855 }
1856
1857 FEM_Mesh::~FEM_Mesh() {
1858 }
1859
1860 FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) {
1861         FEM_Entity *e=NULL;
1862         if (entity>=FEM_ENTITY_FIRST && entity<FEM_ENTITY_LAST) 
1863         { //It's in the right range for an entity code:
1864                 bool isGhost=false;
1865                 if (entity-FEM_ENTITY_FIRST>=FEM_GHOST) {
1866                         entity-=FEM_GHOST;
1867                         isGhost=true;
1868                 }
1869                 if (entity==FEM_NODE)
1870                         e=&node;
1871                 else if (entity>=FEM_ELEM && entity<FEM_ELEM+100) 
1872                 { //It's a kind of element:
1873                         int elType=entity-FEM_ELEM;
1874                         e=&elem.set(elType);
1875                 }
1876                 else if (entity>=FEM_SPARSE && entity<FEM_SPARSE+100) 
1877                 { //It's a kind of sparse:
1878                         int sID=entity-FEM_SPARSE;
1879                         e=&sparse.set(sID);
1880                 }
1881                 
1882                 if (isGhost) //Move from the real to the ghost entity
1883                         e=e->getGhost();
1884         }
1885         
1886         if (e==NULL) //We didn't find an entity!
1887                 FEM_Abort(caller,"Expected an entity type (FEM_NODE, FEM_ELEM, etc.) but got %d",entity);
1888         return e;
1889 }
1890 const FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) const {
1891         /// FIXME: the const version is quite similar to the above, 
1892         /// but it should *not* create new Entity types...
1893         return ((FEM_Mesh *)this)->lookup(entity,caller);
1894 }
1895
1896 void FEM_Mesh::setFemMeshModify(femMeshModify *m){
1897   fmMM = m;
1898 }
1899
1900 void FEM_Mesh::setParfumSA(ParFUMShadowArray *m){
1901   parfumSA = m;
1902 }
1903
1904
1905 femMeshModify *FEM_Mesh::getfmMM(){
1906   return fmMM;
1907 }
1908
1909 void FEM_Mesh::pup(PUP::er &p)  //For migration
1910 {
1911         p.comment(" ------------- Node Data ---------- ");
1912         node.pup(p);
1913
1914         p.comment(" ------------- Element Types ---------- ");
1915         elem.pup(p);
1916         
1917         p.comment("-------------- Sparse Types ------------");
1918         sparse.pup(p);
1919         
1920         p.comment("-------------- Symmetries ------------");
1921         symList.pup(p);
1922         
1923         p|m_isSetting;
1924         p|lastLayerSet;
1925         if(lastLayerSet) {
1926           if(p.isUnpacking()) {
1927             if(lastElemAdjLayer==NULL) lastElemAdjLayer = new FEM_ElemAdj_Layer();
1928           }
1929           lastElemAdjLayer->pup(p);
1930         }
1931         p.comment("-------------- Mesh data --------------");
1932         udata.pup(p);
1933
1934
1935 /* NOTE: for backward file compatability (fem_mesh_vp files),
1936    be sure to add new stuff at the *end* of this routine--
1937    it will be read as zeros for old files. */
1938 }
1939
1940 int FEM_Mesh::chkET(int elType) const {
1941         if ((elType<0)||(elType>=elem.size())) {
1942                 CkError("FEM Error! Bad element type %d!\n",elType);
1943                 CkAbort("FEM Error! Bad element type used!\n");
1944         }
1945         return elType;
1946 }
1947
1948 int FEM_Mesh::nElems(int t_max) const //Return total number of elements before type t_max
1949 {
1950 #ifndef CMK_OPTIMIZE
1951         if (t_max<0 || t_max>elem.size()) {
1952                 CkPrintf("FEM> Invalid element type %d used!\n");
1953                 CkAbort("FEM> Invalid element type");
1954         }
1955 #endif
1956         int ret=0;
1957         for (int t=0;t<t_max;t++){ 
1958                 if (elem.has(t)){
1959                         ret+=elem.get(t).size();
1960                 }
1961         }       
1962         return ret;
1963 }
1964
1965 int FEM_Mesh::getGlobalElem(int elType,int elNo) const
1966 {
1967         int base=nElems(elType); //Global number of first element of this type
1968 #ifndef CMK_OPTIMIZE
1969         if (elNo<0 || elNo>=elem[elType].size()) {
1970                 CkPrintf("FEM> Element number %d is invalid-- element type %d has only %d elements\n",
1971                         elNo,elType,elem[elType].size());
1972                 CkAbort("FEM> Invalid element number, probably passed via FEM_Set_Sparse_elem");
1973         }
1974 #endif
1975         return base+elNo;
1976 }
1977
1978 /// Set our global numbers as 0...n-1 for nodes, elements, and sparse
1979 void FEM_Mesh::setAscendingGlobalno(void) {
1980         node.setAscendingGlobalno();
1981         for (int e=0;e<elem.size();e++)
1982                 if (elem.has(e)) elem[e].setAscendingGlobalno();
1983         for (int s=0;s<sparse.size();s++)
1984                 if (sparse.has(s)) sparse[s].setAscendingGlobalno();
1985 }
1986 void FEM_Mesh::setAbsoluteGlobalno(){
1987         node.setAscendingGlobalno();
1988         for (int e=0;e<elem.size();e++){
1989                 if (elem.has(e)) elem[e].setAscendingGlobalno(nElems(e));
1990         }       
1991 }
1992
1993 void FEM_Mesh::copyOldGlobalno(const FEM_Mesh &m) {
1994         node.copyOldGlobalno(m.node);
1995         for (int e=0;e<m.elem.size();e++)
1996                 if (m.elem.has(e) && e<elem.size() && elem.has(e)) 
1997                         elem[e].copyOldGlobalno(m.elem[e]);
1998         for (int s=0;s<m.sparse.size();s++)
1999                 if (m.sparse.has(s) && s<sparse.size() && sparse.has(s)) 
2000                         sparse[s].copyOldGlobalno(m.sparse[s]);
2001 }
2002
2003 void FEM_Mesh::clearSharedNodes(){
2004         node.shared.clear();
2005 }
2006
2007 void FEM_Mesh::clearGhostNodes(){
2008         node.clearGhost();
2009 }
2010
2011 void FEM_Mesh::clearGhostElems(){
2012         for(int i=0;i<elem.size();i++){
2013                 if(elem.has(i)){
2014                         elem[i].clearGhost();
2015                 }
2016         }
2017 }
2018
2019 void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType) {
2020         if (type<0 || type>maxType) {
2021                 char msg[1024];
2022                 sprintf(msg,"%s %d is not a valid entity type (it must be between %d and %d)",
2023                         entityType,type, 0, maxType-1);
2024                 FEM_Abort(caller,msg);
2025         }
2026 }
2027 void FEM_Is_NULL(const char *caller,const char *entityType,int type) {
2028         char msg[1024];
2029         sprintf(msg,"%s %d was never set--it cannot now be read",entityType,type);
2030         FEM_Abort(caller,msg);
2031 }
2032
2033 void FEM_Mesh::copyShape(const FEM_Mesh &src)
2034 {
2035         node.copyShape(src.node);
2036         for (int t=0;t<src.elem.size();t++) 
2037                 if (src.elem.has(t)) elem.set(t).copyShape(src.elem.get(t));
2038         
2039         for (int s=0;s<src.sparse.size();s++)
2040                 if (src.sparse.has(s)) sparse.set(s).copyShape(src.sparse.get(s));
2041         
2042         setSymList(src.getSymList());
2043 }
2044
2045 /// Extract a list of our entities:
2046 int FEM_Mesh::getEntities(int *entities) {
2047         int len=0;
2048         entities[len++]=FEM_NODE;
2049         for (int t=0;t<elem.size();t++) 
2050                 if (elem.has(t)) entities[len++]=FEM_ELEM+t;
2051         for (int s=0;s<sparse.size();s++)
2052                 if (sparse.has(s)) entities[len++]=FEM_SPARSE+s;
2053         return len;
2054 }
2055
2056
2057 FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead)
2058 {
2059     char fname[256];
2060     static const char *meshFileNames="%s_vp%d_%d.dat";
2061     sprintf(fname, meshFileNames, prefix, nchunks, chunkNo);
2062     FILE *fp = fopen(fname, forRead?"r":"w");
2063     CkPrintf("FEM> %s %s...\n",forRead?"Reading":"Writing",fname);  
2064     if(fp==0) {
2065       FEM_Abort(forRead?"FEM: unable to open input file"
2066         :"FEM: unable to create output file.\n");
2067     }
2068     return fp;
2069 }
2070
2071 static inline void read_version()
2072 {
2073     FILE *f = fopen("FEMVERSION", "r");
2074     if (f!=NULL)  {
2075         fscanf(f, "%d", &femVersion);
2076         if (CkMyPe()==0) CkPrintf("FEM> femVersion detected: %d\n", femVersion);
2077         fclose(f);
2078     }
2079 }
2080
2081 static inline void write_version()
2082 {
2083     FILE *f = fopen("FEMVERSION", "w");
2084     if (f!=NULL)  {
2085                 fprintf(f, "%d", femVersion);
2086                 fclose(f);
2087     }
2088 }
2089
2090 FEM_Mesh *FEM_readMesh(const char *prefix,int chunkNo,int nChunks)
2091 {
2092         // find FEM file version number
2093         static int version_checked = 0;
2094         if (!version_checked) {
2095             version_checked=1;
2096             read_version();
2097         }
2098
2099         FEM_Mesh *ret=new FEM_Mesh;
2100         ret->becomeGetting();
2101         FILE *fp = FEM_openMeshFile(prefix,chunkNo,nChunks,true);
2102         PUP::fromTextFile p(fp);
2103         ret->pup(p);
2104         fclose(fp);
2105         return ret;
2106 }
2107
2108 void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks)
2109 {
2110         if (chunkNo == 0) {
2111            write_version();
2112         }
2113         FILE *fp = FEM_openMeshFile(prefix,chunkNo,nChunks,false);
2114         PUP::toTextFile p(fp);
2115         m->pup(p);
2116         fclose(fp);
2117 }
2118
2119
2120 // Setup the entity FEM_IS_VALID tables
2121 CDECL void FEM_Mesh_allocate_valid_attr(int fem_mesh, int entity_type){  
2122   FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_create_valid_elem");
2123   FEM_Entity *entity = m->lookup(entity_type,"FEM_Mesh_allocate_valid_attr");
2124   entity->allocateValid();
2125 }
2126 FORTRAN_AS_C(FEM_MESH_ALLOCATE_VALID_ATTR,
2127              FEM_Mesh_allocate_valid_attr,
2128              fem_mesh_allocate_valid_attr, 
2129              (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )
2130
2131
2132 // mark an entity as valid
2133 CDECL void FEM_set_entity_valid(int mesh, int entityType, int entityIdx){
2134   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_");
2135   FEM_Entity *entity = m->lookup(entityType,"FEM_");
2136   entity->set_valid(entityIdx,true);
2137 }
2138 FORTRAN_AS_C(FEM_SET_ENTITY_VALID, 
2139              FEM_set_entity_valid, 
2140              fem_set_entity_valid,  
2141              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) ) 
2142
2143
2144 // mark an entity as invalid
2145 CDECL void FEM_set_entity_invalid(int mesh, int entityType, int entityIdx){
2146   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2147   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2148   return entity->set_invalid(entityIdx,true);
2149 }
2150 FORTRAN_AS_C(FEM_SET_ENTITY_INVALID,  
2151              FEM_set_entity_invalid,  
2152              fem_set_entity_invalid,   
2153              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )  
2154
2155
2156 // Determine if an entity is valid
2157 CDECL int FEM_is_valid(int mesh, int entityType, int entityIdx){
2158   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2159   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2160   return (entity->is_valid(entityIdx) != 0);
2161 }
2162 FORTRAN_AS_C(FEM_IS_VALID,   
2163              FEM_is_valid,   
2164              fem_is_valid,    
2165              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )   
2166
2167
2168 // Count number of valid items for this entity type
2169 CDECL int FEM_count_valid(int mesh, int entityType){
2170   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2171   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2172   return entity->count_valid();
2173 }
2174 FORTRAN_AS_C(FEM_COUNT_VALID,    
2175              FEM_count_valid,    
2176              fem_count_valid,     
2177              (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )    
2178  
2179
2180 // Set coordinates for some entity's item number idx 
2181 void FEM_set_entity_coord2(int mesh, int entityType, int idx, double x, double y){
2182   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2183   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2184   entity->set_coord(idx,x,y);
2185 }
2186 void FEM_set_entity_coord3(int mesh, int entityType, int idx, double x, double y, double z){
2187   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2188   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2189   entity->set_coord(idx,x,y,z);
2190 }