Home Research Papers Posters Manuals Talks Download People Help Internal

libs/ck-libs/ParFUM/mesh.C

Go to the documentation of this file.

/*
Finite Element Method (FEM) Framework for Charm++
Parallel Programming Lab, Univ. of Illinois 2002
Orion Sky Lawlor, olawlor@acm.org, 12/20/2002

FEM Implementation file: mesh creation and user-data manipulation.
*/

#include "ParFUM_internals.h"

extern int femVersion;

/*** IDXL Interface ****/
FEM_Comm_Holder::FEM_Comm_Holder(FEM_Comm *sendComm, FEM_Comm *recvComm)
    :comm(sendComm,recvComm)
{
    registered=false;
    idx=-1;
}
void FEM_Comm_Holder::registerIdx(IDXL_Chunk *c) {
    assert(!registered);
    IDXL_Chunk *owner=c;
    if (idx!=-1) // had an old index: try to get it back
        idx=owner->addStatic(&comm,idx);
    else //No old index:
        idx=owner->addStatic(&comm);
}
void FEM_Comm_Holder::pup(PUP::er &p) {
    p|idx;
    if (p.isUnpacking() && idx!=-1)
    { // Try to grab the same index we had on our old processor:
        registerIdx(IDXL_Chunk::get("FEM_Comm_Holder::pup"));
    }
}

FEM_Comm_Holder::~FEM_Comm_Holder(void)
{
    if (registered)
    { // Try to unregister from IDXL:
        const char *caller="FEM_Comm_Holder::~FEM_Comm_Holder";
        IDXL_Chunk *owner=IDXL_Chunk::getNULL();
        if (owner) owner->destroy(idx,caller);
    }
}

/******* FEM_Mesh API ******/

static void checkIsSet(int fem_mesh,bool wantSet,const char *caller) {
    if (FEM_Mesh_is_set(fem_mesh)!=wantSet) {
        const char *msg=wantSet?"This mesh (%d) is not a setting mesh":
            "This mesh (%d) is not a getting mesh";
        FEM_Abort(caller,msg,fem_mesh);
    }
}

/* Connectivity: Map calls to appropriate version of FEM_Mesh_data */
CLINKAGE void
FEM_Mesh_conn(int fem_mesh,int entity,
    int *conn, int firstItem, int length, int width)
{
    FEM_Mesh_data(fem_mesh,entity,FEM_CONN, conn, firstItem,length, FEM_INDEX_0, width);
}
FLINKAGE void
FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(int *fem_mesh,int *entity,
    int *conn, int *firstItem,int *length, int *width)
{
    //Can't just call the C version of this routine, because we use 1-based indices:
    FEM_Mesh_data(*fem_mesh,*entity,FEM_CONN, conn, *firstItem-1,*length, FEM_INDEX_1, *width);
}

CLINKAGE void
FEM_Mesh_set_conn(int fem_mesh,int entity,
    const int *conn, int firstItem,int length, int width)
{
    checkIsSet(fem_mesh,true,"FEM_Mesh_set_conn");
    FEM_Mesh_conn(fem_mesh,entity,(int *)conn,firstItem,length,width);
}
FLINKAGE void
FTN_NAME(FEM_MESH_SET_CONN,fem_mesh_set_conn)(int *fem_mesh,int *entity,
    const int *conn, int *firstItem,int *length, int *width)
{
    checkIsSet(*fem_mesh,true,"fem_mesh_set_conn");
    FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(fem_mesh,entity,(int *)conn,firstItem,length,width);
}

CLINKAGE void
FEM_Mesh_get_conn(int fem_mesh,int entity,
    int *conn, int firstItem,int length, int width)
{
    checkIsSet(fem_mesh,false,"FEM_Mesh_get_conn");
    FEM_Mesh_conn(fem_mesh,entity,conn,firstItem,length,width);
}
FLINKAGE void
FTN_NAME(FEM_MESH_GET_CONN,fem_mesh_get_conn)(int *fem_mesh,int *entity,
    int *conn, int *firstItem,int *length, int *width)
{
    checkIsSet(*fem_mesh,false,"fem_mesh_get_conn");
    FTN_NAME(FEM_MESH_CONN,fem_mesh_conn)(fem_mesh,entity,conn,firstItem,length,width);
}


/* Data: map to FEM_Mesh_offset */
CLINKAGE void
FEM_Mesh_data(int fem_mesh,int entity,int attr,
    void *data, int firstItem,int length, int datatype,int width)
{
    IDXL_Layout lo(datatype,width);
    FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,lo);
}

FORTRAN_AS_C(FEM_MESH_DATA,FEM_Mesh_data,fem_mesh_data,
    (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
    (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
)


CLINKAGE void
FEM_Mesh_set_data(int fem_mesh,int entity,int attr,
    const void *data, int firstItem,int length, int datatype,int width)
{
    checkIsSet(fem_mesh,true,"FEM_Mesh_set_data");
    FEM_Mesh_data(fem_mesh,entity,attr,(void *)data,firstItem,length,datatype,width);
}
FORTRAN_AS_C(FEM_MESH_SET_DATA,FEM_Mesh_set_data,fem_mesh_set_data,
    (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
    (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
)

CLINKAGE void
FEM_Mesh_get_data(int fem_mesh,int entity,int attr,
    void *data, int firstItem,int length, int datatype,int width)
{
    checkIsSet(fem_mesh,false,"FEM_Mesh_get_data");
    FEM_Mesh_data(fem_mesh,entity,attr,data,firstItem,length,datatype,width);
}
FORTRAN_AS_C(FEM_MESH_GET_DATA,FEM_Mesh_get_data,fem_mesh_get_data,
    (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *datatype,int *width),
    (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*datatype,*width)
)

CLINKAGE void
FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,
    void *data, int firstItem,int length, IDXL_Layout_t layout)
{
    const char *caller="FEM_Mesh_data_layout";
    FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,
        IDXL_Layout_List::get().get(layout,caller));
}

FORTRAN_AS_C(FEM_MESH_DATA_LAYOUT,FEM_Mesh_data_layout,fem_mesh_data_layout,
    (int *fem_mesh,int *entity,int *attr,void *data,int *firstItem,int *length,int *layout),
    (*fem_mesh,*entity,*attr,data,*firstItem-1,*length,*layout)
)

CLINKAGE void
FEM_Mesh_data_offset(int fem_mesh,int entity,int attr,
    void *data, int firstItem,int length,
    int type,int width, int offsetBytes,int distanceBytes,int skewBytes)
{
    const char *caller="FEM_Mesh_data_offset";
    FEM_Mesh_data_layout(fem_mesh,entity,attr,data,firstItem,length,
        IDXL_Layout(type,width,offsetBytes,distanceBytes,skewBytes));
}
FORTRAN_AS_C(FEM_MESH_DATA_OFFSET,FEM_Mesh_data_offset,fem_mesh_data_offset,
    (int *fem_mesh,int *entity,int *attr,
     void *data,int *firstItem,int *length,
     int *type,int *width,int *offset,int *distance,int *skew),
    (*fem_mesh,*entity,*attr,
     data,*firstItem-1,*length,
     *type,*width,*offset,*distance,*skew)
)


void FEM_Register_array(int fem_mesh,int entity,int attr,
    void *data, int datatype,int width,int firstItem){
    IDXL_Layout lo(datatype,width);
/*  if(attr == FEM_CONN){
        printf("CONN width %d \n",width);
        int len = FEM_Mesh_get_length(fem_mesh,entity);
        int *connd = (int *)data;
        for(int i=0;i<len;i++){
            printf("%d -> (%d %d %d) \n",i+1,connd[3*i],connd[3*i+1],connd[3*i+2]);
        }
    }
    printf("firstItem %d \n",firstItem);*/
        FEM_Register_array_layout(fem_mesh,entity,attr,data,firstItem,lo);
}

void FEM_Register_array_layout(int fem_mesh,int entity,int attr,
    void *data, IDXL_Layout_t layout,int firstItem){
    const char *caller="FEM_Register_array_layout";
    FEM_Register_array_layout(fem_mesh,entity,attr,data,firstItem,
        IDXL_Layout_List::get().get(layout,caller));

}

/*registration api */
CLINKAGE void
FEM_Register_array(int fem_mesh,int entity,int attr,
    void *data, int datatype,int width)
{
    FEM_Register_array(fem_mesh,entity,attr,data,datatype,width,0);
}

CLINKAGE void
FEM_Register_array_layout(int fem_mesh,int entity,int attr,
    void *data, IDXL_Layout_t layout){
    FEM_Register_array_layout(fem_mesh,entity,attr,data,layout,0);
}


CLINKAGE void
FEM_Register_entity(int fem_mesh,int entity,void *data,
        int len,int max,FEM_Mesh_alloc_fn fn) {
        FEM_Register_entity_impl(fem_mesh,entity,data,len,max,fn);
}

FORTRAN_AS_C(FEM_REGISTER_ARRAY,FEM_Register_array,fem_register_array,
    (int *fem_mesh,int *entity,int *attr,void *data,int *datatype,int *width),(*fem_mesh,*entity,*attr,data,*datatype,*width,0))


FORTRAN_AS_C(FEM_REGISTER_ARRAY_LAYOUT,FEM_Register_array_layout,fem_register_array_layout,
    (int *fem_mesh,int *entity,int *attr,void *data,int *layout),(*fem_mesh,*entity,*attr,data,*layout,0))

FORTRAN_AS_C(FEM_REGISTER_ENTITY,FEM_Register_entity,fem_register_entity,
    (int *fem_mesh,int *entity,void *data,int *len,int *max,FEM_Mesh_alloc_fn fn),(*fem_mesh,*entity,data,*len,*max,fn))


// User data API:
CLINKAGE void
FEM_Mesh_pup(int fem_mesh,int dataTag,FEM_Userdata_fn fn,void *data) {
    const char *caller="FEM_Mesh_pup"; FEMAPI(caller);
    FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
    FEM_Userdata_item &i=m->udata.find(dataTag);
    FEM_Userdata_pupfn f(fn,data);
    if (m->isSetting()) i.store(f);
    else /* m->isGetting() */ {
        if (!i.hasStored())
            FEM_Abort(caller,"Never stored any user data at tag %d",dataTag);
        i.restore(f);
    }
}
FORTRAN_AS_C(FEM_MESH_PUP,FEM_Mesh_pup,fem_mesh_pup,
    (int *m,int *t,FEM_Userdata_fn fn,void *data), (*m,*t,fn,data))

// Accessor API:
CLINKAGE void
FEM_Mesh_set_length(int fem_mesh,int entity,int newLength) {
    const char *caller="FEM_Mesh_set_length"; FEMAPI(caller);
    checkIsSet(fem_mesh,true,caller);
    FEM_Entity_lookup(fem_mesh,entity,caller)->setLength(newLength);
}
FORTRAN_AS_C(FEM_MESH_SET_LENGTH,FEM_Mesh_set_length,fem_mesh_set_length,
    (int *fem_mesh,int *entity,int *newLength),
    (*fem_mesh,*entity,*newLength)
)


CLINKAGE int
FEM_Mesh_get_length(int fem_mesh,int entity) {
    const char *caller="FEM_Mesh_get_length"; FEMAPI(caller);
    int len=FEM_Entity_lookup(fem_mesh,entity,caller)->size();
    return len;
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_GET_LENGTH,FEM_Mesh_get_length,fem_mesh_get_length,
    (int *fem_mesh,int *entity),(*fem_mesh,*entity)
)


CLINKAGE void
FEM_Mesh_set_width(int fem_mesh,int entity,int attr,int newWidth) {
    const char *caller="FEM_Mesh_set_width";
    FEMAPI(caller);
    checkIsSet(fem_mesh,true,caller);
    FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->setWidth(newWidth,caller);
}
FORTRAN_AS_C(FEM_MESH_SET_WIDTH,FEM_Mesh_set_width,fem_mesh_set_width,
    (int *fem_mesh,int *entity,int *attr,int *newWidth),
    (*fem_mesh,*entity,*attr,*newWidth)
)

CLINKAGE int
FEM_Mesh_get_width(int fem_mesh,int entity,int attr) {
    const char *caller="FEM_Mesh_get_width";
    FEMAPI(caller);
    return FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->getWidth();
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_GET_WIDTH,FEM_Mesh_get_width,fem_mesh_get_width,
    (int *fem_mesh,int *entity,int *attr),(*fem_mesh,*entity,*attr)
)

CLINKAGE int
FEM_Mesh_get_datatype(int fem_mesh,int entity,int attr) {
    const char *caller="FEM_Mesh_get_datatype";
    FEMAPI(caller);
    return FEM_Attribute_lookup(fem_mesh,entity,attr,caller)->getDatatype();
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_GET_DATATYPE,FEM_Mesh_get_datatype,fem_mesh_get_datatype,
    (int *fem_mesh,int *entity,int *attr),(*fem_mesh,*entity,*attr)
)

CLINKAGE int
FEM_Mesh_is_set(int fem_mesh) /* return 1 if this is a writing mesh */
{
    return (FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_is_get")->isSetting())?1:0;
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_IS_SET,FEM_Mesh_is_set,fem_mesh_is_set,
    (int *fem_mesh),(*fem_mesh)
)

CLINKAGE int
FEM_Mesh_is_get(int fem_mesh) /* return 1 if this is a readable mesh */
{
    return (!FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_is_get")->isSetting())?1:0;
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_IS_GET,FEM_Mesh_is_get,fem_mesh_is_get,
    (int *fem_mesh),(*fem_mesh)
)

CLINKAGE void
FEM_Mesh_become_get(int fem_mesh) /* Make this a readable mesh */
{ FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_become_get")->becomeGetting(); }
FORTRAN_AS_C(FEM_MESH_BECOME_GET,FEM_Mesh_become_get,fem_mesh_become_get, (int *m),(*m))

CLINKAGE void
FEM_Mesh_become_set(int fem_mesh)
{ FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_become_get")->becomeSetting(); }
FORTRAN_AS_C(FEM_MESH_BECOME_SET,FEM_Mesh_become_set,fem_mesh_become_set, (int *m),(*m))


CLINKAGE IDXL_t
FEM_Comm_shared(int fem_mesh,int entity) {
    const char *caller="FEM_Comm_shared";
    FEMAPI(caller);
    if (entity!=FEM_NODE) FEM_Abort(caller,"Only shared nodes supported");
    return FEM_Mesh_lookup(fem_mesh,caller)->node.
        sharedIDXL.getIndex(IDXL_Chunk::get(caller));
}
FORTRAN_AS_C_RETURN(int,
    FEM_COMM_SHARED,FEM_Comm_shared,fem_comm_shared,
    (int *fem_mesh,int *entity),(*fem_mesh,*entity)
)

CLINKAGE IDXL_t
FEM_Comm_ghost(int fem_mesh,int entity) {
    const char *caller="FEM_Comm_ghost";
    FEMAPI(caller);
    FEM_Entity *e=FEM_Entity_lookup(fem_mesh,entity,caller);
    if (e->isGhost()) FEM_Abort(caller,"Can only call FEM_Comm_ghost on real entity type");
    return e->ghostIDXL.getIndex(IDXL_Chunk::get(caller));
}
FORTRAN_AS_C_RETURN(int,
    FEM_COMM_GHOST,FEM_Comm_ghost,fem_comm_ghost,
    (int *fem_mesh,int *entity),(*fem_mesh,*entity)
)


// Internal API:
void FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,
    void *data, int firstItem,int length, const IDXL_Layout &layout)
{
    if (femVersion == 0 && length==0) return;
    const char *caller="FEM_Mesh_data";
    FEMAPI(caller);
    FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
    FEM_Attribute *a=m->lookup(entity,caller)->lookup(attr,caller);

    if (m->isSetting())
        a->set(data,firstItem,length,layout,caller);
    else /* m->isGetting()*/
        a->get(data,firstItem,length,layout,caller);
}

void FEM_Register_array_layout(int fem_mesh,int entity,int attr,void *data,int firstItem,const IDXL_Layout &layout){
    const char *caller="FEM_Register_array";
    FEMAPI(caller);
    FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
    FEM_Entity *e = m->lookup(entity,caller);
    int length = e->size();
    //should actually be a call on the entity
    int max = e->getMax();
    FEM_Attribute *a = e->lookup(attr,caller);
    if(a->getWidth() == 0){
    //  a->setDatatype(layout.type);
        a->setWidth(layout.width);
    }


    if(m->isSetting()){
      a->register_data(data,length,max,layout,caller);
    }else{
      a->get(data,firstItem,length,layout,caller);
      //replace the attribute's data array with the user's data
      a->register_data(data,max/*length*/,max,layout,caller);
    }
}
void FEM_Register_entity_impl(int fem_mesh,int entity,void *args,int len,int max,FEM_Mesh_alloc_fn fn){
    const char *caller = "FEM_Register_entity";
    FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
/*  if(!m->isSetting()){
        CmiAbort("Register entity called on mesh that can't be written into");
    }
*/
    FEM_Entity *e = m->lookup(entity,caller);
    e->setMaxLength(len,max,args,fn);
}

FEM_Entity *FEM_Entity_lookup(int fem_mesh,int entity,const char *caller) {
    return FEM_Mesh_lookup(fem_mesh,caller)->lookup(entity,caller);
}
FEM_Attribute *FEM_Attribute_lookup(int fem_mesh,int entity,int attr,const char *caller) {
    return FEM_Entity_lookup(fem_mesh,entity,caller)->lookup(attr,caller);
}

CLINKAGE int FEM_Mesh_get_entities(int fem_mesh, int *entities) {
    const char *caller="FEM_Mesh_get_entities";
    FEMAPI(caller);
    return FEM_Mesh_lookup(fem_mesh,caller)->getEntities(entities);
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_GET_ENTITIES,FEM_Mesh_get_entities,fem_mesh_get_entities,
    (int *mesh,int *ent), (*mesh,ent)
)

CLINKAGE int FEM_Mesh_get_attributes(int fem_mesh,int entity,int *attributes) {
    const char *caller="FEM_Mesh_get_attributes";
    FEMAPI(caller);
    return FEM_Entity_lookup(fem_mesh,entity,caller)->getAttrs(attributes);
}
FORTRAN_AS_C_RETURN(int,
    FEM_MESH_GET_ATTRIBUTES,FEM_Mesh_get_attributes,fem_mesh_get_attributes,
    (int *mesh,int *ent,int *attrs), (*mesh,*ent,attrs)
)

/************** FEM_Attribute ****************/

CLINKAGE const char *FEM_Get_datatype_name(int datatype,char *storage) {
    switch(datatype) {
    case FEM_BYTE: return "FEM_BYTE";
    case FEM_INT: return "FEM_INT";
    case FEM_FLOAT: return "FEM_FLOAT";
    case FEM_DOUBLE: return "FEM_DOUBLE";
    case FEM_INDEX_0: return "FEM_INDEX_0";
    case FEM_INDEX_1: return "FEM_INDEX_1";
    };
    sprintf(storage,"unknown datatype code (%d)",datatype);
    return storage;
}

CLINKAGE const char *FEM_Get_attr_name(int attr,char *storage)
{
    if (attr<FEM_ATTRIB_TAG_MAX)
    { //It's a user tag:
        sprintf(storage,"FEM_DATA+%d",attr-FEM_DATA);
        return storage;
    }
    switch(attr) {
    case FEM_CONN: return "FEM_CONN";
    case FEM_SPARSE_ELEM: return "FEM_SPARSE_ELEM";
    case FEM_COOR: return "FEM_COOR";
    case FEM_GLOBALNO: return "FEM_GLOBALNO";
    case FEM_PARTITION: return "FEM_PARTITION";
    case FEM_SYMMETRIES: return "FEM_SYMMETRIES";
    case FEM_NODE_PRIMARY: return "FEM_NODE_PRIMARY";
    case FEM_NODE_ELEM_ADJACENCY: return "FEM_NODE_ELEM_ADJACENCY";
    case FEM_NODE_NODE_ADJACENCY: return "FEM_NODE_NODE_ADJACENCY";
    case FEM_ELEM_ELEM_ADJACENCY: return "FEM_ELEM_ELEM_ADJACENCY";
    case FEM_ELEM_ELEM_ADJ_TYPES: return "FEM_ELEM_ELEM_ADJ_TYPES";
    case FEM_IS_VALID_ATTR: return "FEM_IS_VALID_ATTR";
    case FEM_MESH_SIZING: return "FEM_MESH_SIZING";

    default: break;
    };
    sprintf(storage,"unknown attribute code (%d)",attr);
    return storage;
}

//Abort with a nice error message saying:
// Our <field> was previously set to <cur>; it cannot now be <operation> <next>
void FEM_Attribute::bad(const char *field,bool forRead,int cur,int next,const char *caller) const
{
    char nameStorage[256];
    const char *name=FEM_Get_attr_name(attr,nameStorage);
    char errBuf[1024];
    const char *cannotBe=NULL;
    if (forRead) {
        if (cur==-1) {
            sprintf(errBuf,"The %s %s %s was never set-- it cannot now be read",
                e->getName(),name,field);
        }
        else /* already had a value */
            cannotBe="read as";
    }
    else /* for write */ {
        cannotBe="set to";
    }
    if (cannotBe!=NULL) /* Use standard ... <something> cannot be <something>'d... error message */
        sprintf(errBuf,"The %s %s %s was previously set to %d; it cannot now be %s %d",
            e->getName(),name,field,cur,cannotBe,next);

    FEM_Abort(caller,errBuf);
}


FEM_Attribute::FEM_Attribute(FEM_Entity *e_,int attr_)
        :e(e_),ghost(0),attr(attr_),width(0),datatype(-1), allocated(false)
{
    tryAllocate();
    if (femVersion == 0) width=-1;
}
void FEM_Attribute::pup(PUP::er &p) {
    // e, attr, and ghost are always set by the constructor
    p|width;
    if (p.isUnpacking() && femVersion > 0 && width<0)  width=0;
    p|datatype;
    if (p.isUnpacking()) tryAllocate();
}
void FEM_Attribute::pupSingle(PUP::er &p, int pupindx) {
  return;
}
FEM_Attribute::~FEM_Attribute() {}

void FEM_Attribute::setLength(int next,const char *caller) {
    int cur=getLength();
    if (next==cur) return; //Already set--nothing to do
    if (cur>0) bad("length",false,cur,next, caller);
    e->setLength(next);
    tryAllocate();
}

void FEM_Attribute::setWidth(int next,const char *caller) {
    int cur=getWidth();
    if (next==cur) return; //Already set--nothing to do
    if (cur>0) bad("width",false,cur,next, caller);
    width=next;
    tryAllocate();
    if (ghost) ghost->setWidth(width,caller);
}

void FEM_Attribute::setDatatype(int next,const char *caller) {
    int cur=getDatatype();
    if (next==cur) return; //Already set--nothing to do
    if (cur!=-1) bad("datatype",false,cur,next, caller);
    datatype=next;
    tryAllocate();
    if (ghost) ghost->setDatatype(datatype,caller);
}

void FEM_Attribute::copyShape(const FEM_Attribute &src) {
    setWidth(src.getWidth());
    if (src.getDatatype()!=-1)
      setDatatype(src.getDatatype()); //Automatically calls tryAllocate
}
void FEM_Attribute::set(const void *src, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller)
{
    if (firstItem!=0) { /* If this isn't the start... */
        if (length!=1) /* And we're not setting one at a time */
            CmiAbort("FEM_Mesh_data: unexpected firstItem");
    }

    if (femVersion == 0 && getRealLength() == -1) setLength(length);
    else if (getLength()==0) setLength(length);
    else if (length!=1 && length!=getLength())
        bad("length",false,getLength(),length, caller);

    int width=layout.width;
    if (femVersion==0 && getRealWidth()==-1) setWidth(width);
    else if (getWidth()==0) setWidth(width);
    else if (width!=getWidth())
        bad("width",false,getWidth(),width, caller);

    int datatype=layout.type;
    if (getDatatype()==-1) setDatatype(datatype);
    else if (datatype!=getDatatype())
        bad("datatype",false,getDatatype(),datatype, caller);

    /* Assert: our storage should be allocated now.
       Our subclass will actually copy user data */
}

void FEM_Attribute::get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller)  const
{
    if (length==0) return; //Nothing to get
    if (length!=1 && length!=getLength())
        bad("length",true,getLength(),length, caller);

    int width=layout.width;
    if (width!=getWidth())
        bad("width",true,getWidth(),width, caller);

    int datatype=layout.type;
    if (datatype!=getDatatype())
        bad("datatype",true,getDatatype(),datatype, caller);

    /* our subclass will actually copy into user data */
}

/*check if the layout is the same as earlier */

void FEM_Attribute::register_data(void *user, int length,int max,
    const IDXL_Layout &layout, const char *caller){

        int width=layout.width;
        if (femVersion == 0 && getRealWidth()==-1) setWidth(width);
        else if (getWidth()==0){
            setWidth(width);
        }else{
            if (width!=getWidth()){
                bad("width",false,getWidth(),width, caller);
            }
        }

        int datatype=layout.type;
        if (getDatatype()==-1){
            setDatatype(datatype);
        }else{
            if (datatype!=getDatatype()){
                bad("datatype",false,getDatatype(),datatype, caller);
            }
        }

}

//Check if all three of length, width, and datatype are set.
// If so, call allocate.
void FEM_Attribute::tryAllocate(void) {
    int lenNull, widthNull;
        if (femVersion == 0) {
      // version 0 takes -1 as empty
      lenNull = (getRealLength()==-1);
      widthNull = (getRealWidth()==-1);
    }
    else {
      lenNull = (getLength()==0);
      widthNull = (getWidth()==0);
    }
    if ((!allocated) && !lenNull && !widthNull && getDatatype()!=-1) {
      allocated=true;
        allocate(getMax(),getWidth(),getDatatype());
    }
}

/*********************** DataAttribute *******************/
FEM_DataAttribute::FEM_DataAttribute(FEM_Entity *e,int myAttr)
    :FEM_Attribute(e,myAttr),
     char_data(0),int_data(0),float_data(0),double_data(0)
{
}
void FEM_DataAttribute::pup(PUP::er &p) {
    super::pup(p);
    switch(getDatatype()) {
    case -1: /* not allocated yet */ break;
    case FEM_BYTE:
      if (char_data) {
        /*if(p.isSizing()) {
          char_data->setRowLen(getEntity()->size());
          }*/
        char_data->pup(p);
      }
      break;
    case FEM_INT:
      if (int_data) {
        /*if(p.isSizing()) {
          int_data->setRowLen(getEntity()->size());
          }*/
        int_data->pup(p);
      }
      break;
    case FEM_FLOAT:
      if (float_data) {
        /*if(p.isSizing()) {
          float_data->setRowLen(getEntity()->size());
          }*/
        float_data->pup(p);
      }
      break;
    case FEM_DOUBLE:
      if (double_data) {
        /*if(p.isSizing()) {
          double_data->setRowLen(getEntity()->size());
          }*/
        double_data->pup(p);
      }
      break;
    default: CkAbort("Invalid datatype in FEM_DataAttribute::pup");
    }
}
void FEM_DataAttribute::pupSingle(PUP::er &p, int pupindx) {
    super::pupSingle(p,pupindx);
    switch(getDatatype()) {
    case -1: /* not allocated yet */ break;
    case FEM_BYTE:   if (char_data) char_data->pupSingle(p,pupindx); break;
    case FEM_INT:    if (int_data) int_data->pupSingle(p,pupindx); break;
    case FEM_FLOAT:  if (float_data) float_data->pupSingle(p,pupindx); break;
    case FEM_DOUBLE: if (double_data) double_data->pupSingle(p,pupindx); break;
    default: CkAbort("Invalid datatype in FEM_DataAttribute::pupSingle");
    }
}
FEM_DataAttribute::~FEM_DataAttribute() {
    if (char_data) delete char_data;
    if (int_data) delete int_data;
    if (float_data) delete float_data;
    if (double_data) delete double_data;

}

template <class T>
inline void setTableData(const void *user, int firstItem, int length,
    IDXL_LAYOUT_PARAM, AllocTable2d<T> *table)
{
    for (int r=0;r<length;r++) {
        T *tableRow=table->getRow(firstItem+r);
        for (int c=0;c<width;c++)
            tableRow[c]=IDXL_LAYOUT_DEREF(T,user,r,c);
    }
}

template <class T>
inline void getTableData(void *user, int firstItem, int length,
    IDXL_LAYOUT_PARAM, const AllocTable2d<T> *table)
{
    for (int r=0;r<length;r++) {
        const T *tableRow=table->getRow(firstItem+r);
        for (int c=0;c<width;c++)
            IDXL_LAYOUT_DEREF(T,user,r,c)=tableRow[c];
    }

}

void FEM_DataAttribute::set(const void *u, int f,int l,
        const IDXL_Layout &layout, const char *caller)
{
    super::set(u,f,l,layout,caller);
    switch(getDatatype()) {
    case FEM_BYTE:  setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),char_data); break;
    case FEM_INT: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),int_data); break;
    case FEM_FLOAT: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),float_data); break;
    case FEM_DOUBLE: setTableData(u,f,l,IDXL_LAYOUT_CALL(layout),double_data); break;
    }
}

void FEM_DataAttribute::get(void *u, int f,int l,
        const IDXL_Layout &layout, const char *caller) const
{
    super::get(u,f,l,layout,caller);
    switch(getDatatype()) {
    case FEM_BYTE:  getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),char_data); break;
    case FEM_INT: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),int_data); break;
    case FEM_FLOAT: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),float_data); break;
    case FEM_DOUBLE: getTableData(u,f,l,IDXL_LAYOUT_CALL(layout),double_data); break;
    }
}

void FEM_DataAttribute::register_data(void *u,int l,int max,
        const IDXL_Layout &layout, const char *caller)
{
    super::register_data(u,l,max,layout,caller);
    switch(getDatatype()){
        case FEM_BYTE: char_data->register_data((unsigned char *)u,l,max); break;
        case FEM_INT:    int_data->register_data((int *)u,l,max); break;
        case FEM_FLOAT: float_data->register_data((float *)u,l,max);break;
        case FEM_DOUBLE: double_data->register_data((double *)u,l,max);break;
    }
}

template<class T>
inline AllocTable2d<T> *allocTablePtr(AllocTable2d<T> *src,int len,int wid) {
    if (src==NULL) src=new AllocTable2d<T>;
    src->allocate(wid,len);
    return src;
}
void FEM_DataAttribute::allocate(int l,int w,int datatype)
{
    switch(datatype) {
    case FEM_BYTE:  char_data=allocTablePtr(char_data,l,w); break;
    case FEM_INT: int_data=allocTablePtr(int_data,l,w); break;
    case FEM_FLOAT: float_data=allocTablePtr(float_data,l,w); break;
    case FEM_DOUBLE: double_data=allocTablePtr(double_data,l,w); break;
    default: CkAbort("Invalid datatype in FEM_DataAttribute::allocate");
    };
}
void FEM_DataAttribute::copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity)
{
    const FEM_DataAttribute *dsrc=(const FEM_DataAttribute *)&src;
    switch(getDatatype()) {
    case FEM_BYTE:  char_data->setRow(dstEntity,dsrc->char_data->getRow(srcEntity)); break;
    case FEM_INT:
            int_data->setRow(dstEntity,dsrc->int_data->getRow(srcEntity)); break;
    case FEM_FLOAT: float_data->setRow(dstEntity,dsrc->float_data->getRow(srcEntity)); break;
    case FEM_DOUBLE: double_data->setRow(dstEntity,dsrc->double_data->getRow(srcEntity)); break;
    }
}

template<class T>
inline void interpolateAttrs(AllocTable2d<T> *data,int A,int B,int D,double frac,int width){
  T *rowA = data->getRow(A);
  T *rowB = data->getRow(B);
  T *rowD = data->getRow(D);
  for(int i=0;i<width;i++){
    double val = (double )rowA[i];
    val *= (frac);
    val += (1-frac) *((double )rowB[i]);
    rowD[i] = (T )val;
  }
}

template<class T>
inline void interpolateAttrs(AllocTable2d<T> *data,int *iNodes,int rNode,int k,int width){
  T *row[8];
  for (int i=0; i<k; i++) {
    row[i] = data->getRow(iNodes[i]);
  }
  T *rowR = data->getRow(rNode);
  for(int i=0;i<width;i++){
    double val = 0.0;
    for (int j=0; j<k; j++) {
      val += (double)row[j][i];
    }
    val = val/k;
    rowR[i] = (T )val;
  }
}

template<class T>
inline void minAttrs(AllocTable2d<T> *data,int A,int B,int D,double frac,int width){
  T *rowA = data->getRow(A);
  T *rowB = data->getRow(B);
  T *rowD = data->getRow(D);
  for(int i=0;i<width;i++){
    if(rowA[i] < rowB[i]){
      rowD[i] = rowA[i];
    }else{
      rowD[i] = rowB[i];
    }
  }
}

template<class T>
inline void minAttrs(AllocTable2d<T> *data,int *iNodes,int rNode,int k,int width){
  T *row[8];
  for (int i=0; i<k; i++) {
    row[i] = data->getRow(iNodes[i]);
  }
  T *rowR = data->getRow(rNode);
  for(int i=0;i<width;i++){
    rowR[i] = row[0][i];
    for (int j=1; j<k; j++) {
      if (row[j][i] < rowR[i]) {
    rowR[i] = row[j][i];
      }
    }
  }
}

void FEM_DataAttribute::interpolate(int A,int B,int D,double frac){
  switch(getDatatype()){
  case FEM_BYTE:
    minAttrs(char_data,A,B,D,frac,getWidth());
    break;
  case FEM_INT:
    minAttrs(int_data,A,B,D,frac,getWidth());
    break;
  case FEM_FLOAT:
    interpolateAttrs(float_data,A,B,D,frac,getWidth());
    break;
  case FEM_DOUBLE:
    interpolateAttrs(double_data,A,B,D,frac,getWidth());
    break;
  }
}

void FEM_DataAttribute::interpolate(int *iNodes,int rNode,int k){
  switch(getDatatype()){
  case FEM_BYTE:
    minAttrs(char_data,iNodes,rNode,k,getWidth());
    break;
  case FEM_INT:
    minAttrs(int_data,iNodes,rNode,k,getWidth());
    break;
  case FEM_FLOAT:
    interpolateAttrs(float_data,iNodes,rNode,k,getWidth());
    break;
  case FEM_DOUBLE:
    interpolateAttrs(double_data,iNodes,rNode,k,getWidth());
    break;
  }
}

/*********************** FEM_IndexAttribute *******************/
FEM_IndexAttribute::Checker::~Checker() {}

FEM_IndexAttribute::FEM_IndexAttribute(FEM_Entity *e,int myAttr,FEM_IndexAttribute::Checker *checker_)
    :FEM_Attribute(e,myAttr), idx(0,0,-1), checker(checker_)
{
    setDatatype(FEM_INT);
}
void FEM_IndexAttribute::pup(PUP::er &p) {
    super::pup(p);
    /*if(p.isSizing()) {
      idx.setRowLen(getEntity()->size());
      }*/
    p|idx;
}
void FEM_IndexAttribute::pupSingle(PUP::er &p, int pupindx) {
    super::pupSingle(p,pupindx);
    idx.pupSingle(p,pupindx);
}
FEM_IndexAttribute::~FEM_IndexAttribute() {
    if (checker) delete checker;
}

void FEM_IndexAttribute::allocate(int length,int width,int datatype)
{
    idx.allocate(width,length);
}

static int type2base(int base_type,const char *caller) {
    if (base_type==FEM_INDEX_0) return 0;
    if (base_type==FEM_INDEX_1) return 1;
    FEM_Abort(caller,"You must use the datatype FEM_INDEX_0 or FEM_INDEX_1 with FEM_CONN, not %d",
        base_type);
    return 0; //< for whining compilers
}

void setIndexTableData(const void *user, int firstItem, int length,
    IDXL_LAYOUT_PARAM, AllocTable2d<int> *table,int indexBase)
{
    for (int r=0;r<length;r++) {
        int *tableRow=table->getRow(firstItem+r);
        for (int c=0;c<width;c++)
            tableRow[c]=IDXL_LAYOUT_DEREF(int,user,r,c)-indexBase;
    }
}

void getIndexTableData(void *user, int firstItem, int length,
    IDXL_LAYOUT_PARAM, const AllocTable2d<int> *table,int indexBase)
{
    for (int r=0;r<length;r++) {
        const int *tableRow=table->getRow(firstItem+r);
        for (int c=0;c<width;c++)
            IDXL_LAYOUT_DEREF(int,user,r,c)=tableRow[c]+indexBase;
    }
}

void FEM_IndexAttribute::set(const void *src, int firstItem,int length,
        const IDXL_Layout &layout,const char *caller)
{
    IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
    super::set(src,firstItem,length,lo,caller);

    int indexBase=type2base(layout.type,caller);
    setIndexTableData(src,firstItem,length,IDXL_LAYOUT_CALL(layout),&idx,indexBase);

    if (checker)
        for (int r=0;r<length;r++)
            checker->check(firstItem+r,idx,caller);
}

void FEM_IndexAttribute::get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout,const char *caller) const
{
    IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
    super::get(dest,firstItem,length,lo,caller);

    int indexBase=type2base(layout.type,caller);
    getIndexTableData(dest,firstItem,length,IDXL_LAYOUT_CALL(layout),&idx,indexBase);
}

void FEM_IndexAttribute::register_data(void *user, int length,int max,
        const IDXL_Layout &layout, const char *caller){
    IDXL_Layout lo=layout; lo.type=FEM_INT; //Pretend it's always int data, not INDEX
    super::register_data(user,length,max,lo,caller);

    idx.register_data((int *)user,length,max);
}

void FEM_IndexAttribute::copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity)
{
    const FEM_IndexAttribute *csrc=(const FEM_IndexAttribute *)&src;
    idx.setRow(dstEntity,csrc->idx.getRow(srcEntity));
}

/*************FEM_VarIndexAttribute***************/

FEM_VarIndexAttribute::FEM_VarIndexAttribute(FEM_Entity *e,int myAttr)
    :FEM_Attribute(e,myAttr)
{
  oldlength = 0;
    allocate(getMax(),getWidth(),getDatatype());
    setDatatype(FEM_INT);
}

void FEM_VarIndexAttribute::pup(PUP::er &p){
    super::pup(p);
    p | idx;
    p|oldlength;
}

void FEM_VarIndexAttribute::pupSingle(PUP::er &p, int pupindx){
    super::pupSingle(p,pupindx);
    p|idx[pupindx];
}

void FEM_VarIndexAttribute::set(const void *src,int firstItem,int length,
        const IDXL_Layout &layout,const char *caller){
        printf("set not yet implemented for FEM_VarIndexAttribute \n");
}

void FEM_VarIndexAttribute::get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller) const{
     printf("get not yet implemented for FEM_VarIndexAttribute \n");

}

void FEM_VarIndexAttribute::copyEntity(int dstEntity,const FEM_Attribute &_src,int srcEntity){
    FEM_VarIndexAttribute &src = (FEM_VarIndexAttribute &)_src;
    const CkVec<CkVec<ElemID> > &srcTable = src.get();
        CkVec<ElemID> temp = srcTable[srcEntity];
    idx.insert(dstEntity, temp);
}

void FEM_VarIndexAttribute::print(){
    for(int i=0;i<idx.size();i++){
        printf("%d -> ",i);
        for(int j=0;j<idx[i].size();j++){
            printf("(%d) ",((idx[i])[j]).id);
        }
        printf("\n");
    }
}

int FEM_VarIndexAttribute::findInRow(int row, const ElemID &data){
    if(row >= idx.length()){
        return -1;
    }
    CkVec<ElemID> &rowVec = idx[row];
    for(int i=0;i<rowVec.length();i++){
        if(data == rowVec[i]){
            return i;
        }
    }
    return -1;
}

/********************** Entity **************************/

CLINKAGE const char *FEM_Get_entity_name(int entity,char *storage)
{
    char *dest=storage;
    if (entity<FEM_ENTITY_FIRST || entity>=FEM_ENTITY_LAST) {
        sprintf(dest,"unknown entity code (%d)",entity);
    }
    else {
        if (entity>FEM_ENTITY_FIRST+FEM_GHOST) {
            sprintf(dest,"FEM_GHOST+");
            dest+=strlen(dest); /* we want "FEM_GHOST+foo" */
            entity-=FEM_GHOST;
        }
        if (entity==FEM_NODE)
            sprintf(dest,"FEM_NODE");
        else if (entity>=FEM_SPARSE)
            sprintf(dest,"FEM_SPARSE+%d",entity-FEM_SPARSE);
        else /* entity>=FEM_ELEM */
            sprintf(dest,"FEM_ELEM+%d",entity-FEM_ELEM);
    }
    return storage;
}

FEM_Entity::FEM_Entity(FEM_Entity *ghost_) //Default constructor
  :length(0), max(0),ghost(ghost_), coord(0), sym(0), globalno(0), valid(0), meshSizing(0),
   ghostIDXL(ghost?&ghostSend:NULL, ghost?&ghost->ghostRecv:NULL),resize(NULL),
   invalidList(NULL),first_invalid(0),last_invalid(0),invalidListLen(0),invalidListAllLen(0)
{
    //No attributes initially
    if (femVersion == 0) {
        length=-1;
        max=-1;
    }
}
void FEM_Entity::pup(PUP::er &p) {
    p|length;
    p|max;
    if (p.isUnpacking() && femVersion > 0 && length<0)  length=0;
    p.comment(" Ghosts to send out: ");
    ghostSend.pup(p);
    p.comment(" Ghosts to recv: ");
    ghostRecv.pup(p);
    p.comment(" Ghost IDXL tag: ");
    ghostIDXL.pup(p);

    int nAttributes=attributes.size();
    p|nAttributes;
    for (int a=0;a<nAttributes;a++)
    {
    /* Beautiful hack: the FEM_Attribute objects are normally hideously cross-linked
       with their owning classes.  Thus instead of trying to rebuild the FEM_Attributes
       from scratch here, we just use the existing "lookup" method to demand-create those
       that need it, or just *find* those that are already there.

       This is a much better fit for this situation than using the general PUP::able.
     */
      int attr=0; // The attribute type we're pupping
      FEM_Attribute *r=NULL;
      if (!p.isUnpacking()) { //Send side: we already know the source
        r=attributes[a];
        attr=r->getAttr();
      }
      p|attr;
      if (p.isUnpacking()) { //Recv side: create (or recycle) the destination
        r=lookup(attr,"FEM_Entity::pup");
      }

      { //Put the human-readable attribute name in the output file:
        char attrNameStorage[256];
        p.comment(FEM_Get_attr_name(attr,attrNameStorage));
      }
      r->pup(p);
    }
    p|first_invalid;
    p|last_invalid;
    p|invalidListLen;
    p|invalidListAllLen;
    if(p.isUnpacking()) {
      delete[] invalidList; //was just created in allocateValid
      if(invalidListAllLen>0) invalidList = new int[invalidListAllLen];
    }
    if(invalidListAllLen>0) {
      PUParray(p, (int*)invalidList, invalidListAllLen);
    }
    if(p.isDeleting()) {
      delete[] invalidList;
    }

    if (ghost!=NULL) {
        p.comment(" ---- Ghost attributes ---- ");
        ghost->pup(p);
    }
}

FEM_Entity::~FEM_Entity()
{
    delete ghost;
    for (int a=0;a<attributes.size();a++)
        delete attributes[a];
}

void FEM_Entity::copyShape(const FEM_Entity &src) {
    for (int a=0;a<src.attributes.size();a++)
    { // We need each of his attributes:
        const FEM_Attribute *Asrc=src.attributes[a];
        FEM_Attribute *Adst=lookup(Asrc->getAttr(),"FEM_Entity::copyShape");
        Adst->copyShape(*Asrc);
    }
    if (ghost) ghost->copyShape(*src.ghost);
}

void FEM_Entity::setLength(int newlen, bool f)
{
    if (!resize) {
        if (size() != newlen) {
            length = newlen;
            // Each of our attributes need to be expanded for our new length:
            for (int a=0; a<attributes.size(); a++) {
                attributes[a]->reallocate();
            }
        }
    }
    else {
      int otherlen = length;
      if(!f) {
    otherlen = newlen;
      }
      length = newlen;
        if (length > max) {
            if (max > 4) {
                max = max + (max >> 2);
            } else {
                max = max + 10;
            }
            for (int a=0;a<attributes.size();a++){
                int code = attributes[a]->getAttr();
                if (!(code <= FEM_ATTRIB_TAG_MAX ||
                            code == FEM_CONN || code == FEM_COORD ||
              code == FEM_BOUNDARY)) { //user can store bpundary & connectivity also
                    attributes[a]->reallocate();
                }
            }
            // call resize with args max n;
            //CkPrintf("Resize called \n");
            resize(args,&otherlen,&max); //resets length to otherlen
            //resize(args,&length,&max);
        length = newlen;
        }
    }
}

void FEM_Entity::reserveLength(int newlen)
{
  if (newlen > max) {
    max = newlen;
    for (int a=0;a<attributes.size();a++){
      attributes[a]->reallocate();
    }
  }
}

void FEM_Entity::allocateValid(void) {
  if (!valid){
    valid=new FEM_DataAttribute(this,FEM_IS_VALID_ATTR);
    add(valid);
    valid->setWidth(1); //Only 1 flag per node
    valid->setLength(size());
    valid->setDatatype(FEM_BYTE);
    valid->reallocate();

    // Set all to valid initially
    for(int i=0;i<size();i++) {
      valid->getChar()(i,0)=1;
    }
    if(true) { //maintains a list of invalid elements. FASTER
      invalidListLen=0; invalidListAllLen=16;
      invalidList = new int[invalidListAllLen];
      //its ok to waste a few bytes for an entity, especially when it can
      //bring up the performance of the average case
      for(int i=0; i<invalidListAllLen; i++) {
    invalidList[i] = -1; //means this position is empty
      }
    }
    else {
      first_invalid = last_invalid = 0;
    }
  }
}

void FEM_Entity::set_valid(int idx, bool isNode){
  if(true) { //maintains a list of invalid elements. FASTER
    CkAssert(idx < size() && idx >=0);
    valid->getChar()(idx,0)=1;
    if(invalidListLen>0) {
      if(invalidList[invalidListLen-1]==idx) {
    invalidListLen--; //pull out the last entry from the invalid list
      }
    }
  }
  else {
    int size1 = size();
    if(size1==0) {
      valid->getChar()(idx,0)=1;
    }
    else {
      CkAssert(idx < size1 && idx >=0 && first_invalid<=last_invalid);
      valid->getChar()(idx,0)=1;
      if(idx == first_invalid) {
    // Move first_invalid to the next invalid entry
    while((first_invalid<last_invalid) && is_valid(first_invalid)){
      first_invalid++;
    }
      }
      else if(idx == last_invalid) {
    // Move last_invalid to the previous invalid entry
    while((first_invalid<last_invalid) && is_valid(last_invalid)) {
      last_invalid--;
    }
      }
    }
    // If we have no invalid elements left, then put both pointers to 0
    if( first_invalid == last_invalid && is_valid(first_invalid) )
      first_invalid = last_invalid = 0;
  }
}

void FEM_Entity::set_invalid(int idx, bool isNode){
  if(true) { //maintains a list of invalid elements. FASTER
    CkAssert(idx < size() && idx >=0);
    valid->getChar()(idx,0)=0;
    if(invalidListLen==invalidListAllLen) {
      invalidListAllLen = invalidListAllLen<<1;
      int* tmp = invalidList;
      invalidList = new int[invalidListAllLen];
      memcpy(invalidList,tmp,invalidListLen*sizeof(int));
      delete[] tmp;
    }
    invalidList[invalidListLen] = idx;
    invalidListLen++;
  }
  else {
    CkAssert(idx < size() && idx >=0 && first_invalid<=last_invalid);
    valid->getChar()(idx,0)=0;

    // If there are currently no invalid entities
    if(first_invalid==0 && last_invalid==0 && is_valid(0)){
      first_invalid = last_invalid = idx;
      return;
    }
    if(idx < first_invalid){
      first_invalid = idx;
    }
    if(idx > last_invalid){
      last_invalid = idx;
    }
    // TODO:
    // We should probably have an algorithm for shrinking the entire attribute
    // array if we invalidate the final element. In this case we should scan backwards
    // to find the largest indexed valid entity and resize down to it.
    //
    // It may be necessary to modify the idxl lists if we do this type of shrinking.
    // Someone needs to confirm whether that is necessary. If not, then it should be
    // simple to allow shrinking of the number of nodes or elements.
  }
}


int FEM_Entity::is_valid(int idx){
    CkAssert(idx < size() && idx >=0);
    return valid->getChar()(idx,0);
}

int FEM_Entity::is_valid_any_idx(int idx){
  if(idx == -1 || idx >= size())
    return false;
  else if(idx < 0)
    return ghost->is_valid_any_idx(FEM_To_ghost_index(idx));
  else
    return valid->getChar()(idx,0);
}

int FEM_Entity::is_valid_nonghost_idx(int idx){
  if(idx < 0 || idx >= size())
    return false;
  else
    return valid->getChar()(idx,0);
}

int FEM_Entity::count_valid(){
  int count=0;
  for(int i=0;i<size();i++)
    if(is_valid(i)) count++;
  return count;
}

void FEM_Entity::set_all_invalid(){
    // Set all to invalid
    for(int i=0;i<size();i++) {
      valid->getChar()(i,0)=0;
    }
}


int FEM_Entity::get_next_invalid(FEM_Mesh *m, bool isNode, bool isGhost){
  int retval=0;
//  if(true) { //maintains a list of invalid elements. FASTER
    if(invalidListLen>0) {
      retval = invalidList[invalidListLen-1];
    }
    else {
      retval = size();
      setLength(retval+1,true);
    }
//  }
/*  else {
    int size1 = size();
    if(size1==0) { //special case because for size=0, first_invalid=last_invalid=0
      retval= 0;
      setLength(1);
    }
    else {
      CkAssert(!is_valid(first_invalid) || first_invalid==0);
      // if we have an invalid entity to return
      bool flag1 = false;
      if(!is_valid(first_invalid)){
    retval = first_invalid;
    if(isNode && !isGhost) { //it is a node & the entity is not a ghost entity
      while(retval <= last_invalid) {
        if(!is_valid(retval)) {
          if(m->getfmMM()->fmLockN[retval].haslocks()) {
        retval++;
          }
          else if(hasConn(retval)) { //has some connectivity
        retval++;
          }
          else {
        flag1 = true;
        break;
          }
        }
        else retval++;
      }
    }
    else if(isNode) {
      while(retval <= last_invalid) {
        if(!is_valid(retval)) {
          if(hasConn(retval)) { //has some connectivity
        retval++;
          }
          else {
        flag1 = true;
        break;
          }
        }
        else retval++;
      }
    }
    else{
      // resize array and return new entity
      flag1 = true;
    }
      }
      if(!flag1) {
    retval = size1;
    setLength(retval+1);
      }
    }
  }*/
  set_valid(retval,isNode);
  return retval;
}

void FEM_Entity::setMaxLength(int newLen,int newMaxLen,void *pargs,FEM_Mesh_alloc_fn fn){
  //CkPrintf("resize fn %p \n",fn);
    max = newMaxLen;
    resize = fn;
    args = pargs;
    setLength(newLen);
}


void FEM_Entity::copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity) {
    FEM_Entity *dp=this; //Destination entity
    const FEM_Entity *sp=&src;
    for (int a=0;a<sp->attributes.size();a++)
    { //We copy each of his attributes:
        const FEM_Attribute *Asrc=sp->attributes[a];
        FEM_Attribute *Adst=dp->lookup(Asrc->getAttr(),"FEM_Entity::copyEntity");
        Adst->copyEntity(dstEntity,*Asrc,srcEntity);
    }
}

int FEM_Entity::push_back(const FEM_Entity &src,int srcEntity) {
    int dstEntity=size();
    setLength(dstEntity+1);
    copyEntity(dstEntity,src,srcEntity);
    return dstEntity;
}

void FEM_Entity::add(FEM_Attribute *attribute) {
    if (ghost!=NULL)
    { //Look up (or create) the ghost attribute, too:
        attribute->setGhost(ghost->lookup(attribute->getAttr(),"FEM_Entity::add"));
    }
    attributes.push_back(attribute);
}

FEM_Attribute *FEM_Entity::lookup(int attr,const char *caller) {
    //Try to find an existing attribute (FIXME: keep attributes in a map, to speed this up)
    for (int a=0;a<attributes.size();a++) {
        if (attributes[a]->getAttr()==attr)
            return attributes[a];
    }

    //If we get here, no existing attribute fits the bill: create one
    create(attr,caller);

    // If create did its job, the next lookup should succeed:
    return lookup(attr,caller);
}

void FEM_Entity::create(int attr,const char *caller) {
  if (attr<=FEM_ATTRIB_TAG_MAX) //It's a valid user data tag
    add(new FEM_DataAttribute(this,attr));
  else if (attr==FEM_COORD)
    allocateCoord();
  else if (attr==FEM_SYMMETRIES)
    allocateSym();
  else if (attr==FEM_GLOBALNO)
    allocateGlobalno();
  else if (attr==FEM_IS_VALID_ATTR)
    allocateValid();
  else if (attr==FEM_MESH_SIZING)
    allocateMeshSizing();
  else if(attr == FEM_CHUNK){
    FEM_IndexAttribute *chunkNo= new FEM_IndexAttribute(this,FEM_CHUNK,NULL);
    add(chunkNo);
    chunkNo->setWidth(1);
  }
  else if(attr == FEM_BOUNDARY){
      allocateBoundary();
  }else if(attr == FEM_ADAPT_FACE_ADJ){
      FEM_DataAttribute *adaptAdj = new FEM_DataAttribute(this,attr);
      adaptAdj->setDatatype(FEM_BYTE);
      add(adaptAdj);
  }else if(attr == FEM_ADAPT_EDGE_ADJ){
      FEM_DataAttribute *adaptAdj = new FEM_DataAttribute(this,attr);
      adaptAdj->setDatatype(FEM_BYTE);
      add(adaptAdj);
  }else if(attr == FEM_ADAPT_LOCK){
      FEM_DataAttribute *adaptLock = new FEM_DataAttribute(this,attr);
      adaptLock->setDatatype(FEM_INT);
      adaptLock->setWidth(2);
      add(adaptLock);
  }else if(attr == FEM_ADAPT_LOCK_PRIO){
      FEM_DataAttribute *adaptLockPrio = new FEM_DataAttribute(this,attr);
      adaptLockPrio->setDatatype(FEM_DOUBLE);
      adaptLockPrio->setWidth(1);
      add(adaptLockPrio);
  } else if (attr == FEM_PARTITION) {
      FEM_DataAttribute* partition = new FEM_DataAttribute(this, attr);
      partition->setDatatype(FEM_INT);
      partition->setWidth(1);
      add(partition);
  }
  else {
    //It's an unrecognized tag: abort
    char attrNameStorage[256], msg[1024];
    sprintf(msg,"Could not locate the attribute %s for entity %s",
            FEM_Get_attr_name(attr,attrNameStorage), getName());
    FEM_Abort(caller,msg);
  }
}

void FEM_Entity::allocateCoord(void) {
    if (coord) CkAbort("FEM_Entity::allocateCoord called, but already allocated");
    coord=new FEM_DataAttribute(this,FEM_COORD);
    add(coord); // coord will be deleted by FEM_Entity now
    coord->setDatatype(FEM_DOUBLE);
}

void FEM_Entity::allocateSym(void) {
    if (sym) CkAbort("FEM_Entity::allocateSym called, but already allocated");
    sym=new FEM_DataAttribute(this,FEM_SYMMETRIES);
    add(sym); // sym will be deleted via attributes list now
    sym->setWidth(1);
    sym->setDatatype(FEM_BYTE); //Same as FEM_Symmetries_t
}

void FEM_Entity::setSymmetries(int r,FEM_Symmetries_t s)
{
    if (!sym) {
        if (s==0) return; //Don't bother allocating just for 0
        allocateSym();
    }
    sym->getChar()(r,0)=s;
}


/*
 * Set the coordinates for a node or other entity.
 * Use the appropriate 2d or 3d version.
 *
 * Note that first the function attempts to find coordinates in
 * the FEM_COORD attribute's array "*coord". If this fails, then
 * it will use the user's FEM_DATA field. Little error checking is
 * done, so the functions may crash if used inappropriately.
 */
inline void FEM_Entity::set_coord(int idx, double x, double y){
  if(coord){
    coord->getDouble()(idx,0)=x;
    coord->getDouble()(idx,1)=y;
  }
  else {
    FEM_DataAttribute* attr =   (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
    attr->getDouble()(idx,0)=x;
    attr->getDouble()(idx,1)=y;
  }
}

inline void FEM_Entity::set_coord(int idx, double x, double y, double z){
  if(coord){
    coord->getDouble()(idx,0)=x;
    coord->getDouble()(idx,1)=y;
    coord->getDouble()(idx,2)=z;
  }
  else {
    FEM_DataAttribute* attr =   (FEM_DataAttribute*)lookup(FEM_DATA,"set_coord");
    attr->getDouble()(idx,0)=x;
    attr->getDouble()(idx,1)=y;
    attr->getDouble()(idx,2)=z;
  }
}


void FEM_Entity::get_coord(int idx, double &x, double &y){
  if(coord){
      if(idx<0){ // ghost
        x = ghost->coord->getDouble()(FEM_From_ghost_index(idx),0);
        y = ghost->coord->getDouble()(FEM_From_ghost_index(idx),1);
      } else { // non-ghost
        x = coord->getDouble()(idx,0);
        y = coord->getDouble()(idx,1);
      }
      
    }
  else {
  
    if(idx<0){ // ghost
      FEM_DataAttribute* attr = (FEM_DataAttribute*)ghost->lookup(FEM_DATA,"get_coord");
      x = attr->getDouble()(FEM_From_ghost_index(idx),0);
      y = attr->getDouble()(FEM_From_ghost_index(idx),1);
    } else  { // non-ghost
      FEM_DataAttribute* attr = (FEM_DataAttribute*)lookup(FEM_DATA,"get_coord");
      x = attr->getDouble()(idx,0);
      y = attr->getDouble()(idx,1);
    }
    
  }
  
}

void FEM_Entity::get_coord(int idx, double &x, double &y, double &z){
  if(coord){
    if(idx<0){ // ghost
      x = ghost->coord->getDouble()(FEM_From_ghost_index(idx),0);
      y = ghost->coord->getDouble()(FEM_From_ghost_index(idx),1);
      z = ghost->coord->getDouble()(FEM_From_ghost_index(idx),2);
    } else { // non-ghost
      x = coord->getDouble()(idx,0);
      y = coord->getDouble()(idx,1);
      z = coord->getDouble()(idx,2);
    }
    
  }
  else {
     
    if(idx<0){ // ghost
      FEM_DataAttribute* attr = (FEM_DataAttribute*)ghost->lookup(FEM_DATA,"get_coord");
      x = attr->getDouble()(FEM_From_ghost_index(idx),0);
      y = attr->getDouble()(FEM_From_ghost_index(idx),1);
      z = attr->getDouble()(FEM_From_ghost_index(idx),2);
    } else { // non-ghost
      FEM_DataAttribute* attr = (FEM_DataAttribute*)lookup(FEM_DATA,"get_coord");
      x = attr->getDouble()(idx,0);
      y = attr->getDouble()(idx,1);
      z = attr->getDouble()(idx,2);
    }
    
  }
}





void FEM_Entity::allocateGlobalno(void) {
    if (globalno) CkAbort("FEM_Entity::allocateGlobalno called, but already allocated");
    globalno=new FEM_IndexAttribute(this,FEM_GLOBALNO,NULL);
    add(globalno); // globalno will be deleted via attributes list now
    globalno->setWidth(1);
}

void FEM_Entity::allocateMeshSizing(void) {
  if (meshSizing)
    CkAbort("FEM_Entity::allocateMeshSizing called, but already allocated");
  meshSizing=new FEM_DataAttribute(this,FEM_MESH_SIZING);
  add(meshSizing); // meshSizing will be deleted via attributes list now
  meshSizing->setWidth(1);
  meshSizing->setDatatype(FEM_DOUBLE);
}

double FEM_Entity::getMeshSizing(int r) {
  if (!meshSizing) {
    allocateMeshSizing();
    return -1.0;
  }
  if (r >= 0)  return meshSizing->getDouble()(r,0);
  else  return ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0);
}

void FEM_Entity::setMeshSizing(int r,double s)
{
  if (!meshSizing) allocateMeshSizing();
  if (s <= 0.0) return;
  if (r >= 0)  meshSizing->getDouble()(r,0)=s;
  else ghost->meshSizing->getDouble()(FEM_To_ghost_index(r),0)=s;
}

void FEM_Entity::setMeshSizing(double *sf)
{
  if (!meshSizing) allocateMeshSizing();
  int len = size();
  for (int i=0; i<len; i++)
    meshSizing->getDouble()(i,0)=sf[i];
}

void FEM_Entity::allocateBoundary(){
    boundary = new FEM_DataAttribute(this,FEM_BOUNDARY);
    add(boundary);
    boundary->setWidth(1);
    boundary->setDatatype(FEM_INT);
}


int FEM_Entity::isBoundary(int idx){
    return boundary->getInt()(idx,0);
}

void FEM_Entity::setGlobalno(int r,int g) {
    if (!globalno) allocateGlobalno();
    globalno->get()(r,0)=g;
}
void FEM_Entity::setAscendingGlobalno(void) {
    if (!globalno) {
        allocateGlobalno();
        int len=size();
        for (int i=0;i<len;i++) globalno->get()(i,0)=i;
    }
}
void FEM_Entity::setAscendingGlobalno(int base) {
    if (!globalno) {
        allocateGlobalno();
        int len=size();
        for (int i=0;i<len;i++) globalno->get()(i,0)=i+base;
    }
}
void FEM_Entity::copyOldGlobalno(const FEM_Entity &e) {
    if ((!hasGlobalno()) && e.hasGlobalno() && size()>=e.size()) {
        for (int i=0;i<size();i++)
            setGlobalno(i,e.getGlobalno(i));
    }
}

/*
 * method to clear ghosts associated with this entity
 * */

void FEM_Entity::clearGhost(){
    CkAssert(ghost != NULL);
    ghostSend.clear();
    ghost->setLength(0);
    ghost->ghostRecv.clear();
}

/********************** Node *****************/
FEM_Node::FEM_Node(FEM_Node *ghost_)
  :FEM_Entity(ghost_), primary(0), sharedIDXL(&shared,&shared),
   elemAdjacency(0),nodeAdjacency(0)
{}

void FEM_Node::allocatePrimary(void) {
  if (primary) CkAbort("FEM_Node::allocatePrimary called, but already allocated");
  primary=new FEM_DataAttribute(this,FEM_NODE_PRIMARY);
  add(primary); // primary will be deleted by FEM_Entity now
  primary->setWidth(1); //Only 1 flag per node
  primary->setDatatype(FEM_BYTE);
}

bool FEM_Node::hasConn(int idx) {
  if((elemAdjacency->get()[idx].length() > 0)||(nodeAdjacency->get()[idx].length() > 0))
    return true;
  else return false;
}

void FEM_Node::pup(PUP::er &p) {
    p.comment(" ---------------- Nodes ------------------ ");
    super::pup(p);
    p.comment(" ---- Shared nodes ----- ");
    shared.pup(p);
    p.comment(" shared nodes IDXL ");
    sharedIDXL.pup(p);
}
FEM_Node::~FEM_Node() {
}


const char *FEM_Node::getName(void) const {return "FEM_NODE";}

void FEM_Node::create(int attr,const char *caller) {
  if (attr==FEM_NODE_PRIMARY)
    allocatePrimary();
  else if(attr == FEM_NODE_ELEM_ADJACENCY)
    allocateElemAdjacency();
  else if(attr == FEM_NODE_NODE_ADJACENCY)
    allocateNodeAdjacency();
  else
    super::create(attr,caller);
}


/********************** Elem *****************/
class FEM_Elem_Conn_Checker : public FEM_IndexAttribute::Checker {
    const FEM_Entity &sizeSrc;
    const FEM_Entity *sizeSrc2;
public:
    FEM_Elem_Conn_Checker(const FEM_Entity &sizeSrc_,const FEM_Entity *sizeSrc2_)
        :sizeSrc(sizeSrc_), sizeSrc2(sizeSrc2_) {}

    void check(int row,const BasicTable2d<int> &table,const char *caller) const {
        const int *idx=table.getRow(row);
        int n=table.width();
        int max=sizeSrc.size();
        if (sizeSrc2) max+=sizeSrc2->size();
        for (int i=0;i<n;i++)
            if ((idx[i]<0) || (idx[i]>=max))
            { /* This index is out of bounds: */
                if (idx[i]<0)
                    FEM_Abort(caller,"Connectivity entry %d's value, %d, is negative",row,idx[i]);
                else /* (idx[i]>=max) */
                    FEM_Abort(caller,
                        "Connectivity entry %d's value, %d, should be less than the number of nodes, %d",
                        row,idx[i],max);
            }
    }
};

FEM_Elem::FEM_Elem(const FEM_Mesh &mesh, FEM_Elem *ghost_)
  :FEM_Entity(ghost_), elemAdjacency(0), elemAdjacencyTypes(0)
{
    FEM_IndexAttribute::Checker *c;
    if (isGhost()) // Ghost elements can point to both real as well as ghost nodes
        c=new FEM_Elem_Conn_Checker(mesh.node, mesh.node.getGhost());
    else /* is real */ //Real elements only point to real nodes
        c=new FEM_Elem_Conn_Checker(mesh.node, NULL);
    conn=new FEM_IndexAttribute(this,FEM_CONN,c);
    add(conn); // conn will be deleted by FEM_Entity now
}
void FEM_Elem::pup(PUP::er &p) {
    p.comment(" ------------- Element data ---------- ");
    FEM_Entity::pup(p);
}
FEM_Elem::~FEM_Elem() {
}


void FEM_Elem::create(int attr,const char *caller) {
  // We need to catch both FEM_ELEM_ELEM_ADJACENCY and FEM_ELEM_ELEM_ADJ_TYPES,
  // since if either one falls through to
  // the super::create(), it will not know what to do, and will fail
  //
  // Note: allocateElemAdjacency() will create both attribute fields since they
  //       should always be used together.

  if(attr == FEM_CONN) ;
    //who allocates the conn?
  else if(attr == FEM_ELEM_ELEM_ADJACENCY)
    allocateElemAdjacency();
  else if(attr == FEM_ELEM_ELEM_ADJ_TYPES)
    allocateElemAdjacency();
  else
    super::create(attr,caller);
}


const char *FEM_Elem::getName(void) const {
    return "FEM_ELEM";
}

bool FEM_Elem::hasConn(int idx) {
  return false;
}

/********************* Sparse ******************/
class FEM_Sparse_Elem_Checker : public FEM_IndexAttribute::Checker {
    const FEM_Mesh &mesh;
public:
    FEM_Sparse_Elem_Checker(const FEM_Mesh &mesh_) :mesh(mesh_) {}

    void check(int row,const BasicTable2d<int> &table,const char *caller) const {
        //assert: table.getWidth==2
        const int *elem=table.getRow(row);
        int maxT=mesh.elem.size();
        if ((elem[0]<0) || (elem[1]<0))
            FEM_Abort(caller,"Sparse element entry %d's values, %d and %d, are negative",
                row,elem[0],elem[1]);
        int t=elem[0];
        if (t>=maxT)
            FEM_Abort(caller,"Sparse element entry %d's element type, %d, is too big",
                row,elem[0]);
        if (elem[1]>=mesh.elem[t].size())
            FEM_Abort(caller,"Sparse element entry %d's element index, %d, is too big",
                row,elem[1]);
    }
};

FEM_Sparse::FEM_Sparse(const FEM_Mesh &mesh_,FEM_Sparse *ghost_)
    :FEM_Elem(mesh_,ghost_), elem(0), mesh(mesh_)
{
}
void FEM_Sparse::allocateElem(void) {
    if (elem) CkAbort("FEM_Sparse::allocateElem called, but already allocated");
    FEM_IndexAttribute::Checker *checker=new FEM_Sparse_Elem_Checker(mesh);
    elem=new FEM_IndexAttribute(this,FEM_SPARSE_ELEM,checker);
    add(elem); //FEM_Entity will delete elem now
    elem->setWidth(2); //SPARSE_ELEM consists of pairs: element type, element number
}
void FEM_Sparse::pup(PUP::er &p) {
    p.comment(" ------------- Sparse Element ---------- ");
    super::pup(p);
}
FEM_Sparse::~FEM_Sparse() {
}

const char *FEM_Sparse::getName(void) const { return "FEM_SPARSE"; }

void FEM_Sparse::create(int attr,const char *caller) {
    if (attr==FEM_SPARSE_ELEM)
        allocateElem();
    else /*super*/ FEM_Entity::create(attr,caller);
}


/******************* Mesh *********************/
FEM_Mesh::FEM_Mesh()
    :node(new FEM_Node(NULL)),
     elem(*this,"FEM_ELEM"),
     sparse(*this,"FEM_SPARSE"),
     lastElemAdjLayer(NULL),
     fmMM(NULL)
{
    m_isSetting=true; //Meshes start out setting
    lastLayerSet = false;
}

FEM_Mesh::~FEM_Mesh() {
}

FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) {
    FEM_Entity *e=NULL;
    if (entity>=FEM_ENTITY_FIRST && entity<FEM_ENTITY_LAST)
    { //It's in the right range for an entity code:
        bool isGhost=false;
        if (entity-FEM_ENTITY_FIRST>=FEM_GHOST) {
            entity-=FEM_GHOST;
            isGhost=true;
        }
        if (entity==FEM_NODE)
            e=&node;
        else if (entity>=FEM_ELEM && entity<FEM_ELEM+100)
        { //It's a kind of element:
            int elType=entity-FEM_ELEM;
            e=&elem.set(elType);
        }
        else if (entity>=FEM_SPARSE && entity<FEM_SPARSE+100)
        { //It's a kind of sparse:
            int sID=entity-FEM_SPARSE;
            e=&sparse.set(sID);
        }

        if (isGhost) //Move from the real to the ghost entity
            e=e->getGhost();
    }

    if (e==NULL) //We didn't find an entity!
        FEM_Abort(caller,"Expected an entity type (FEM_NODE, FEM_ELEM, etc.) but got %d",entity);
    return e;
}
const FEM_Entity *FEM_Mesh::lookup(int entity,const char *caller) const {
    return ((FEM_Mesh *)this)->lookup(entity,caller);
}

void FEM_Mesh::setFemMeshModify(femMeshModify *m){
  fmMM = m;
}

void FEM_Mesh::setParfumSA(ParFUMShadowArray *m){
  parfumSA = m;
}


femMeshModify *FEM_Mesh::getfmMM(){
  return fmMM;
}

void FEM_Mesh::pup(PUP::er &p)  //For migration
{
    p.comment(" ------------- Node Data ---------- ");
    node.pup(p);

    p.comment(" ------------- Element Types ---------- ");
    elem.pup(p);

    p.comment("-------------- Sparse Types ------------");
    sparse.pup(p);

    p.comment("-------------- Symmetries ------------");
    symList.pup(p);

    p|m_isSetting;
    p|lastLayerSet;
    if(lastLayerSet) {
      if(p.isUnpacking()) {
        if(lastElemAdjLayer==NULL) lastElemAdjLayer = new FEM_ElemAdj_Layer();
      }
      lastElemAdjLayer->pup(p);
    }
    p.comment("-------------- Mesh data --------------");
    udata.pup(p);


/* NOTE: for backward file compatability (fem_mesh_vp files),
   be sure to add new stuff at the *end* of this routine--
   it will be read as zeros for old files. */
}

int FEM_Mesh::chkET(int elType) const {
    if ((elType<0)||(elType>=elem.size())) {
        CkError("FEM Error! Bad element type %d!\n",elType);
        CkAbort("FEM Error! Bad element type used!\n");
    }
    return elType;
}

int FEM_Mesh::nElems(int t_max) const //Return total number of elements before type t_max
{
#if CMK_ERROR_CHECKING
    if (t_max<0 || t_max>elem.size()) {
        CkPrintf("FEM> Invalid element type %d used!\n");
        CkAbort("FEM> Invalid element type");
    }
#endif
    int ret=0;
    for (int t=0;t<t_max;t++){
        if (elem.has(t)){
            ret+=elem.get(t).size();
        }
    }
    return ret;
}

int FEM_Mesh::getGlobalElem(int elType,int elNo) const
{
    int base=nElems(elType); //Global number of first element of this type
#if CMK_ERROR_CHECKING
    if (elNo<0 || elNo>=elem[elType].size()) {
        CkPrintf("FEM> Element number %d is invalid-- element type %d has only %d elements\n",
            elNo,elType,elem[elType].size());
        CkAbort("FEM> Invalid element number, probably passed via FEM_Set_Sparse_elem");
    }
#endif
    return base+elNo;
}

void FEM_Mesh::setAscendingGlobalno(void) {
    node.setAscendingGlobalno();
    for (int e=0;e<elem.size();e++)
        if (elem.has(e)) elem[e].setAscendingGlobalno();
    for (int s=0;s<sparse.size();s++)
        if (sparse.has(s)) sparse[s].setAscendingGlobalno();
}
void FEM_Mesh::setAbsoluteGlobalno(){
    node.setAscendingGlobalno();
    for (int e=0;e<elem.size();e++){
        if (elem.has(e)) elem[e].setAscendingGlobalno(nElems(e));
    }
}

void FEM_Mesh::copyOldGlobalno(const FEM_Mesh &m) {
    node.copyOldGlobalno(m.node);
    for (int e=0;e<m.elem.size();e++)
        if (m.elem.has(e) && e<elem.size() && elem.has(e))
            elem[e].copyOldGlobalno(m.elem[e]);
    for (int s=0;s<m.sparse.size();s++)
        if (m.sparse.has(s) && s<sparse.size() && sparse.has(s))
            sparse[s].copyOldGlobalno(m.sparse[s]);
}

void FEM_Mesh::clearSharedNodes(){
    node.shared.clear();
}

void FEM_Mesh::clearGhostNodes(){
    node.clearGhost();
}

void FEM_Mesh::clearGhostElems(){
    for(int i=0;i<elem.size();i++){
        if(elem.has(i)){
            elem[i].clearGhost();
        }
    }
}

void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType) {
    if (type<0 || type>maxType) {
        char msg[1024];
        sprintf(msg,"%s %d is not a valid entity type (it must be between %d and %d)",
            entityType,type, 0, maxType-1);
        FEM_Abort(caller,msg);
    }
}
void FEM_Is_NULL(const char *caller,const char *entityType,int type) {
    char msg[1024];
    sprintf(msg,"%s %d was never set--it cannot now be read",entityType,type);
    FEM_Abort(caller,msg);
}

void FEM_Mesh::copyShape(const FEM_Mesh &src)
{
    node.copyShape(src.node);
    for (int t=0;t<src.elem.size();t++)
        if (src.elem.has(t)) elem.set(t).copyShape(src.elem.get(t));

    for (int s=0;s<src.sparse.size();s++)
        if (src.sparse.has(s)) sparse.set(s).copyShape(src.sparse.get(s));

    setSymList(src.getSymList());
}

int FEM_Mesh::getEntities(int *entities) {
    int len=0;
    entities[len++]=FEM_NODE;
    for (int t=0;t<elem.size();t++)
        if (elem.has(t)) entities[len++]=FEM_ELEM+t;
    for (int s=0;s<sparse.size();s++)
        if (sparse.has(s)) entities[len++]=FEM_SPARSE+s;
    return len;
}


FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead)
{
    char fname[256];
    static const char *meshFileNames="%s_vp%d_%d.dat";
    sprintf(fname, meshFileNames, prefix, nchunks, chunkNo);
    FILE *fp = fopen(fname, forRead?"r":"w");
    CkPrintf("FEM> %s %s...\n",forRead?"Reading":"Writing",fname);
    if(fp==0) {
      FEM_Abort(forRead?"FEM: unable to open input file"
        :"FEM: unable to create output file.\n");
    }
    return fp;
}

static inline void read_version()
{
    FILE *f = fopen("FEMVERSION", "r");
    if (f!=NULL)  {
    fscanf(f, "%d", &femVersion);
    if (CkMyPe()==0) CkPrintf("FEM> femVersion detected: %d\n", femVersion);
    fclose(f);
    }
}

static inline void write_version()
{
    FILE *f = fopen("FEMVERSION", "w");
    if (f!=NULL)  {
        fprintf(f, "%d", femVersion);
        fclose(f);
    }
}

FEM_Mesh *FEM_readMesh(const char *prefix, int chunkNo, int nChunks)
{
    // find FEM file version number
    static int version_checked = 0;
    if (!version_checked) {
        version_checked = 1;
        read_version();
    }

    FEM_Mesh *ret = new FEM_Mesh;
    ret->becomeGetting();
        FILE *fp = FEM_openMeshFile(prefix, chunkNo, nChunks, true);
    PUP::fromTextFile p(fp);
    ret->pup(p);
    fclose(fp);

#ifdef PRINT_SHARED_NODE_INFO
        // For Abhinav, print out the neighbors for this vp
    // CkPrintf("%d: Finding Neighbors for VP\n", chunkNo);
        
        FEM_Comm &shared = ret->node.shared;    // Shared nodes: The type is really an IDXL_Side
        int numNeighborVPs = shared.size();
        // CkPrintf("%d: communicates with %d neighbors\n", chunkNo, numNeighborVPs);
        
        for(int i=0; i<numNeighborVPs; i++){
          IDXL_List list = shared.getLocalList(i);

      double x, y, sumx = 0.0, sumy = 0.0;
      int nnodes = ret->node.size();
      for(int n=0; n<nnodes; n++){
        ret->node.get_coord(n, x, y);
        sumx += x;
        sumy += y;
      }

          CkPrintf("chunk %d at %f %f communicates with chunk %d through %d shared nodes\n", chunkNo, sumx/(double)nnodes, sumy/(double)nnodes, list.getDest(), list.size()); 
        }
#endif

    return ret;
}

void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks)
{
    if (chunkNo == 0) {
       write_version();
    }
        FILE *fp = FEM_openMeshFile(prefix,chunkNo,nChunks,false);
    PUP::toTextFile p(fp);
    m->pup(p);
    fclose(fp);
}

// A subroutine that can enable projections tracing
CLINKAGE void FEM_traceBegin(){
  traceBegin();
}
FORTRAN_AS_C(FEM_TRACEBEGIN,
             FEM_traceBegin,
             fem_tracebegin,
             (),  () )

// A subroutine that can disable projections tracing                                                                                
CLINKAGE void FEM_traceEnd(){
  traceEnd();
}
FORTRAN_AS_C(FEM_TRACEEND,
             FEM_traceEnd,
             fem_traceend,
             (),  () )




// Setup the entity FEM_IS_VALID tables
CLINKAGE void FEM_Mesh_allocate_valid_attr(int fem_mesh, int entity_type){
  FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entity_type,"FEM_Mesh_allocate_valid_attr");
  entity->allocateValid();
}
FORTRAN_AS_C(FEM_MESH_ALLOCATE_VALID_ATTR,
             FEM_Mesh_allocate_valid_attr,
             fem_mesh_allocate_valid_attr,
             (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )


// mark an entity as valid
CLINKAGE void FEM_set_entity_valid(int mesh, int entityType, int entityIdx){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_");
  FEM_Entity *entity = m->lookup(entityType,"FEM_");
  entity->set_valid(entityIdx,true);
}
FORTRAN_AS_C(FEM_SET_ENTITY_VALID,
             FEM_set_entity_valid,
             fem_set_entity_valid,
             (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )


// mark an entity as invalid
CLINKAGE void FEM_set_entity_invalid(int mesh, int entityType, int entityIdx){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
  return entity->set_invalid(entityIdx,true);
}
FORTRAN_AS_C(FEM_SET_ENTITY_INVALID,
             FEM_set_entity_invalid,
             fem_set_entity_invalid,
             (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )


// Determine if an entity is valid
CLINKAGE int FEM_is_valid(int mesh, int entityType, int entityIdx){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
  return (entity->is_valid(entityIdx) != 0);
}
FORTRAN_AS_C(FEM_IS_VALID,
             FEM_is_valid,
             fem_is_valid,
             (int *fem_mesh, int *entity_type, int *entityIdx),  (*fem_mesh, *entity_type, *entityIdx) )


// Count number of valid items for this entity type
CLINKAGE int FEM_count_valid(int mesh, int entityType){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
  return entity->count_valid();
}
FORTRAN_AS_C(FEM_COUNT_VALID,
             FEM_count_valid,
             fem_count_valid,
             (int *fem_mesh, int *entity_type),  (*fem_mesh, *entity_type) )


CLINKAGE int FEM_is_node_shared(int mesh_id, int idx) {
    FEM_Mesh* mesh = FEM_Mesh_lookup(mesh_id,"FEM_is_node_shared");
    const IDXL_Rec *rec = mesh->node.shared.getRec(idx);
    return (int)(rec != NULL);
}
FORTRAN_AS_C(FEM_IS_NODE_SHARED,
        FEM_is_node_shared,
        fem_is_node_shared,
        (int* fem_mesh, int* node_idx), (*fem_mesh, *node_idx))


CLINKAGE int FEM_find_node_owner(int mesh_id, int idx) {
    FEM_Mesh* mesh = FEM_Mesh_lookup(mesh_id,"FEM_find_node_owner");
    const IDXL_Rec *rec = mesh->node.shared.getRec(idx);
    int min = FEM_My_partition();
    if (!rec) {
        return min;
    }
    for (int i=0; i<rec->getShared(); i++) {
        int chk = rec->getChk(i);
        if (chk < min) {
            min = chk;
        }
    }
    return min;
}
FORTRAN_AS_C(FEM_FIND_NODE_OWNER,
        FEM_find_node_owner,
        fem_find_node_owner,
        (int* fem_mesh, int* node_idx), (*fem_mesh, *node_idx))


// Set coordinates for some entity's item number idx
void FEM_set_entity_coord2(int mesh, int entityType, int idx, double x, double y){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
  entity->set_coord(idx,x,y);
}
void FEM_set_entity_coord3(int mesh, int entityType, int idx, double x, double y, double z){
  FEM_Mesh *m=FEM_Mesh_lookup(mesh,"FEM_Mesh_create_valid_elem");
  FEM_Entity *entity = m->lookup(entityType,"FEM_Mesh_allocate_valid_attr");
  entity->set_coord(idx,x,y,z);
}

---