Add AirgapTorque : PostProcessorElement

    Torque calculation from body forces is complicated by discontinuities in field gradients across material boundaries

Fix for Mesh::insert_internal_boundaries with small airgaps
    This method is probably not 100% fixed, but it is more stable than before (SalientPoleSynchrel regression tests run with CylindricalAirgap<2>)

src/Mesh/src/BoundaryConstraint.cpp

@@ -7,7 +7,7 @@ bool BoundaryConstraint::is_uniform(Mesh const &m) const {
     double_t delta_0{dc.delta()};
 
     for (size_t i : DartConstraints) {
-        double delta_i{m.dart_constraint(i).delta()};
+        double_t delta_i{m.dart_constraint(i).delta()};
         if (std::abs(delta_0 - delta_i) > FLT_EPSILON) {
             return false;
         }
@@ -16,7 +16,7 @@ bool BoundaryConstraint::is_uniform(Mesh const &m) const {
     return true;
 }
 
-double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
+double_t BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
     double_t delta_min{DBL_MAX};
 
     for (size_t i : DartConstraints) {
@@ -27,16 +27,16 @@ double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
     return delta_min;
 }
 
-double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
-    std::cout << "//TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued" << std::endl;
-    std::cout << "//TODO: Constrained dart enqueuing needs to be revisited, and need to make sure enqueued edges are unique" << std::endl;
-    double delta_max{delta_min};
+double_t BoundaryConstraint::queue_uniform(Mesh &m, double_t delta_min) const {
+    // TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued
+    // TODO: Constrained dart enqueuing needs to be revisited, and need to make sure enqueued edges are unique
+
+    double_t delta_max{delta_min};
 
     for (size_t i : DartConstraints) {
         DartConstraint const &dc = m.dart_constraint(i);
         double_t delta = std::abs(dc.S1 - dc.S0);
         if (delta > delta_min) {
-            //m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
             m.add_dart_to_queue(dc.dart());
             delta_max = std::max(delta_max, delta);
         }
@@ -46,12 +46,7 @@ double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
 }
 
 void BoundaryConstraint::add_to_queue(Mesh &m, size_t dart) const {
-    if (UniformDiscretization) {
-        queue_uniform(m, smallest_parametric_edge_length(m) / 2.0);
-    } else {
-        //m.add_dart_to_queue(dart);
-        m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
-    }
+    m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
 }
 
 void BoundaryConstraint::insert_parameters(Mesh &mesh, std::vector<double_t> s) {
@@ -63,7 +58,7 @@ void BoundaryConstraint::insert_parameters(Mesh &mesh, std::vector<double_t> s)
 
 void BoundaryConstraint::make_uniform(Mesh &m) { make_uniform(m, smallest_parametric_edge_length(m)); };
 
-void BoundaryConstraint::make_uniform(Mesh &m, double delta_min) {
+void BoundaryConstraint::make_uniform(Mesh &m, double_t delta_min) {
     delta_min *= (1 + FLT_EPSILON);
     double_t delta_max{DBL_MAX};
     while (delta_max > delta_min) {
@@ -119,12 +114,14 @@ void BoundaryConstraint::refine(Mesh &m, double_t tol) {
 
             double_t len_c = ConstraintCurve->length(dc.S0, dc.S1);
             double_t len_e = m.edge(dc.dart()).length(m);
-            double_t err_i = std::abs(len_c - len_e) / len_c;
+            double_t err_i = std::abs(len_c - len_e);
 
-            if (err_i > tol) {
+            if (len_c == 0) {
+                throw std::runtime_error("Detected a zero-length boundary constraint. Does your model have a small aspect ratio region?");
+            } else if (err_i > tol * len_c) {
                 m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
 
-                err_max = std::max(err_max, err_i);
+                err_max = std::max(err_max, err_i / len_c);
             }
         }
 

src/Mesh/src/Mesh.h

@@ -184,6 +184,8 @@ public:
         return std::vector<std::shared_ptr<MappedBoundaryPair>>(MappedBoundaries.end() - cbp_vec.size(), MappedBoundaries.end());
     }
 
+    std::shared_ptr<MappedBoundaryPair> mapped_boundary(size_t i) { return MappedBoundaries[i]; };
+
     DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
 
     DartTriangle const &triangle(size_t i) const { return Triangles[i]; };
@@ -384,6 +386,8 @@ private:
 
     void enforce_boundary_mappings();
 
+    void enforce_uniform_boundary_discretizations();
+
     void get_triangles();
 
     void insert_from_queue();

src/Physics/src/PostProcessorInterface.cpp

@@ -29,7 +29,7 @@ void IntegralAverage::build() {
     OperatorMatrix /= area;
 }
 
-double_t LortenzTorque::operator()(Oe::VectorXd const &v) {
+double_t LorentzTorque::operator()(Oe::VectorXd const &v) {
     auto pptr = Physics.lock();
 
     auto d1mg = pptr->derivatives();
@@ -68,5 +68,42 @@ double_t LortenzTorque::operator()(Oe::VectorXd const &v) {
         torque += (Jz * wts(q) * (x.col(q) * Hx + y.col(q) * Hy)).sum();
     }
 
+    return torque;
+}
+
+double_t AirgapTorque::operator()(Oe::VectorXd const &v) {
+    auto pptr = Physics.lock();
+
+    auto dmg = pptr->derivatives();
+    auto wts = pptr->weights();
+
+    Oe::ArrayXXd x = Oe::ArrayXXd::Zero(pptr->domain()->elements_size(),pptr->domain()->quadrature_size());
+    Oe::ArrayXXd y = Oe::ArrayXXd::Zero(pptr->domain()->elements_size(),pptr->domain()->quadrature_size());
+
+    for (size_t const &j : Region->elements()) {
+        auto qpts = pptr->domain()->quadrature_points(j);
+        for (size_t q = 0; q != qpts.size(); ++q) {
+            x(j, q) = qpts[q].x();
+            y(j, q) = qpts[q].y();
+        }
+    }
+
+    double_t torque{0.0};
+    for (size_t q = 0; q != pptr->domain()->quadrature_size(); ++q) {
+        Oe::ArrayXd By = dmg.dx(q).transpose() * v;
+        Oe::ArrayXd Bx = -dmg.dy(q).transpose() * v;
+
+        for (size_t const &j : Region->elements()) {
+            double_t r = hypot(x(j,q),y(j,q));
+            double_t a = atan2(y(j,q),x(j,q));
+
+            double_t Br = Bx(j) * cos(a) + By(j) * sin(a);
+            double_t Bt = By(j) * cos(a) - Bx(j) * sin(a);
+
+            torque += wts(q)(j) * Br * Bt * r;
+        }
+    }
+    torque *= (Scaling * Length) / (mu_0 * (OuterRadius - InnerRadius));
+
     return torque;
 }
\ No newline at end of file

src/Physics/src/PostProcessorInterface.h

@@ -36,9 +36,9 @@ protected:
     Oe::VectorXd OperatorMatrix;
 };
 
-class LortenzTorque : public PostProcessorElement {
+class LorentzTorque : public PostProcessorElement {
 public:
-    LortenzTorque(std::shared_ptr<PhysicsInterface> p, std::vector<std::shared_ptr<DiscreteRegion>> r) : PostProcessorElement{p}, Regions{r} {};
+    LorentzTorque(std::shared_ptr<PhysicsInterface> p, std::vector<std::shared_ptr<DiscreteRegion>> r) : PostProcessorElement{p}, Regions{r} {};
 
     void build() override {}; // No-op for now...
 
@@ -49,6 +49,23 @@ protected:
 
 };
 
+class AirgapTorque : public PostProcessorElement {
+public:
+    AirgapTorque(std::shared_ptr<PhysicsInterface> p, std::shared_ptr<DiscreteRegion> r, double_t ri, double_t ro, double_t ls, double_t s)
+            : PostProcessorElement{p}, Region{r}, InnerRadius{ri}, OuterRadius{ro}, Length{ls}, Scaling{s} {};
+
+    void build() override {}; // No-op
+
+    double_t operator()(Oe::VectorXd const &v) override;
+
+protected:
+    std::shared_ptr<DiscreteRegion> Region;
+    double_t InnerRadius;
+    double_t OuterRadius;
+    double_t Length;
+    double_t Scaling;
+};
+
 class PostProcessorInterface {
 public:
     double_t operator()(std::string key, Oe::VectorXd const &val) { return (*Elements[key])(val); };

src/Sketch/src/CylindricalAirgap.h

@@ -9,6 +9,7 @@
 #include "Angle.h"
 
 // TODO: Move to a library for higher order functions
+// TODO: Need a better meshing strategy for small gaps
 
 template <size_t Nl>
 class CylindricalAirgap {
@@ -40,48 +41,39 @@ public:
         double_t r = ri;
         double_t dr = (ro - ri) / (Nl + 1);
 
-        auto v0 = vi0;
-        auto v1 = vi1;
-        for (size_t i = 0; i != Nl; ++i) {
+        r += dr;
+
+        // create airgap Vertices
+        Vertex0[0] = sketch.new_element<Vertex>(r * std::cos(aa), r * std::sin(aa));
+        Vertex1[0] = sketch.new_element<Vertex>(r * std::cos(aa + da), r * std::sin(aa + da));
+
+        // create airgap edge LineSegments
+        Line0[0] = sketch.new_element<LineSegment>(vi0, Vertex0[0]);
+        Line1[0] = sketch.new_element<LineSegment>(vi1, Vertex1[0]);
+
+        // create airgap interface CircularArc
+        Interface[0] = sketch.new_element<CircularArc>(Vertex0[0], Vertex1[0], vc, r);
+
+        for (size_t i = 1; i != Nl; ++i) {
             r += dr;
 
-            // Create Airgap Vertices
+            // create airgap Vertices
             Vertex0[i] = sketch.new_element<Vertex>(r * std::cos(aa), r * std::sin(aa));
             Vertex1[i] = sketch.new_element<Vertex>(r * std::cos(aa + da), r * std::sin(aa + da));
 
-            // Create Airgap Edge LineSegments
-            Line0[i] = sketch.new_element<LineSegment>(v0, Vertex0[i]);
-            Line1[i] = sketch.new_element<LineSegment>(v1, Vertex1[i]);
+            // create airgap edge LineSegments
+            Line0[i] = sketch.new_element<LineSegment>(Vertex0[i-1], Vertex0[i]);
+            Line1[i] = sketch.new_element<LineSegment>(Vertex1[i-1], Vertex1[i]);
 
-            // Create Airgap CircularArc Interface
+            // create airgap interface CircularArc interface
             Interface[i] = sketch.new_element<CircularArc>(Vertex0[i], Vertex1[i], vc, r);
-
-            // Next
-            v0 = Vertex0[i];
-            v1 = Vertex1[i];
         }
 
-        // Last line segments
-        Line0[Nl+1] = sketch.new_element<LineSegment>(v0,vo0);
-        Line1[Nl+1] = sketch.new_element<LineSegment>(v1,vo1);
-
-            /*
-            // Create Airgap Vertices
-            Vertex0 = sketch.new_element<Vertex>(ra * std::cos(aa), ra * std::sin(aa));
-            Vertex1 = sketch.new_element<Vertex>(ra * std::cos(aa + da), ra*std::sin(aa + da));
-
-            // Create Airgap Edge LineSegments
-            Line0ia = sketch.new_element<LineSegment>(vi0, Vertex0);
-            Line0ao = sketch.new_element<LineSegment>(Vertex0, vo0);
-
-            Line1ia = sketch.new_element<LineSegment>(vi1, Vertex1);
-            Line1ao = sketch.new_element<LineSegment>(Vertex1, vo1);
-
-            // Create Airgap CircularArc Interface
-            AirgapInterface = sketch.new_element<CircularArc>(Vertex0, Vertex1, vc, ra);
+        // last line segments
+        Line0[Nl] = sketch.new_element<LineSegment>(Vertex0[Nl-1],vo0);
+        Line1[Nl] = sketch.new_element<LineSegment>(Vertex1[Nl-1],vo1);
 
-            // TODO: Should constraints be added?
-             */
+        // TODO: Should constraints be added?
     }
 
     std::shared_ptr<Vertex> vertex0(size_t i) const { return Vertex0[i]; };

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -240,8 +240,8 @@ public:
 
         me.create();
 
-        me.MaximumElementSize = 0.1;
-        me.MinimumElementSize = 0.01;
+        me.MaximumElementSize = 0.05;
+        me.MinimumElementSize = 0.005;
         me.MinimumElementQuality = 0.5;
 
         for (auto &pb : periodic_boundary) {
@@ -982,30 +982,16 @@ public:
         EXPECT_LE(residual, FLT_EPSILON * rro);
 
         // Rotor part of airgap
-        std::cout << vr1->x() << "," << vr1->y() << std::endl;
-        std::cout << vs0->x() << "," << vs0->y() << std::endl;
-
-        CylindricalAirgap<1> cyag{sk,origin,vr1,vs0,90.0};
-
-        //auto vra0 = sk.new_element<Vertex>(rro + airgap / 2.0, 0.0);
-        //auto vra1 = sk.new_element<Vertex>(0.0, rro + airgap / 2.0);
-
-        auto fra0 = sk.new_element<Fixation>(cyag.vertex0(0));
-        auto fra1 = sk.new_element<Fixation>(cyag.vertex1(0));
-
-        //auto vr1p = sk.select_periodic_vertex(vr1, origin, 90.0);
-        //auto lra0 = sk.new_element<LineSegment>(vr1, vra0);
-        //auto lra1 = sk.new_element<LineSegment>(vr1p, vra1);
+        size_t const Nl{2};
+        CylindricalAirgap<Nl> cyag{sk, origin, vr1, vs0, 90.0};
+        //CylindricalAirgap cyag{sk, origin, vr1, vs0, 90.0};
 
-        //continuous_airgap = sk.new_element<CircularArc>(vra0, vra1, origin, rro + airgap / 2.0);
-        continuous_airgap = cyag.interface(0);
-
-        // Stator part of airgap
-        //auto vsa0 = vs0;
-        //auto vsa1 = sk.select_periodic_vertex(vsa0, origin, 90.0);
+        for (size_t i = 0; i != Nl; ++i) {
+            sk.new_element<Fixation>(cyag.vertex0(i));
+            sk.new_element<Fixation>(cyag.vertex1(i));
+        }
 
-        //auto lsa0 = sk.new_element<LineSegment>(vra0, vsa0);
-        //auto lsa1 = sk.new_element<LineSegment>(vra1, vsa1);
+        continuous_airgap = cyag.interface(1);
 
         // Build and Save
         residual = sk.solve();
@@ -1039,21 +1025,27 @@ public:
         auto discrete_airgap = me.boundary_constraint(continuous_airgap);
         discrete_airgap->uniform_discretization(true);
 
-        // Create initial mesh
-        me.create();
-        me.MaximumElementSize = 5.0e-3;
+        // Create initial mesh;
         me.MinimumElementSize = airgap;
+
+        me.create();
+
+        EXPECT_TRUE(me.edges_are_valid());
+        EXPECT_TRUE(me.edges_are_optimal());
+
+        me.save_as(SAVE_DIR, std::string("salient_pole_synchrel_mesh"));
+
+        //me.MaximumElementSize = 5.0e-3;
         me.MinimumElementQuality = 0.5;
 
         // Refine and save
         me.refine();
+
         me.save_as(SAVE_DIR, std::string("salient_pole_synchrel_mesh"));
 
         // Create FiniteElementMesh
         fem = std::make_shared<FiniteElementMesh<2, 2>>(me); // TODO: Look at speeding this up somehow, try row major sparse matrix operators, sparse vectors for evaluation of magnetization
 
-        //TODO: Force density is F = curl(B) X H
-
         for (auto &pb : periodic_boundary) {
             EXPECT_EQ(fem->boundary(pb.curve0()).size(), 1);
             EXPECT_EQ(fem->boundary(pb.curve1()).size(), 1);
@@ -1108,6 +1100,10 @@ public:
 
         stator_iron = me.select_contour(Point{r_stator_iron * cos(M_PI / 4.0), r_stator_iron * cos(M_PI / 4.0)});
         ASSERT_TRUE(stator_iron);
+
+        double_t r_agc = rsi*1.0/2.0 + rro*1.0/2.0;
+        airgap_contour = me.select_contour(Point{r_agc*cos(M_PI/4.0),r_agc*sin(M_PI/4.0)});
+        ASSERT_TRUE(airgap_contour);
     }
 
     double_t poles;
@@ -1132,6 +1128,7 @@ public:
     std::shared_ptr<Contour const> phase_c_contour;
     std::shared_ptr<Contour const> rotor_iron;
     std::shared_ptr<Contour const> stator_iron;
+    std::shared_ptr<Contour const> airgap_contour;
 
     decltype(sk.select_periodic_boundary_pairs(origin, 90.0)) periodic_boundary;
     decltype(sk.select_radial_boundary(origin, rso)) outer_boundary;
@@ -1142,6 +1139,7 @@ TEST_F(Salient_Pole_Synchrel, Test) {
 
     // Set material properties
     MaterialProperties electrical_steel = JFE_35JN210();
+    //MaterialProperties electrical_steel = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
     fem->region(rotor_iron)->material(electrical_steel);
     fem->region(stator_iron)->material(electrical_steel);
 
@@ -1168,6 +1166,8 @@ TEST_F(Salient_Pole_Synchrel, Test) {
     // Post Processor
     PostProcessorInterface pp;
 
+    pp.add<AirgapTorque>(msph,"Torque",fem->region(airgap_contour), rro*2.0/3.0 + rsi*1.0/3.0, rro*1.0/3.0 + rsi*2.0/3.0, 1.0, 4.0);
+
     // Solve
     Oe::VectorXd v_prev;
     msph->assemble();
@@ -1191,11 +1191,13 @@ TEST_F(Salient_Pole_Synchrel, Test) {
         fem->write_vector(Bx, By, SAVE_DIR, std::string("salient_pole_synchrel_B_") + std::to_string(iter));
 
         // Increment Position
-        *position += 2;
-        t += 2 * dt;
-        fargs["t"] = t;
+        *position += 32;
+        //t += 32 * dt;
+        //fargs["t"] = t;
 
         msph->assemble();
+
+        std::cout << iter << "," << pp("Torque",solution->v) << std::endl;
     }
 }
 
@@ -1206,8 +1208,8 @@ TEST_F(Salient_Pole_Synchrel, Test_OpenMP) {
 
 
     #ifdef _OPENMP
-    omp_set_num_threads(omp_get_max_threads() / 4);
-    size_t N = (size_t)(std::sqrt(omp_get_max_threads() / 4));
+    size_t N = (size_t)sqrt((omp_get_max_threads() / 4)); // NxN matrix of t/J values
+    omp_set_num_threads(N * N);
 
     // Input Arguments
     std::vector <double_t> t;
@@ -1223,6 +1225,7 @@ TEST_F(Salient_Pole_Synchrel, Test_OpenMP) {
 
     // Materials
     MaterialProperties electrical_steel = JFE_35JN210();
+    //MaterialProperties electrical_steel = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
     fem->region(rotor_iron)->material(electrical_steel);
     fem->region(stator_iron)->material(electrical_steel);
 

test/LibraryIntegration/test_Annuli.cpp

@@ -281,8 +281,8 @@ public:
         ms->add_magnetic_insulation({outer_boundaries[0],outer_boundaries[1],outer_boundaries[2],outer_boundaries[3]});
 
         PostProcessorInterface pp;
-        pp.add<LortenzTorque>(ms, "Inner Torque", std::vector{fem->region(inner_region)});
-        pp.add<LortenzTorque>(ms, "Outer Torque", std::vector{fem->region(outer_region)});
+        pp.add<LorentzTorque>(ms, "Inner Torque", std::vector{fem->region(inner_region)});
+        pp.add<LorentzTorque>(ms, "Outer Torque", std::vector{fem->region(outer_region)});
 
         ms->assemble();
 

test/LibraryIntegration/test_HalfAnnuli.cpp

@@ -51,14 +51,14 @@ public:
     void SetUp_Mesh() {
         me = Mesh(sk);
 
-        me.add_mapped_boundary_pair(sk.select_periodic_boundary_pairs(origin, 90.0));
+        me.add_mapped_boundary_pair(sk.select_periodic_boundary_pairs(origin, 180.0));
 
         me.boundary_constraint(inner_boundary)->uniform_discretization(true);
 
         me.create();
 
-        me.MaximumElementSize = 0.125;
-        me.MinimumElementSize = me.MaximumElementSize / 8.0;
+        me.MaximumElementSize = 1.0 / 16.0;
+        me.MinimumElementSize = me.MaximumElementSize / 20.0;
         me.MinimumElementQuality = 0.5;
 
         bool result = me.refine();
@@ -115,11 +115,11 @@ public:
         auto pb = ms->add_periodic_boundary(sk.select_periodic_boundary_pairs(origin, 180.0), false);
         for (size_t i = 0; i != pb->map().size(); ++i) {
             auto vm = pb->map()[i];
-            /*
-            std::cout << vm << std::endl;
-            std::cout << fem->node(vm.first()) << std::endl;
-            std::cout << fem->node(vm.second()) << std::endl;
-             */
+
+            //std::cout << vm << std::endl;
+            //std::cout << fem->node(vm.first()) << std::endl;
+            //std::cout << fem->node(vm.second()) << std::endl;
+
             auto n0 = fem->node(vm.first());
             auto n1 = fem->node(vm.second());
 
@@ -237,8 +237,8 @@ public:
         //double_t dt = 1.0 / (2.0 * position->size());
 
         PostProcessorInterface pp;
-        pp.add<LortenzTorque>(ms, "Inner Torque", std::vector{fem->region(inner_region)});
-        pp.add<LortenzTorque>(ms, "Outer Torque", std::vector{fem->region(outer_region)});
+        pp.add<LorentzTorque>(ms, "Inner Torque", std::vector{fem->region(inner_region)});
+        pp.add<LorentzTorque>(ms, "Outer Torque", std::vector{fem->region(outer_region)});
 
         ms->assemble();
 

test/Mesh/Mesh_Test.cpp

@@ -431,6 +431,79 @@ TEST(Mesh, create_horseshoe_domain) {
     m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
 }
 
+TEST(Mesh, create_three_layer_airgap) {
+    Sketch sketch;
+
+    auto v0 = sketch.new_element<Vertex>(0.0,0.0);
+
+    auto v1 = sketch.new_element<Vertex>(80.00,0.0);
+    auto v2 = sketch.new_element<Vertex>(80.25,0.0);
+    auto v3 = sketch.new_element<Vertex>(80.50,0.0);
+    auto v4 = sketch.new_element<Vertex>(80.75,0.0);
+
+    auto v5 = sketch.new_element<Vertex>(0.0,80.00);
+    auto v6 = sketch.new_element<Vertex>(0.0,80.25);
+    auto v7 = sketch.new_element<Vertex>(0.0,80.50);
+    auto v8 = sketch.new_element<Vertex>(0.0,80.75);
+
+    auto l0 = sketch.new_element<LineSegment>(v0,v1);
+    auto l1 = sketch.new_element<LineSegment>(v1,v2);
+    auto l2 = sketch.new_element<LineSegment>(v2,v3);
+    auto l3 = sketch.new_element<LineSegment>(v3,v4);
+
+    auto l4 = sketch.new_element<LineSegment>(v0,v5);
+    auto l5 = sketch.new_element<LineSegment>(v5,v6);
+    auto l6 = sketch.new_element<LineSegment>(v6,v7);
+    auto l7 = sketch.new_element<LineSegment>(v7,v8);
+
+    auto c0 = sketch.new_element<CircularArc>(v1,v5,v0,80.00);
+    //auto c1 = sketch.new_element<CircularArc>(v2,v6,v0,80.25);
+    auto c2 = sketch.new_element<CircularArc>(v3,v7,v0,80.50);
+    //auto c3 = sketch.new_element<CircularArc>(v4,v8,v0,80.75);
+
+    auto v9 = sketch.new_element<Vertex>(132.0,0.0);
+    auto v10 = sketch.new_element<Vertex>(0.0,132.0);
+    auto l8 = sketch.new_element<LineSegment>(v4,v9);
+    auto l9 = sketch.new_element<LineSegment>(v8,v10);
+    auto c4 = sketch.new_element<CircularArc>(v9,v10,v0,132.0);
+
+    double_t r = 80.25;
+    double_t da = M_PI / 2.0 / 4.0;
+    auto v11 = sketch.new_element<Vertex>(r * cos(da), r * sin(da));
+    auto v12 = sketch.new_element<Vertex>(r * cos(2 * da), r * sin(2 * da));
+    auto v13 = sketch.new_element<Vertex>(r * cos(3 * da), r * sin(3 * da));
+    auto c5 = sketch.new_element<CircularArc>(v2,v11,v0,r);
+    auto c6 = sketch.new_element<CircularArc>(v11,v12,v0,r);
+    auto c7 = sketch.new_element<CircularArc>(v12,v13,v0,r);
+    auto c8 = sketch.new_element<CircularArc>(v13,v6,v0,r);
+
+    size_t N = 4;
+    r = 80.75;
+    da = M_PI / 2.0 / N;
+    std::vector<std::shared_ptr<Vertex>> v8075;
+    for (size_t i = 1; i != N; ++i) {
+        v8075.push_back(sketch.new_element<Vertex>(r * cos(i * da), r * sin(i * da)));
+    }
+
+    auto c9 = sketch.new_element<CircularArc>(v4,v8075[0],v0,r);
+    auto c10 = sketch.new_element<CircularArc>(v8075.back(),v8,v0,r);
+
+    for (size_t i = 0; i != v8075.size() - 1; ++i) {
+        sketch.new_element<CircularArc>(v8075[i],v8075[i+1],v0,r);
+    }
+
+    EXPECT_LE(sketch.solve(), 132.0 * FLT_EPSILON);
+
+    sketch.build();
+
+    sketch.save_as<SaveMethod::Rasterize>(SAVE_DIR, "three_layer_airgap");
+
+    Mesh mesh{sketch};
+
+    EXPECT_TRUE(mesh.edges_are_valid());
+    EXPECT_TRUE(mesh.edges_are_optimal());
+}
+
 TEST(Mesh, create_I_shaped_domain) {
     std::string test_name{"i_domain"};