Liebe Freundinnen und Freunde des Neuen Jahres,
natürlich kann, sollte und - geradezu - muss der dreiKantSpatial[] so
verfasst werden, dass {0,0,0} kein besonderer Punkt mehr ist.
Grüsse
Udo.
Clear[dreiKantRaw]
dreiKantRaw[s_, k1_, k2_, k3_] :=
Block[{k1s, k2s, k3s, x1, x2, x3, n1, n2, n3, pep12, pep13, pep21,
pep23, pep31, pep32, targetR, targetM, t, t3, s1, s2, s3, u1, u2,
u3, res, p, q, g, bCommon},
If[MatrixRank[Subtract[#, s] & /@ {k1, k2, k3}] != 3,
Print["Singular. Bye."];
Return[$Failed]
];
{k1s, k2s, k3s} = Normalize /@ rightHand[k1 - s, k2 - s, k3 - s];
(* Die Spitze des Dreikants sei s, der Spitzenpunkt *)
{x1, x2, x3} = N[s + #] & /@ {k1s, k2s, k3s};
n1 = Normalize /@ dreiBein[x1 - s, x2 - s];
pep12 = Parallelepiped[s, n1];
pep13 =
Parallelepiped[s,
Times[{1, -1, 1}, Normalize /@ dreiBein[x1 - s, x3 - s]]];
n2 = Normalize /@ dreiBein[x2 - s, x3 - s];
pep23 = Parallelepiped[s, n2];
pep21 =
Parallelepiped[s,
Times[{1, -1, 1}, Normalize /@ dreiBein[x2 - s, x1 - s]]];
n3 = Normalize /@ dreiBein[x3 - s, x1 - s];
pep31 = Parallelepiped[s, n3];
pep32 =
Parallelepiped[s,
Times[{1, -1, 1}, Normalize /@ dreiBein[x3 - s, x2 - s]]];
targetR =
RegionIntersection[Region[pep12], Region[pep13], Region[pep21],
Region[pep23], Region[pep31], Region[pep32]];
targetM = If[ (* Is the condition too strong? *)
Head[targetR[[1]]] === EmptyRegion \[Or]
Head[targetR[[1]]] === BooleanRegion,
{}, (* else *)
If[Head[targetR[[1]]] === Parallelepiped,
Check[
ConvexHullMesh[
Plus[targetR[[1, 1]], #] & /@
Join[{{0, 0, 0}}, targetR[[1, 2]]]],
{}, {BoundaryMeshRegion::bsuncl}
], (* else *)
Check[
ConvexHullMesh[targetR[[1, 1]]],
{}, {BoundaryMeshRegion::bsuncl}
]
]
];
res = If[targetM === {},
{}, (* else *)
Check[
ConicOptimization[-t,
{VectorGreaterEqual[{{s1, s2, t3}, 0}, {"PowerCone", 1/2}],
VectorGreaterEqual[{{t3, s3, t}, 0}, {"PowerCone", 2/3}],
t3 >= 0, s1 >= 0, s2 >= 0, s3 >= 0,
Map[({u1, u2, u3} + # {s1, s2, s3} \[Element] targetM) &,
Tuples[{0, 1}, 3]]},
{t, t3, s1, s2, s3, u1, u2, u3}],
{}, {ConicOptimization::tcnstr}
]
];
If[res === {},
recordOutlier[{s, k1, k2, k3}];
{Null} ,(* else *)
p = Plus @@ ({{u1, u2, u3}, {s1, s2, s3}/2.} /. res);
q = LinearSolve[Transpose[n1], p - s];
{n1[[1]], n1[[3]]} = q[[1 ;; 3 ;; 2]] {n1[[1]], n1[[3]]};
x1 = s + n1[[1]];
q = LinearSolve[Transpose[n2], p - s];
{n2[[1]], n2[[3]]} = q[[1 ;; 3 ;; 2]] {n2[[1]], n2[[3]]};
x2 = s + n2[[1]];
q = LinearSolve[Transpose[n3], p - s];
{n3[[1]], n3[[3]]} = q[[1 ;; 3 ;; 2]] {n3[[1]], n3[[3]]};
x3 = s + n3[[1]];
{
{
Opacity[1/E],
Polygon[{
{s, x1, x1 + n1[[3]], x2}, {s, x2, x2 + n2[[3]], x3}, {s, x3,
x3 + n3[[3]], x1},
{p, x3 + n3[[3]], x1, x1 + n1[[3]]}, {p, x1 + n1[[3]], x2,
x2 + n2[[3]]}, {p, x2 + n2[[3]], x3, x3 + n3[[3]]}
}]
},
{AbsolutePointSize[9], Black, Point[s]},
{AbsolutePointSize[9], Gray, Point[{x1, x2, x3}]},
{AbsolutePointSize[9], Pink,
Point[{x1 + n1[[3]], x2 + n2[[3]], x3 + n3[[3]]}]},
{AbsolutePointSize[9], Red, Point[p]}
}
]
] /; MatrixQ[{s, k1, k2, k3}, NumericQ] &&
Dimensions[{s, k1, k2, k3}][[2]] ==
3 && (Alternatives @@ Join[s, k1, k2, k3]) \[Element] Reals
Clear[dreiKantSpatial]
(* ---- dreiKantSpatial: background free, i.e. {0,0,0} is not a \
special point ------------ *)
dreiKantSpatial[s_, k1_, k2_, k3_] :=
Graphics3D[
Join[
dreiKantRaw[s, k1, k2, k3],
dreiKantRaw[s, s - (k1 - s), k2, k3],
dreiKantRaw[s, k1, s - (k2 - s), k3],
dreiKantRaw[s, k1, k2, s - (k3 - s)],
dreiKantRaw[s, k1, s - (k2 - s), s - (k3 - s)],
dreiKantRaw[s, s - (k1 - s), k2, s - (k3 - s)],
dreiKantRaw[s, s - (k1 - s), s - (k2 - s), k3],
dreiKantRaw[s, s - (k1 - s), s - (k2 - s), s - (k3 - s)]
], Boxed -> False
] /; MatrixQ[{s, k1, k2, k3}, NumericQ] &&
Dimensions[{s, k1, k2, k3}][[2]] ==
3 && (Alternatives @@ Join[s, k1, k2, k3]) \[Element] Reals
_______________________________________________
DMUG Deutschsprachiges Mathematica-Forum demug@XXXXXXX.ch
http://www.mathematica.ch/mailman/listinfo/demug
Archiv: http://www.mathematica.ch/archiv.html
DMUG-Archiv, http://www.mathematica.ch/archiv.html