Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -192,7 +192,7 @@ updateCoeffs()
     scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
     scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
-    const vector& myRayId = dom.IRay(rayId).d();
+    const vector& myRayId = dom.IRay(rayId).dAve();
 
     // Use updated Ir while iterating over rays
     // avoids to used lagged qin

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -87,14 +87,6 @@ void Foam::radiation::fvDOM::initialise()
     // 2D
     else if (mesh_.nSolutionD() == 2)
     {
-        // Currently 2D solution is limited to the x-y plane
-        if (mesh_.solutionD()[vector::Z] != -1)
-        {
-            FatalErrorInFunction
-                << "Currently 2D solution is limited to the x-y plane"
-                << exit(FatalError);
-        }
-
         scalar thetai = piByTwo;
         scalar deltaTheta = pi;
         nRay_ = 4*nPhi_;
@@ -127,14 +119,6 @@ void Foam::radiation::fvDOM::initialise()
     // 1D
     else
     {
-        // Currently 1D solution is limited to the x-direction
-        if (mesh_.solutionD()[vector::X] != 1)
-        {
-            FatalErrorInFunction
-                << "Currently 1D solution is limited to the x-direction"
-                << exit(FatalError);
-        }
-
         scalar thetai = piByTwo;
         scalar deltaTheta = pi;
         nRay_ = 2;
@@ -187,20 +171,28 @@ void Foam::radiation::fvDOM::initialise()
         );
     }
 
-    Info<< "fvDOM : Allocated " << IRay_.size()
-        << " rays with average orientation:" << nl;
 
+    // Calculate the maximum solid angle
     forAll(IRay_, rayId)
     {
         if (omegaMax_ <  IRay_[rayId].omega())
         {
             omegaMax_ = IRay_[rayId].omega();
         }
-        Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : "
-            << '\t' << IRay_[rayId].omega() << nl;
     }
 
-    Info<< endl;
+
+    Info<< typeName << ": Created " << IRay_.size() << " rays with average "
+        << "directions (dAve) and solid angles (omega)" << endl;
+    Info<< incrIndent;
+    forAll(IRay_, rayId)
+    {
+        Info<< indent
+            << "Ray " << IRay_[rayId].I().name() << ": "
+            << "dAve = " << IRay_[rayId].dAve() << ", "
+            << "omega = " << IRay_[rayId].omega() << endl;
+    }
+    Info<< decrIndent << endl;
 }
 
 
@@ -242,8 +234,8 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
             "qem",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh_,
         dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -49,11 +49,14 @@ Description
         solverFreq   1; // Number of flow iterations per radiation iteration
     \endverbatim
 
-    The total number of solid angles is  4*nPhi*nTheta.
+    In 1-D the ray directions are bound to one of the X, Y or Z directions. The
+    total number of solid angles is 2. nPhi and nTheta are ignored.
 
-    In 1D the direction of the rays is X (nPhi and nTheta are ignored)
+    In 2-D the ray directions are within one of the X-Y, X-Z or Y-Z planes. The
-    In 2D the direction of the rays is on X-Y plane (only nPhi is considered)
+    total number of solid angles is 4*nPhi. nTheta is ignored.
-    In 3D (nPhi and nTheta are considered)
+
+    In 3D the rays span all directions. The total number of solid angles is
+    4*nPhi*nTheta.
 
 SourceFiles
     fvDOM.C

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -137,6 +137,42 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
     );
 
+    // Transform directions so that they fall inside the bounds of reduced
+    // dimension cases
+    if (mesh_.nSolutionD() == 2)
+    {
+        vector meshDir(vector::zero);
+        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+        {
+            if (mesh_.geometricD()[cmpt] == -1)
+            {
+                meshDir[cmpt] = 1;
+            }
+        }
+        const vector normal(vector(0, 0, 1));
+
+        const tensor coordRot = rotationTensor(normal, meshDir);
+
+        dAve_ = coordRot & dAve_;
+        d_ = coordRot & d_;
+    }
+    else if (mesh_.nSolutionD() == 1)
+    {
+        vector meshDir(vector::zero);
+        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+        {
+            if (mesh_.geometricD()[cmpt] == 1)
+            {
+                meshDir[cmpt] = 1;
+            }
+        }
+        const vector normal(vector(1, 0, 0));
+
+        dAve_ = (dAve_ & normal)*meshDir;
+        d_ = (d_ & normal)*meshDir;
+    }
+
+
     autoPtr<volScalarField> IDefaultPtr;
 
     forAll(ILambda_, lambdaI)
@@ -210,12 +246,12 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 
     scalar maxResidual = -GREAT;
 
+    const surfaceScalarField Ji(dAve_ & mesh_.Sf());
+
     forAll(ILambda_, lambdaI)
     {
         const volScalarField& k = dom_.aLambda(lambdaI);
 
-        const surfaceScalarField Ji(dAve_ & mesh_.Sf());
-
         fvScalarMatrix IiEq
         (
             fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")