diff --git a/libs/Allwmake b/libs/Allwmake
index a494487..16f1269 100755
--- a/libs/Allwmake
+++ b/libs/Allwmake
@@ -6,5 +6,6 @@ targetType=libso
wmake $targetType chemistryModel_POSTECH
wmake $targetType combustionModels_POSTECH
+wmake $targetType thermos
#------------------------------------------------------------------------------
diff --git a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
index 8ed9abd..8a49cbc 100644
--- a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
+++ b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C
@@ -48,6 +48,9 @@ Contact
#include "IFstream.H"
#include "OFstream.H"
#include "Math.H" //Mathematical functions for CMC (AMC, gammaln, TDMA)
+
+#include "LinearMesh.H"
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@@ -73,6 +76,8 @@ int main(int argc, char *argv[])
#include "startSummary.H" //Make logSummary file for CMC
#include "readCMCProperties.H" //Read and set fields for CMC calculation
+ #include "createLinearMesh.H"
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@@ -87,6 +92,8 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
+ #include "QiEqn.H"
+
#include "rhoEqn.H"
#include "UEqn.H"
if(init_start_CMC == true)
diff --git a/solvers_post/LagrangianCMCFoam/LinearMesh.H b/solvers_post/LagrangianCMCFoam/LinearMesh.H
new file mode 100644
index 0000000..862eae3
--- /dev/null
+++ b/solvers_post/LagrangianCMCFoam/LinearMesh.H
@@ -0,0 +1,337 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Class
+ Foam::LinearMesh
+
+Description
+
+SourceFiles
+ LinearMeshI.H
+ LinearMesh.C
+ LinearMeshIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LinearMesh_H
+#define LinearMesh_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class LinearMesh;
+Istream& operator>>(Istream&, LinearMesh&);
+Ostream& operator<<(Ostream&, const LinearMesh&);
+
+
+/*---------------------------------------------------------------------------*\
+ Class LinearMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class LinearMesh
+{
+ // Private data
+
+ //- Lower end value of 1-D domain
+ scalar x0_;
+
+ //- Upper end value of 1-D domain
+ scalar x1_;
+
+ //- Total number of cells
+ label nCells_;
+
+ //- Total number of cells
+ pointField points_;
+
+ //- Total number of cells
+ faceList faces_;
+
+ //- Total number of cells
+ labelList owner_;
+
+ //- Total number of cells
+ labelList neighbour_;
+
+ //- Total number of cells
+ label nPatches_;
+
+ //- Total number of cells
+ wordList patchNames_;
+
+ //- Total number of cells
+ labelList patchStartIndices_;
+
+ //- Total number of cells
+ labelList patchSizes_;
+
+
+ // Private Member Functions
+
+ //- Disallow Construct null
+ LinearMesh();
+
+ //- Disallow default bitwise copy construct
+ LinearMesh(const LinearMesh&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const LinearMesh&);
+
+ //- Generate points and faces
+ void init(const scalarField &etaValue)
+ {
+ for(label i = 0; i < (nCells_+1); i++)
+ {
+ points_[4*i+0] = vector(etaValue[i], 0.0, 0.0);
+ points_[4*i+1] = vector(etaValue[i], 0.1, 0.0);
+ points_[4*i+2] = vector(etaValue[i], 0.1, 0.1);
+ points_[4*i+3] = vector(etaValue[i], 0.0, 0.1);
+ }
+
+ const label nBoundaryFaces = 2;
+ const label nInternalFaces = nCells_-1;
+ const label nEmptyFaces = 4*nCells_;
+
+ for(label i = 0; i < nInternalFaces; i++)
+ {
+ IStringStream faceFormatStream ("4(4 5 6 7)");
+ labelField faceFormat (faceFormatStream);
+ faces_[i] = face(faceFormat + 4*i);
+ }
+
+ {
+ IStringStream faceFormatStream ("4(0 3 2 1)");
+ faces_[nInternalFaces] = face(faceFormatStream);
+ }
+
+ {
+ IStringStream faceFormatStream ("4(4 5 6 7)");
+ labelField faceFormat (faceFormatStream);
+ faces_[nInternalFaces+1] = face(faceFormat + 4*nInternalFaces);
+ }
+
+ for(label i = 0; i < nCells_; i++)
+ {
+ IStringStream faceFormatStream ("4(4(0 1 5 4)4(1 2 6 5)4(2 3 7 6)4(3 0 4 7))");
+ List faceFormats (faceFormatStream);
+
+ label fi = nInternalFaces + nBoundaryFaces + 4*i;
+ faces_[fi+0] = face(faceFormats[0] + 4*i);
+ faces_[fi+1] = face(faceFormats[1] + 4*i);
+ faces_[fi+2] = face(faceFormats[2] + 4*i);
+ faces_[fi+3] = face(faceFormats[3] + 4*i);
+ }
+
+
+ for(label i = 0; i < nInternalFaces; i++)
+ {
+ owner_[i] = i;
+ neighbour_[i] = i+1;
+ }
+
+ {
+ owner_[nInternalFaces] = 0;
+ }
+
+ {
+ owner_[nInternalFaces+1] = nCells_ - 1;
+ }
+
+ for(label i = 0; i < nCells_; i++)
+ {
+ label fi = nInternalFaces + nBoundaryFaces + 4*i;
+ owner_[fi+0] = i;
+ owner_[fi+1] = i;
+ owner_[fi+2] = i;
+ owner_[fi+3] = i;
+ }
+
+
+ nPatches_ = 3;
+
+ IStringStream sInd
+ (
+ word("3(")
+ + name(nInternalFaces) + " "
+ + name(nInternalFaces+1) + " "
+ + name(nInternalFaces+2) + " "
+ + word(")")
+ );
+ patchStartIndices_ = labelList(sInd);
+
+ IStringStream sSize(word("3(1 1 ")+name(4*nCells_)+word(")"));
+ patchSizes_ = labelList(sSize);
+
+ IStringStream sName("3(lowerEnd upperEnd Sides)");
+ patchNames_ = wordList(sName);
+ }
+
+public:
+
+ // Static data members
+
+ //- Static data staticData
+ // static const dataType staticData;
+
+
+ // Constructors
+
+ //- Construct from components
+ LinearMesh(const scalar x0, const scalar x1, const label nCells)
+ :
+ x0_(x0),
+ x1_(x1),
+ nCells_(nCells),
+ points_(4*(nCells+1)),
+ faces_(5*nCells+1),
+ owner_(faces_.size()),
+ neighbour_(nCells-1)
+ {
+
+ scalarField etaValue(nCells+1, 0.0);
+ scalar delta = (x1 - x0) / scalar(nCells);
+ forAll(etaValue, i)
+ {
+ etaValue[i] = x0 + delta*i;
+ }
+
+
+ init(etaValue);
+
+ }
+
+ //- Construct from components
+ LinearMesh(const scalarField &x)
+ :
+ x0_(x[0]),
+ x1_(x.last()),
+ nCells_(x.size()-1),
+ points_(4*(nCells_+1)),
+ faces_(5*nCells_+1),
+ owner_(faces_.size()),
+ neighbour_(nCells_-1)
+ {
+
+ init(x);
+
+ }
+
+ //- Construct from Istream
+ LinearMesh(Istream&);
+
+ //- Construct as copy
+ // LinearMesh(const LinearMesh&);
+
+
+ // Selectors
+
+ //- Select null constructed
+ // static autoPtr New();
+
+
+ //- Destructor
+ ~LinearMesh() {}
+
+
+ // Member Functions
+
+ // Access
+ //- Lower end value of 1-D domain
+ scalar x0() const {return x0_;};
+
+ //- Upper end value of 1-D domain
+ scalar x1() const {return x1_;};
+
+ //- Total number of cells
+ label nCells() const {return nCells_;};
+
+ //- Total number of cells
+ const pointField &points() const {return points_;};
+
+ //- Total number of cells
+ const faceList &faces() const {return faces_;};
+
+ //- Total number of cells
+ const labelList &owner() const {return owner_;};
+
+ //- Total number of cells
+ const labelList &neighbour() const {return neighbour_;};
+
+ //- Total number of cells
+ label nPatches () const {return nPatches_;};
+
+ //- Total number of cells
+ word patchName (const label i) const {return patchNames_[i];};
+
+ //- Total number of cells
+ label patchStartIndex (const label i) const {return patchStartIndices_[i];};
+
+ //- Total number of cells
+ label patchSize (const label i) const {return patchSizes_[i];};
+
+
+ // Check
+
+ // Edit
+
+ // Write
+
+
+ // Member Operators
+
+ // void operator=(const LinearMesh&);
+
+
+ // Friend Functions
+
+ // Friend Operators
+
+ // IOstream Operators
+
+ friend Istream& operator>>(Istream&, LinearMesh&);
+ friend Ostream& operator<<(Ostream&, const LinearMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// #include "LinearMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/solvers_post/LagrangianCMCFoam/QiEqn.H b/solvers_post/LagrangianCMCFoam/QiEqn.H
new file mode 100644
index 0000000..c8e84c8
--- /dev/null
+++ b/solvers_post/LagrangianCMCFoam/QiEqn.H
@@ -0,0 +1,29 @@
+{
+ fvScalarMatrix testSource(testEtaMesh, testEtaMesh.dimensions()*dimVol/dimTime);
+
+ testSource.source() = 1.0;
+
+ fvScalarMatrix QiEqn
+ (
+ fvm::ddt(testEtaMesh)
+// + mvConvection->fvmDiv(phi, he)
+// + fvc::ddt(rho, K) + fvc::div(phi, K)
+// + (
+// he.name() == "ea"
+// ? fvc::div
+// (
+// fvc::absolute(phi/fvc::interpolate(rho), U),
+// p,
+// "div(phiv,p)"
+// )
+// : -dpdt
+// )
+ - fvm::laplacian(DQi, testEtaMesh)
+ ==
+ testSource
+ );
+
+ QiEqn.relax();
+
+ QiEqn.solve();
+}
diff --git a/solvers_post/LagrangianCMCFoam/createLinearMesh.H b/solvers_post/LagrangianCMCFoam/createLinearMesh.H
new file mode 100644
index 0000000..7749bcf
--- /dev/null
+++ b/solvers_post/LagrangianCMCFoam/createLinearMesh.H
@@ -0,0 +1,145 @@
+Info<< "Constructing single cell mesh" << nl << endl;
+
+LinearMesh lMesh(0, 1, 100);
+
+pointField CMCPoints(lMesh.points());
+faceList CMCFaces(lMesh.faces());
+labelList CMCOwners(lMesh.owner());
+labelList CMCNeighbors(lMesh.neighbour());
+
+fvMesh mfMesh
+(
+ IOobject
+ (
+ "mfSpace",
+ runTime.timeName(),
+ runTime,
+ IOobject::READ_IF_PRESENT
+ ),
+ xferMove>(CMCPoints),
+ CMCFaces.xfer(),
+ CMCOwners.xfer(),
+ CMCNeighbors.xfer()
+);
+
+if (mfMesh.boundaryMesh().empty())
+{
+ List mfPatches(3);
+
+ mfPatches[0] = new polyPatch
+ (
+ lMesh.patchName(0),
+ lMesh.patchSize(0),
+ lMesh.patchStartIndex(0),
+ 0,
+ mfMesh.boundaryMesh(),
+ polyPatch::typeName
+ );
+
+ mfPatches[1] = new polyPatch
+ (
+ lMesh.patchName(1),
+ lMesh.patchSize(1),
+ lMesh.patchStartIndex(1),
+ 1,
+ mfMesh.boundaryMesh(),
+ polyPatch::typeName
+ );
+
+ mfPatches[2] = new emptyPolyPatch
+ (
+ lMesh.patchName(2),
+ lMesh.patchSize(2),
+ lMesh.patchStartIndex(2),
+ 2,
+ mfMesh.boundaryMesh(),
+ emptyPolyPatch::typeName
+ );
+
+ mfMesh.addFvPatches(mfPatches);
+}
+
+IStringStream tempStream("3(fixedValue fixedValue empty)");
+wordList testEtaPatchTypeNames(tempStream);
+
+volScalarField testEtaMesh
+(
+ IOobject
+ (
+ "testEtaMesh",
+ runTime.timeName(),
+ mfMesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mfMesh,
+ dimensionedScalar("testEtaMesh", dimless, 1.0),
+ testEtaPatchTypeNames
+);
+
+const volVectorField &cEtaMesh(mfMesh.C());
+forAll(cEtaMesh, celli)
+{
+ const scalar x = cEtaMesh[celli].x();
+ if (x > 0.5)
+ {
+ testEtaMesh[celli] = 1.0;
+ }
+ else
+ {
+ testEtaMesh[celli] = 0.0;
+ }
+
+}
+
+volScalarField DQi
+(
+ IOobject
+ (
+ "DQi",
+ runTime.timeName(),
+ mfMesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mfMesh,
+ dimensionedScalar("DQi", dimArea/dimTime, 0.01),
+ testEtaPatchTypeNames
+);
+
+PtrList > FStruct (nf);
+for(label k = 0 ; k < nf ; k++)
+{
+ FStruct.set (k, new PtrList (Ysize+3));
+
+ PtrList &Fk(FStruct[k]);
+ Info << Fk.size() << endl;
+ for (label i = 0; i < Ysize + 3; i++)
+ {
+ Fk.set (i,
+ new volScalarField
+ (
+ IOobject
+ (
+ "test",
+ runTime.timeName(),
+ mfMesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mfMesh,
+ dimensionedScalar("test", dimless, 0.0),
+ testEtaPatchTypeNames
+ )
+ );
+
+ }
+}
+
+for(label k = 0 ; k < nf ; k++)
+{
+ for (label i = 0; i < Ysize + 3; i++)
+ {
+ FStruct[k][i].rename("test."+name(k)+"."+name(i));
+ }
+}