Home Research Papers Posters Manuals Talks Download People Help Internal

libs/ck-libs/fem/fem_mesh.h

Go to the documentation of this file.

/*
Charm++ Finite Element Framework:
C++ implementation file: Mesh representation and manipulation

This file lists the classes used to represent and manipulate 
Finite Element Meshes inside the FEM framework.

Orion Sky Lawlor, olawlor@acm.org, 1/3/2003
*/
#ifndef __CHARM_FEM_MESH_H
#define __CHARM_FEM_MESH_H

#include "charm-api.h"

#ifndef __CHARM_IDXL_COMM_H
#  include "idxl_comm.h" /*for IDXL_Side */
#endif
#ifndef __CHARM_IDXL_LAYOUT_H
#  include "idxl_layout.h" /*for IDXL_Layout */
#endif

#if defined(_WIN32) && defined(max)
#undef max
#endif

// Map the IDXL names to the old FEM names (FIXME: change all references, too)
typedef IDXL_Side FEM_Comm;
typedef IDXL_List FEM_Comm_List;
typedef IDXL_Rec FEM_Comm_Rec;
class IDXL_Chunk;

class FEM_Comm_Holder {
    IDXL comm; //Our communication lists.
    bool registered; //We're registered with IDXL
    IDXL_Comm_t idx; //Our IDXL index, or -1 if we're unregistered
    void registerIdx(IDXL_Chunk *c);
public:
    FEM_Comm_Holder(FEM_Comm *sendComm, FEM_Comm *recvComm);
    void pup(PUP::er &p);
    ~FEM_Comm_Holder(void);
    
    inline IDXL_Comm_t getIndex(IDXL_Chunk *c) {
        if (idx==-1) registerIdx(c);
        return idx;
    }
};


/*\@{*/


typedef unsigned char FEM_Symmetries_t;


class FEM_Sym_Desc : public PUP::able {
public:
    virtual ~FEM_Sym_Desc();

    virtual CkVector3d applyLoc(const CkVector3d &loc) const =0;
    
    virtual CkVector3d applyVec(const CkVector3d &vec) const =0;
    
    friend inline void operator|(PUP::er &p,FEM_Sym_Desc &a) {a.pup(p);}
    friend inline void operator|(PUP::er &p,FEM_Sym_Desc* &a) {
        PUP::able *pa=a;  p(&pa);  a=(FEM_Sym_Desc *)pa;
    }
};

class FEM_Sym_Linear : public FEM_Sym_Desc {
    CkVector3d shift; //Offset to add to locations
public:
    FEM_Sym_Linear(const CkVector3d &shift_) :shift(shift_) {}
    FEM_Sym_Linear(CkMigrateMessage *m) {}
    
    CkVector3d applyLoc(const CkVector3d &loc) const {return loc+shift;}
    
    virtual CkVector3d applyVec(const CkVector3d &vec) const {return vec;}
    
    virtual void pup(PUP::er &p);
    PUPable_decl(FEM_Sym_Linear);
};

class FEM_Sym_List {
    //This lists the different kinds of symmetry
    CkPupAblePtrVec<FEM_Sym_Desc> sym; 
    
    FEM_Sym_List(const FEM_Sym_List &src); //NOT DEFINED: copy constructor
public:
    FEM_Sym_List();
    void operator=(const FEM_Sym_List &src); //Assignment operator
    ~FEM_Sym_List();
    
    FEM_Symmetries_t add(FEM_Sym_Desc *desc);
    
    void applyLoc(CkVector3d *loc,FEM_Symmetries_t sym) const;
    
    void applyVec(CkVector3d *vec,FEM_Symmetries_t sym) const;
    
    void pup(PUP::er &p);
};


template <class T>
class BasicTable2d : public CkNoncopyable {
protected:
    int rows; //Number of entries in table
    int cols; //Size of each entry in table
    T *table; //Data in table [rows * cols]
public:
    BasicTable2d(T *src,int cols_,int rows_) 
        :rows(rows_), cols(cols_), table(src) {}
    
    inline int size(void) const {return rows;}
    inline int width(void) const {return cols;}
    
    T *getData(void) {return table;}
    const T *getData(void) const {return table;}
    
    T operator() (int r,int c) const {return table[c+r*cols];}
    T &operator() (int r,int c) {return table[c+r*cols];}
    
    //Get a pointer to a row of the table:
    inline T *getRow(int r) {return &table[r*cols];}
    inline const T *getRow(int r) const {return &table[r*cols];}
    inline void getRow(int r,T *dest,T idxBase=0) const {
        const T *src=getRow(r);
        for (int c=0;c<cols;c++) dest[c]=src[c]+idxBase;
    }
    inline T *operator[](int r) {return getRow(r);}
    inline const T *operator[](int r) const {return getRow(r);}
    inline void setRow(int r,const T *src,T idxBase=0) {
        T *dest=getRow(r);
        for (int c=0;c<cols;c++) dest[c]=src[c]-idxBase;
    }
    inline void setRow(int r,T value) {
        T *dest=getRow(r);
        for (int c=0;c<cols;c++) dest[c]=value;
    }
    
    void set(const T *src,T idxBase=0) {
        for (int r=0;r<rows;r++) 
        for (int c=0;c<cols;c++)
            table[c+r*cols]=src[c+r*cols]-idxBase;
    }
    void setTranspose(const T *srcT,int idxBase=0) {
        for (int r=0;r<rows;r++) 
        for (int c=0;c<cols;c++)
            table[c+r*cols]=srcT[r+c*rows]-idxBase;
    }
    void get(T *dest,T idxBase=0) const {
        for (int r=0;r<rows;r++) 
        for (int c=0;c<cols;c++)
            dest[c+r*cols]=table[c+r*cols]+idxBase;
    }
    void getTranspose(T *destT,int idxBase=0) const {
        for (int r=0;r<rows;r++) 
        for (int c=0;c<cols;c++)
            destT[r+c*rows]=table[c+r*cols]+idxBase;
    }
    void set(T value) {
        for (int r=0;r<rows;r++) setRow(r,value);
    }
};

template <class T>
class AllocTable2d : public BasicTable2d<T> {
    int max; //Maximum number of rows that can be used without reallocation
    T fill; //Value to fill uninitialized regions with
    T *allocTable; // the table that I allocated
public:
    AllocTable2d(int cols_=0,int rows_=0,T fill_=0) 
        :BasicTable2d<T>(NULL,cols_,rows_), max(0), fill(fill_),allocTable(NULL)
    {
        if (this->rows>0) allocate(this->rows);
    }
    ~AllocTable2d() {delete[] allocTable;}
    void allocate(int rows_) { 
        allocate(this->width(),rows_);
    }
    void allocate(int cols_,int rows_,int max_=0) { 
        if (cols_==this->cols && rows_<max) {
            //We have room--just update the size:
            this->rows=rows_;
            return;
        }
        if (max_==0) { //They gave no suggested size-- pick one:
            if (rows_==this->rows+1) //Growing slowly: grab a little extra
                max_=10+rows_+(rows_>>2); //  125% plus 10
            else // for a big change, just go with the minimum needed: 
                max_=rows_;
        }

        int oldRows=this->rows;
        this->cols=cols_;
        this->rows=rows_;
        this->max=max_;
        this->table=new T[max*this->cols];
        //Preserve old table entries (FIXME: assumes old cols is unchanged)
        int copyRows=0;
        if (allocTable!=NULL) { 
            copyRows=oldRows;
            if (copyRows>max) copyRows=max;
            memcpy(this->table,allocTable,sizeof(T)*this->cols*copyRows);
            delete[] allocTable;
        }else{
            for (int r=copyRows;r<max;r++)
                this->setRow(r,fill);
        }
        allocTable = this->table;
    }
    
    void pup(PUP::er &p) {
        p|this->rows; p|this->cols;
        if (this->table==NULL) allocate(this->rows);
        p(this->table,this->rows*this->cols); //T better be a basic type, or this won't compile!
    }

    void pupSingle(PUP::er &p, int pupindx) {
      p|this->table[pupindx];
    }

    friend void operator|(PUP::er &p,AllocTable2d<T> &t) {t.pup(p);}

    T *push_back(void) {
        if (this->rows>=max) 
        { //Not already enough room for the new row:
            int newMax=max+(max/4)+16; //Grow 25% longer
            allocate(this->cols,this->rows,newMax);
        }
        this->rows++;
        return getRow(this->rows-1);
    }

    void register_data(T *user,int len,int max_){
        if(allocTable != NULL){
            delete [] allocTable;
            allocTable = NULL;
        }   
        this->table = user;
        this->rows = len;
        max = max_;
    }
};


class FEM_Entity; // Forward declaration
class FEM_Mesh; 


CLINKAGE const char *FEM_Get_entity_name(int entity,char *storage);

CLINKAGE const char *FEM_Get_attr_name(int attr,char *storage);



class FEM_Attribute {
    FEM_Entity *e; //Owning entity (to get length, etc.)
    FEM_Attribute *ghost; // Ghost attribute, which has the same width and datatype as us (or 0)
    int attr; // My attribute code (e.g., FEM_DATA+7, FEM_CONN, etc.)
    
    int width; //Number of columns in our table of data (initially unknown)
    int datatype; //Datatype of entries (initially unknown)
    bool allocated; //True if subclass allocate has been called.
    
    //Abort with a nice error message saying: 
    // Our <field> was previously set to <cur>; it cannot now be <next>
    void bad(const char *field,bool forRead,int cur,int next,const char *caller) const;
    
protected:
    virtual void allocate(int length,int width,int datatype) =0;
public:
    FEM_Attribute(FEM_Entity *owner_,int myAttr_);
    virtual void pup(PUP::er &p);
    virtual void pupSingle(PUP::er &p, int pupindx);
    virtual ~FEM_Attribute();
    
    inline void setGhost(FEM_Attribute *ghost_) { ghost=ghost_;}
    
    inline bool isGhost(void) const { return ghost!=NULL; }
    
    inline int getAttr(void) const {return attr;}

    inline FEM_Entity *getEntity(void) {return e;}
    
    int getLength(void) const;
    int getRealLength(void) const;

    int getMax();
    
    inline int getWidth(void) const { return width<0?0:width; }
    inline int getRealWidth(void) const { return width; }
    
    inline int getDatatype(void) const { return datatype; }
    
    void setLength(int l,const char *caller="");
    
    void setWidth(int w,const char *caller="");
    
    void setDatatype(int dt,const char *caller="");
    
    virtual void copyShape(const FEM_Attribute &src);
    
    void tryAllocate(void);
    
    inline void reallocate(void) { allocated=false; tryAllocate(); }
    
    inline bool isAllocated(void) const {return allocated;}
    
    virtual void set(const void *src, int firstItem,int length, 
        const IDXL_Layout &layout, const char *caller);
    
    virtual void get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller) const;
    
    virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity) =0;

    virtual void register_data(void *dest, int length,int max,
        const IDXL_Layout &layout, const char *caller);
    
};
PUPmarshall(FEM_Attribute)



class FEM_DataAttribute : public FEM_Attribute {
    typedef FEM_Attribute super;
    AllocTable2d<unsigned char> *char_data; //Non-NULL for getDatatype==FEM_BYTE
    AllocTable2d<int> *int_data; //Non-NULL for getDatatype==FEM_INT
    AllocTable2d<float> *float_data; //Non-NULL for getDatatype==FEM_FLOAT
    AllocTable2d<double> *double_data; //Non-NULL for getDatatype==FEM_DOUBLE
protected:
    virtual void allocate(int length,int width,int datatype);
public:
    FEM_DataAttribute(FEM_Entity *owner,int myAttr);
    virtual void pup(PUP::er &p);
    virtual void pupSingle(PUP::er &p, int pupindx);
    ~FEM_DataAttribute();
    
    AllocTable2d<unsigned char> &getChar(void) {return *char_data;}
    const AllocTable2d<unsigned char> &getChar(void) const {return *char_data;}
    
    AllocTable2d<int> &getInt(void) {return *int_data;}
    const AllocTable2d<int> &getInt(void) const {return *int_data;}
    
    AllocTable2d<double> &getDouble(void) {return *double_data;}

    
    virtual void set(const void *src, int firstItem,int length, 
        const IDXL_Layout &layout, const char *caller);
    
    virtual void get(void *dest, int firstItem,int length, 
        const IDXL_Layout &layout, const char *caller) const;
    
    virtual void register_data(void *dest, int length,int max,
        const IDXL_Layout &layout, const char *caller);

    
    virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);

    /*used during refining to extrapolate the values of a node */
    void interpolate(int A,int B,int D,double frac);
    void interpolate(int *iNodes,int rNode,int k);
};
PUPmarshall(FEM_DataAttribute)


class FEM_IndexAttribute : public FEM_Attribute {
public:
    class Checker {
    public:
        virtual ~Checker();
        
        virtual void check(int row,const BasicTable2d<int> &table,
            const char *caller) const =0;
    };
private:
    typedef FEM_Attribute super;
    AllocTable2d<int> idx;
    Checker *checker; //Checks indices (or NULL). This attribute will destroy this object
protected:
    virtual void allocate(int length,int width,int datatype);
public:
    FEM_IndexAttribute(FEM_Entity *owner,int myAttr, Checker *checker_=NULL);
    virtual void pup(PUP::er &p);
    virtual void pupSingle(PUP::er &p, int pupindx);
    ~FEM_IndexAttribute();
    
    AllocTable2d<int> &get(void) {return idx;}
    const AllocTable2d<int> &get(void) const {return idx;}
    
    virtual void set(const void *src, int firstItem,int length, 
        const IDXL_Layout &layout, const char *caller);
    
    virtual void get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller) const;

    virtual void register_data(void *dest, int length, int max,
        const IDXL_Layout &layout, const char *caller);
    
    virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);
};
PUPmarshall(FEM_IndexAttribute)

/*
    This table maps an entity to a list of integer indices 
    of unknown length.
    The node to element adjacency array is an example, where a node
    is mapped to a list of element indices of unknown length.
*/
class FEM_VarIndexAttribute : public FEM_Attribute{
public:
 class ID{
    //type is negative for ghost elements
    //id refers to the index in the entity list
    public:
        int type;
        int id;
        ID(){
            type=-1;
            id = -1;
        };
        ID(int _type,int _id){
          if(_id < 0) {
            type = -(_type+1);
            id = FEM_To_ghost_index(_id);
          }
          else {
            type = _type;
            id = _id;
          }
        };
        bool operator ==(const ID &rhs)const {
            return (type == rhs.type) && (id == rhs.id);
        }
        virtual void pup(PUP::er &p){
            p | type;
            p | id;
        };

        static ID createNodeID(int type,int node){
          ID temp(type, node);
          return temp;
        }
        int getSignedId() {
          if(type<0){
            return FEM_From_ghost_index(id);
          }
          else return id;
        }
 };
private:
    typedef FEM_Attribute super;
    CkVec<CkVec<ID> > idx;
    int oldlength;
protected:
    virtual void allocate(int _length,int _width,int _datatype){
      if(_length > oldlength){
        setWidth(1,"allocate"); //there is 1 vector per entity 
        oldlength = _length*2;
        idx.reserve(oldlength);
        for(int i=idx.size();i<oldlength;i++){
          CkVec<ID> tempVec;
          idx.insert(i,tempVec);
        }
      }
    };
public:
    FEM_VarIndexAttribute(FEM_Entity *owner,int myAttr);
    ~FEM_VarIndexAttribute(){};
    virtual void pup(PUP::er &p);
    virtual void pupSingle(PUP::er &p, int pupindx);
    CkVec<CkVec<ID> > &get(){return idx;};
    const CkVec<CkVec<ID> > &get() const {return idx;};
    
    virtual void set(const void *src,int firstItem,int length,
        const IDXL_Layout &layout,const char *caller);

    virtual void get(void *dest, int firstItem,int length,
        const IDXL_Layout &layout, const char *caller) const;
    
    virtual void copyEntity(int dstEntity,const FEM_Attribute &src,int srcEntity);

    int findInRow(int row,const ID &data);
    
    void print();
};


class l2g_t;
class FEM_Entity {
    typedef CkVec<FEM_Symmetries_t> sym_t;
    int length; // Number of entities of our type currently in the mesh
    int max;    // Maximum number of entities of our type in the mesh that will be allowed
    FEM_Mesh_alloc_fn resize; // should be a function pointer to the actual resize function later
    void *args; // arguments to the resize function
    
    CkVec<FEM_Attribute *> attributes;
    
    FEM_DataAttribute *coord; 
    void allocateCoord(void);
    
    FEM_DataAttribute *sym; 
    void allocateSym(void);
    
    FEM_IndexAttribute *globalno;
    void allocateGlobalno(void);
    
    /*
        used to allocate the integer array for storing the boundary
        values associated with an entity. 
    */
    void allocateBoundary();


    FEM_DataAttribute *meshSizing;
    void allocateMeshSizing(void);

    /*
        used to allocate the char array for storing whether each entity is valid
        When a node/element is deleted the flag in the valid table is set to 0.
        Additionally, we keep track of the first and last occurence in the array of
        invalid indices. This should make searching for slots to reuse quicker.
    */
    FEM_DataAttribute *valid;
    unsigned int first_invalid, last_invalid;
    
protected:
    virtual void create(int attr,const char *caller);
    
    void add(FEM_Attribute *attribute);
 public:

    FEM_Entity *ghost; // Our ghost entity type, or NULL if we're the ghost
        
    FEM_Comm ghostSend; //Non-ghosts we send out (only set for real entities)
    FEM_Comm ghostRecv; //Ghosts we recv into (only set for ghost entities)

    FEM_Entity(FEM_Entity *ghost_); //Default constructor
    void pup(PUP::er &p);
    virtual ~FEM_Entity();
    
    bool isGhost(void) const {return ghost==NULL;}
    
    FEM_Entity *getGhost(void) {return ghost;}
    const FEM_Entity *getGhost(void) const {return ghost;}
    
    inline int size(void) const {return length==-1?0:length;}
    inline int realsize(void) const {return length;}

    // return the maximum size 
    inline int getMax() { if(max > 0) return max; else return length;}
    
    virtual const char *getName(void) const =0;
    
    void copyShape(const FEM_Entity &src);
    
    void setLength(int newlen);

    void setMaxLength(int newLen,int newMaxLen,void *args,FEM_Mesh_alloc_fn fn);
    
    void copyEntity(int dstEntity,const FEM_Entity &src,int srcEntity);

    int push_back(const FEM_Entity &src,int srcEntity);
    
    FEM_Attribute *lookup(int attr,const char *caller);

    
    int getAttrs(int *attrs) const {
        int len=attributes.size();
        for (int i=0;i<len;i++) 
            attrs[i]=attributes[i]->getAttr();
        return len;
    }
    
    void allocateValid();
    void set_valid(unsigned int idx, bool isNode);
    void set_invalid(unsigned int idx, bool isNode);
    int is_valid(unsigned int idx);
    unsigned int count_valid();
    unsigned int get_next_invalid(FEM_Mesh *m, bool isNode, bool isGhost);
    virtual bool hasConn(unsigned int idx)=0;
    void set_coord(int idx, double x, double y);
    void set_coord(int idx, double x, double y, double z);

    CkVec<FEM_Attribute *>* getAttrVec(){
        return &attributes;
    }
    
    // Stupidest possible coordinate access
    inline FEM_DataAttribute *getCoord(void) {return coord;}
    inline const FEM_DataAttribute *getCoord(void) const {return coord;}
    
    //Symmetry array access:
    const FEM_Symmetries_t *getSymmetries(void) const {
        if (sym==NULL) return NULL;
        else return (const FEM_Symmetries_t *)sym->getChar()[0];
    }
    FEM_Symmetries_t getSymmetries(int r) const { 
        if (sym==NULL) return FEM_Symmetries_t(0);
        else return sym->getChar()(r,0);
    }
    void setSymmetries(int r,FEM_Symmetries_t s);
    
    //Global numbering array access
    bool hasGlobalno(void) const {return globalno!=0;}
    int getGlobalno(int r) const {
        if (globalno==0) return -1; // Unknown global number
        return globalno->get()(r,0);
    }
    void setGlobalno(int r,int g);
    void setAscendingGlobalno(void);
    void setAscendingGlobalno(int base);
    void copyOldGlobalno(const FEM_Entity &e);

    // Mesh sizing array access
    bool hasMeshSizing(void) const {return meshSizing!=0;}
    double getMeshSizing(int r); 
    void setMeshSizing(int r,double s);
    void setMeshSizing(double *sf);
    
    //Ghost comm. list access
    FEM_Comm &setGhostSend(void) { return ghostSend; }
    const FEM_Comm &getGhostSend(void) const { return ghostSend; }
    FEM_Comm &setGhostRecv(void) { 
        if (ghost==NULL) return ghostRecv;
        else return ghost->ghostRecv; 
    }
    const FEM_Comm &getGhostRecv(void) const { return ghost->ghostRecv; }
    FEM_Comm_Holder ghostIDXL; //IDXL interface
    
    void addVarIndexAttribute(int code){
        FEM_VarIndexAttribute *varAttribute = new FEM_VarIndexAttribute(this,code);
        add(varAttribute);
    }
    
    void print(const char *type,const IDXL_Print_Map &map);
};
PUPmarshall(FEM_Entity)

// Now that we have FEM_Entity, we can define attribute lenth, as entity length
inline int FEM_Attribute::getLength(void) const { return e->size(); }
inline int FEM_Attribute::getRealLength(void) const { return e->realsize(); }
inline int FEM_Attribute::getMax(){ return e->getMax();}

class FEM_Elem;
class FEM_Node : public FEM_Entity {
    typedef FEM_Entity super;

    FEM_DataAttribute *primary; 
    void allocatePrimary(void);
    
    void allocateElemAdjacency();
    void allocateNodeAdjacency();

    FEM_VarIndexAttribute *elemAdjacency; // stores the node to element adjacency vector 
    FEM_VarIndexAttribute *nodeAdjacency; // stores the node to node adjacency vector 
    typedef FEM_VarIndexAttribute::ID var_id;
protected:
    virtual void create(int attr,const char *caller);
public:
    FEM_Comm shared; //Shared nodes
    FEM_Comm_Holder sharedIDXL; //IDXL interface to shared nodes
    
    FEM_Node(FEM_Node *ghost_);
    void pup(PUP::er &p);
    ~FEM_Node();
    
    virtual const char *getName(void) const;
    
    inline bool getPrimary(int nodeNo) const {
        if (primary==NULL) return true; //Everything must be primary
        else return primary->getChar()(nodeNo,0);
    }
    inline void setPrimary(int nodeNo,bool isPrimary) {
        if (primary==NULL) allocatePrimary();
        primary->getChar()(nodeNo,0)=isPrimary;
    }
    void fillElemAdjacencyTable(int type,const FEM_Elem &elem);
    void setElemAdjacency(int type,const FEM_Elem &elem);
    void fillNodeAdjacency(const FEM_Elem &elem);
    void setNodeAdjacency(const FEM_Elem &elem);
    void fillNodeAdjacencyForElement(int node,int nodesPerElem,const int *conn,FEM_VarIndexAttribute *adjacencyAttr);
    bool hasConn(unsigned int idx);
    void print(const char *type,const IDXL_Print_Map &map);
};
PUPmarshall(FEM_Node)


class FEM_Elem:public FEM_Entity {
    typedef FEM_Entity super;
protected:
    typedef AllocTable2d<int> conn_t;

    int tuplesPerElem;

    // The following are attributes that will commonly be used:
    FEM_IndexAttribute *conn;                // FEM_CONN attribute: element-to-node mapping 
    FEM_IndexAttribute *elemAdjacency;       // FEM_ELEM_ELEM_ADJACENCY attribute
    FEM_IndexAttribute *elemAdjacencyTypes;  // FEM_ELEM_ELEM_ADJ_TYPES attribute

public:

    FEM_Elem(const FEM_Mesh &mesh_, FEM_Elem *ghost_);
    void pup(PUP::er &p);
    ~FEM_Elem();
    
    virtual const char *getName(void) const;
    
    inline conn_t &setConn(void) {return conn->get();}
    inline const conn_t &getConn(void) const {return conn->get();}
    
    void print(const char *type,const IDXL_Print_Map &map);

    void create(int attr,const char *caller);

    void allocateElemAdjacency();

    FEM_IndexAttribute *getElemAdjacency(){return elemAdjacency;}

    // Backward compatability routines:
    int getConn(int elem,int nodeNo) const {return conn->get()(elem,nodeNo);}
    int getNodesPer(void) const {return conn->get().width();}
    int *connFor(int i) {return conn->get().getRow(i);}
    const int *connFor(int i) const {return conn->get().getRow(i);}
    void connIs(int i,const int *src) {conn->get().setRow(i,src);}
    bool hasConn(unsigned int idx);
};
PUPmarshall(FEM_Elem)




class FEM_Sparse : public FEM_Elem {
    typedef FEM_Elem super;
    typedef AllocTable2d<int> elem_t;
    
    FEM_IndexAttribute *elem; //FEM_SPARSE_ELEM attribute
    void allocateElem(void);
    const FEM_Mesh &mesh; //Reference to enclosing mesh, for error checking
protected:
    virtual void create(int attr,const char *caller);
public:
    FEM_Sparse(const FEM_Mesh &mesh_, FEM_Sparse *ghost_);
    void pup(PUP::er &p);
    virtual ~FEM_Sparse();
    
    virtual const char *getName(void) const;
    
    bool hasElements(void) const {return elem!=NULL;}
    
    inline elem_t &setElem(void) {return elem->get();}
    inline const elem_t &getElem(void) const {return elem->get();}
};
PUPmarshall(FEM_Sparse)


class FEM_Userdata_pupfn {
    FEM_Userdata_fn fn;
    void *data;
public:
    FEM_Userdata_pupfn(FEM_Userdata_fn fn_,void *data_)
        :fn(fn_), data(data_) {}
    void pup(PUP::er &p) {
        (fn)((pup_er)&p,data);
    }
};

class FEM_Userdata_item {
    CkVec<char> data; 
public:
    int tag; //User-assigned identifier
    FEM_Userdata_item(int tag_=-1) {tag=tag_;}
    
    bool hasStored(void) const {return data.size()!=0;}
    
    void store(FEM_Userdata_pupfn &f) {
        data.resize(PUP::size(f));
        PUP::toMemBuf(f,&data[0],data.size());
    }
    void restore(FEM_Userdata_pupfn &f) {
        PUP::fromMemBuf(f,&data[0],data.size());
    }
    
    void pup(PUP::er &p) {
        p|tag;
        p|data;
    }
};

class FEM_Userdata_list {
    CkVec<FEM_Userdata_item> list;
public:
    FEM_Userdata_item &find(int tag) {
        for (int i=0;i<list.size();i++)
            if (list[i].tag==tag)
                return list[i];
        // It's not in the list-- add it.
        list.push_back(FEM_Userdata_item(tag));
        return list[list.size()-1];
    }
    int size(){
        return list.size();
    }
    void pup(PUP::er &p) {p|list;}
};


void FEM_Index_Check(const char *caller,const char *entityType,int type,int maxType);
void FEM_Is_NULL(const char *caller,const char *entityType,int type);

template <class T>
class FEM_Entity_Types {
    CkVec<T *> types; // Our main storage for different entity types
    const FEM_Mesh &mesh;
    const char *name; //FEM_SPARSE or FEM_ELEM, or some such.

public:
    FEM_Entity_Types(const FEM_Mesh &mesh_,const char *name_) 
        :mesh(mesh_), name(name_) {}
    void pup(PUP::er &p) { 
        // Can't just use p|types, because T has a funky constructor:
        int n=types.size();
        p|n;
        for (int i=0;i<n;i++) {
            int isNULL=0;
            if (!p.isUnpacking()) isNULL=(types[i]==NULL);
            p|isNULL;
            if (!isNULL) set(i,"pup").pup(p);
        }
    }
    ~FEM_Entity_Types() {
        for (int i=0;i<types.size();i++)
            if (types[i]) delete types[i];
    }
    
    inline int size(void) const {return types.size();}
    
    const T &get(int type,const char *caller="") const {
        FEM_Index_Check(caller,name,type,types.size());
        const T *ret=types[type];
        if (ret==NULL) FEM_Is_NULL(caller,name,type);
        return *ret;
    }
    
    bool has(int type) const { 
        if (type>=types.size()) return false;
        return types[type]!=NULL; 
    }
    
    T &set(int type,const char *caller="") {
        if (type<0) FEM_Index_Check(caller,name,type,types.size());
        while (types.size()<=type) types.push_back(NULL); //Make room for new type
        if (types[type]==NULL) { //Have to allocate a new T:
            T *ghost=new T(mesh,NULL);
            types[type]=new T(mesh,ghost);
        }
        return *types[type];
    }
    
    inline T &operator[] (int type) { return set(type); }
    inline const T &operator[] (int type) const { return get(type); }

};

// Map fortran element (>=1) or node (0) marker to C version (>=1, -1)
inline int zeroToMinusOne(int i) {
    if (i==0) return -1;
    return i;
}




class FEM_ElemAdj_Layer;
class femMeshModify;
class FEM_Mesh : public CkNoncopyable {
  FEM_Sym_List symList;
  bool m_isSetting;
  femMeshModify *fmMM;
  
  void checkElemType(int elType,const char *caller) const;
  void checkSparseType(int uniqueID,const char *caller) const; 
  
  FEM_ElemAdj_Layer* lastElemAdjLayer;

 public:
  void setFemMeshModify(femMeshModify *m);
  
  FEM_Mesh();
  void pup(PUP::er &p); //For migration
  ~FEM_Mesh();
  
  FEM_Node node; 
  
  FEM_Entity_Types<FEM_Elem> elem;
  
  FEM_Entity_Types<FEM_Sparse> sparse;
  
  FEM_Userdata_list udata;
  
  void setSymList(const FEM_Sym_List &src) {symList=src;}
  const FEM_Sym_List &getSymList(void) const {return symList;}
  
  void copyShape(const FEM_Mesh &src);

  // Get the fem mesh modification object associated with this mesh or partition  
  femMeshModify *getfmMM();

  //Return this type of element, given an element type
  FEM_Entity &setCount(int elTypeOrMinusOne) {
    if (elTypeOrMinusOne==-1) return node;
    else return elem[chkET(elTypeOrMinusOne)];
  }
  const FEM_Entity &getCount(int elTypeOrMinusOne) const {
    if (elTypeOrMinusOne==-1) return node;
    else return elem[chkET(elTypeOrMinusOne)];
  }
  FEM_Elem &setElem(int elType) {return elem[chkET(elType)];}
  const FEM_Elem &getElem(int elType) const {return elem[chkET(elType)];}
  int chkET(int elType) const; //Check this element type-- abort if it's bad
  
  FEM_Entity *lookup(int entity,const char *caller);
  const FEM_Entity *lookup(int entity,const char *caller) const;
  
  inline bool isSetting(void) const {return m_isSetting;}
  void becomeSetting(void) {m_isSetting=true;}
  void becomeGetting(void) {m_isSetting=false;}
  
  int nElems() const //Return total number of elements (of all types)
    {return nElems(elem.size());}
  int nElems(int t) const;
  int getGlobalElem(int elType,int elNo) const;
  void setAscendingGlobalno(void);
  void setAbsoluteGlobalno();
  
  void copyOldGlobalno(const FEM_Mesh &m);
  void print(int idxBase);//Write a human-readable description to CkPrintf
  int getEntities(int *entites);
  
  
  
  /********** New methods ***********/
  /*
    This method creates the mapping from a node to all the elements that are 
    incident on it . It assumes the presence of one layer of ghost nodes that
    share a node.
  */
  void createNodeElemAdj();
  void createNodeNodeAdj();
  void createElemElemAdj();
  
  FEM_ElemAdj_Layer *getElemAdjLayer(void);
  
  // Terry's adjacency accessors & modifiers

  //  ------- Element-to-element: preserve initial ordering relative to nodes
  void e2e_getAll(int e, int *neighborss, int etype=0);
  int e2e_getNbr(int e, short idx, int etype=0);
  int e2e_getIndex(int e, int nbr, int etype=0);
  void e2e_setAll(int e, int *neighbors, int etype=0);
  void e2e_setIndex(int e, short idx, int newElem, int etype=0);
  void e2e_replace(int e, int oldNbr, int newNbr, int etype=0);
  void e2e_removeAll(int e, int etype=0);


  //  ------- Element-to-node: preserve initial ordering
  void e2n_getAll(int e, int *adjnodes, int etype=0);
  int e2n_getNode(int e, short idx, int etype=0);
  short e2n_getIndex(int e, int n, int etype=0);
  void e2n_setAll(int e, int *adjnodes, int etype=0);
  void e2n_setIndex(int e, short idx, int newNode, int etype=0);
  void e2n_replace(int e, int oldNode, int newNode, int etype=0);
  void e2n_removeAll(int e, int etype=0);


  //  ------- Node-to-node
  void n2n_getAll(int n, int **adjnodes, int *sz);
  void n2n_add(int n, int newNode);
  void n2n_remove(int n, int oldNode);
  void n2n_replace(int n, int oldNode, int newNode);
  void n2n_removeAll(int n);
  int n2n_exists(int n, int queryNode);

  //  ------- Node-to-element
  void n2e_getAll(int n, int **adjelements, int *sz);
  void n2e_add(int n, int newElem);
  void n2e_remove(int n, int oldElem);
  void n2e_replace(int n, int oldElem, int newElem);
  void n2e_removeAll(int n);

  int getElementOnEdge(int n1, int n2);

  void get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) ;





}; 
PUPmarshall(FEM_Mesh)
FEM_Mesh *FEM_Mesh_lookup(int fem_mesh,const char *caller);
FEM_Entity *FEM_Entity_lookup(int fem_mesh,int entity,const char *caller);
FEM_Attribute *FEM_Attribute_lookup(int fem_mesh,int entity,int attr,const char *caller);

void FEM_Mesh_data_layout(int fem_mesh,int entity,int attr,     
    void *data, int firstItem,int length, const IDXL_Layout &layout);

//registration internal api
void FEM_Register_array_layout(int fem_mesh,int entity,int attr,void *data,int firstItem,const IDXL_Layout &layout);
void FEM_Register_entity_impl(int fem_mesh,int entity,void *args,int len,int max,FEM_Mesh_alloc_fn fn);
FEM_Mesh *FEM_Mesh_assemble(int nchunks,FEM_Mesh **chunks);

FILE *FEM_openMeshFile(const char *prefix,int chunkNo,int nchunks,bool forRead);
FEM_Mesh *FEM_readMesh(const char *prefix,int chunkNo,int nChunks);
void FEM_writeMesh(FEM_Mesh *m,const char *prefix,int chunkNo,int nChunks);

/*\@}*/


#endif

---