From 944eb1164e4a1d3151b618b2fc4d7e0a8a761960 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 18 Jul 2013 13:06:36 +0100
Subject: [PATCH] ENH: pointConstraint: calculate unconstrained directions

---
 .../pointConstraint/pointConstraint.H         |  4 ++
 .../pointConstraint/pointConstraintI.H        | 43 +++++++++++++++++++
 2 files changed, 47 insertions(+)

diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H
index 5b48463befa..04de1932202 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H
@@ -83,6 +83,10 @@ public:
 
         //- Return the accumulated constraint transformation tensor
         inline tensor constraintTransformation() const;
+
+        //- Return the accumulated unconstrained directions. Directions
+        //  coded as first n rows of tensor.
+        inline void unconstrainedDirections(label& n, tensor& vecs) const;
 };
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H
index 0915a301d4d..82c1957c3b2 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H
@@ -136,6 +136,49 @@ Foam::tensor Foam::pointConstraint::constraintTransformation() const
 }
 
 
+void Foam::pointConstraint::unconstrainedDirections(label& n, tensor& tt)
+const
+{
+    n = 3-first();
+
+    FixedList<vector, 3> vecs;
+
+    if (first() == 0)
+    {
+        vecs[0] = vector(1, 0, 0);
+        vecs[1] = vector(0, 1, 0);
+        vecs[2] = vector(0, 0, 1);
+    }
+    else if (first() == 1)
+    {
+        const vector& planeDir = second();
+
+        vecs[0] = vector(1, 0, 0) - planeDir.x()*planeDir;
+
+        if (mag(vecs[0].x()) < 1e-3)
+        {
+            vecs[0] = vector(0, 1, 0) - planeDir.y()*planeDir;
+        }
+
+        vecs[0] /= mag(vecs[0]);
+        vecs[1] = vecs[0] ^ planeDir;
+        vecs[1] /= mag(vecs[1]);
+    }
+    else if (first() == 2)
+    {
+        vecs[0] = second();
+    }
+
+    // Knock out remaining vectors
+    for (direction dir = n; dir < vecs.size(); dir++)
+    {
+        vecs[dir] = vector::zero;
+    }
+
+    tt = tensor(vecs[0], vecs[1], vecs[2]);
+}
+
+
 void Foam::combineConstraintsEqOp::operator()
 (
     pointConstraint& x,
-- 
GitLab