me/Tetrahedrons
Archived
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

a bit faster multivector calculations

Browse Source
This commit is contained in:
David Allemang
2018-04-01 20:52:57 -04:00
parent 9cb8b1ab09
commit a823b33f0b
2 changed files with 138 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

164
Tetrahedrons/MVec3d.cs
View File

@@ -44,45 +44,14 @@ namespace Tetrahedrons
{
}
public static MVec3d Scalar(double e)
{
return new MVec3d(e, 0, 0, 0, 0, 0, 0, 0);
}
public static MVec3d Vector(double e1, double e2, double e3)
{
return new MVec3d(0, e1, e2, e3, 0, 0, 0, 0);
}
public static MVec3d Vector(Vector3d v)
{
return new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0);
}
public static MVec3d Bivector(double e23, double e31, double e12)
{
return new MVec3d(0, 0, 0, 0, e23, e31, e12, 0);
}
public static MVec3d Trivector(double e123)
{
return new MVec3d(0, 0, 0, 0, 0, 0, 0, e123);
}
public static MVec3d Complex(double real, double imag)
{
return new MVec3d(real, 0, 0, 0, 0, 0, 0, imag);
}
public static MVec3d Complex(double w, double i, double j, double k)
{
return new MVec3d(w, 0, 0, 0, i, j, k, 0);
}
public static MVec3d Complex(double angle, MVec3d plane)
{
return Math.Cos(angle) + Math.Sin(angle) * plane;
}
public static MVec3d Scalar(double e) => new MVec3d(e, 0, 0, 0, 0, 0, 0, 0);
public static MVec3d Vector(double e1, double e2, double e3) => new MVec3d(0, e1, e2, e3, 0, 0, 0, 0);
public static MVec3d Vector(Vector3d v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0);
public static MVec3d Bivector(double e23, double e31, double e12) => new MVec3d(0, 0, 0, 0, e23, e31, e12, 0);
public static MVec3d Trivector(double e123) => new MVec3d(0, 0, 0, 0, 0, 0, 0, e123);
public static MVec3d Complex(double real, double imag) => new MVec3d(real, 0, 0, 0, 0, 0, 0, imag);
public static MVec3d Complex(double w, double i, double j, double k) => new MVec3d(w, 0, 0, 0, i, j, k, 0);
public static MVec3d Complex(double angle, MVec3d plane) => Math.Cos(angle) + Math.Sin(angle) * plane;
public static MVec3d Add(MVec3d m, MVec3d n)
{
@@ -102,16 +71,105 @@ namespace Tetrahedrons
public static MVec3d Mult(MVec3d m, MVec3d n)
{
return new MVec3d(
m.E * n.E + m.E1 * n.E1 + m.E2 * n.E2 + m.E3 * n.E3 - m.E23 * n.E23 - m.E31 * n.E31 - m.E12 * n.E12 - m.E123 * n.E123,
m.E * n.E1 + m.E1 * n.E - m.E2 * n.E12 + m.E3 * n.E31 - m.E23 * n.E123 - m.E31 * n.E3 + m.E12 * n.E2 - m.E123 * n.E23,
m.E * n.E2 + m.E1 * n.E12 + m.E2 * n.E - m.E3 * n.E23 + m.E23 * n.E3 - m.E31 * n.E123 - m.E12 * n.E1 - m.E123 * n.E31,
m.E * n.E3 - m.E1 * n.E31 + m.E2 * n.E23 + m.E3 * n.E - m.E23 * n.E2 + m.E31 * n.E1 - m.E12 * n.E123 - m.E123 * n.E12,
m.E * n.E23 + m.E1 * n.E123 + m.E2 * n.E3 - m.E3 * n.E2 + m.E23 * n.E - m.E31 * n.E12 + m.E12 * n.E31 + m.E123 * n.E1,
m.E * n.E31 - m.E1 * n.E3 + m.E2 * n.E123 + m.E3 * n.E1 + m.E23 * n.E12 + m.E31 * n.E - m.E12 * n.E23 + m.E123 * n.E2,
m.E * n.E12 + m.E1 * n.E2 - m.E2 * n.E1 + m.E3 * n.E123 - m.E23 * n.E31 + m.E31 * n.E23 + m.E12 * n.E + m.E123 * n.E3,
m.E * n.E123 + m.E1 * n.E23 + m.E2 * n.E31 + m.E3 * n.E12 + m.E23 * n.E1 + m.E31 * n.E2 + m.E12 * n.E3 + m.E123 * n.E
);
var r = new MVec3d();
if (Math.Abs(m.E) > 1e-8)
{
r.E += m.E * n.E;
r.E1 += m.E * n.E1;
r.E2 += m.E * n.E2;
r.E3 += m.E * n.E3;
r.E23 += m.E * n.E23;
r.E31 += m.E * n.E31;
r.E12 += m.E * n.E12;
r.E123 += m.E * n.E123;
}
if (Math.Abs(m.E1) > 1e-8)
{
r.E += +m.E1 * n.E1;
r.E1 += +m.E1 * n.E;
r.E2 += +m.E1 * n.E12;
r.E3 += -m.E1 * n.E31;
r.E23 += +m.E1 * n.E123;
r.E31 += -m.E1 * n.E3;
r.E12 += +m.E1 * n.E2;
r.E123 += +m.E1 * n.E23;
}
if (Math.Abs(m.E2) > 1e-8)
{
r.E += +m.E2 * n.E2;
r.E1 += -m.E2 * n.E12;
r.E2 += +m.E2 * n.E;
r.E3 += +m.E2 * n.E23;
r.E23 += +m.E2 * n.E3;
r.E31 += +m.E2 * n.E123;
r.E12 += -m.E2 * n.E1;
r.E123 += +m.E2 * n.E31;
}
if (Math.Abs(m.E3) > 1e-8)
{
r.E += +m.E3 * n.E3;
r.E1 += +m.E3 * n.E31;
r.E2 += -m.E3 * n.E23;
r.E3 += +m.E3 * n.E;
r.E23 += -m.E3 * n.E2;
r.E31 += +m.E3 * n.E1;
r.E12 += +m.E3 * n.E123;
r.E123 += +m.E3 * n.E12;
}
if (Math.Abs(m.E23) > 1e-8)
{
r.E += -m.E23 * n.E23;
r.E1 += -m.E23 * n.E123;
r.E2 += +m.E23 * n.E3;
r.E3 += -m.E23 * n.E2;
r.E23 += +m.E23 * n.E;
r.E31 += +m.E23 * n.E12;
r.E12 += -m.E23 * n.E31;
r.E123 += +m.E23 * n.E1;
}
if (Math.Abs(m.E31) > 1e-8)
{
r.E += -m.E31 * n.E31;
r.E1 += -m.E31 * n.E3;
r.E2 += -m.E31 * n.E123;
r.E3 += +m.E31 * n.E1;
r.E23 += -m.E31 * n.E12;
r.E31 += +m.E31 * n.E;
r.E12 += +m.E31 * n.E23;
r.E123 += +m.E31 * n.E2;
}
if (Math.Abs(m.E12) > 1e-8)
{
r.E += -m.E12 * n.E12;
r.E1 += +m.E12 * n.E2;
r.E2 += -m.E12 * n.E1;
r.E3 += -m.E12 * n.E123;
r.E23 += +m.E12 * n.E31;
r.E31 += -m.E12 * n.E23;
r.E12 += +m.E12 * n.E;
r.E123 += +m.E12 * n.E3;
}
if (Math.Abs(m.E123) > 1e-8)
{
r.E += -m.E123 * n.E123;
r.E1 += -m.E123 * n.E23;
r.E2 += -m.E123 * n.E31;
r.E3 += -m.E123 * n.E12;
r.E23 += +m.E123 * n.E1;
r.E31 += +m.E123 * n.E2;
r.E12 += +m.E123 * n.E3;
r.E123 += +m.E123 * n.E;
}
return r;
}
public double Norm2 => E * E + E1 * E1 + E2 * E2 + E3 * E3 + E23 * E23 + E31 * E31 + E12 * E12 + E123 * E123;
@@ -193,12 +251,18 @@ namespace Tetrahedrons
public static MVec3d operator &(MVec3d m, MVec3d n) => Inner(m, n);
public static MVec3d operator ^(MVec3d m, MVec3d n) => Outer(m, n);
public Vector3d VectorPart => new Vector3d(E1, E2, E3);
public MVec3d ScalarPart => Scalar(E);
public MVec3d VectorPart => Vector(E1, E2, E3);
public MVec3d BivectorPart => Bivector(E23, E31, E12);
public MVec3d TrivectorPart => Trivector(E123);
public static implicit operator MVec3d(Vector3 v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0);
public static implicit operator MVec3d(Vector3d v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0);
public static implicit operator MVec3d(double s) => new MVec3d(s, 0, 0, 0, 0, 0, 0, 0);
public static explicit operator double(MVec3d m) => m.E;
public static explicit operator Vector3d(MVec3d m) => new Vector3d(m.E1, m.E2, m.E3);
public override string ToString()
{
return $"({E:0.00} {E1:0.00}e1 {E2:0.00}e2 {E3:0.00}e3 {E23:0.00}e23 {E31:0.00}e31 {E12:0.00}e12 {E123:0.00}e123)";

57
Tetrahedrons/TetrahedronWindow.cs
View File

@@ -38,7 +38,7 @@ namespace Tetrahedrons
{
var pts = new Vector3d[p.Points.Length];
for (var i = 0; i < p.Points.Length; i++)
pts[i] = (p.Points[i] + m).VectorPart;
pts[i] = (Vector3d) (p.Points[i] + m);
return new Polygon(pts, p.Color);
}
@@ -46,7 +46,7 @@ namespace Tetrahedrons
{
var pts = new Vector3d[p.Points.Length];
for (var i = 0; i < p.Points.Length; i++)
pts[i] = (p.Points[i] - m).VectorPart;
pts[i] = (Vector3d) (p.Points[i] - m);
return new Polygon(pts, p.Color);
}
@@ -54,7 +54,7 @@ namespace Tetrahedrons
{
var pts = new Vector3d[p.Points.Length];
for (var i = 0; i < p.Points.Length; i++)
pts[i] = (p.Points[i] * m).VectorPart;
pts[i] = (Vector3d) (p.Points[i] * m);
return new Polygon(pts, p.Color);
}
@@ -62,7 +62,7 @@ namespace Tetrahedrons
{
var pts = new Vector3d[p.Points.Length];
for (var i = 0; i < p.Points.Length; i++)
pts[i] = (p.Points[i] ^ m).VectorPart;
pts[i] = (Vector3d) (p.Points[i] ^ m);
return new Polygon(pts, p.Color);
}
@@ -70,7 +70,7 @@ namespace Tetrahedrons
{
var pts = new Vector3d[p.Points.Length];
for (var i = 0; i < p.Points.Length; i++)
pts[i] = (p.Points[i] & m).VectorPart;
pts[i] = (Vector3d) (p.Points[i] & m);
return new Polygon(pts, p.Color);
}
}
@@ -78,21 +78,12 @@ namespace Tetrahedrons
public class Tetrahedron
{
public readonly Vector3d[] Points;
public Color Color;
public Color WireColor;
public Color FillColor = Color.DodgerBlue;
public Color SurfaceColor = Color.GhostWhite;
public Color WireColor = Color.Maroon;
public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4) : this(p1, p2, p3, p4, Color.GhostWhite)
public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4)
{
}
public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4, Color color) : this(p1, p2, p3, p4, color, Color.Maroon)
{
}
public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4, Color color, Color wireColor)
{
Color = color;
WireColor = wireColor;
Points = new[] {p1, p2, p3, p4};
}
@@ -115,7 +106,7 @@ namespace Tetrahedrons
else
{
GL.Begin(PrimitiveType.Triangles);
GL.Color3(Color);
GL.Color3(SurfaceColor);
for (var i = 0; i < 4; i++)
for (var j = i + 1; j < 4; j++)
@@ -130,12 +121,12 @@ namespace Tetrahedrons
}
}
public Polygon Intersection(MVec3d blade, Vector3d pivot)
public Polygon Intersection(MVec3d blade, MVec3d pivot)
{
var pts = new Vector3d[4];
for (var i = 0; i < 4; i++)
{
pts[i] = Points[i] - pivot;
pts[i] = (Vector3d) (Points[i] - pivot);
}
var a = blade * MVec3d.Unit123;
@@ -155,18 +146,18 @@ namespace Tetrahedrons
var t = ((blade ^ pts[i]) * ~(blade ^ (pts[j] - pts[i]))).E;
var point = pts[i] + t * (pts[j] - pts[i]);
vecs.Add(point + pivot);
vecs.Add((Vector3d) (point + pivot));
}
}
return new Polygon(vecs.ToArray());
return new Polygon(vecs.ToArray()) {Color = FillColor};
}
}
public class TetrahedronWindow : GameWindow
{
private Matrix4 _proj3d;
private Matrix4 _proj2d;
private Matrix4 _proj3D;
private Matrix4 _proj2D;
private Matrix4 _view;
private double _t;
@@ -208,14 +199,14 @@ namespace Tetrahedrons
if (_div)
{
GL.MatrixMode(MatrixMode.Projection);
GL.LoadMatrix(ref _proj2d);
GL.LoadMatrix(ref _proj2D);
_polyg.Draw();
}
else
{
GL.MatrixMode(MatrixMode.Projection);
GL.LoadMatrix(ref _proj3d);
GL.LoadMatrix(ref _proj3D);
GL.MatrixMode(MatrixMode.Modelview);
GL.LoadMatrix(ref _view);
@@ -227,12 +218,12 @@ namespace Tetrahedrons
GL.Begin(PrimitiveType.Points);
GL.Color3(Color.DarkBlue);
GL.Vertex3(_pivot.VectorPart);
GL.Vertex3((Vector3d) _pivot);
GL.End();
GL.Begin(PrimitiveType.Lines);
GL.Vertex3(_pivot.VectorPart);
GL.Vertex3((_pivot + (_b2 ^ _b1) * MVec3d.Unit123 * .5).VectorPart);
GL.Vertex3((Vector3d) _pivot);
GL.Vertex3((Vector3d) (_pivot + (_b2 ^ _b1) * MVec3d.Unit123 * .5));
GL.End();
GL.Enable(EnableCap.DepthTest);
@@ -248,8 +239,8 @@ namespace Tetrahedrons
{
base.OnUpdateFrame(e);
_proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2);
_proj2d = Matrix4.CreateOrthographic(4, 4f * Height / Width, -1, 1);
_proj3D = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2);
_proj2D = Matrix4.CreateOrthographic(4, 4f * Height / Width, -1, 1);
_t += e.Time;
@@ -301,7 +292,7 @@ namespace Tetrahedrons
// var blade = pln;
// var pivot = Vector3d.Zero;
_polyg = _tetra.Intersection(_b1 ^ _b2, _pivot.VectorPart);
_polyg = _tetra.Intersection(_b1 ^ _b2, _pivot);
if (_div)
{