First cut of GradateMesh, moved SetReferenceMesh into fem_adapt_algs
authorAaron Becker <abecker3@illinois.edu>
Tue, 25 Oct 2005 18:59:02 +0000 (18:59 +0000)
committerAaron Becker <abecker3@illinois.edu>
Tue, 25 Oct 2005 18:59:02 +0000 (18:59 +0000)
src/libs/ck-libs/fem/fem_adapt_algs.C
src/libs/ck-libs/fem/fem_adapt_algs.h
src/libs/ck-libs/fem/fem_adapt_if.h

index 7601a840d2f32a98cd3c0d37171e0bf7b80d84f9..a4dead0dc2179c2b1739e40d467f65fbc3d7bffc 100644 (file)
@@ -246,16 +246,108 @@ void FEM_Adapt_Algs::SetMeshSize(int method, double factor, double *sizes)
       theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
     }
   }
       theMesh->elem[0].setMeshSizing(i, factor*avgEdgeLength);
     }
   }
-  else if (method == 3) { // scale existing mesh size
-    for (int i=0; i<numElements; i++) {
-      theMesh->elem[0].setMeshSizing(i, factor*theMesh->elem[0].getMeshSizing(i));
-    }
-  }
-  else if (method == 4) { // mesh sizing has been set independently; use as is
+  else if (method == 3) { // mesh sizing has been set independently; use as is
     return;
   }
 }
 
     return;
   }
 }
 
+
+void FEM_Adapt_Algs::SetReferenceMesh()
+{
+    // for each element, set its size to its average edge length
+    // TODO: do we need to run this loop for element types other than 0?
+    
+    double avgLength;
+    int width = theMesh->elem[0].getConn().width();
+    int* eConn = (int*)malloc(width*sizeof(int));
+    int numElements = theMesh->elem[0].size();
+    
+    for (int i=0; i<numElements; ++i, avgLength=0) {
+        theMesh->e2n_getAll(i, eConn);
+        
+        for (int j=0; j<width-1; ++j) {
+            avgLength += length(eConn[j], eConn[j+1]);
+        }
+        avgLength += length(eConn[0], eConn[width-1]);
+        avgLength /= width;
+        theMesh->elem[0].setMeshSizing(i, avgLength);      
+    }
+    free(eConn);
+}
+
+
+void FEM_Adapt_Algs::GradateMesh(double smoothness)
+{
+    // Resize mesh elements to avoid jumps in element size
+    // Algorithm based on h-shock correction, described in
+    // Mesh Gradation Control, Borouchaki et al
+    // IJNME43 1998 www.ann.jussieu.fr/~frey/publications/ijnme4398.pdf
+
+    const double beta = smoothness;
+
+    double maxShock = 0;
+    double minShock = 1e10;
+    int iteration = 0;
+
+    int* adjNodes, *boundNodes;
+    int nadjNodes, nnodes;
+    int meshNum = FEM_Mesh_default_read();
+
+    nnodes = theMesh->node.size();
+    boundNodes = new int[nnodes];
+    FEM_Mesh_data(meshNum, FEM_NODE, FEM_BOUNDARY, 
+            boundNodes, 0, nnodes, FEM_INT, 1);
+
+
+    printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
+    while (maxShock > beta) {
+        
+        for (int node=0; node<nnodes; ++node) {
+            // only do internal nodes (?)
+            if (boundNodes[node] < 0 || !FEM_is_valid(meshNum, FEM_NODE, node))
+                continue;
+            
+            theMesh->n2n_getAll(node, &adjNodes, &nadjNodes);
+            for (int adjNode=0; adjNode<nadjNodes; ++adjNode) {
+                // get length of edge
+                // (do we need to worry about ghost/non-ghost?)
+                double edgelen = length(node, adjNode);
+                
+                // get adjacent elemnents and their sizes
+                int e1, e2;
+                theMesh->get2ElementsOnEdge(node, adjNode, &e1, &e2);
+                double s1, s2;
+                s1 = theMesh->elem[0].getMeshSizing(e1);
+                s2 = theMesh->elem[0].getMeshSizing(e2);
+                
+                // h-shock=max(size ratio)^(1/edge length)
+                double ratio = (s1 > s2) ? s1/s2 : s2/s1;
+                double hs = pow(ratio, 1.0/edgelen);
+                if (hs > maxShock) maxShock = hs;
+                if (hs < minShock) minShock = hs;
+                
+                // if hs > beta, resize the larger elt:
+                // new size = old size / eta^2
+                // eta = (beta / h-shock)^(edge length)
+                //     = beta^(edge length) / size ratio
+                if (hs > beta) {
+                    double etasq = pow(pow(beta, edgelen) / ratio, 2.0);
+                    if (s1 > s2) {
+                        theMesh->elem[0].setMeshSizing(e1, etasq);
+                    } else {
+                        theMesh->elem[0].setMeshSizing(e2, etasq);
+                    }
+                }
+            }
+        }
+        printf("Finished iteration %d\n", iteration++);
+        printf("Max shock:%8.3f\n", maxShock);
+        printf("Min shock:%8.3f\n", minShock);
+        printf("Target:%8.3f\n", beta);
+    }
+}
+
+
 int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
   int noEle = theMesh->elem[0].size();
   int *con = (int*)malloc(3*sizeof(int));
 int FEM_Adapt_Algs::simple_refine(double targetA, double xmin, double ymin, double xmax, double ymax) {
   int noEle = theMesh->elem[0].size();
   int *con = (int*)malloc(3*sizeof(int));
index 38208c19739280f12f901b1120a4fdc68a963d16..8a20849bf2a656f2fcffff10c152ed1ea785c3fe 100644 (file)
@@ -100,6 +100,11 @@ class FEM_Adapt_Algs {
       if method = 1, keep regional mesh sizes, and scale by factor
       if method = 2, uses sizes to size mesh by regions */
   void FEM_Remesh(int qm, int method, double factor, double *sizes);
       if method = 1, keep regional mesh sizes, and scale by factor
       if method = 2, uses sizes to size mesh by regions */
   void FEM_Remesh(int qm, int method, double factor, double *sizes);
+  
+  /// Set sizes on mesh elements based on their average edge length
+  void SetReferenceMesh();
+  /// Adjust sizes on mesh elements to avoid sharp discontinuities
+  void GradateMesh(double smoothness);
  private:
   // Helper methods
   /// Performs refinement; returns number of modifications
  private:
   // Helper methods
   /// Performs refinement; returns number of modifications
index 3512c690e3039005a9cdf1ea3bb26a204699e9e6..de061371e595510cff5d7814f39bb454d2b73da8 100644 (file)
@@ -64,30 +64,22 @@ FDECL void FTN_NAME(FEM_ADAPT_SETELEMENTSSIZEFIELD,fem_adapt_setelementssizefiel
 
 void FEM_ADAPT_SetReferenceMesh(int meshID) {
     FEM_Mesh* mesh = FEM_Mesh_lookup(meshID, "FEM_ADAPT_Init");
 
 void FEM_ADAPT_SetReferenceMesh(int meshID) {
     FEM_Mesh* mesh = FEM_Mesh_lookup(meshID, "FEM_ADAPT_Init");
-
-    // for each element, set its size to its average edge length
-    // TODO: do we need to run this loop for element types other than 0?
-    
-    double avgLength;
-    int width = mesh->elem[0].getConn().width();
-    int* eConn = (int*)malloc(width*sizeof(int));
-    int numElements = mesh->elem[0].size();
-    
-    for (int i=0; i<numElements; ++i, avgLength=0) {
-        mesh->e2n_getAll(i, eConn);
-        
-        for (int j=0; j<width-1; ++j) {
-            avgLength += CtvAccess(_adaptAlgs)->length(eConn[j], eConn[j+1]);
-        }
-        avgLength += CtvAccess(_adaptAlgs)->length(eConn[0], eConn[width-1]);
-        avgLength /= width;
-        mesh->elem[0].setMeshSizing(i, avgLength);      
-    }
-    free(eConn);
+    mesh->getfmMM()->getfmAdaptAlgs()->SetReferenceMesh();
 }
 FDECL void FTN_NAME(FEM_ADAPT_SETREFERENCEMESH, fem_adapt_setreferencemesh)(int* meshID)
 {
     FEM_ADAPT_SetReferenceMesh(*meshID);
 }
 
 }
 FDECL void FTN_NAME(FEM_ADAPT_SETREFERENCEMESH, fem_adapt_setreferencemesh)(int* meshID)
 {
     FEM_ADAPT_SetReferenceMesh(*meshID);
 }
 
+
+void FEM_ADAPT_GradateMesh(int meshID, double smoothness)
+{
+    FEM_Mesh* mesh = FEM_Mesh_lookup(meshID, "FEM_ADAPT_Init");
+    mesh->getfmMM()->getfmAdaptAlgs()->GradateMesh(smoothness);
+}
+FDECL void FTN_NAME(FEM_ADAPT_GRADATEMESH, fem_adapt_gradatemesh)(int* meshID, double* smoothness)
+{
+    FEM_ADAPT_GradateMesh(*meshID, *smoothness);
+}
+
 #endif
 #endif