Added support for feature detection for 2-d triangular meshes, along with some auxili...
authorIsaac Dooley <idooley2@illinois.edu>
Fri, 9 Nov 2007 03:49:41 +0000 (03:49 +0000)
committerIsaac Dooley <idooley2@illinois.edu>
Fri, 9 Nov 2007 03:49:41 +0000 (03:49 +0000)
src/libs/ck-libs/ParFUM/Makefile
src/libs/ck-libs/ParFUM/ParFUM.h
src/libs/ck-libs/ParFUM/ParFUM_internals.h
src/libs/ck-libs/ParFUM/mesh_adjacency.C
src/libs/ck-libs/ParFUM/mesh_feature_detect.C [new file with mode: 0644]

index 51aa999ab1dc80e7fe4a3c5ef830d36a03296c19..d8066f431dfabd2c2afa4b03077ac978c6bb04ba 100644 (file)
@@ -12,7 +12,7 @@ GENHEADERS= ParFUM.def.h ParFUM.decl.h ParFUM_Adapt.decl.h ParFUM_Adapt.def.h Pa
 
 HEADDEP= $(GENHEADERS) $(HEADERS) $(INTERNALHEADERS)
 
-OBJS=ParFUM.o mesh.o symmetries.o partition.o map.o compat.o call_init.o parallel_part.o mesh_modify.o mesh_adjacency.o adapt.o adapt_algs.o adapt_if.o interpolate.o lock.o util.o lock_node.o adapt_lock.o ParFUM_SA.o global_numbering.o import.o adapt_adj.o bulk_adapt_ops.o
+OBJS=ParFUM.o mesh.o symmetries.o partition.o map.o compat.o call_init.o parallel_part.o mesh_modify.o mesh_adjacency.o adapt.o adapt_algs.o adapt_if.o interpolate.o lock.o util.o lock_node.o adapt_lock.o ParFUM_SA.o global_numbering.o import.o adapt_adj.o bulk_adapt_ops.o mesh_feature_detect.o
 
 COMPAT=compat_init.o compat_finit.o compat_driver.o compat_fdriver.o 
 LIB=libmoduleParFUM
@@ -101,6 +101,9 @@ parallel_part.o: parallel_part.C $(HEADDEP)
 mesh_adjacency.o: mesh_adjacency.C $(HEADDEP)
        $(CHARMC) -c mesh_adjacency.C $(INCS)
 
+mesh_feature_detect.o: mesh_feature_detect.C $(HEADDEP)
+       $(CHARMC) -c mesh_feature_detect.C $(INCS)
+
 adapt.o: adapt.C $(HEADDEP)
        $(CHARMC) -c adapt.C $(INCS)
 
index d9d3d7d01f0a2fdb2660f4490e19c36b270728af..400b9cf9cc68b63483b51c4cf45453e4d6dccb96 100644 (file)
@@ -275,6 +275,10 @@ extern "C" {
   int FEM_is_valid(int mesh, int entityType, int entityIdx);
   int FEM_count_valid(int mesh, int entityType);
 
+  /* Create and modify the Detected Topological and Geometric Mesh Features */
+  void FEM_Mesh_detect_features(int fem_mesh);
+
+
   /* Easy set method for coordinates, that may be helpful when creating a mesh */
   void FEM_set_entity_coord2(int mesh, int entityType, int entityIdx, double x, double y);
   void FEM_set_entity_coord3(int mesh, int entityType, int entityIdx, double x, double y, double z);
index 724167a9a97ff9cd678aa11a8ede686e5a6f7258..1d2da6174345b73dc0018cc471e7be2aba59187c 100644 (file)
@@ -22,6 +22,7 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <vector>
+#include <set>
 
 #include "charm-api.h" /* for Fortran name mangling FTN_NAME */
 #include "ckvector3d.h"
@@ -889,6 +890,9 @@ class FEM_Entity {
   void set_coord(int idx, double x, double y);
   void set_coord(int idx, double x, double y, double z);
 
+  void get_coord(int idx, double &x, double &y);
+  void get_coord(int idx, double &x, double &y, double &z);
+
   /** Expose the attribute vector for refining. breaks modularity but more efficient
   */
   CkVec<FEM_Attribute *>* getAttrVec(){
@@ -1346,7 +1350,7 @@ class FEM_Mesh : public CkNoncopyable {
 
   FEM_ElemAdj_Layer *getElemAdjLayer(void);
 
-  // Terry's adjacency accessors & modifiers
+  // Adjacency accessors & modifiers
 
   //  ------- Element-to-element: preserve initial ordering relative to nodes
   /// Place all of element e's adjacent elements in neighbors; assumes
@@ -1426,6 +1430,32 @@ class FEM_Mesh : public CkNoncopyable {
 
   /// Get two elements adjacent to both n1 and n2
   void get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2) ;
+
+  /// Count the number of elements adjacent to both n1 and n2
+  int countElementsOnEdge(int n1, int n2);
+  
+  
+  
+  // The experimental new code for feature detection below has not yet been widely tested, but it does work for some example
+  
+  /// list of edges on the boundary(adjacent to at least one local node)
+  std::set<std::pair<int,int> > edgesOnBoundary;
+  
+  /// list of the nodes found in edgesOnBoundary
+  std::set<int> verticesOnBoundary;
+  
+  /** list of corners on the mesh. These are nodes in edgesOnBoundary that have small angles between their adjacent boundary edges
+  @note the criterion angle is specified in the preprocessor define CORNER_ANGLE_CUTOFF in mesh_feature_detect.C
+  */
+  std::set<int> cornersOnBoundary;
+    
+  /** Detect features of the mesh, storing results in verticesOnBoundary, edgesOnBoundary, cornersOnBoundary
+      @note currently only 2-d triangular meshes(with triangles in FEM_ELEM+0) are supported
+      @note We require at least a single node-adjacent ghost layer for this to work correctly
+      @note We do not yet pup the structures that store the results, so don't use these with migration yet
+  */
+  void detectFeatures();
+
 };
 PUPmarshall(FEM_Mesh);
 
index 01c8178878c29f306d7aa19010841515b2d10e10..154969961a1b4166b987358fce85ddb0d2ee82ef 100644 (file)
@@ -1233,3 +1233,38 @@ void FEM_Mesh::get2ElementsOnEdge(int n1, int n2, int *result_e1, int *result_e2
   delete[] n2AdjElems;
 }
 
+
+/** Count the number of elements on edge (n1, n2) */
+int FEM_Mesh::countElementsOnEdge(int n1, int n2) {
+  if (n1==n2) {
+    CkPrintf("ERROR: You called countElementsOnEdge() with two identical nodes %d, and %d \n", n1, n2);
+    CkExit();
+  }
+  
+  int *n1AdjElems=0, *n2AdjElems=0;
+  int n1NumElems, n2NumElems;
+
+  CkAssert(node.is_valid_any_idx(n1));
+  CkAssert(node.is_valid_any_idx(n2));
+
+  n2e_getAll(n1, n1AdjElems, n1NumElems);
+  n2e_getAll(n2, n2AdjElems, n2NumElems);
+  CkAssert(n1AdjElems!=0);
+  CkAssert(n2AdjElems!=0);
+  int count=0;
+
+
+  for (int i=0; i<n1NumElems; i++) {
+    for (int j=0; j<n2NumElems; j++) {
+      if (n1AdjElems[i] == n2AdjElems[j]) {
+        count++;
+      }
+    }
+  }
+  delete[] n1AdjElems;
+  delete[] n2AdjElems;
+
+  return count;
+}
+
+
diff --git a/src/libs/ck-libs/ParFUM/mesh_feature_detect.C b/src/libs/ck-libs/ParFUM/mesh_feature_detect.C
new file mode 100644 (file)
index 0000000..6bea277
--- /dev/null
@@ -0,0 +1,299 @@
+/***************************************************
+ * fem_feature_detect.C
+ *
+ * Feature Detection
+ *
+ * Authors include: Isaac
+ *
+ * The experimental new code for feature detection below has not yet been widely tested, but it does work for some example
+ */
+
+#include "ParFUM.h"
+#include "ParFUM_internals.h"
+#include "import.h"
+
+#include <map>
+#include <set>
+#include <utility>
+#include <cmath>
+#define CORNER_ANGLE_CUTOFF  (3.0/4.0*3.14159)
+
+
+using namespace std;
+
+CDECL void
+FEM_Mesh_detect_features(int fem_mesh) {
+  const char *caller="FEM_Mesh_detect_features"; FEMAPI(caller);
+  FEM_Mesh *m=FEM_Mesh_lookup(fem_mesh,caller);
+  m->detectFeatures();
+}
+FORTRAN_AS_C(FEM_MESH_DETECT_FEATURES,
+             FEM_Mesh_detect_features,
+             fem_mesh_detect_features,
+             (int *fem_mesh), (*fem_mesh))
+
+    
+/**
+  A feature detection routine that will eventually work on triangular and tetrahedral meshes. Currently it just works on oriented triangle meshes in serial.
+
+  @note The oriented restriction is here because I am considering the following two edges to be distinct (n1,n2) and (n2,n1). To fix this I could use a class which representst the edges, so that when I put the edges in their container things would work out, or I could order the nodes in the edge by node id, or I could search for (n2,n1) before inserting (n1,n2).
+
+*/
+void FEM_Mesh::detectFeatures() {
+  CkPrintf("detectFeatures() has been called\n");
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+     
+  
+  /** The vertices that look like they are on a boundary of {local mesh unioned with its ghosts} 
+      This is used as an intermediate result.
+  */
+  std::set<std::pair<int,int> > edgesOnBoundaryWithGhosts;
+  
+  
+  // Assume for now we have an oriented triangular mesh in FEM_ELEM+0
+  int nodesPerElem = elem[0].getNodesPer();
+  CkAssert(nodesPerElem==3);
+  
+  // Make sure the valid attributes have been set for the elements
+  // This is actually not sufficient. We should also verify that the adjacencies have been computed
+  elem[0].allocateValid();
+  node.allocateValid();
+  
+ CkPrintf("%d:  %d Nodes, %d Ghost Nodes\n", rank, node.size(), node.ghost->size());
+ CkPrintf("%d:  %d Elem, %d Ghost Elem\n", rank, elem[0].size(), elem[0].ghost->size());
+   
+  
+  for(int n=0; n<node.size(); n++){
+  // Get the coordinates of the 3 nodes
+    double x,y;
+      
+    node.get_coord(n,x,y);
+    CkPrintf("Node %d is at %.1lf,%.1lf\n", n, x, y);    
+  }
+
+  // Find topological edges of this mesh
+
+  for (int ele=0;ele<elem[0].size();ele++) {
+    // if this is a valid element
+    if (elem[0].is_valid(ele)) {
+      
+      const int *conn = elem[0].connFor(ele);
+        
+      int n1 = conn[0];
+      int n2 = conn[1];
+      int n3 = conn[2];
+
+      int numElementsOnEdge;
+      
+      numElementsOnEdge = countElementsOnEdge(n1,n2);
+      if (numElementsOnEdge==1) {
+        edgesOnBoundaryWithGhosts.insert(std::make_pair(n1,n2));
+      }
+      
+      numElementsOnEdge = countElementsOnEdge(n2,n3);
+      if (numElementsOnEdge==1) {
+        edgesOnBoundaryWithGhosts.insert(std::make_pair(n2,n3));
+      }
+      
+      numElementsOnEdge = countElementsOnEdge(n3,n1);
+      if (numElementsOnEdge==1) {
+        edgesOnBoundaryWithGhosts.insert(std::make_pair(n3,n1));
+      }
+         
+    }
+    }
+    
+    // And find anything that looks like an edge around the ghost layer
+    for (int ele=0;ele<elem[0].ghost->size();ele++) {
+    // if this is a valid element
+      if (elem[0].ghost->is_valid(ele)) {
+      
+        FEM_Elem *ghosts = (FEM_Elem*)elem[0].ghost;
+        const int *conn = ghosts->connFor(ele);
+        
+        int n1 = conn[0];
+        int n2 = conn[1];
+        int n3 = conn[2];
+
+        int numElementsOnEdge;
+      
+        numElementsOnEdge = countElementsOnEdge(n1,n2);
+        if (numElementsOnEdge==1) {
+          edgesOnBoundaryWithGhosts.insert(std::make_pair(n1,n2));
+        }
+      
+        numElementsOnEdge = countElementsOnEdge(n2,n3);
+        if (numElementsOnEdge==1) {
+          edgesOnBoundaryWithGhosts.insert(std::make_pair(n2,n3));
+        }
+      
+        numElementsOnEdge = countElementsOnEdge(n3,n1);
+        if (numElementsOnEdge==1) {
+          edgesOnBoundaryWithGhosts.insert(std::make_pair(n3,n1));
+        }
+         
+      }
+    }
+    
+    
+    
+    // Produce the final list of edges on the boundary(pruning things on the boundary that we can't verify are on the real boundary)
+  
+    for(std::set<std::pair<int,int> >::iterator iter=edgesOnBoundaryWithGhosts.begin(); iter!=edgesOnBoundaryWithGhosts.end(); ++iter){
+      int n1 = (*iter).first;  // The experimental new code for feature detection below has not yet been widely tested, but it does work for some example
+  
+      int n2 = (*iter).second;
+  
+      // Include any edge that adjacent to at least one local node. With a 1-deep node-based ghost layer, all such edges will be guaranteed to be in the true boundary
+      if(n1>=0 || n2>=0){
+        edgesOnBoundary.insert(make_pair(n1,n2));
+        verticesOnBoundary.insert(n1);
+        verticesOnBoundary.insert(n2);
+      }
+      
+    }
+   
+    
+  CkPrintf("%d: Found a total of %d edges(with at least one local node) on boundary of mesh: ", rank, edgesOnBoundary.size());
+  for(std::set<std::pair<int,int> >::iterator iter=edgesOnBoundary.begin();iter!=edgesOnBoundary.end();++iter){
+    int n1 = (*iter).first;
+    int n2 = (*iter).second;
+    double x1,y1,x2,y2;
+    node.get_coord(n1,x1,y1);
+    node.get_coord(n2,x2,y2);
+    CkPrintf("    %.2lf,%.2lf---%.2lf,%.2lf", x1,y1, x2,y2);
+  }
+  CkPrintf("\n");
+  
+  
+  
+  CkPrintf("%d: Found a total of %d nodes on boundary described earlier:", rank, verticesOnBoundary.size());
+  for(std::set<int>::iterator iter=verticesOnBoundary.begin();iter!=verticesOnBoundary.end();++iter){
+    int n = (*iter);
+    double x,y;
+    node.get_coord(n,x,y);
+    CkPrintf("    %.2lf,%.2lf  ", x,y);
+  }
+  CkPrintf("\n");
+    
+  
+  
+  // Find geometric corners of this mesh
+
+  double COS_CORNER_ANGLE_CUTOFF = cos(CORNER_ANGLE_CUTOFF);
+  
+  
+// build map from boundary vertex to adjacent boundary vertices
+//   std::map<int,int> verticesOnBoundary;
+//   std::map<std::pair<int,int>,int> edgesOnBoundary;
+//   
+  std::map<int,std::pair<int,int> >  vertexToAdjacentBoundaryVertices;
+  
+  for(std::set<std::pair<int,int> >::iterator iter=edgesOnBoundaryWithGhosts.begin();iter!=edgesOnBoundaryWithGhosts.end();++iter){
+    int n1 = (*iter).first;
+    int n2 = (*iter).second;
+    
+    // --------------- add n1's neighbor if n1 is local ---------------
+    if(n1>=0){
+       // see if we have an existing entry
+      std::map<int,std::pair<int,int> >::iterator v = vertexToAdjacentBoundaryVertices.find(n1);
+      if( v == vertexToAdjacentBoundaryVertices.end() ){
+        // no entry yet, create one
+        vertexToAdjacentBoundaryVertices[n1] = make_pair( n2, -1 );
+      } 
+      else if((*v).first == n2){
+        // do nothing, we've seen this one already
+      }
+      else{
+        // we have recorded a different adjacent vertex 
+        // already, but now we add this one
+        CkAssert((*v).second.second==-1);
+        (*v).second.second = n2;
+      }
+    }
+    
+    // --------------- add n2's neighbor if n2 is local---------------
+    // see if we have an existing entry
+    if(n2>=0){
+      std::map<int,std::pair<int,int> >::iterator v = vertexToAdjacentBoundaryVertices.find(n1);
+      v = vertexToAdjacentBoundaryVertices.find(n2);
+      if( v == vertexToAdjacentBoundaryVertices.end() ){
+        // no entry yet, create one
+        vertexToAdjacentBoundaryVertices[n2] = make_pair( n1, -1 );
+        
+      } 
+      else if((*v).first == n1){
+        // do nothing, we've seen this one already
+      }
+      else{
+        // we have recorded a different adjacent vertex 
+        // already, but now we add this one
+        CkAssert((*v).second.second==-1);
+        (*v).second.second = n1;
+      }
+    }
+    
+  }
+  
+  
+  // Iterate over all vertices on the boundary(these will include some ghost nodes adjacent to local nodes)
+  for(std::set<int>::iterator iter=verticesOnBoundary.begin();iter!=verticesOnBoundary.end();iter++){
+    int n = *iter;  
+    
+    
+    // Only consider the local nodes
+    if(n>=0){
+      
+      // Find the adjacent 2 edges on boundary
+      int n1 = vertexToAdjacentBoundaryVertices[n].first;
+      int n2 = vertexToAdjacentBoundaryVertices[n].second;
+      
+      // Get the coordinates of the 3 nodes
+      double x1,y1,x2,y2,x,y;
+      
+      node.get_coord(n,x,y);
+      node.get_coord(n1,x1,y1);
+      node.get_coord(n2,x2,y2);
+      
+      // Compute the angle between the two edges
+      
+      double v1x = x1-x; // a vector from node n to n1
+      double v1y = y1-y;
+      
+      double v2x = x2-x; // a vector from node n to n2
+      double v2y = y2-y;
+      
+      double cos_theta = (v1x*v2x+v1y*v2y)/(sqrt(v1x*v1x+v1y*v1y) * sqrt(v2x*v2x+v2y*v2y) );
+      
+          
+      CkPrintf("nodes (%.2lf,%.2lf)   (%.2lf,%.2lf)   (%.2lf,%.2lf)  form angle cos=%.3lf\n", x,y,x1,y1,x2,y2, cos_theta);
+      
+      // Apply cuttoff criterion and save 
+        
+      if(cos_theta > COS_CORNER_ANGLE_CUTOFF) {
+        // Any angle less than the cutoff should be considered a corner
+        cornersOnBoundary.insert(n);
+      }
+    }
+    
+  }
+  
+  CkPrintf("Found a total of %d corner nodes on boundary of mesh:\n", cornersOnBoundary.size());
+  
+  for(std::set<int>::iterator iter=cornersOnBoundary.begin();iter!=cornersOnBoundary.end();++iter){
+    int n = (*iter);
+    double x,y;
+    node.get_coord(n,x,y);
+    CkPrintf(" %d=(%.2lf,%.2lf)", n, x,y);
+  }
+  CkPrintf("\n");
+  
+  
+}
+
+
+
+