Home Research Papers Posters Manuals Talks Download People Help Internal

libs/ck-libs/ParFUM/bulk_adapt.h

Go to the documentation of this file.

/* File: bulk_adapt.h
 * Author: Terry Wilmarth
 */
#ifndef __ParFUM_Bulk_Adapt_H
#define __ParFUM_Bulk_Adapt_H

#define SLIVERAREA 1.0e-18
#define REFINE_TOL 1.1  // Refine elements with average edge length > 
                        // REFINE_TOL*desiredEdgeLength
#define COARSEN_TOL 0.8 // Coarsen element with average edge length <
                        // COARSEN_TOL*desiredEdgeLength
#define QUALITY_MIN 0.7

// Refine Methods
#define SCALED_SIZING 0
#define ABSOLUTE_SIZING 1

class Bulk_Adapt;
CtvExtern(Bulk_Adapt *, _bulkAdapt);


class Element_Bucket {
 public:
  typedef struct {
    int elID;     // the element ID
    double len;   // the sorting criteria
    int leftIdx;  // lesser elements
    int rightIdx; // greater elements
  } elemEntry;
  elemEntry *elements;
  int numElements;
  int bucketSz;
  int root;
  Element_Bucket() { 
    root = -1;
    numElements = bucketSz = 0;
    elements = NULL;
  }
  void Alloc(int n) {
    elements = new elemEntry[n];
    for (int i=0; i<n; i++) {
      elements[i].elID = -1;
      elements[i].len = -1.0;
      elements[i].leftIdx = -1;
      elements[i].rightIdx = -1;
    }
    bucketSz = n;
    numElements = 0;
    root = -1;
  }
  void Resize(int n) {
    elemEntry *newBkt = new elemEntry[n];
    for (int i=0; i<n; i++) {
      newBkt[i].elID = -1;
      newBkt[i].len = -1.0;
      newBkt[i].leftIdx = -1;
      newBkt[i].rightIdx = -1;
    }
    for (int i=0; i<bucketSz; i++) {
      newBkt[i].elID = elements[i].elID;
      newBkt[i].len = elements[i].len;
      newBkt[i].leftIdx = elements[i].leftIdx;
      newBkt[i].rightIdx = elements[i].rightIdx;
    }
    bucketSz = n;
    if (elements) delete[] elements;
    elements = newBkt;
  }
  void Clear() {
    if (elements) delete[] elements;
    bucketSz = numElements = 0;
    root = -1;
  }
  bool IsBucketEmpty() { return(numElements == 0); }
  void Insert(int elID, double len) {
    CkAssert(numElements <= bucketSz);
    if (numElements == bucketSz) {
      Resize(bucketSz+100);
    }
    CkAssert(elID < bucketSz);
    if (root==-1) {
      root = 0;
      elements[0].elID = elID;
      elements[0].len = len;
      elements[0].leftIdx = elements[0].rightIdx = -1;
      numElements++;
    }
    else {
      int pos = findEmptySlot();
      elements[pos].elID = elID;
      elements[pos].len = len;
      elements[pos].leftIdx = elements[pos].rightIdx = -1;
      if ((len < elements[root].len) || ((len == elements[root].len) && (elID < elements[root].elID))) { 
    // go left
    if (elements[root].leftIdx == -1) {
      elements[root].leftIdx = pos;
    }
    else {
      Insert_help(elements[root].leftIdx, pos);
    }
    numElements++;
      }
      else if ((len > elements[root].len) || ((len == elements[root].len) && (elID >= elements[root].elID))) { 
    // go right
    if (elements[root].rightIdx == -1) {
      elements[root].rightIdx = pos;
    }
    else {
      Insert_help(elements[root].rightIdx, pos);
    }
    numElements++;
      }
      else { // invalid case
    CkAbort("ERROR: Element_Bucket::Insert: unreachable case!\n");
      }
    }
    sanity_check();
  }

  void Insert_help(int subtree, int pos) {
    if ((elements[pos].len < elements[subtree].len) || 
    ((elements[pos].len == elements[subtree].len) && (elements[pos].elID < elements[subtree].elID))) { 
      // go left
      if (elements[subtree].leftIdx == -1) {
    elements[subtree].leftIdx = pos;
      }
      else {
    Insert_help(elements[subtree].leftIdx, pos);
      }
    }
    else if ((elements[pos].len > elements[subtree].len) || 
         ((elements[pos].len == elements[subtree].len) && (elements[pos].elID >=elements[subtree].elID))) { 
      // go right
      if (elements[subtree].rightIdx == -1) {
    elements[subtree].rightIdx = pos;
      }
      else {
    Insert_help(elements[subtree].rightIdx, pos);
      }
    }
    else { // invalid case
      CkAbort("ERROR: Element_Bucket::Insert_help: unreachable case!\n");
    }
  }

  int Remove(double *len) {
    // ASSERT: bucket is not empty; however, will return -1 in this case
    int elID, leftIdx, rightIdx;
    if (numElements == 0) {
      return -1;
    }
    int idx = root;
    int parent = root;
    while (elements[idx].leftIdx != -1) {
      if (parent != idx) parent = idx;
      idx = elements[idx].leftIdx;
    }
    elID = elements[idx].elID;
    CkAssert((elID < bucketSz) && (elID >=0));
    (*len) = elements[idx].len;
    if (elements[idx].rightIdx == -1) { // no children
      elements[idx].elID = -1;
      elements[idx].len = -1.0;
      elements[parent].leftIdx = -1;
    }
    else { // there is a right child
      elements[idx].elID = -1;
      elements[idx].len = -1.0;
      if (parent == idx) {
    root = elements[idx].rightIdx;
      }
      else {
    elements[parent].leftIdx = elements[idx].rightIdx;
      }
      elements[idx].rightIdx = -1;
    }
    numElements--;
    if (numElements==0) {
      root=-1;
    }
    CkPrintf("In Element_Bucket::Remove: the elId to be returned is: %d\n", elID);
    sanity_check();
    return elID;
  }

  int findEmptySlot() {
    for (int i=0; i<numElements+1; i++) {
      if (elements[i].elID == -1) {
    return i;
      }
    }
    CkAbort("ERROR: Element_Bucket::findEmptySlot: more than numElements solts appear to be full!\n");
  }

  void sanity_check() {
    int idx = root;
    int count = 0;
    CkAssert((root == -1) || ((root < bucketSz) && (root >=0)));
    if (root == -1) {
      CkAssert(numElements == 0);
    }
    else {
      sanity_check_helper(idx, &count);
      CkAssert(count == numElements);
    }
  }
  
  void sanity_check_helper(int idx, int *count) {
    (*count)++;
    CkAssert((elements[idx].elID < bucketSz) && (elements[idx].elID >= 0));
    CkAssert(elements[idx].len > 0.0);
    int leftIdx, rightIdx;
    leftIdx = elements[idx].leftIdx;
    CkAssert((leftIdx == -1) || ((leftIdx < bucketSz) && (leftIdx >=0)));
    rightIdx = elements[idx].rightIdx;
    CkAssert((rightIdx == -1) || ((rightIdx < bucketSz) && (rightIdx >=0)));
    if (leftIdx != -1) {
      sanity_check_helper(leftIdx, count);
    }
    if (rightIdx != -1) {
      sanity_check_helper(rightIdx, count);
    }
  }
};



class Bulk_Adapt {
  friend class FEM_Interpolate;
 protected: 
  FEM_Mesh *theMesh; 
  int theMeshID;
    int dim, elemType;
 public:
  Bulk_Adapt() { theMesh = NULL; }
  Bulk_Adapt(FEM_Mesh *m, int d) { theMesh = m; dim=d; }
  void Bulk_Adapt_SetMesh(FEM_Mesh *m) { theMesh = m; }
  void Bulk_Adapt_SetDim(int dimension) { dim = dimension; }
  void Bulk_Adapt_SetElemType(int et) { elemType = et; }
  ~Bulk_Adapt() {}
  void pup(PUP::er &p) {
    p|dim;
  }

  void Bulk_Adapt_Init(int mid, FEM_Mesh *mp, int d) { 
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    theMesh = mp;
    theMeshID = mid;
    dim = d;
    CkPrintf("[%d] Adapt init...\n", rank);
    FEM_Mesh_allocate_valid_attr(theMeshID,FEM_NODE);
    FEM_Mesh_allocate_valid_attr(theMeshID,FEM_ELEM);
    MPI_Barrier(MPI_COMM_WORLD);
    CkPrintf("[%d] Creating adapt adjacencies...\n", rank);
    CreateAdaptAdjacencies(theMeshID, 0);
    ParFUM_SA_Init(theMeshID);  
    MPI_Barrier(MPI_COMM_WORLD);
    CkPrintf("[%d] End Adapt init...\n", rank);
  }

  void ParFUM_Refine(int qm, int method, double factor, double *sizes);
  void ParFUM_Coarsen(int qm, int method, double factor, double *sizes);
  void ParFUM_AdaptMesh(int qm, int method, double factor, double *sizes);
  void ParFUM_Smooth(int qm, int method);
  void ParFUM_Repair(int qm);
  void ParFUM_Remesh(int qm, int method, double factor, double *sizes);
  
  void ParFUM_SetReferenceMesh();
  void ParFUM_GradateMesh(double smoothness);
  void findEdgeLengths(int elemID, double *avgEdgeLength, double *maxEdgeLength, int *maxEdge, 
               double *minEdgeLength, int *minEdge);

 private:
  // Helper methods
  int Refine_h(int qm, int method, double factor, double *sizes);
  int Coarsen_h(int qm, int method, double factor, double *sizes);

  void SetMeshSize(int method, double factor, double *sizes);

  double length(int n1, int n2);
  double length3D(int n1, int n2);
  double length(double *n1_coord, double *n2_coord);
  double length3D(double *n1_coord, double *n2_coord);

  double getArea(int n1, int n2, int n3);
  double getArea(double *n1_coord, double *n2_coord, double *n3_coord);
  double getSignedArea(int n1, int n2, int n3);
  double getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord);
  int getCoord(int n1, double *crds);
  int getShortestEdge(int n1, int n2, int n3, int* shortestEdge);

  double getAreaQuality(int elem);
  void ensureQuality(int n1, int n2, int n3);
  bool controlQualityF(int n1, int n2, int n3, int n4);
  bool controlQualityR(int n1, int n2, int n3, int n4);
  bool controlQualityR(double *n1_coord, double *n2_coord, double *n3_coord);
  bool controlQualityC(int n1, int n2, int n3, double *n4_coord);
  bool controlQualityC(double *n1_coord, double *n2_coord, double *n3_coord, double *n4_coord);
  bool flipOrBisect(int elId, int n1, int n2, int maxEdgeIdx, double maxlen);
};

#endif

---