Import manifold code from fieldviz

2 files changed, 343 insertions, 0 deletions

diff --git a/src/curve.cpp b/src/curve.cpp
new file mode 100644
index 0000000..0b49bef
--- /dev/null
+++ b/src/curve.cpp
@@ -0,0 +1,192 @@
+#include <iostream>
+#include <cmath>
+#include <fml/curve.h>
+
+using namespace std;
+
+namespace fml {
+    vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+    {
+        vec3 diff = this->b - this->a, sum = 0;
+
+        scalar len = diff.magnitude();
+
+        vec3 diffnorm = diff / len, s = this->a, ds = diffnorm * dl;
+
+        scalar l;
+
+        for(l = 0; l < len; l += dl, s += ds)
+            sum += integrand(s, ds);
+
+        if(l < len)
+            sum += integrand(s, diffnorm * (len - l));
+
+        return sum;
+    }
+
+    vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+    {
+        //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+        scalar r = this->radius.magnitude(), dtheta = dl / r;
+
+        //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+        quat rot = quat::from_angleaxis(dtheta, normal);
+
+        //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+        scalar len = this->angle * r, l;
+
+        vec3 sum = 0;
+
+        quat radius = this->radius;
+
+        for(l = 0; l < len; l += dl)
+        {
+            vec3 ds = this->normal.cross(radius).normalize() * dl;
+            sum += integrand(this->center + radius, ds);
+
+            radius = radius.rotateby(rot);
+        }
+
+        if(l < len)
+        {
+            dl = len - l;
+            dtheta = dl / r;
+
+            vec3 ds = this->normal.cross(radius).normalize() * dl;
+            sum += integrand(this->center + radius, ds);
+        }
+
+        return sum;
+    }
+
+    vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+    {
+        //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+        scalar r = this->radius.magnitude(), dtheta = dl / r;
+
+        //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+        quat rot = quat::from_angleaxis(dtheta, normal);
+
+        //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+        /* only flat (doesn't include stretched length) */
+        scalar len = this->angle * r, l;
+
+        vec3 sum = 0;
+
+        quat radius = this->radius;
+
+        /* how far along the axis we've moved */
+        vec3 axis = 0;
+
+        /* we add this axial component each iteration */
+        vec3 dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+        for(l = 0; l < len; l += dl, axis += dp)
+        {
+            vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+            sum += integrand(this->origin + axis + radius, ds);
+
+            /* rotate by dtheta about normal */
+            radius = radius.rotateby(rot);
+        }
+
+        if(l < len)
+        {
+            dl = len - l;
+            dtheta = dl / r;
+            dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+            vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+            sum += integrand(this->origin + axis + radius, ds);
+        }
+
+        return sum;
+    }
+
+    vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+    {
+        /* this applies to the minor rotation */
+        scalar minor_r = this->minor_r;
+
+        scalar major_r = this->major_radius.magnitude();
+
+        /* how far parallel to the major median each dl takes us */
+        scalar pitch_length = this->pitch * major_r;
+
+        scalar dpar = this->pitch * dl / (2 * M_PI * sqrt(minor_r * minor_r + pitch_length * pitch_length)) * major_r;
+
+        /* how far each dl rotates the major vector */
+        scalar major_dtheta = dpar / major_r;
+
+        /* how long perpendicular to the major median each dl takes us (from
+         * Pythagoras) */
+        scalar dper = sqrt(dl * dl - dpar * dpar);
+
+        scalar minor_dtheta = dper / minor_r;
+
+        /* minor normal is binormal to major normal and radius */
+        vec3 minor_normal = this->major_normal.cross(this->major_radius).normalize();
+
+        /* quaternion describing rotation of minor vector */
+        quat minor_rot = quat::from_angleaxis(minor_dtheta, minor_normal);
+
+        quat major_rot = quat::from_angleaxis(major_dtheta, major_normal);
+
+        vec3 sum = 0;
+
+        quat major_radius = this->major_radius;
+        quat minor_radius = this->major_radius.normalize() * this->minor_r;
+
+        //cout << "minor_dtheta = " << minor_dtheta << endl;
+        //cout << "minor_dper = " << dper << endl;
+        //cout << "dl = " << dl << endl;
+        //cout << "dpar = " << dpar << endl;
+        //cout << "pitch = " << this->pitch << endl;
+
+        for(scalar theta = 0; theta < this->major_angle; theta += major_dtheta)
+        {
+            vec3 dpar_vec = minor_normal * dpar;
+            vec3 ds = minor_normal.cross(minor_radius).normalize() * dper + dpar_vec;
+            sum += integrand(this->origin + major_radius + minor_radius, ds);
+
+            /* we need to rotate some vectors: first the major radius and
+             * minor normal about the major normal, and then the minor
+             * radius */
+            /* we also need to rotate the minor radius rotation quaternion
+             * about the origin */
+            major_radius = major_radius.rotateby(major_rot);
+            minor_normal = minor_normal.rotateby(major_rot);
+
+            minor_radius = minor_radius.rotateby(minor_rot);
+
+            minor_rot = quat::from_angleaxis(minor_dtheta, minor_normal);
+            //minor_rot = major_rot * minor_rot;
+            //minor_crot = minor_rot.conjugate();
+
+            //cout << "Minor radius is " << minor_radius << endl;
+            //cout << "Minor rot quat is " << minor_rot << endl;
+        }
+
+#if 0
+        if(l < len)
+        {
+            dl = len - l;
+            dtheta = dl / r;
+            rot = quat::from_angleaxis(dtheta, normal);
+            crot = rot.conjugate();
+            dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+            vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+            sum += integrand(this->origin + axis + radius, ds);
+        }
+#endif
+
+        return sum;
+    }
+}
diff --git a/src/surface.cpp b/src/surface.cpp
new file mode 100644
index 0000000..e3fb061
--- /dev/null
+++ b/src/surface.cpp
@@ -0,0 +1,151 @@
+#include <iostream>
+#include <cmath>
+#include <fml/surface.h>
+
+namespace fml {
+    vec3 Plane::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+    {
+        vec3 sum = 0;
+
+        vec3 dA = (d * this->v1).cross(d * this->v2);
+
+        for(scalar s = 0; s < 1; s += d)
+            for(scalar t = 0; t < 1; t += d)
+            {
+                vec3 p = this->p0 + s * this->v1 + t * this->v2;
+                sum += integrand(p, dA);
+            }
+        return sum;
+    }
+
+    vec3 Disk::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar dr) const
+    {
+        vec3 sum = 0;
+
+        scalar radius = this->radius.magnitude();
+        vec3 radnorm = this->radius.normalize();
+
+        /* chosen so that the outermost ring will consist of square area
+         * elements */
+        scalar dtheta = dr / radius;
+
+        quat rot = quat::from_angleaxis(dtheta, this->normal);
+
+        for(scalar r = 0; r < radius; r += dr)
+        {
+            vec3 s = this->center + radnorm * r;
+
+            /* area element is constant for given r */
+            vec3 dA = this->normal * r * dr * dtheta;
+
+            for(scalar theta = 0; theta < this->angle; theta += dtheta)
+            {
+                sum += integrand(s, dA);
+
+                s = s.rotateby(rot);
+            }
+        }
+
+        return sum;
+    }
+
+    vec3 Sphere::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+    {
+        vec3 sum = 0;
+        /*
+         * Coordinate reference (right-handed):
+         *
+         *  ^ z
+         *  |
+         *  |    ^ y
+         *  |   /
+         *  |  /
+         *  | /
+         *  |/
+         *  O---------> x
+         *
+         */
+
+        /* we will rotate this unit vector clockwise in the X-Z plane by d
+         * each outer loop (scale and offset later) */
+        vec3 rad = vec3(0, 0, 1.0);
+
+        scalar r_sq = this->radius * this->radius;
+
+        quat roty = quat::from_angleaxis(d, vec3(0, 1, 0));
+        quat rotz = quat::from_angleaxis(d, vec3(0, 0, 1));
+
+        for(scalar phi = 0; phi < M_PI; phi += d)
+        {
+            /* operate on a copy to avoid accumulating roundoff error */
+            vec3 rad2 = rad;
+
+            /* from Jacobian: dA = r^2 * sin(phi) * dphi * dtheta */
+            scalar dA = r_sq * sin(phi) * d * d;
+
+            for(scalar theta = 0; theta < 2*M_PI; theta += d)
+            {
+                sum += integrand(this->center + this->radius * rad2, dA * rad2);
+
+                /* rotate rad2 around the z axis */
+                rad2 = rad2.rotateby(rotz);
+            }
+
+            /* rotate radius clockwise around y */
+            rad = rad.rotateby(roty);
+        }
+
+        return sum;
+    }
+
+    vec3 OpenCylinder::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+    {
+        /*
+         * We loop along the axis length, rotating a vector around it as
+         * we go.
+         *
+         * Offset is trivial.
+         *
+         *        rad/v
+         *          ^       x
+         *         .|x     . x
+         *        . | x   .   x
+         *    O-------x-------------> axis
+         *    x   .   x   .   x   .
+         *     x .     x .     x .
+         *      x       x       x
+         *
+         */
+
+        /* TODO: normalize d so all integrals run in 1/d^n time? */
+
+        vec3 v = vec3::any_unit_normal(this->axis);
+        quat rot = quat::from_angleaxis(d, this->axis);
+
+        scalar axis_len = this->axis.magnitude();
+        vec3 norm_axis = this->axis.normalize();
+
+        vec3 sum = 0;
+
+        for(scalar l = 0; l < axis_len; l += d)
+        {
+            vec3 rad = v;
+            for(scalar theta = 0; theta < 2*M_PI; theta += d)
+            {
+                sum += integrand(this->origin + l * norm_axis + this->radius * rad, rad);
+                rad = rad.rotateby(rot);
+            }
+        }
+
+        return sum;
+    }
+
+    vec3 ClosedCylinder::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+    {
+        vec3 sum = OpenCylinder::integrate(integrand, d);
+
+        sum += cap1.integrate(integrand, d) + cap2.integrate(integrand, d);
+
+        return sum;
+    }
+}