aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CMakeLists.txt2
-rw-r--r--include/fml/curve.h113
-rw-r--r--include/fml/fml.h9
-rw-r--r--include/fml/manifold.h20
-rw-r--r--include/fml/surface.h121
-rw-r--r--src/curve.cpp192
-rw-r--r--src/surface.cpp151
7 files changed, 608 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 53de8b9..86f33d4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,9 @@ cmake_minimum_required(VERSION 3.9)
project(libfml VERSION 1.0 DESCRIPTION "Frank's Math Library")
add_library(fml SHARED
+ src/curve.cpp
src/quat.cpp
+ src/surface.cpp
src/vec2.cpp
src/vec3.cpp)
diff --git a/include/fml/curve.h b/include/fml/curve.h
new file mode 100644
index 0000000..8a705cd
--- /dev/null
+++ b/include/fml/curve.h
@@ -0,0 +1,113 @@
+#ifndef CURVE_H
+#define CURVE_H
+
+#include <cmath>
+#include <iostream>
+
+#include "fml.h"
+
+#include "manifold.h"
+
+namespace fml {
+
+/* All curves inherit this class (which is empty because Manifold has
+ * everything we need). */
+ class Curve : public Manifold {
+ public:
+ const int dimension() const { return 1; }
+ };
+
+ class LineSegment : public Curve {
+ private:
+ vec3 a, b;
+ public:
+ LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+
+ const char *name() const { return "LineSegment"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+ };
+
+ std::istream operator>>(std::istream &is, LineSegment &ls);
+
+ class Arc : public Curve {
+ private:
+ vec3 center;
+
+ /* these are relative to the center (direction will be determined
+ * by RHR of normal), and should be orthonormal */
+ vec3 radius, normal;
+
+ /* how many radians the arc extends for (can be greater than 2pi) */
+ scalar angle;
+ public:
+ Arc(vec3 c_, vec3 r_, vec3 n_, scalar th) : center(c_), radius(r_), normal(n_), angle(th) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+
+ const char *name() const { return "Arc"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+ };
+
+ std::istream operator>>(std::istream &is, LineSegment &ls);
+
+ class Spiral : public Curve {
+ private:
+ vec3 origin;
+
+ /* these are relative to the center (direction will be determined
+ * by RHR of normal), and should be orthonormal */
+ vec3 radius, normal;
+
+ /* how many radians the arc extends for (can be greater than 2pi) */
+ scalar angle;
+
+ /* linear distance between turns (2pi) */
+ scalar pitch;
+ public:
+ Spiral(vec3 c_, vec3 r_, vec3 n_, scalar th, scalar p) : origin(c_), radius(r_), normal(n_), angle(th), pitch(p) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+
+ const char *name() const { return "Solenoid"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+ };
+
+ std::istream operator>>(std::istream &is, LineSegment &ls);
+
+ class Toroid : public Curve {
+ private:
+ vec3 origin;
+
+ /* these are relative to the center (direction will be determined
+ * by RHR of normal), and should be orthonormal */
+ vec3 major_radius, major_normal;
+
+ /* "thickness" of toroid */
+ scalar minor_r;
+
+ /* how many radians (about the center) the toroid extends for
+ * (can be greater than 2pi) */
+ scalar major_angle;
+
+ /* central angle between successive turns (2pi rotation of small
+ * radius vector) */
+ scalar pitch;
+ public:
+ Toroid() {};
+ Toroid(vec3 o, vec3 maj_r, vec3 maj_n, scalar ang, scalar min_r, scalar p) : origin(o), major_radius(maj_r), major_normal(maj_n), major_angle(ang), minor_r(min_r), pitch(p) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+
+ const char *name() const { return "Toroid"; };
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+ };
+
+ std::istream operator>>(std::istream &is, LineSegment &ls);
+}
+#endif
diff --git a/include/fml/fml.h b/include/fml/fml.h
index 046db5a..bef1851 100644
--- a/include/fml/fml.h
+++ b/include/fml/fml.h
@@ -4,7 +4,16 @@
* Copyright (C) 2019 Franklin Wei
*/
+#ifndef FML_H
+#define FML_H
+
typedef float scalar;
+#include "curve.h"
+#include "manifold.h"
#include "quat.h"
+#include "surface.h"
+#include "vec2.h"
#include "vec3.h"
+
+#endif
diff --git a/include/fml/manifold.h b/include/fml/manifold.h
new file mode 100644
index 0000000..4f3664f
--- /dev/null
+++ b/include/fml/manifold.h
@@ -0,0 +1,20 @@
+#ifndef MANIFOLD_H
+#define MANIFOLD_H
+
+#include <cmath>
+#include <iostream>
+
+#include "fml.h"
+#include "vec3.h"
+
+namespace fml {
+ /* All manifolds inherit this class */
+ class Manifold {
+ public:
+ virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const = 0;
+ virtual const char *name() const = 0;
+ virtual const int dimension() const = 0; // 0 = point, 1 = curve, 2 = surface, 3 = solid
+ };
+};
+
+#endif
diff --git a/include/fml/surface.h b/include/fml/surface.h
new file mode 100644
index 0000000..e42efe0
--- /dev/null
+++ b/include/fml/surface.h
@@ -0,0 +1,121 @@
+#ifndef SURFACE_H
+#define SURFACE_H
+
+#include "fml.h"
+#include "manifold.h"
+
+namespace fml {
+ /* All surfaces inherit this class */
+ /* The exact meaning of d is surface-dependent (see class
+ * definitions below); the limit of the calculated integral must
+ * approach its true value as d->0. d must be > 0. */
+ class Surface : public Manifold {
+ public:
+ const int dimension() const { return 2; }
+ };
+
+ class Plane : public Surface {
+ private:
+ /* The surface is specified by all points p = p0 + s v1 + t v2,
+ * such that 0 <= {s, t} < 1
+ *
+ * v1 and v2 must NOT be parallel.
+ *
+ * v1 and v2 should (but do not have to be) be normal (this is
+ * motivated primarily by usability; if the two vectors are indeed
+ * normal, then m1 and m2 have the nice geometric meaning of side
+ * length).
+ *
+ * d = ds = dt.
+ * dA will be in the direction of v1 x v2.
+ */
+ vec3 p0, v1, v2;
+
+ public:
+ Plane(vec3 p, vec3 _v1, vec3 _v2) : p0(p), v1(_v1), v2(_v2) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Plane"; }
+ };
+
+/* Flat, circular disk (of varying extent) */
+ class Disk : public Surface {
+ private:
+ /* This represents a (possibly incomplete) circular disk in space,
+ * with a center, radius, and normal vector as specified.
+ *
+ * `angle' specifies the extent of the disk; angle=2pi for a full
+ * circle.
+ *
+ * We use:
+ * d = dr.
+ * d_theta = d / r.
+ *
+ * This makes it so that differential area elements on the edge of
+ * the disk are square.
+ *
+ * dA will be in the direction of normal.
+ */
+ vec3 center, radius, normal;
+ scalar angle;
+
+ public:
+ Disk(vec3 c, vec3 r, vec3 n, scalar a) : center(c), radius(r), normal(n), angle(a) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Disk"; }
+ };
+
+/* Hollow, spherical shell */
+ class Sphere : public Surface {
+ private:
+ /* d = dtheta = dphi */
+
+ vec3 center;
+ scalar radius;
+
+ public:
+ Sphere(vec3 c, scalar r) : center(c), radius(r) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Sphere"; }
+ };
+
+/* Cylinder without end caps */
+ class OpenCylinder : public Surface {
+ private:
+ /* Like this:
+ *
+ * ___________________
+ * / \ \
+ * origin ---------------> axis
+ * \_/_________________/
+ *
+ */
+
+ vec3 origin, axis;
+ scalar radius;
+
+ public:
+ OpenCylinder(vec3 o, vec3 a, scalar r) : origin(o), axis(a), radius(r) {}
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "OpenCylinder"; }
+ };
+
+/* Capped cylinder */
+ class ClosedCylinder : public OpenCylinder {
+ private:
+ Disk cap1, cap2;
+
+ public:
+ ClosedCylinder(vec3 o, vec3 a, scalar r) : OpenCylinder(o, a, r),
+ cap1(o, vec3::any_unit_normal(a), -a.normalize(), 2*M_PI),
+ cap2(o + a, vec3::any_unit_normal(a), a.normalize(), 2*M_PI) {}
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "ClosedCylinder"; }
+ };
+}
+
+#endif
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;
+ }
+}
generated by cgit v1.1 at 2026-04-24 19:44:31 +0000