add: Additional test cases, method to choose level set spacing

Author: cgoates

scripts/clevis_quad.jou

@@ -0,0 +1,27 @@
+reset
+import step "/Users/caleb/Downloads/bracket(1).step" noheal color
+surf all copy
+delete vol 1
+webcut body 3 4 with plane yplane offset -0.0127
+webcut body 18 19 with plane zplane offset 0.0254
+create volume surface all heal
+
+surf 37 38 39 40 scheme polyhedron
+composite create surface 17 25 24 23 20  keep angle 15
+
+surf all size 0.005
+
+#curve 101 71 70 97 72 81 interval 1
+#curve 101 71 70 97 72 81 scheme equal
+#mesh curve 101 71 70 97 72 81
+
+mesh surface 37 33 38 26 27 28 39 40 35
+
+delete mesh surface 37  propagate
+copy mesh surface 38  onto  surface 37  source curve 146  source vertex 85  target curve 145  target vertex 58  
+
+delete mesh surface 40  propagate
+copy mesh surface 39  onto  surface 40  source curve 151  source vertex 49  target curve 150  target vertex 69 
+
+#vol 21 scheme sweep source surf 37 33 38 26 27 28 39 40 35 target surf 41
+#mesh vol all

scripts/clevis_tet.jou

@@ -0,0 +1,17 @@
+reset
+import step "/Users/caleb/Downloads/bracket(1).step" noheal color
+surf all copy
+delete vol 1
+webcut body 3 4 with plane yplane offset -0.0127
+webcut body 2 3 4 5 14 15 16 17 18 19 with plane zplane offset 0.0254
+create volume surface all heal
+
+vol all scheme tet
+mesh vol all
+
+nodeset 1 add surface 26 27 28 39 40 41 42 53 54 55 56
+nodeset 1 name "start"
+nodeset 2 add surface 23 24 25 37 38 43 44
+nodeset 2 name "end"
+
+export abaqus "/Users/caleb/sweeps/attempt-sweep/test/data/clevis.inp"  mesh_only overwrite  everything 

scripts/flange_imprinted.jou

@@ -0,0 +1,34 @@
+reset
+import step "/Users/caleb/Downloads/flange(1).step" noheal color 
+create surface rectangle width 8 height 9 xplane 
+move body 2 y 4 z 3.5
+imprint body all
+
+create surface rectangle width 9 height 2 xplane 
+move body 3 y -7 z 3.5
+imprint body 3 1
+rotate Surface 25 angle 45  about Z include_merged 
+imprint body 3 1
+rotate Surface 25 angle 45  about Z include_merged 
+imprint body 3 1
+rotate Surface 25 angle 45  about Z include_merged 
+imprint body 3 1
+rotate Surface 25 angle 90  about Z include_merged 
+imprint body 3 1
+rotate Surface 25 angle 45  about Z include_merged 
+imprint body 3 1
+rotate Surface 25 angle 45  about Z include_merged 
+imprint body 3 1
+
+delete body all except 1
+
+vol 1 scheme tet
+mesh vol 1
+
+nodeset 1 add surface 5 4 3 2 1  
+nodeset 1 name "end"
+
+nodeset 2 add surface 7  
+nodeset 2 name "start"
+
+export abaqus "/Users/caleb/sweeps/attempt-sweep/test/data/flange_imprinted.inp"  mesh_only overwrite  everything 

src/main.cpp

@@ -157,6 +157,83 @@ std::vector<int> verifyDelaunayTetrahedralization( const SimplicialComplex& mesh
     return failed_tetrahedra;
 }
 
+std::vector<double> computeEquallySpacedLevelSetValues( const topology::TetMeshCombinatorialMap& mesh,
+                                                        const Eigen::VectorXd& vertex_values,
+                                                        const double max_spacing )
+{
+    // 1. Compute gradient magnitude for each triangle (constant per triangle)
+    std::vector<double> gradient_magnitudes;
+    std::vector<std::pair<double, double>> triangle_value_ranges;
+
+    const topology::CombinatorialMapBoundary bdry( mesh );
+
+    const auto vert_ids = indexingOrError( bdry, 0 );
+
+    iterateCellsWhile( bdry, 2, [&]( const topology::Face& face ) {
+        const auto tri = triangleOfFace( mesh, face );
+
+        const double f0 = vertex_values( vert_ids( topology::Vertex( face.dart() ) ) );
+        const double f1 = vertex_values( vert_ids( topology::Vertex( phi( bdry, 1, face.dart() ).value() ) ) );
+        const double f2 = vertex_values( vert_ids( topology::Vertex( phi( bdry, -1, face.dart() ).value() ) ) );
+
+        const Eigen::Vector3d field_vals( f0, f1, f2 );
+        const Eigen::Vector3d grad = gradient( tri, field_vals );
+        const double grad_magnitude = grad.norm();
+
+        if( grad_magnitude < 1e-12 ) return true; // Skip degenerate/flat triangles
+
+        gradient_magnitudes.push_back( grad_magnitude );
+
+        const double f_min = std::min( { f0, f1, f2 } );
+        const double f_max = std::max( { f0, f1, f2 } );
+        triangle_value_ranges.push_back( { f_min, f_max } );
+        return true;
+    } );
+
+    // 2. Step through [0,1], at each step finding the minimum gradient
+    //    in the "active" triangles (those intersected by current level set)
+    std::vector<double> levels = { 0.0 };
+    double current_value = 0.0;
+
+    while( current_value < 1.0 - 1e-10 )
+    {
+        // Find min gradient among triangles that contain current_value
+        double min_grad = std::numeric_limits<double>::infinity();
+        for( size_t i = 0; i < gradient_magnitudes.size(); ++i )
+        {
+            if( triangle_value_ranges.at( i ).first <= current_value &&
+                current_value <= triangle_value_ranges.at( i ).second )
+            {
+                min_grad = std::min( min_grad, gradient_magnitudes.at( i ) );
+            }
+        }
+
+        if( std::isinf( min_grad ) || min_grad < 1e-12 )
+        {
+            // No active triangles or flat region
+            current_value = 1.0;
+        }
+        else
+        {
+            // To achieve max spacing d: dv = d * ||∇f||_min
+            const double dv = max_spacing * min_grad;
+            current_value += dv;
+        }
+
+        if( current_value < 1.0 )
+        {
+            levels.push_back( std::min( current_value, 1.0 ) );
+        }
+    }
+
+    if( levels.back() < 1.0 - 1e-10 )
+    {
+        levels.push_back( 1.0 );
+    }
+
+    return levels;
+}
+
 int main( int argc, char* argv[] )
 {
     const std::vector<std::string> input_args(argv + 1, argv + argc);
@@ -172,6 +249,10 @@ int main( int argc, char* argv[] )
                 return SweepInputTestCases::ventricle();
             else if( std::find( input_args.begin(), input_args.end(), "flange" ) != input_args.end() )
                 return SweepInputTestCases::flange();
+            else if( std::find( input_args.begin(), input_args.end(), "flange_imprinted" ) != input_args.end() )
+                return SweepInputTestCases::flange_imprinted();
+            else if( std::find( input_args.begin(), input_args.end(), "clevis" ) != input_args.end() )
+                return SweepInputTestCases::clevis();
             else if( std::find( input_args.begin(), input_args.end(), "cube" ) != input_args.end() )
                 return SweepInputTestCases::twelveTetCube();
             else if( std::find( input_args.begin(), input_args.end(), "bunny" ) != input_args.end() )
@@ -356,6 +437,9 @@ int main( int argc, char* argv[] )
             sweepEmbedding( map, sweep_input.zero_bcs, sweep_input.one_bcs, normals );
         std::cout << "Laplace time 2: " << t.stop( 4 ) << std::endl;
 
+        std::cout << "Level set values: " << computeEquallySpacedLevelSetValues( map, ans, 0.7 ) << std::endl;
+
+
         const Eigen::Matrix3Xd grad = gradientsWithBoundaryCorrection( map, sides, ans, normals );
 
         if( std::find( input_args.begin(), input_args.end(), "output-laplace" ) != input_args.end() )

test/SimplicialComplexTestCases.hpp

@@ -229,6 +229,16 @@ class SweepInputTestCases
         return sweep_input;
     }
 
+    static SweepInput flange_imprinted()
+    {
+        return io::loadINPFile( SRC_HOME "/test/data/flange_imprinted.inp", "start", "end" );
+    }
+
+    static SweepInput clevis()
+    {
+        return io::loadINPFile( SRC_HOME "/test/data/clevis.inp", "start", "end" );
+    }
+
     static SweepInput bunny()
     {
         SweepInput sweep_input = io::loadINPFile( SRC_HOME "/test/data/stanford_bunny.inp", "SS1_S3", "SS2_S3" );