Split Constraint.h/Constraint.cpp into multiple h/cpp files

.idea/Oersted.iml

@@ -10,28 +10,50 @@
       <sourceFolder url="file://$MODULE_DIR$/lib/GoogleTest/googletest/cmake/internal_utils.cmake" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/CMakeLists.txt" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/CMakeLists.txt" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/CircularArc.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Curve.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Sketch.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Angle.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Sketch.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Coincident.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Vertical.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Angle.h" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Vertex.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Pattern.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Contour.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Curve.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Pattern.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Radius.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Coincident.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constellation.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Radius.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Horizontal.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Vertex.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/LineSegment.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Curve.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Tangency.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Rotation.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Branch.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Branch.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Length.h" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constraint.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/CircularArc.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Vertex.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Sketch.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Star.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Symmetry.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constraint.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Length.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constellation.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/LineSegment.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Contour.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Fixation.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Vertical.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/CircularArc.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Distance.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Symmetry.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Distance.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Star.h" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/LineSegment.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constellation.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Fixation.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Sketch.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Tangency.h" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Contour.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Star.h" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Star.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Pattern.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Constraint.cpp" isTestSource="false" />
-      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Branch.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Pattern.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Rotation.h" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/Horizontal.cpp" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/src/Sketch/src/CircularArc.h" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/include/Sketch.hpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/src/Sketch/CMakeLists.txt" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/main.cpp" isTestSource="false" />

src/Sketch/CMakeLists.txt

@@ -10,13 +10,17 @@ set(SOURCE_FILES
             ./src/CircularArc.h ./src/CircularArc.cpp
 
         ./src/Constraint.h ./src/Constraint.cpp
+            ./src/Angle.h           ./src/Angle.cpp
+            ./src/Coincident.cpp    ./src/Coincident.h
+            ./src/Distance.cpp      ./src/Distance.h
+            ./src/Fixation.cpp      ./src/Fixation.h
         ./src/Pattern.h ./src/Pattern.cpp
 
         ./src/Branch.h ./src/Branch.cpp
         ./src/Star.h ./src/Star.cpp
         ./src/Constellation.h ./src/Constellation.cpp
 
-        ./src/Contour.h ./src/Contour.cpp)
+        ./src/Contour.h ./src/Contour.cpp src/Horizontal.cpp src/Horizontal.h src/Length.cpp src/Length.h src/Radius.cpp src/Radius.h src/Rotation.cpp src/Rotation.h src/Symmetry.cpp src/Symmetry.h src/Tangency.cpp src/Tangency.h src/Vertical.cpp src/Vertical.h)
 
 add_library(sketch STATIC ${SOURCE_FILES})
 

src/Sketch/include/Sketch.hpp

@@ -6,10 +6,9 @@
 #include "../src/Vertex.h"
 
 #include "../src/Curve.h"
-    #include "../src/LineSegment.h"
-    #include "../src/CircularArc.h"
 
 #include "../src/Constraint.h"
+
 #include "../src/Pattern.h"
 
 #include "../src/Branch.h"

src/Sketch/src/Angle.cpp

+#include "Sketch.hpp"
+
+void Angle::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) {
+    /*
+        The dot product of two unit vectors equals the cosine of the angle between them.
+        The cross product of two unit vectros equals the sine of the angle between them.
+        To solve the constraint, use the smaller of the dot produce and cross product equations.
+        Do this because the derivative of the dot (resp. cross) product with respect to the
+        parameters is zero when the cosine (resp. sine) is equal to one.
+    */
+
+    const double x00 = Line0->Start->x();
+    const double x01 = Line0->End->x();
+    const double x10 = Line1->Start->x();
+    const double x11 = Line1->End->x();
+
+    const double y00 = Line0->Start->y();
+    const double y01 = Line0->End->y();
+    const double y10 = Line1->Start->y();
+    const double y11 = Line1->End->y();
+
+    double vx0 = x01 - x00;
+    double vy0 = y01 - y00;
+    double vx1 = x11 - x10;
+    double vy1 = y11 - y10;
+
+    double d0 = sqrt(vx0 * vx0 + vy0 * vy0);
+    double d1 = sqrt(vx1 * vx1 + vy1 * vy1);
+
+    vx0 /= d0;
+    vy0 /= d0;
+    vx1 /= d1;
+    vy1 /= d1;
+
+    double dot = (vx0 * vx1 + vy0 * vy1);
+    double cross = (vx0 * vy1 - vy0 * vx1);
+    double scale = std::fmax(d0, d1);
+    double f;
+
+    if (abs(dot) < abs(cross)) {
+        // Solve (dot product) = cos(Dim)
+        r(EquationIndex) = scale * (dot - cos(M_PI * Dim / 180.0));
+
+        f = scale * (vx1 - dot * vx0) / d0;
+        J(EquationIndex, Line0->Start->X->get_index()) -= f;
+        J(EquationIndex, Line0->End->X->get_index()) += f;
+
+        f = scale * (vx0 - dot * vx1) / d1;
+        J(EquationIndex, Line1->Start->X->get_index()) -= f;
+        J(EquationIndex, Line1->End->X->get_index()) += f;
+
+        f = scale * (vy1 - dot * vy0) / d0;
+        J(EquationIndex, Line0->Start->Y->get_index()) -= f;
+        J(EquationIndex, Line0->End->Y->get_index()) += f;
+
+        f = scale * (vy0 - dot * vy1) / d1;
+        J(EquationIndex, Line1->Start->Y->get_index()) -= f;
+        J(EquationIndex, Line1->End->Y->get_index()) += f;
+    } else {
+        // Solve (cross product) = sin(Dim)
+        r(EquationIndex) = scale * (cross - sin(M_PI * Dim / 180.0));
+
+        f = scale * (vy1 - cross * vx0) / d0;
+        J(EquationIndex, Line0->Start->X->get_index()) -= f;
+        J(EquationIndex, Line0->End->X->get_index()) += f;
+
+        f = scale * (vy0 + cross * vx1) / d1;
+        J(EquationIndex, Line1->Start->X->get_index()) += f;
+        J(EquationIndex, Line1->End->X->get_index()) -= f;
+
+        f = scale * (vx1 + cross * vy0) / d0;
+        J(EquationIndex, Line0->Start->Y->get_index()) += f;
+        J(EquationIndex, Line0->End->Y->get_index()) -= f;
+
+        f = scale * (vx0 - cross * vy1) / d1;
+        J(EquationIndex, Line1->Start->Y->get_index()) -= f;
+        J(EquationIndex, Line1->End->Y->get_index()) += f;
+    }
+}
\ No newline at end of file

src/Sketch/src/Angle.h

+#ifndef OERSTED_ANGLE_H
+#define OERSTED_ANGLE_H
+
+class Angle : public Constraint {
+public:
+    LineSegment *Line0, *Line1;
+
+    double Dim;
+
+    // Constructors
+    Angle(LineSegment &l0, LineSegment &l1, double angle) : Line0(&l0), Line1(&l1), Dim(angle) {};
+
+    // Public Member Functions
+    size_t set_equation_index(size_t i) override {
+        EquationIndex = i;
+        return 1;
+    };
+
+    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
+};
+
+
+#endif //OERSTED_ANGLE_H

src/Sketch/src/Coincident.cpp

+#include "Sketch.hpp"
+
+template<>
+void Coincident<CircularArc>::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) {
+    const double rc = Element->radius();
+    const double xc = Element->Center->x();
+    const double yc = Element->Center->y();
+    const double xp = Point->x();
+    const double yp = Point->y();
+
+    double dx = xp - xc;
+    double dy = yp - yc;
+    double dr = sqrt(dx * dx + dy * dy);
+
+    dx /= dr;
+    dy /= dr;
+
+    r(EquationIndex) = dr - rc;
+
+    J(EquationIndex, Element->Radius->get_index()) = -dr / rc;
+
+    J(EquationIndex, Element->Center->X->get_index()) -= dx;
+    J(EquationIndex, Element->Center->Y->get_index()) -= dy;
+
+    J(EquationIndex, Point->X->get_index()) += dx;
+    J(EquationIndex, Point->Y->get_index()) += dy;
+}
+
+template
+class Coincident<CircularArc>;
+
+template<>
+void Coincident<LineSegment>::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) {
+    const double xp = Point->x();
+    const double yp = Point->y();
+
+    const double x0 = Element->Start->x();
+    const double y0 = Element->Start->y();
+    const double x1 = Element->End->x();
+    const double y1 = Element->End->y();
+
+    double dx0 = x0 - xp;
+    double dy0 = y0 - yp;
+    double dr0 = sqrt(dx0 * dx0 + dy0 * dy0);
+
+    double dx1 = x1 - xp;
+    double dy1 = y1 - yp;
+    double dr1 = sqrt(dx1 * dx1 + dy1 * dy1);
+
+    dx0 /= dr0;
+    dy0 /= dr0;
+    dx1 /= dr1;
+    dy1 /= dr1;
+
+    double cross = (dx0 * dy1 - dy0 * dx1);
+    double scale = std::fmax(dr0, dr1);
+
+    r(EquationIndex) = scale * cross;
+
+    double f0, f1;
+
+    f0 = scale * (dy1 - cross * dx0) / dr0;
+    f1 = scale * (dy0 + cross * dx1) / dr1;
+
+    J(EquationIndex, Point->X->get_index()) -= f0 - f1;
+    J(EquationIndex, Element->Start->X->get_index()) += f0;
+    J(EquationIndex, Element->End->X->get_index()) -= f1;
+
+    f0 = scale * (dx1 + cross * dy0) / dr0;
+    f1 = scale * (dx0 - cross * dy1) / dr1;
+
+    J(EquationIndex, Point->Y->get_index()) += f0 - f1;
+    J(EquationIndex, Element->Start->Y->get_index()) -= f0;
+    J(EquationIndex, Element->End->Y->get_index()) += f1;
+}
+
+template
+class Coincident<LineSegment>;
\ No newline at end of file

src/Sketch/src/Coincident.h

+#ifndef OERSTED_COINCIDENT_H
+#define OERSTED_COINCIDENT_H
+
+template<class T>
+class Coincident : public Constraint {
+public:
+    Vertex *Point;
+    T *Element;
+
+    // Constructors
+    Coincident(Vertex &p, T &e) : Point(&p), Element(&e) {};
+
+    // Public Member Functions
+    size_t set_equation_index(size_t i) override {
+        EquationIndex = i;
+        return 1;
+    };
+
+    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
+};
+
+#endif //OERSTED_COINCIDENT_H

src/Sketch/src/Constraint.h

@@ -8,192 +8,18 @@ public:
     void register_parameters(Sketch *s) override {};
 };
 
-class Angle : public Constraint {
-public:
-    LineSegment *Line0, *Line1;
-
-    double Dim;
-
-    // Constructors
-    Angle(LineSegment &l0, LineSegment &l1, double angle) : Line0(&l0), Line1(&l1), Dim(angle) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-template<class T>
-class Coincident : public Constraint {
-public:
-    Vertex *Point;
-    T *Element;
-
-    // Constructors
-    Coincident(Vertex &p, T &e) : Point(&p), Element(&e) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-template<class T>
-class Distance : public Constraint {
-public:
-    T *Element0;
-    T *Element1;
-
-    double Dim;
-
-    // Constructors
-    Distance(T &e0, T &e1, double d) : Element0(&e0), Element1(&e1), Dim(d) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override;
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Fixation : public Constraint {
-public:
-    Vertex *Point;
-    Vertex *Dim;
-
-    // Constructors
-    Fixation(Vertex &v);
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 2;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Horizontal : public Constraint {
-public:
-    LineSegment *Line;
-
-    // Constructors
-    Horizontal(LineSegment &l) : Line(&l) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Length : public Constraint {
-public:
-    LineSegment *Line;
-
-    double Dim;
+#include "Angle.h"
+#include "Coincident.h"
+#include "Distance.h"
+#include "Constraint.h"
+#include "Fixation.h"
+#include "Horizontal.h"
+#include "Length.h"
+#include "Radius.h"
+#include "Rotation.h"
+#include "Symmetry.h"
+#include "Tangency.h"
+#include "Vertical.h"
 
-    // Constructors
-    Length(LineSegment &c, double length) : Line(&c), Dim(length) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Radius : public Constraint {
-public:
-    CircularArc *Arc;
-
-    double Dim;
-
-    // Constructors
-    Radius(CircularArc &c, double r) : Arc(&c), Dim(r) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Rotation : public Constraint {
-public:
-    Vertex *V0;
-    Vertex *V1;
-    Vertex *Origin;
-    double Angle;
-
-    Rotation(Vertex &v0, Vertex &v1, Vertex &origin, double a) : V0(&v0), V1(&v1), Origin(&origin), Angle(a) {};
-
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 2;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r);
-};
-
-class Symmetry : public Constraint {
-public:
-    Vertex *V0;
-    Vertex *V1;
-    LineSegment *SymmetryLine;
-
-    Symmetry(Vertex &v0, Vertex &v1, LineSegment &line) : V0(&v0), V1(&v1), SymmetryLine(&line) {};
-
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 2;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Tangency : public Constraint {
-public:
-    CircularArc *Arc;
-    LineSegment *Line;
-
-    // Constructors
-    Tangency(CircularArc &ca, LineSegment &ls) : Arc(&ca), Line(&ls) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
-
-class Vertical : public Constraint {
-public:
-    LineSegment *Line;
-
-    // Constructors
-    Vertical(LineSegment &l) : Line(&l) {};
-
-    // Public Member Functions
-    size_t set_equation_index(size_t i) override {
-        EquationIndex = i;
-        return 1;
-    };
-
-    void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) override;
-};
 
 #endif //OERSTED_CONSTRAINT_H
\ No newline at end of file

src/Sketch/src/Curve.h

@@ -45,7 +45,7 @@ public:
     virtual bool on_segment(const Vertex *v) const final; // true if vertex is on curve segment
     virtual bool on_segment(const Vertex *v, const Vertex *origin, const double angle) const final;
 
-    // Curve-Curve Comparions
+    // Curve-Curve Comparison
     virtual bool is_identical(const Curve *c) const = 0; // true if (input curve) XOR (object curve) is a set with measure < tol
     virtual bool is_identical(const Curve *c, const Vertex *origin, const double angle) const = 0;
 
@@ -72,4 +72,7 @@ protected:
     virtual bool on_segment(const double x, const double y) const = 0;
 };
 
+#include "LineSegment.h"
+#include "CircularArc.h"
+
 #endif //OERSTED_CURVE_H
\ No newline at end of file

src/Sketch/src/Distance.cpp

+#include "Sketch.hpp"
+
+template<>
+size_t Distance<CircularArc>::set_equation_index(size_t i) {
+    EquationIndex = i;
+    return 1;
+};
+
+template<>
+void Distance<CircularArc>::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) {
+    /*
+        There are two different cases depending on whether or not the smaller circle is intended to be
+        exterior or interior to the larger circle. At present, the intent is discriminated by the position
+        of the smaller circle origin relative to the larger circle. The interior behavior is unstable in
+        the sense that for large enough distances, the smaller circle will get pushed outside the larger
+        circle and find a new equilibrium there.
+    */
+
+    const double r0 = Element0->radius();
+    const double x0 = Element0->Center->x();
+    const double y0 = Element0->Center->y();
+
+    const double r1 = Element1->radius();
+    const double x1 = Element1->Center->x();
+    const double y1 = Element1->Center->y();
+
+    double dx = x1 - x0;
+    double dy = y1 - y0;
+    double dr = sqrt(dx * dx + dy * dy);
+
+    if (dr >= std::fmax(r0, r1)) {
+        // Exterior
+        r(EquationIndex) = dr - r0 - r1 - Dim;
+        J(EquationIndex, Element0->Radius->get_index()) -= (dr - Dim) / (r0 + r1);
+        J(EquationIndex, Element1->Radius->get_index()) -= (dr - Dim) / (r0 + r1);
+    } else {
+        // Interior
+        if (r0 > r1) {
+            r(EquationIndex) = Dim - (r0 - r1 - dr);
+            J(EquationIndex, Element0->Radius->get_index()) -= (Dim + dr) / (r0 - r1);
+            J(EquationIndex, Element1->Radius->get_index()) += (Dim + dr) / (r0 - r1);
+        } else {
+            r(EquationIndex) = Dim - (r1 - r0 - dr);
+            J(EquationIndex, Element0->Radius->get_index()) += (Dim + dr) / (r1 - r0);
+            J(EquationIndex, Element1->Radius->get_index()) -= (Dim + dr) / (r1 - r0);
+        }
+    }
+
+    dx /= dr;
+    J(EquationIndex, Element0->Center->X->get_index()) -= dx;
+    J(EquationIndex, Element1->Center->X->get_index()) += dx;
+
+    dy /= dr;
+    J(EquationIndex, Element0->Center->Y->get_index()) -= dy;
+    J(EquationIndex, Element1->Center->Y->get_index()) += dy;
+}
+
+template
+class Distance<CircularArc>;
+
+template<>
+size_t Distance<LineSegment>::set_equation_index(size_t i) {
+    EquationIndex = i;
+    return 2;
+};
+
+template<>
+void Distance<LineSegment>::update(Eigen::MatrixXd &J, Eigen::VectorXd &r) {
+    /*
+        Pick one line, draw vectors from the center of the line to each of the remaining line's end points.
+        The area of the triangle formed by the vectors and line should equal half of the length of the line times
+        the distance beween the lines (one half base times height). The signed area of the triangle is equal to