Merge branch 'master' of /home/hunt2/OpenFOAM/OpenFOAM-dev

applications/solvers/incompressible/oodles/oodles.C

@@ -66,6 +66,9 @@ int main(int argc, char *argv[])
           + sgsModel->divDevBeff(U)
         );
 
+        // Optionally ensure diagonal-dominance of the momentum matrix
+        UEqn.relax();
+
         if (momentumPredictor)
         {
             solve(UEqn == -fvc::grad(p));

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C

@@ -206,7 +206,7 @@ Foam::label Foam::quadraticFitData::calcFit
     // calculate the matrix of the polynomial components
     scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
 
-    for(label ip = 0; ip < C.size(); ip++)
+        for(label ip = 0; ip < C.size(); ip++)
     {
         const point& p = C[ip];
 
@@ -228,9 +228,8 @@ Foam::label Foam::quadraticFitData::calcFit
         pz /= scale;
 
         label is = 0;
-        B[ip][is++] = wts[ip];
-
-        B[ip][is++] = wts[ip]*px;
+        B[ip][is++] = wts[0]*wts[ip];
+        B[ip][is++] = wts[0]*wts[ip]*px;
         B[ip][is++] = wts[ip]*sqr(px);
 
         if (dim_ >= 2)
@@ -255,14 +254,14 @@ Foam::label Foam::quadraticFitData::calcFit
     label nSVDzeros = 0;
 
     bool goodFit = false;
-    for(scalar tol = SMALL; tol < 0.1 && !goodFit; tol *= 10)
+    for(int iIt = 0; iIt < 35 && !goodFit; iIt++)
     {
-        SVD svd(B, tol);
+        SVD svd(B, SMALL);
 
-        scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
-        scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
+        scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
+        scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
 
-        goodFit = sign(fit0) == sign(fit1);
+        goodFit = sign(fit0) == sign(fit1) && fit0 < 1 && fit1 < 1;
 
         if (goodFit)
         {
@@ -270,12 +269,27 @@ Foam::label Foam::quadraticFitData::calcFit
             fit_[faci][1] = fit1;
             for(label i = 2; i < stencilSize; i++)
             {
-                fit_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
+                fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
             }
 
             singVals = svd.S();
             nSVDzeros = svd.nZeros();
         }
+        else // (not good fit so increase weight in the centre and for linear)
+        {
+            wts[0] *= 10;
+            wts[1] *= 10;
+            for(label i = 0; i < B.n(); i++)
+            {
+                B[i][0] *= 10;
+                B[i][1] *= 10;
+            }
+            for(label j = 0; j < B.m(); j++)
+            {
+                B[0][j] *= 10;
+                B[1][j] *= 10;
+            }
+        }
     }
     if (!goodFit)
     {
@@ -283,7 +297,7 @@ Foam::label Foam::quadraticFitData::calcFit
         (
             "quadraticFitData::calcFit(const pointField&, const label)"
             ) << "For face " << faci << endl
-            << "Fit not good even once tol >= 0.1"
+            << "Cannot find good fit"
             << exit(FatalError);
     }