diff options
-rw-r--r-- | CMakeLists.txt | 2 | ||||
-rw-r--r-- | src/curve.cpp | 190 | ||||
-rw-r--r-- | src/curve.h | 111 | ||||
-rw-r--r-- | src/main.cpp | 2 | ||||
-rw-r--r-- | src/manifold.h | 17 | ||||
-rw-r--r-- | src/surface.cpp | 149 | ||||
-rw-r--r-- | src/surface.h | 121 |
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 |