aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt2
-rw-r--r--src/curve.cpp190
-rw-r--r--src/curve.h111
-rw-r--r--src/main.cpp2
-rw-r--r--src/manifold.h17
-rw-r--r--src/surface.cpp149
-rw-r--r--src/surface.h121
7 files changed, 1 insertions, 591 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a99d693..ba50e10 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
cmake_minimum_required (VERSION 2.6)
project (fieldviz)
-add_executable(fieldviz src/main.cpp src/curve.cpp src/surface.cpp)
+add_executable(fieldviz src/main.cpp)
add_definitions(-std=c++14 -O2 -g)
diff --git a/src/curve.cpp b/src/curve.cpp
deleted file mode 100644
index 83bec06..0000000
--- a/src/curve.cpp
+++ /dev/null
@@ -1,190 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include "curve.h"
-
-using namespace std;
-
-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/curve.h b/src/curve.h
deleted file mode 100644
index a3d743a..0000000
--- a/src/curve.h
+++ /dev/null
@@ -1,111 +0,0 @@
-#ifndef CURVE_H
-#define CURVE_H
-
-#include <cmath>
-#include <iostream>
-#include <fml/fml.h>
-
-#include "manifold.h"
-
-using 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/src/main.cpp b/src/main.cpp
index f56209f..17eb0e5 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -13,8 +13,6 @@
#include "gnuplot_i.hpp"
-#include "curve.h"
-#include "surface.h"
#include <fml/fml.h>
// under $HOME
diff --git a/src/manifold.h b/src/manifold.h
deleted file mode 100644
index f7fad5b..0000000
--- a/src/manifold.h
+++ /dev/null
@@ -1,17 +0,0 @@
-#ifndef MANIFOLD_H
-#define MANIFOLD_H
-
-#include <cmath>
-#include <iostream>
-#include <fml/fml.h>
-using 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/src/surface.cpp b/src/surface.cpp
deleted file mode 100644
index 8ab52d4..0000000
--- a/src/surface.cpp
+++ /dev/null
@@ -1,149 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include "surface.h"
-
-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;
-}
diff --git a/src/surface.h b/src/surface.h
deleted file mode 100644
index 4726521..0000000
--- a/src/surface.h
+++ /dev/null
@@ -1,121 +0,0 @@
-#ifndef SURFACE_H
-#define SURFACE_H
-
-#include <fml/fml.h>
-#include "manifold.h"
-
-using 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