GradateMesh bug fixes
authorAaron Becker <abecker3@illinois.edu>
Wed, 26 Oct 2005 16:27:44 +0000 (16:27 +0000)
committerAaron Becker <abecker3@illinois.edu>
Wed, 26 Oct 2005 16:27:44 +0000 (16:27 +0000)
src/libs/ck-libs/fem/fem_adapt_algs.C

index a4dead0dc2179c2b1739e40d467f65fbc3d7bffc..5281de3e4600fc948ac22749b9c702550f936a6c 100644 (file)
@@ -285,8 +285,7 @@ void FEM_Adapt_Algs::GradateMesh(double smoothness)
 
     const double beta = smoothness;
 
-    double maxShock = 0;
-    double minShock = 1e10;
+    double maxShock, minShock;
     int iteration = 0;
 
     int* adjNodes, *boundNodes;
@@ -300,7 +299,9 @@ void FEM_Adapt_Algs::GradateMesh(double smoothness)
 
 
     printf("Running h-shock mesh gradation with beta=%.3f\n", beta);
-    while (maxShock > beta) {
+    do {
+        maxShock = 0;
+        minShock = 1e10;
         
         for (int node=0; node<nnodes; ++node) {
             // only do internal nodes (?)
@@ -311,40 +312,55 @@ void FEM_Adapt_Algs::GradateMesh(double smoothness)
             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);
+                double edgelen = length(node, adjNodes[adjNode]);
+                if (edgelen > 1e6) {
+                    printf("Edge (%d, %d) has length %8.3f\n", node, adjNodes[adjNode], edgelen);
+                }
                 
                 // get adjacent elemnents and their sizes
                 int e1, e2;
-                theMesh->get2ElementsOnEdge(node, adjNode, &e1, &e2);
+                theMesh->get2ElementsOnEdge(node, adjNodes[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)
+                CkAssert(s1 > 1e-6 && s2 > 1e-6 && "Size 0 element");
+                CkAssert(edgelen > 1e-6 && "Length 0 edge");
+                CkAssert(edgelen < 1e6 && "Length inf edge");
+
                 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 > 100) {
+                    printf("hs: %8.5f\ns1: %8.5f\ns2: %8.5f\nedgelen: %8.5f\nratio: %8.5f\n", hs, s1, s2, edgelen, ratio);
+                }
                 
                 // 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);
+                    double etasq = pow(beta, edgelen) / ratio;
+                    etasq *= etasq;
                     if (s1 > s2) {
-                        theMesh->elem[0].setMeshSizing(e1, etasq);
+                        theMesh->elem[0].setMeshSizing(e1, s1 / etasq);
                     } else {
-                        theMesh->elem[0].setMeshSizing(e2, etasq);
+                        theMesh->elem[0].setMeshSizing(e2, s2 / etasq);
                     }
                 }
             }
-        }
+            delete[] adjNodes;
+        } 
         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);
-    }
+    } while (maxShock > beta);
+
+    delete[] boundNodes;
 }