From 3e49bb14e5e0167d6d33ced18bfa3d904b771b48 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 11 Oct 2010 16:24:07 +0100
Subject: [PATCH] BUG: plane : fix 3-plane intersection

---
 .../meshes/primitiveShapes/plane/plane.C      | 31 +++++++++----------
 .../meshes/primitiveShapes/plane/plane.H      |  2 +-
 2 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
index d4407a7d68..8ced550982 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
@@ -272,10 +272,10 @@ const Foam::point& Foam::plane::refPoint() const
 }
 
 
-// Return coefficcients for plane equation: ax + by + cz + d = 0
-Foam::scalarList Foam::plane::planeCoeffs() const
+// Return coefficients for plane equation: ax + by + cz + d = 0
+Foam::FixedList<Foam::scalar, 4> Foam::plane::planeCoeffs() const
 {
-    scalarList C(4);
+    FixedList<scalar, 4> C(4);
 
     scalar magX = mag(unitVector_.x());
     scalar magY = mag(unitVector_.y());
@@ -291,8 +291,8 @@ Foam::scalarList Foam::plane::planeCoeffs() const
         }
         else
         {
-            C[0] = 0;
-            C[1] = 0;
+            C[0] = unitVector_.x()/unitVector_.z();
+            C[1] = unitVector_.y()/unitVector_.z();
             C[2] = 1;
         }
     }
@@ -300,14 +300,14 @@ Foam::scalarList Foam::plane::planeCoeffs() const
     {
         if (magY > magZ)
         {
-            C[0] = 0;
+            C[0] = unitVector_.x()/unitVector_.y();
             C[1] = 1;
             C[2] = unitVector_.z()/unitVector_.y();
         }
         else
         {
-            C[0] = 0;
-            C[1] = 0;
+            C[0] = unitVector_.x()/unitVector_.z();
+            C[1] = unitVector_.y()/unitVector_.z();
             C[2] = 1;
         }
     }
@@ -422,19 +422,18 @@ Foam::point Foam::plane::planePlaneIntersect
     const plane& plane3
 ) const
 {
-    List<scalarList> pcs(3);
-    pcs[0]= planeCoeffs();
-    pcs[1]= plane2.planeCoeffs();
-    pcs[2]= plane3.planeCoeffs();
+    FixedList<scalar, 4> coeffs1(planeCoeffs());
+    FixedList<scalar, 4> coeffs2(plane2.planeCoeffs());
+    FixedList<scalar, 4> coeffs3(plane3.planeCoeffs());
 
     tensor a
     (
-        pcs[0][0],pcs[0][1],pcs[0][2],
-        pcs[1][0],pcs[1][1],pcs[1][2],
-        pcs[2][0],pcs[2][1],pcs[2][2]
+        coeffs1[0],coeffs1[1],coeffs1[2],
+        coeffs2[0],coeffs2[1],coeffs2[2],
+        coeffs3[0],coeffs3[1],coeffs3[2]
     );
 
-    vector b(pcs[0][3],pcs[1][3],pcs[2][3]);
+    vector b(coeffs1[3],coeffs2[3],coeffs3[3]);
 
     return (inv(a) & (-b));
 }
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
index 44eb20cbc0..89fb6afd28 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
@@ -149,7 +149,7 @@ public:
 
         //- Return coefficients for the
         //  plane equation: ax + by + cz + d = 0
-        scalarList planeCoeffs() const;
+        FixedList<scalar, 4> planeCoeffs() const;
 
         //- Return nearest point in the plane for the given point
         point nearestPoint(const point& p) const;