Liebe Freundinnen und Freunde von Korrekturen,zunächst einmal war die Singularitätskontrolle falsch; nach der Korrektur sieht man Fälle wie den folgenden (die Koordinaten können aus dem File dreikant-outlier.txt ausgelesen werden), bei dem nach Massgabe der aktuellen Oktantenbedingungen kein Polarpunkt existiert (kein-polarpunkt.jpg).
Grüsse
Udo.
Clear[dreiBein]
(* dreiBein gives left-handed triads *)
dreiBein[x1_, x2_] := FoldList[Cross, x1, {x2, x1}]
Clear[recordOutlier]
recordOutlier[l_List] :=
Block[{fileN =
FileNameJoin[{$TemporaryDirectory, "dreikant-outlier.txt"}]},
If[FileExistsQ[fileN],
PutAppend[l, fileN], (* else *)
Put[l, fileN]
]
] /; Length[l] > 0
Clear[rightHand]
rightHand[v1_, v2_, v3_] := Block[{a = VectorAngle[Cross[v1, v2], v3]},
If[a > \[Pi]/2,
a = VectorAngle[Cross[v2, v1], v3];
If[a > \[Pi]/2,
Print["Numerical complanar. Fail."];
{v1, v2, v3},(* else *)
{v2, v1, v3}
], (* else *)
{v1, v2, v3}
]
] /; MatrixQ[{v1, v2, v3}] && Length[v1] == 3
Clear[octantControl]
octantControl::nedge = "Empty Intersection on edge `1` found.";
octantControl::nface = "Empty Intersection on face `1` found.";
octantControl::nwtf = "No Empty Intersection found.";
octantControl[pep12_Parallelepiped, pep13_Parallelepiped, (* Line[{s,
x1}] *)
pep21_Parallelepiped, pep23_Parallelepiped, (* Line[{s,x2}] *)
pep31_Parallelepiped, pep32_Parallelepiped (* Line[{s,x3}] *)] :=
Block[{bContinue = True, lfd = 0, res},
While[bContinue,
++lfd;
Which[
lfd == 1,
(* check the edges: internal error *)
res =
Check[RegionIntersection[Region[pep12],
Region[pep13]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == \!\(\*
TagBox[
StyleBox[
RowBox[{"Region", "[",
RowBox[{"EmptyRegion", "[", "3", "]"}], "]"}],
ShowSpecialCharacters->False,
ShowStringCharacters->True,
NumberMarks->True],
FullForm]\),
bContinue = False;
Message[octantControl::nedge, 1];
res = {Opacity[0.2], LightGreen, pep12, pep13}
],
lfd == 2,
res =
Check[RegionIntersection[Region[pep21],
Region[pep23]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == \!\(\*
TagBox[
StyleBox[
RowBox[{"Region", "[",
RowBox[{"EmptyRegion", "[", "3", "]"}], "]"}],
ShowSpecialCharacters->False,
ShowStringCharacters->True,
NumberMarks->True],
FullForm]\),
bContinue = False;
Message[octantControl::nedge, 2];
res = {Opacity[0.2], LightGreen, pep21, pep23}
],
lfd == 3,
res =
Check[RegionIntersection[Region[pep31],
Region[pep32]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == Region[EmptyRegion[3]],
bContinue = False;
Message[octantControl::nedge, 3];
res = {Opacity[0.2], LightGreen, pep31, pep32}
],
(* check the 4 octants to a face with its two edges *)
lfd == 4,
res =
Check[RegionIntersection[Region[pep12], Region[pep13],
Region[pep21],
Region[pep23]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == Region[EmptyRegion[3]],
bContinue = False;
Message[octantControl::nface, 1, 2];
res = {Opacity[0.2], LightGreen, pep12, pep13, pep21, pep23}
],
lfd == 5,
res =
Check[RegionIntersection[Region[pep21], Region[pep23],
Region[pep31],
Region[pep32]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == Region[EmptyRegion[3]],
bContinue = False;
Message[octantControl::nface, 2, 3];
res = {Opacity[0.2], LightGreen, pep21, pep23, pep31, pep32}
],
lfd == 6,
res =
Check[RegionIntersection[Region[pep31], Region[pep32],
Region[pep12],
Region[pep13]], {}, {BoundaryMeshRegion::bsuncl}];
If[res === {} \[Or] res == Region[EmptyRegion[3]],
bContinue = False;
Message[octantControl::nface, 3, 1];
res = {Opacity[0.2], LightGreen, pep31, pep32, pep12, pep13}
],
True,
bContinue = False;
Message[octantControl::nwtf];
res = {Opacity[0.2], LightGreen, pep12, pep13, pep21, pep23,
pep31, pep32}
]
];
res
]
Clear[dreiKant]
dreiKant[s_, k1_, k2_, k3_, bRegion_: False] :=
Block[{k1s, k2s, k3s, x1, x2, x3, n1, n2, n3, pep12, pep13, pep21,
pep23, pep31, pep32, targetR, targetM, res, t, t3, s1, s2, s3, u1,
u2, u3, 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 === {},
Print["targetR = ", FullForm[targetR]];
Print["targetM = ", FullForm[targetM]];
recordOutlier[{s, k1, k2, k3}];
(* Graphics to show *)
g = {Join[{Blue, Thick , Line[{s, x1}], Line[{s, x2}],
Line[{s, x3}]},
octantControl[pep12, pep13, pep21, pep23, pep31, pep32]
], Ticks -> Automatic, Axes -> True,
AxesLabel -> {"X", "Y", "Z"}
}, (* else *)
(* Polarpunkt p *)
p = Plus @@ ({{u1, u2, u3}, {s1, s2, s3}/2.} /. res);
(* Graphcis to show *)
g = If[bRegion,
{
{
{Green, Opacity[0.2], targetM, Yellow, Opacity[0.8],
Cuboid[{u1, u2, u3}, {u1 + s1, u2 + s2, u3 + s3}] /. res},
Polygon[{
{s, x1, x1 + n1[[3]], x2}, {s, x2, x2 + n2[[3]], x3}, {s,
x3, x3 + n3[[3]], x1}}],
{PointSize[0.03], Black, Point[s]},
{PointSize[0.03], Gray, Point[{x1, x2, x3}]},
{PointSize[0.03], Pink,
Point[{x1 + n1[[3]], x2 + n2[[3]], x3 + n3[[3]]}]}
},
Ticks -> Automatic, Axes -> True, AxesLabel -> {"X", "Y", "Z"}
}, (* else *)
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]];
bCommon =
If[Sign[((p - s).n1[[2]])/Norm[p - s]] !=
Sign[((x3 - s).n1[[2]])/Norm[x3 - s]]
||
Sign[((p - s).n2[[2]])/Norm[p - s]] !=
Sign[((x1 - s).n2[[2]])/Norm[x1 - s]]
||
Sign[((p - s).n3[[2]])/Norm[p - s]] !=
Sign[((x2 - s).n3[[2]])/Norm[x2 - s]]
(* --------- *)
||
Sign[((x3 + n3[[3]] - s).n1[[2]])/Norm[x3 + n3[[3]] - s]] !=
Sign[((x3 - s).n1[[2]])/Norm[x3 - s]]
||
Sign[((x3 + n3[[3]] - s).n2[[2]])/Norm[x3 + n3[[3]] - s]] !=
Sign[((x1 - s).n2[[2]])/Norm[x1 - s]]
(* --------- *)
||
Sign[((x1 + n1[[3]] - s).n2[[2]])/Norm[x1 + n1[[3]] - s]] !=
Sign[((x1 - s).n2[[2]])/Norm[x1 - s]]
||
Sign[((x1 + n1[[3]] - s).n3[[2]])/Norm[x1 + n1[[3]] - s]] !=
Sign[((x2 - s).n3[[2]])/Norm[x2 - s]]
(* --------- *)
||
Sign[((x2 + n2[[3]] - s).n3[[2]])/Norm[x2 + n2[[3]] - s]] !=
Sign[((x2 - s).n3[[2]])/Norm[x2 - s]]
||
Sign[((x2 + n2[[3]] - s).n1[[2]])/Norm[x2 + n2[[3]] - s]] !=
Sign[((x3 - s).n1[[2]])/Norm[x3 - s]],
recordOutlier[{s, k1, k2, k3}];
False, (* else *)
True
];
{
{
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]]}
}],
{PointSize[0.03], Black, Point[s]},
{PointSize[0.03], Gray, Point[{x1, x2, x3}]},
{PointSize[0.03], Pink,
Point[{x1 + n1[[3]], x2 + n2[[3]], x3 + n3[[3]]}]},
{PointSize[0.03], Red, Point[p]}
}, Ticks -> Automatic, Axes -> True,
AxesLabel -> {"X", "Y", "Z"},
PlotLabel -> If[bCommon, "Common", "Outlier"]
}
]
];
Graphics3D @@ g
] /; MatrixQ[{s, k1, k2, k3}, NumericQ] &&
Dimensions[{s, k1, k2, k3}][[2]] ==
3 && (Alternatives @@ Join[s, k1, k2, k3]) \[Element] Reals &&
bRegion \[Element] Booleans
kein-polarpunkt.jpg
Description: JPEG image
_______________________________________________ DMUG Deutschsprachiges Mathematica-Forum demug@XXXXXXX.ch http://www.mathematica.ch/mailman/listinfo/demug Archiv: http://www.mathematica.ch/archiv.html
