dyna3/engine-proto/engine.jl

180 lines
4.7 KiB
Julia
Raw Normal View History

2024-01-24 16:16:24 +00:00
module Engine
export Construction, mprod
2024-01-24 16:16:24 +00:00
import Subscripts
2024-01-24 16:16:24 +00:00
using LinearAlgebra
using AbstractAlgebra
2024-01-24 16:16:24 +00:00
using Groebner
# --- primitve elements ---
2024-01-27 17:28:29 +00:00
abstract type Element{T} end
mutable struct Point{T} <: Element{T}
coords::Vector{MPolyRingElem{T}}
vec::Union{Vector{MPolyRingElem{T}}, Nothing}
2024-01-27 17:28:29 +00:00
rel::Nothing
## [to do] constructor argument never needed?
2024-01-27 17:28:29 +00:00
Point{T}(
coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
2024-01-27 17:28:29 +00:00
vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
) where T = new(coords, vec, nothing)
end
##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ]
function buildvec!(pt::Point)
coordring = parent(pt.coords[1])
pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
end
2024-01-27 17:28:29 +00:00
mutable struct Sphere{T} <: Element{T}
coords::Vector{MPolyRingElem{T}}
vec::Union{Vector{MPolyRingElem{T}}, Nothing}
2024-01-27 17:28:29 +00:00
rel::Union{MPolyRingElem{T}, Nothing}
2024-01-24 16:16:24 +00:00
## [to do] constructor argument never needed?
2024-01-27 17:28:29 +00:00
Sphere{T}(
coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
2024-01-27 17:28:29 +00:00
vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
rel::Union{MPolyRingElem{T}, Nothing} = nothing
) where T = new(coords, vec, rel)
2024-01-24 16:16:24 +00:00
end
##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
function buildvec!(sph::Sphere)
coordring = parent(sph.coords[1])
sph.vec = sph.coords
2024-01-27 17:28:29 +00:00
sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
end
const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
)
coordname(elem::Element, index) = coordnames[nameof(typeof(elem))][index]
function pushcoordname!(coordnamelist, indexed_elem::Tuple{Any, Element}, coordindex)
elemindex, elem = indexed_elem
name = coordname(elem, coordindex)
if !isnothing(name)
subscript = Subscripts.sub(string(elemindex))
push!(coordnamelist, Symbol(name, subscript))
end
end
function takecoord!(coordlist, indexed_elem::Tuple{Any, Element}, coordindex)
elem = indexed_elem[2]
if !isnothing(coordname(elem, coordindex))
push!(elem.coords, popfirst!(coordlist))
end
end
# --- primitive relations ---
abstract type Relation{T} end
mprod(v, w) = v[1]*w[2] + w[1]*v[2] - dot(v[3:end], w[3:end])
2024-01-27 17:28:29 +00:00
# elements: point, sphere
struct LiesOn{T} <: Relation{T}
2024-01-27 17:28:29 +00:00
elements::Vector{Element{T}}
LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
end
2024-01-27 17:28:29 +00:00
equation(rel::LiesOn) = dot(rel.elements[1].vec, rel.elements[2].vec)
# elements: sphere, sphere
struct AlignsWithBy{T} <: Relation{T}
2024-01-27 17:28:29 +00:00
elements::Vector{Element{T}}
cos_angle::T
2024-01-27 17:28:29 +00:00
LiesOn{T}(sph1::Point{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
end
2024-01-27 17:28:29 +00:00
equation(rel::AlignsWithBy) = dot(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
# --- constructions ---
mutable struct Construction{T}
2024-01-27 17:28:29 +00:00
elements::Set{Element{T}}
relations::Set{Relation{T}}
2024-01-27 17:28:29 +00:00
function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
allelements = union(elements, (rel.elements for rel in relations)...)
new{T}(allelements, relations)
end
end
2024-01-27 17:28:29 +00:00
function Base.push!(ctx::Construction{T}, elem::Element{T}) where T
push!(ctx.elements, elem)
end
2024-01-27 17:28:29 +00:00
function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
push!(ctx.relations, rel)
union!(ctx.elements, rel.elements)
end
function realize(ctx::Construction{T}) where T
# collect coordinate names
2024-01-27 17:28:29 +00:00
coordnamelist = Symbol[]
elemenum = enumerate(ctx.elements)
for coordindex in 1:5
for indexed_elem in elemenum
pushcoordname!(coordnamelist, indexed_elem, coordindex)
end
end
display(collect(elemenum))
display(coordnamelist)
println()
# construct coordinate ring
2024-01-27 17:28:29 +00:00
coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
2024-01-24 16:16:24 +00:00
# retrieve coordinates
for (_, elem) in elemenum
empty!(elem.coords)
end
for coordindex in 1:5
for indexed_elem in elemenum
takecoord!(coordqueue, indexed_elem, coordindex)
end
end
# construct coordinate vectors
2024-01-27 17:28:29 +00:00
for (_, elem) in elemenum
buildvec!(elem)
display(elem.coords)
display(elem.vec)
println()
2024-01-24 16:16:24 +00:00
end
2024-01-27 17:28:29 +00:00
# turn relations into equations
vcat(
equation.(ctx.relations),
[elem.rel for elem in ctx.elements if !isnothing(elem.rel)]
)
2024-01-24 16:16:24 +00:00
end
end
# ~~~ sandbox setup ~~~
2024-01-24 16:16:24 +00:00
a = Engine.Point{Rational{Int64}}()
s = Engine.Sphere{Rational{Int64}}()
2024-01-27 17:28:29 +00:00
a_on_s = Engine.LiesOn{Rational{Int64}}(a, s)
ctx = Engine.Construction{Rational{Int64}}(elements = Set([a]), relations= Set([a_on_s]))
eqns_a_s = Engine.realize(ctx)
b = Engine.Point{Rational{Int64}}()
b_on_s = Engine.LiesOn{Rational{Int64}}(b, s)
Engine.push!(ctx, b)
2024-01-27 17:28:29 +00:00
Engine.push!(ctx, s)
Engine.push!(ctx, b_on_s)
eqns_ab_s = Engine.realize(ctx)