Added some functional get_coord() functions that will in all reasonable cases return...
[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
1330 /** Determine if an index refers to a valid non-ghost entity. Index must be between 0 and size(). Changing this assumption will break things. If required use the other functions is_valid_any_idx() or is_valid_nonghost_idx() */
1331 int FEM_Entity::is_valid(int idx){
1332     CkAssert(idx < size() && idx >=0);
1333     return valid->getChar()(idx,0);
1334 }
1335
1336 /** Determine if an index refers to a valid entity, with negatives corresponding to a valid ghost */
1337 int FEM_Entity::is_valid_any_idx(int idx){
1338   if(idx == -1 || idx >= size())
1339         return false;
1340   else if(idx < 0)
1341         return ghost->is_valid_any_idx(FEM_To_ghost_index(idx));
1342   else
1343     return valid->getChar()(idx,0);
1344 }
1345
1346 /** Determine if an index refers to a valid non-ghost entity */
1347 int FEM_Entity::is_valid_nonghost_idx(int idx){
1348   if(idx < 0 || idx >= size())
1349     return false;
1350   else
1351     return valid->getChar()(idx,0);
1352 }
1353
1354 /** Linear scan to count valid entities */
1355 int FEM_Entity::count_valid(){
1356   int count=0;
1357   for(int i=0;i<size();i++)
1358         if(is_valid(i)) count++;
1359   return count;
1360 }
1361
1362 void FEM_Entity::set_all_invalid(){
1363     // Set all to invalid
1364     for(int i=0;i<size();i++) {
1365       valid->getChar()(i,0)=0;
1366     }
1367 }
1368
1369
1370 /// Get an entry(entity index) that corresponds to an invalid entity
1371 /// The invalid slot in the tables can then be reused when "creating" a new element or node
1372 /// We either return an empty slot, or resize the array and return a value at the end
1373 /// If someone has a better name for this function, please change it.
1374 int FEM_Entity::get_next_invalid(FEM_Mesh *m, bool isNode, bool isGhost){
1375   int retval=0;
1376 //  if(true) { //maintains a list of invalid elements. FASTER
1377     if(invalidListLen>0) {
1378       retval = invalidList[invalidListLen-1];
1379     }
1380     else {
1381       retval = size();
1382       setLength(retval+1,true);
1383     }
1384 //  }
1385 /*  else {
1386     int size1 = size();
1387     if(size1==0) { //special case because for size=0, first_invalid=last_invalid=0
1388       retval= 0;
1389       setLength(1);
1390     }
1391     else {
1392       CkAssert(!is_valid(first_invalid) || first_invalid==0);
1393       // if we have an invalid entity to return
1394       bool flag1 = false;
1395       if(!is_valid(first_invalid)){
1396         retval = first_invalid;
1397         if(isNode && !isGhost) { //it is a node & the entity is not a ghost entity
1398           while(retval <= last_invalid) {
1399             if(!is_valid(retval)) {
1400               if(m->getfmMM()->fmLockN[retval].haslocks()) {
1401                 retval++;
1402               }
1403               else if(hasConn(retval)) { //has some connectivity
1404                 retval++;
1405               }
1406               else {
1407                 flag1 = true;
1408                 break;
1409               }
1410             }
1411             else retval++;
1412           }
1413         }
1414         else if(isNode) {
1415           while(retval <= last_invalid) {
1416             if(!is_valid(retval)) {
1417               if(hasConn(retval)) { //has some connectivity
1418                 retval++;
1419               }
1420               else {
1421                 flag1 = true;
1422                 break;
1423               }
1424             }
1425             else retval++;
1426           }
1427         }
1428         else{
1429           // resize array and return new entity
1430           flag1 = true;
1431         }
1432       }
1433       if(!flag1) {
1434         retval = size1;
1435         setLength(retval+1);
1436       }
1437     }
1438   }*/
1439   set_valid(retval,isNode);
1440   return retval;
1441 }
1442
1443 void FEM_Entity::setMaxLength(int newLen,int newMaxLen,void *pargs,FEM_Mesh_alloc_fn fn){
1444   //CkPrintf("resize fn %p \n",fn);
1445         max = newMaxLen;
1446         resize = fn;
1447         args = pargs;
1448         setLength(newLen);
1449 };
1450
1451
1452 /// Copy src[srcEntity] into our dstEntity.
1453 void FEM_Entity::copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity) {
1454         FEM_Entity *dp=this; //Destination entity
1455         const FEM_Entity *sp=&src;
1456         for (int a=0;a<sp->attributes.size();a++)
1457         { //We copy each of his attributes:
1458                 const FEM_Attribute *Asrc=sp->attributes[a];
1459                 FEM_Attribute *Adst=dp->lookup(Asrc->getAttr(),"FEM_Entity::copyEntity");
1460                 Adst->copyEntity(dstEntity,*Asrc,srcEntity);
1461         }
1462 }
1463
1464 /// Add room for one more entity, with initial values from src[srcEntity],
1465 /// and return the new entity's index.
1466 int FEM_Entity::push_back(const FEM_Entity &src,int srcEntity) {
1467         int dstEntity=size();
1468         setLength(dstEntity+1);
1469         copyEntity(dstEntity,src,srcEntity);
1470         return dstEntity;
1471 }
1472
1473 /// Add this attribute to this kind of Entity.
1474 /// This method is normally called by the default lookup method.
1475 void FEM_Entity::add(FEM_Attribute *attribute) {
1476         if (ghost!=NULL)
1477         { //Look up (or create) the ghost attribute, too:
1478                 attribute->setGhost(ghost->lookup(attribute->getAttr(),"FEM_Entity::add"));
1479         }
1480         attributes.push_back(attribute);
1481 }
1482
1483 /**
1484  * Find this attribute (from an FEM_ATTR code) of this entity.
1485  * The default implementation searches the list of userdata attributes;
1486  * subclasses with other attributes should override this routine.
1487  */
1488 FEM_Attribute *FEM_Entity::lookup(int attr,const char *caller) {
1489         //Try to find an existing attribute (FIXME: keep attributes in a map, to speed this up)
1490         for (int a=0;a<attributes.size();a++) {
1491                 if (attributes[a]->getAttr()==attr)
1492                         return attributes[a];
1493         }
1494
1495         //If we get here, no existing attribute fits the bill: create one
1496         create(attr,caller);
1497
1498         // If create did its job, the next lookup should succeed:
1499         return lookup(attr,caller);
1500 }
1501
1502 /**
1503  * Create a new attribute from an FEM_ATTR code.
1504  * The default implementation handles FEM_DATA tags; entity-specific
1505  * attributes (like FEM_CONN) need to be overridden and created
1506  * by subclasses.
1507  */
1508 void FEM_Entity::create(int attr,const char *caller) {
1509   if (attr<=FEM_ATTRIB_TAG_MAX) //It's a valid user data tag
1510         add(new FEM_DataAttribute(this,attr));
1511   else if (attr==FEM_COORD)
1512         allocateCoord();
1513   else if (attr==FEM_SYMMETRIES)
1514         allocateSym();
1515   else if (attr==FEM_GLOBALNO)
1516         allocateGlobalno();
1517   else if (attr==FEM_IS_VALID_ATTR)
1518         allocateValid();
1519   else if (attr==FEM_MESH_SIZING)
1520         allocateMeshSizing();
1521   else if(attr == FEM_CHUNK){
1522         FEM_IndexAttribute *chunkNo= new FEM_IndexAttribute(this,FEM_CHUNK,NULL);
1523         add(chunkNo);
1524         chunkNo->setWidth(1);
1525   }
1526   else if(attr == FEM_BOUNDARY){
1527       allocateBoundary();
1528   }else if(attr == FEM_ADAPT_FACE_ADJ){
1529       FEM_DataAttribute *adaptAdj = new FEM_DataAttribute(this,attr);
1530       adaptAdj->setDatatype(FEM_BYTE);
1531       add(adaptAdj);
1532   }else if(attr == FEM_ADAPT_EDGE_ADJ){
1533       FEM_DataAttribute *adaptAdj = new FEM_DataAttribute(this,attr);
1534       adaptAdj->setDatatype(FEM_BYTE);
1535       add(adaptAdj);
1536   }else if(attr == FEM_ADAPT_LOCK){
1537       FEM_DataAttribute *adaptLock = new FEM_DataAttribute(this,attr);
1538       adaptLock->setDatatype(FEM_INT);
1539       adaptLock->setWidth(2);
1540       add(adaptLock);
1541   }
1542   else {
1543         //It's an unrecognized tag: abort
1544         char attrNameStorage[256], msg[1024];
1545         sprintf(msg,"Could not locate the attribute %s for entity %s",
1546                         FEM_Get_attr_name(attr,attrNameStorage), getName());
1547         FEM_Abort(caller,msg);
1548   }
1549 }
1550
1551 void FEM_Entity::allocateCoord(void) {
1552         if (coord) CkAbort("FEM_Entity::allocateCoord called, but already allocated");
1553         coord=new FEM_DataAttribute(this,FEM_COORD);
1554         add(coord); // coord will be deleted by FEM_Entity now
1555         coord->setDatatype(FEM_DOUBLE);
1556 }
1557
1558 void FEM_Entity::allocateSym(void) {
1559         if (sym) CkAbort("FEM_Entity::allocateSym called, but already allocated");
1560         sym=new FEM_DataAttribute(this,FEM_SYMMETRIES);
1561         add(sym); // sym will be deleted via attributes list now
1562         sym->setWidth(1);
1563         sym->setDatatype(FEM_BYTE); //Same as FEM_Symmetries_t
1564 }
1565
1566 void FEM_Entity::setSymmetries(int r,FEM_Symmetries_t s)
1567 {
1568         if (!sym) {
1569                 if (s==0) return; //Don't bother allocating just for 0
1570                 allocateSym();
1571         }
1572         sym->getChar()(r,0)=s;
1573 }
1574
1575
1576 /*
1577  * Set the coordinates for a node or other entity.
1578  * Use the appropriate 2d or 3d version.
1579  *
1580  * Note that first the function attempts to find coordinates in
1581  * the FEM_COORD attribute's array "*coord". If this fails, then
1582  * it will use the user's FEM_DATA field. Little error checking is
1583  * done, so the functions may crash if used inappropriately.
1584  */
1585 inline void FEM_Entity::set_coord(int idx, double x, double y){
1586   if(coord){
1587         coord->getDouble()(idx,0)=x;
1588         coord->getDouble()(idx,1)=y;
1589   }
1590   else {
1591         FEM_DataAttribute* attr =       (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
1592         attr->getDouble()(idx,0)=x;
1593         attr->getDouble()(idx,1)=y;
1594   }
1595 }
1596
1597 inline void FEM_Entity::set_coord(int idx, double x, double y, double z){
1598   if(coord){
1599         coord->getDouble()(idx,0)=x;
1600         coord->getDouble()(idx,1)=y;
1601         coord->getDouble()(idx,2)=z;
1602   }
1603   else {
1604         FEM_DataAttribute* attr =       (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
1605         attr->getDouble()(idx,0)=x;
1606         attr->getDouble()(idx,1)=y;
1607         attr->getDouble()(idx,2)=z;
1608   }
1609 }
1610
1611
1612 /** Get coordinates for a given 2-D node, if idx<0 the corresponding ghost node is used 
1613 @note there are a couple places coordinates have been historically stored. We can find either of them here.
1614 */
1615 void FEM_Entity::get_coord(int idx, double &x, double &y){
1616   if(coord){
1617       if(idx<0){ // ghost
1618         x = ghost->coord->getDouble()(FEM_From_ghost_index(idx),0);
1619         y = ghost->coord->getDouble()(FEM_From_ghost_index(idx),1);
1620       } else { // non-ghost
1621         x = coord->getDouble()(idx,0);
1622         y = coord->getDouble()(idx,1);
1623       }
1624       
1625     }
1626   else {
1627   
1628     if(idx<0){ // ghost
1629       FEM_DataAttribute* attr = (FEM_DataAttribute*)ghost->lookup(FEM_DATA,"get_coord");
1630       x = attr->getDouble()(FEM_From_ghost_index(idx),0);
1631       y = attr->getDouble()(FEM_From_ghost_index(idx),1);
1632     } else  { // non-ghost
1633       FEM_DataAttribute* attr = (FEM_DataAttribute*)lookup(FEM_DATA,"get_coord");
1634       x = attr->getDouble()(idx,0);
1635       y = attr->getDouble()(idx,1);
1636     }
1637     
1638   }
1639   
1640 }
1641
1642 /** Get coordinates for a given 3-D node, if idx<0 the corresponding ghost node is used 
1643 @note there are a couple places coordinates have been historically stored. We can find either of them here.
1644  */
1645 void FEM_Entity::get_coord(int idx, double &x, double &y, double &z){
1646   if(coord){
1647     if(idx<0){ // ghost
1648       x = ghost->coord->getDouble()(FEM_From_ghost_index(idx),0);
1649       y = ghost->coord->getDouble()(FEM_From_ghost_index(idx),1);
1650       z = ghost->coord->getDouble()(FEM_From_ghost_index(idx),2);
1651     } else { // non-ghost
1652       x = coord->getDouble()(idx,0);
1653       y = coord->getDouble()(idx,1);
1654       z = coord->getDouble()(idx,2);
1655     }
1656     
1657   }
1658   else {
1659      
1660     if(idx<0){ // ghost
1661       FEM_DataAttribute* attr = (FEM_DataAttribute*)ghost->lookup(FEM_DATA,"get_coord");
1662       x = attr->getDouble()(FEM_From_ghost_index(idx),0);
1663       y = attr->getDouble()(FEM_From_ghost_index(idx),1);
1664       z = attr->getDouble()(FEM_From_ghost_index(idx),2);
1665     } else { // non-ghost
1666       FEM_DataAttribute* attr = (FEM_DataAttribute*)lookup(FEM_DATA,"get_coord");
1667       x = attr->getDouble()(idx,0);
1668       y = attr->getDouble()(idx,1);
1669       z = attr->getDouble()(idx,2);
1670     }
1671     
1672   }
1673 }
1674
1675
1676
1677
1678
1679 void FEM_Entity::allocateGlobalno(void) {
1680         if (globalno) CkAbort("FEM_Entity::allocateGlobalno called, but already allocated");
1681         globalno=new FEM_IndexAttribute(this,FEM_GLOBALNO,NULL);
1682         add(globalno); // globalno will be deleted via attributes list now
1683         globalno->setWidth(1);
1684 }
1685
1686 void FEM_Entity::allocateMeshSizing(void) {
1687   if (meshSizing)
1688     CkAbort("FEM_Entity::allocateMeshSizing called, but already allocated");
1689   meshSizing=new FEM_DataAttribute(this,FEM_MESH_SIZING);
1690   add(meshSizing); // meshSizing will be deleted via attributes list now
1691   meshSizing->setWidth(1);
1692   meshSizing->setDatatype(FEM_DOUBLE);
1693 }
1694
1695 double FEM_Entity::getMeshSizing(int r) {
1696   if (!meshSizing) {
1697     allocateMeshSizing();
1698     return -1.0;
1699   }
1700   if (r >= 0)  return meshSizing->getDouble()(r,0);
1701   else  return ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0);
1702 }
1703
1704 void FEM_Entity::setMeshSizing(int r,double s)
1705 {
1706   if (!meshSizing) allocateMeshSizing();
1707   if (s <= 0.0) return;
1708   if (r >= 0)  meshSizing->getDouble()(r,0)=s;
1709   else ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0)=s;
1710 }
1711
1712 void FEM_Entity::setMeshSizing(double *sf)
1713 {
1714   if (!meshSizing) allocateMeshSizing();
1715   int len = size();
1716   for (int i=0; i<len; i++)
1717     meshSizing->getDouble()(i,0)=sf[i];
1718 }
1719
1720 void FEM_Entity::allocateBoundary(){
1721         FEM_DataAttribute *bound = new FEM_DataAttribute(this,FEM_BOUNDARY);
1722         add(bound);
1723         bound->setWidth(1);
1724         bound->setDatatype(FEM_INT);
1725 }
1726
1727 void FEM_Entity::setGlobalno(int r,int g) {
1728         if (!globalno) allocateGlobalno();
1729         globalno->get()(r,0)=g;
1730 }
1731 void FEM_Entity::setAscendingGlobalno(void) {
1732         if (!globalno) {
1733                 allocateGlobalno();
1734                 int len=size();
1735                 for (int i=0;i<len;i++) globalno->get()(i,0)=i;
1736         }
1737 }
1738 void FEM_Entity::setAscendingGlobalno(int base) {
1739         if (!globalno) {
1740                 allocateGlobalno();
1741                 int len=size();
1742                 for (int i=0;i<len;i++) globalno->get()(i,0)=i+base;
1743         }
1744 }
1745 void FEM_Entity::copyOldGlobalno(const FEM_Entity &e) {
1746         if ((!hasGlobalno()) && e.hasGlobalno() && size()>=e.size()) {
1747                 for (int i=0;i<size();i++)
1748                         setGlobalno(i,e.getGlobalno(i));
1749         }
1750 }
1751
1752 /*
1753  * method to clear ghosts associated with this entity
1754  * */
1755
1756 void FEM_Entity::clearGhost(){
1757         CkAssert(ghost != NULL);
1758         ghostSend.clear();
1759         ghost->setLength(0);
1760         ghost->ghostRecv.clear();
1761 };
1762
1763 /********************** Node *****************/
1764 FEM_Node::FEM_Node(FEM_Node *ghost_)
1765   :FEM_Entity(ghost_), primary(0), sharedIDXL(&shared,&shared),
1766    elemAdjacency(0),nodeAdjacency(0)
1767 {}
1768
1769 void FEM_Node::allocatePrimary(void) {
1770   if (primary) CkAbort("FEM_Node::allocatePrimary called, but already allocated");
1771   primary=new FEM_DataAttribute(this,FEM_NODE_PRIMARY);
1772   add(primary); // primary will be deleted by FEM_Entity now
1773   primary->setWidth(1); //Only 1 flag per node
1774   primary->setDatatype(FEM_BYTE);
1775 }
1776
1777 bool FEM_Node::hasConn(int idx) {
1778   if((elemAdjacency->get()[idx].length() > 0)||(nodeAdjacency->get()[idx].length() > 0))
1779     return true;
1780   else return false;
1781 }
1782
1783 void FEM_Node::pup(PUP::er &p) {
1784         p.comment(" ---------------- Nodes ------------------ ");
1785         super::pup(p);
1786         p.comment(" ---- Shared nodes ----- ");
1787         shared.pup(p);
1788         p.comment(" shared nodes IDXL ");
1789         sharedIDXL.pup(p);
1790 }
1791 FEM_Node::~FEM_Node() {
1792 }
1793
1794
1795 const char *FEM_Node::getName(void) const {return "FEM_NODE";}
1796
1797 void FEM_Node::create(int attr,const char *caller) {
1798   if (attr==FEM_NODE_PRIMARY)
1799         allocatePrimary();
1800   else if(attr == FEM_NODE_ELEM_ADJACENCY)
1801         allocateElemAdjacency();
1802   else if(attr == FEM_NODE_NODE_ADJACENCY)
1803         allocateNodeAdjacency();
1804   else
1805         super::create(attr,caller);
1806 }
1807
1808
1809 /********************** Elem *****************/
1810 /// This checker verifies that FEM_Elem::conn's entries are valid node indices.
1811 class FEM_Elem_Conn_Checker : public FEM_IndexAttribute::Checker {
1812         const FEM_Entity &sizeSrc;
1813         const FEM_Entity *sizeSrc2;
1814 public:
1815         FEM_Elem_Conn_Checker(const FEM_Entity &sizeSrc_,const FEM_Entity *sizeSrc2_)
1816                 :sizeSrc(sizeSrc_), sizeSrc2(sizeSrc2_) {}
1817
1818         void check(int row,const BasicTable2d<int> &table,const char *caller) const {
1819                 const int *idx=table.getRow(row);
1820                 int n=table.width();
1821                 int max=sizeSrc.size();
1822                 if (sizeSrc2) max+=sizeSrc2->size();
1823                 for (int i=0;i<n;i++)
1824                         if ((idx[i]<0) || (idx[i]>=max))
1825                         { /* This index is out of bounds: */
1826                                 if (idx[i]<0)
1827                                         FEM_Abort(caller,"Connectivity entry %d's value, %d, is negative",row,idx[i]);
1828                                 else /* (idx[i]>=max) */
1829                                         FEM_Abort(caller,
1830                                                 "Connectivity entry %d's value, %d, should be less than the number of nodes, %d",
1831                                                 row,idx[i],max);
1832                         }
1833         }
1834 };
1835
1836 FEM_Elem::FEM_Elem(const FEM_Mesh &mesh, FEM_Elem *ghost_)
1837   :FEM_Entity(ghost_), elemAdjacency(0), elemAdjacencyTypes(0)
1838 {
1839         FEM_IndexAttribute::Checker *c;
1840         if (isGhost()) // Ghost elements can point to both real as well as ghost nodes
1841                 c=new FEM_Elem_Conn_Checker(mesh.node, mesh.node.getGhost());
1842         else /* is real */ //Real elements only point to real nodes
1843                 c=new FEM_Elem_Conn_Checker(mesh.node, NULL);
1844         conn=new FEM_IndexAttribute(this,FEM_CONN,c);
1845         add(conn); // conn will be deleted by FEM_Entity now
1846 }
1847 void FEM_Elem::pup(PUP::er &p) {
1848         p.comment(" ------------- Element data ---------- ");
1849         FEM_Entity::pup(p);
1850 }
1851 FEM_Elem::~FEM_Elem() {
1852 }
1853
1854
1855 void FEM_Elem::create(int attr,const char *caller) {
1856   // We need to catch both FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES,
1857   // since if either one falls through to
1858   // the super::create(), it will not know what to do, and will fail
1859   //
1860   // Note: allocateElemAdjacency() will create both attribute fields since they
1861   //       should always be used together.
1862
1863   if(attr == FEM_CONN) ;
1864     //who allocates the conn?
1865   else if(attr == FEM_ELEM_ELEM_ADJACENCY)
1866     allocateElemAdjacency();
1867   else if(attr == FEM_ELEM_ELEM_ADJ_TYPES)
1868     allocateElemAdjacency();
1869   else
1870     super::create(attr,caller);
1871 }
1872
1873
1874 const char *FEM_Elem::getName(void) const {
1875         return "FEM_ELEM";
1876 }
1877
1878 bool FEM_Elem::hasConn(int idx) {
1879   return false;
1880 }
1881
1882 /********************* Sparse ******************/
1883 /**
1884  * This checker makes sure FEM_Sparse::elem's two element indices
1885  * (element type, element index) are valid.
1886  */
1887 class FEM_Sparse_Elem_Checker : public FEM_IndexAttribute::Checker {
1888         const FEM_Mesh &mesh;
1889 public:
1890         FEM_Sparse_Elem_Checker(const FEM_Mesh &mesh_) :mesh(mesh_) {}
1891
1892         void check(int row,const BasicTable2d<int> &table,const char *caller) const {
1893                 //assert: table.getWidth==2
1894                 const int *elem=table.getRow(row);
1895                 int maxT=mesh.elem.size();
1896                 if ((elem[0]<0) || (elem[1]<0))
1897                         FEM_Abort(caller,"Sparse element entry %d's values, %d and %d, are negative",
1898                                 row,elem[0],elem[1]);
1899                 int t=elem[0];
1900                 if (t>=maxT)
1901                         FEM_Abort(caller,"Sparse element entry %d's element type, %d, is too big",
1902                                 row,elem[0]);
1903                 if (elem[1]>=mesh.elem[t].size())
1904                         FEM_Abort(caller,"Sparse element entry %d's element index, %d, is too big",
1905                                 row,elem[1]);
1906         }
1907 };
1908
1909 FEM_Sparse::FEM_Sparse(const FEM_Mesh &mesh_,FEM_Sparse *ghost_)
1910         :FEM_Elem(mesh_,ghost_), elem(0), mesh(mesh_)
1911 {
1912 }
1913 void FEM_Sparse::allocateElem(void) {
1914         if (elem) CkAbort("FEM_Sparse::allocateElem called, but already allocated");
1915         FEM_IndexAttribute::Checker *checker=new FEM_Sparse_Elem_Checker(mesh);
1916         elem=new FEM_IndexAttribute(this,FEM_SPARSE_ELEM,checker);
1917         add(elem); //FEM_Entity will delete elem now
1918         elem->setWidth(2); //SPARSE_ELEM consists of pairs: element type, element number
1919 }
1920 void FEM_Sparse::pup(PUP::er &p) {
1921         p.comment(" ------------- Sparse Element ---------- ");
1922         super::pup(p);
1923 }
1924 FEM_Sparse::~FEM_Sparse() {
1925 }
1926
1927 const char *FEM_Sparse::getName(void) const { return "FEM_SPARSE"; }
1928
1929 void FEM_Sparse::create(int attr,const char *caller) {
1930         if (attr==FEM_SPARSE_ELEM)
1931                 allocateElem();
1932         else /*super*/ FEM_Entity::create(attr,caller);
1933 }
1934
1935
1936 /******************* Mesh *********************/
1937 FEM_Mesh::FEM_Mesh()
1938         :node(new FEM_Node(NULL)),
1939          elem(*this,"FEM_ELEM"),
1940          sparse(*this,"FEM_SPARSE"),
1941      lastElemAdjLayer(NULL),
1942          fmMM(NULL)
1943 {
1944         m_isSetting=true; //Meshes start out setting
1945         lastLayerSet = false;
1946 }
1947
1948 FEM_Mesh::~FEM_Mesh() {
1949 }
1950
1951 FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) {
1952         FEM_Entity *e=NULL;
1953         if (entity>=FEM_ENTITY_FIRST && entity<FEM_ENTITY_LAST)
1954         { //It's in the right range for an entity code:
1955                 bool isGhost=false;
1956                 if (entity-FEM_ENTITY_FIRST>=FEM_GHOST) {
1957                         entity-=FEM_GHOST;
1958                         isGhost=true;
1959                 }
1960                 if (entity==FEM_NODE)
1961                         e=&node;
1962                 else if (entity>=FEM_ELEM && entity<FEM_ELEM+100)
1963                 { //It's a kind of element:
1964                         int elType=entity-FEM_ELEM;
1965                         e=&elem.set(elType);
1966                 }
1967                 else if (entity>=FEM_SPARSE && entity<FEM_SPARSE+100)
1968                 { //It's a kind of sparse:
1969                         int sID=entity-FEM_SPARSE;
1970                         e=&sparse.set(sID);
1971                 }
1972
1973                 if (isGhost) //Move from the real to the ghost entity
1974                         e=e->getGhost();
1975         }
1976
1977         if (e==NULL) //We didn't find an entity!
1978                 FEM_Abort(caller,"Expected an entity type (FEM_NODE, FEM_ELEM, etc.) but got %d",entity);
1979         return e;
1980 }
1981 const FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) const {
1982         /// FIXME: the const version is quite similar to the above,
1983         /// but it should *not* create new Entity types...
1984         return ((FEM_Mesh *)this)->lookup(entity,caller);
1985 }
1986
1987 void FEM_Mesh::setFemMeshModify(femMeshModify *m){
1988   fmMM = m;
1989 }
1990
1991 void FEM_Mesh::setParfumSA(ParFUMShadowArray *m){
1992   parfumSA = m;
1993 }
1994
1995
1996 femMeshModify *FEM_Mesh::getfmMM(){
1997   return fmMM;
1998 }
1999
2000 void FEM_Mesh::pup(PUP::er &p)  //For migration
2001 {
2002         p.comment(" ------------- Node Data ---------- ");
2003         node.pup(p);
2004
2005         p.comment(" ------------- Element Types ---------- ");
2006         elem.pup(p);
2007
2008         p.comment("-------------- Sparse Types ------------");
2009         sparse.pup(p);
2010
2011         p.comment("-------------- Symmetries ------------");
2012         symList.pup(p);
2013
2014         p|m_isSetting;
2015         p|lastLayerSet;
2016         if(lastLayerSet) {
2017           if(p.isUnpacking()) {
2018             if(lastElemAdjLayer==NULL) lastElemAdjLayer = new FEM_ElemAdj_Layer();
2019           }
2020           lastElemAdjLayer->pup(p);
2021         }
2022         p.comment("-------------- Mesh data --------------");
2023         udata.pup(p);
2024
2025
2026 /* NOTE: for backward file compatability (fem_mesh_vp files),
2027    be sure to add new stuff at the *end* of this routine--
2028    it will be read as zeros for old files. */
2029 }
2030
2031 int FEM_Mesh::chkET(int elType) const {
2032         if ((elType<0)||(elType>=elem.size())) {
2033                 CkError("FEM Error! Bad element type %d!\n",elType);
2034                 CkAbort("FEM Error! Bad element type used!\n");
2035         }
2036         return elType;
2037 }
2038
2039 int FEM_Mesh::nElems(int t_max) const //Return total number of elements before type t_max
2040 {
2041 #ifndef CMK_OPTIMIZE
2042         if (t_max<0 || t_max>elem.size()) {
2043                 CkPrintf("FEM> Invalid element type %d used!\n");
2044                 CkAbort("FEM> Invalid element type");
2045         }
2046 #endif
2047         int ret=0;
2048         for (int t=0;t<t_max;t++){
2049                 if (elem.has(t)){
2050                         ret+=elem.get(t).size();
2051                 }
2052         }
2053         return ret;
2054 }
2055
2056 int FEM_Mesh::getGlobalElem(int elType,int elNo) const
2057 {
2058         int base=nElems(elType); //Global number of first element of this type
2059 #ifndef CMK_OPTIMIZE
2060         if (elNo<0 || elNo>=elem[elType].size()) {
2061                 CkPrintf("FEM> Element number %d is invalid-- element type %d has only %d elements\n",
2062                         elNo,elType,elem[elType].size());
2063                 CkAbort("FEM> Invalid element number, probably passed via FEM_Set_Sparse_elem");
2064         }
2065 #endif
2066         return base+elNo;
2067 }
2068
2069 /// Set our global numbers as 0...n-1 for nodes, elements, and sparse
2070 void FEM_Mesh::setAscendingGlobalno(void) {
2071         node.setAscendingGlobalno();
2072         for (int e=0;e<elem.size();e++)
2073                 if (elem.has(e)) elem[e].setAscendingGlobalno();
2074         for (int s=0;s<sparse.size();s++)
2075                 if (sparse.has(s)) sparse[s].setAscendingGlobalno();
2076 }
2077 void FEM_Mesh::setAbsoluteGlobalno(){
2078         node.setAscendingGlobalno();
2079         for (int e=0;e<elem.size();e++){
2080                 if (elem.has(e)) elem[e].setAscendingGlobalno(nElems(e));
2081         }
2082 }
2083
2084 void FEM_Mesh::copyOldGlobalno(const FEM_Mesh &m) {
2085         node.copyOldGlobalno(m.node);
2086         for (int e=0;e<m.elem.size();e++)
2087                 if (m.elem.has(e) && e<elem.size() && elem.has(e))
2088                         elem[e].copyOldGlobalno(m.elem[e]);
2089         for (int s=0;s<m.sparse.size();s++)
2090                 if (m.sparse.has(s) && s<sparse.size() && sparse.has(s))
2091                         sparse[s].copyOldGlobalno(m.sparse[s]);
2092 }
2093
2094 void FEM_Mesh::clearSharedNodes(){
2095         node.shared.clear();
2096 }
2097
2098 void FEM_Mesh::clearGhostNodes(){
2099         node.clearGhost();
2100 }
2101
2102 void FEM_Mesh::clearGhostElems(){
2103         for(int i=0;i<elem.size();i++){
2104                 if(elem.has(i)){
2105                         elem[i].clearGhost();
2106                 }
2107         }
2108 }
2109
2110 void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType) {
2111         if (type<0 || type>maxType) {
2112                 char msg[1024];
2113                 sprintf(msg,"%s %d is not a valid entity type (it must be between %d and %d)",
2114                         entityType,type, 0, maxType-1);
2115                 FEM_Abort(caller,msg);
2116         }
2117 }
2118 void FEM_Is_NULL(const char *caller,const char *entityType,int type) {
2119         char msg[1024];
2120         sprintf(msg,"%s %d was never set--it cannot now be read",entityType,type);
2121         FEM_Abort(caller,msg);
2122 }
2123
2124 void FEM_Mesh::copyShape(const FEM_Mesh &src)
2125 {
2126         node.copyShape(src.node);
2127         for (int t=0;t<src.elem.size();t++)
2128                 if (src.elem.has(t)) elem.set(t).copyShape(src.elem.get(t));
2129
2130         for (int s=0;s<src.sparse.size();s++)
2131                 if (src.sparse.has(s)) sparse.set(s).copyShape(src.sparse.get(s));
2132
2133         setSymList(src.getSymList());
2134 }
2135
2136 /// Extract a list of our entities:
2137 int FEM_Mesh::getEntities(int *entities) {
2138         int len=0;
2139         entities[len++]=FEM_NODE;
2140         for (int t=0;t<elem.size();t++)
2141                 if (elem.has(t)) entities[len++]=FEM_ELEM+t;
2142         for (int s=0;s<sparse.size();s++)
2143                 if (sparse.has(s)) entities[len++]=FEM_SPARSE+s;
2144         return len;
2145 }
2146
2147
2148 FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead)
2149 {
2150     char fname[256];
2151     static const char *meshFileNames="%s_vp%d_%d.dat";
2152     sprintf(fname, meshFileNames, prefix, nchunks, chunkNo);
2153     FILE *fp = fopen(fname, forRead?"r":"w");
2154     CkPrintf("FEM> %s %s...\n",forRead?"Reading":"Writing",fname);
2155     if(fp==0) {
2156       FEM_Abort(forRead?"FEM: unable to open input file"
2157         :"FEM: unable to create output file.\n");
2158     }
2159     return fp;
2160 }
2161
2162 static inline void read_version()
2163 {
2164     FILE *f = fopen("FEMVERSION", "r");
2165     if (f!=NULL)  {
2166         fscanf(f, "%d", &femVersion);
2167         if (CkMyPe()==0) CkPrintf("FEM> femVersion detected: %d\n", femVersion);
2168         fclose(f);
2169     }
2170 }
2171
2172 static inline void write_version()
2173 {
2174     FILE *f = fopen("FEMVERSION", "w");
2175     if (f!=NULL)  {
2176                 fprintf(f, "%d", femVersion);
2177                 fclose(f);
2178     }
2179 }
2180
2181 FEM_Mesh *FEM_readMesh(const char *prefix,int chunkNo,int nChunks)
2182 {
2183         // find FEM file version number
2184         static int version_checked = 0;
2185         if (!version_checked) {
2186             version_checked=1;
2187             read_version();
2188         }
2189
2190         FEM_Mesh *ret=new FEM_Mesh;
2191         ret->becomeGetting();
2192         FILE *fp = FEM_openMeshFile(prefix,chunkNo,nChunks,true);
2193         PUP::fromTextFile p(fp);
2194         ret->pup(p);
2195         fclose(fp);
2196
2197
2198
2199 #ifdef PRINT_SHARED_NODE_INFO
2200         /** For Abhinav, print out the neighbors for this vp */
2201         
2202         CkPrintf("%d: Finding Neighbors for VP\n", chunkNo);
2203         
2204         FEM_Comm &shared = ret->node.shared; ///<Shared nodes     The type is really an IDXL_Side
2205         
2206         int numNeighborVPs = shared.size();
2207         
2208         CkPrintf("%d: communicates with %d neighbors\n", chunkNo, numNeighborVPs);
2209         
2210         for(int i=0;i<numNeighborVPs;i++){
2211           IDXL_List list = shared.getLocalList(i);
2212             CkPrintf("%d: communicates %d shared nodes with chunk %d\n", chunkNo , list.size() , list.getDest()); 
2213         }
2214 #endif
2215
2216
2217         return ret;
2218 }
2219
2220 void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks)
2221 {
2222         if (chunkNo == 0) {
2223            write_version();
2224         }
2225         FILE *fp = FEM_openMeshFile(prefix,chunkNo,nChunks,false);
2226         PUP::toTextFile p(fp);
2227         m->pup(p);
2228         fclose(fp);
2229 }
2230
2231
2232 // Setup the entity FEM_IS_VALID tables
2233 CDECL void FEM_Mesh_allocate_valid_attr(int fem_mesh, int entity_type){
2234   FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_create_valid_elem");
2235   FEM_Entity *entity = m->lookup(entity_type,"FEM_Mesh_allocate_valid_attr");
2236   entity->allocateValid();
2237 }
2238 FORTRAN_AS_C(FEM_MESH_ALLOCATE_VALID_ATTR,
2239              FEM_Mesh_allocate_valid_attr,
2240              fem_mesh_allocate_valid_attr,
2241              (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )
2242
2243
2244 // mark an entity as valid
2245 CDECL void FEM_set_entity_valid(int mesh, int entityType, int entityIdx){
2246   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_");
2247   FEM_Entity *entity = m->lookup(entityType,"FEM_");
2248   entity->set_valid(entityIdx,true);
2249 }
2250 FORTRAN_AS_C(FEM_SET_ENTITY_VALID,
2251              FEM_set_entity_valid,
2252              fem_set_entity_valid,
2253              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )
2254
2255
2256 // mark an entity as invalid
2257 CDECL void FEM_set_entity_invalid(int mesh, int entityType, int entityIdx){
2258   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2259   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2260   return entity->set_invalid(entityIdx,true);
2261 }
2262 FORTRAN_AS_C(FEM_SET_ENTITY_INVALID,
2263              FEM_set_entity_invalid,
2264              fem_set_entity_invalid,
2265              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )
2266
2267
2268 // Determine if an entity is valid
2269 CDECL int FEM_is_valid(int mesh, int entityType, int entityIdx){
2270   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2271   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2272   return (entity->is_valid(entityIdx) != 0);
2273 }
2274 FORTRAN_AS_C(FEM_IS_VALID,
2275              FEM_is_valid,
2276              fem_is_valid,
2277              (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )
2278
2279
2280 // Count number of valid items for this entity type
2281 CDECL int FEM_count_valid(int mesh, int entityType){
2282   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2283   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2284   return entity->count_valid();
2285 }
2286 FORTRAN_AS_C(FEM_COUNT_VALID,
2287              FEM_count_valid,
2288              fem_count_valid,
2289              (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )
2290
2291
2292 // Set coordinates for some entity's item number idx
2293 void FEM_set_entity_coord2(int mesh, int entityType, int idx, double x, double y){
2294   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2295   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2296   entity->set_coord(idx,x,y);
2297 }
2298 void FEM_set_entity_coord3(int mesh, int entityType, int idx, double x, double y, double z){
2299   FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
2300   FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
2301   entity->set_coord(idx,x,y,z);
2302 }