IHTP01inverseLaplacian.H

Go to the documentation of this file.

#ifndef IHTP01inverseLaplacian_H
#define IHTP01inverseLaplacian_H

class IHTP01inverseLaplacian_paramBC: public inverseLaplacianProblem_paramBC
{
    public:
        explicit IHTP01inverseLaplacian_paramBC(int argc, char* argv[])
            :
            inverseLaplacianProblem_paramBC(argc, argv),
            T(_T()),
            mesh(_mesh()),
            runTime(_runTime())
        {
            hotSide_ind = mesh.boundaryMesh().findPatchID("hotSide");
            coldSide_ind = mesh.boundaryMesh().findPatchID("coldSide");
            thermocouplesRead = 0;
        }

        volScalarField& T;

        fvMesh& mesh;

        Time& runTime;

        double a;

        double b;

        double c;

        double d;

        List<scalar> heatFlux_gammaEx1;

        List<scalar> heatFlux_gammaEx2;

        List<scalar> heatFlux_gammaEx3;

        List<scalar> heatFlux_gammaEx4;

        label gammaEx1_ind;

        label gammaEx2_ind;

        label gammaEx3_ind;

        label gammaEx4_ind;

        PtrList<volScalarField> gField;
 
 
        //--------------------------------------------------------------------------
        void solveAdditional()
        {
            restart();
            fvMesh& mesh = _mesh();
            simpleControl& simple = _simple();
            volScalarField Tad(_T());
            Foam::Time& runTime = _runTime();
            set_valueFraction();
            List<scalar> RobinBC = - Tf;
            forAll(mesh.boundaryMesh(), patchI)
            {
                if (patchI == mesh.boundaryMesh().findPatchID("coldSide"))
                {
                    ITHACAutilities::assignMixedBC(Tad, patchI, RobinBC, refGrad, valueFraction);
                }
                else if (patchI == mesh.boundaryMesh().findPatchID("hotSide"))
                {
                    ITHACAutilities::assignBC(Tad, patchI, homogeneousBC);
                }
                else if (patchI == mesh.boundaryMesh().findPatchID("gammaEx1"))
                {
                    ITHACAutilities::assignBC(Tad, patchI, -heatFlux_gammaEx1);
                }
                else if (patchI == mesh.boundaryMesh().findPatchID("gammaEx2"))
                {
                    ITHACAutilities::assignBC(Tad, patchI, -heatFlux_gammaEx2);
                }
                else if (patchI == mesh.boundaryMesh().findPatchID("gammaEx3"))
                {
                    ITHACAutilities::assignBC(Tad, patchI, -heatFlux_gammaEx3);
                }
                else if (patchI == mesh.boundaryMesh().findPatchID("gammaEx4"))
                {
                    ITHACAutilities::assignBC(Tad, patchI, -heatFlux_gammaEx4);
                }
            }
 
#if defined(OFVER) && (OFVER == 6)
 
            while (simple.loop(runTime))
#else
            while (simple.loop())
#endif
            {
                while (simple.correctNonOrthogonal())
                {
                    fvScalarMatrix TEqn
                    (
                        fvm::laplacian(DT, Tad)
                    );
                    TEqn.solve();
                }
            }
 
            Tad_base.append(Tad.clone());
            //Reinitializing runTime
            instantList Times = runTime.times();
            runTime.setTime(Times[1], 1);
            addSol = fieldValueAtThermocouples(Tad);
            ITHACAstream::exportSolution(Tad, "1",
                                         folderOffline,
                                         "Tad");
        }
 
        //--------------------------------------------------------------------------
        void postProcess(word folder, word heatFluxFieldName, scalar innerField = 0.0)
        {
#include"postProcess.H"
        }
 
 
        //--------------------------------------------------------------------------
        void assignDirectBC()
        {
#include"directBC.H"
        }
 
        //--------------------------------------------------------------------------
        void solveTrue()
        {
            word outputFolder = "./ITHACAoutput/true/";
#include"solveTrue.H"
        }
 
        //--------------------------------------------------------------------------
        Eigen::VectorXd bestApproximator()
        {
            fvMesh& mesh = _mesh();
            int Nbasis = gBaseFunctions.size();
            Eigen::MatrixXd mass(Nbasis, Nbasis);
            Eigen::VectorXd source(Nbasis);
            forAll(gBaseFunctions, baseI)
            {
                forAll(gBaseFunctions, baseJ)
                {
                    mass(baseI, baseJ) = ITHACAutilities::L2productOnPatch(mesh,
                                         gBaseFunctions[baseI], gBaseFunctions[baseJ], "hotSide");
                }
 
                source(baseI) = ITHACAutilities::L2productOnPatch(mesh, gBaseFunctions[baseI],
                                gTrue, "hotSide");
            }
 
            return mass.fullPivLu().solve(source); //Vector of coefficients
        }
 
 
        //--------------------------------------------------------------------------
        Eigen::VectorXd bestInterpolation()
        {
            int Nbasis = gBaseFunctions.size();
            Eigen::VectorXd coeff(Nbasis);
            Eigen::MatrixXd Phi = Eigen::MatrixXd::Zero(Nbasis, Nbasis);
            forAll(gBaseFunctions, baseI)
            {
                scalar centerX = thermocouplesPos[baseI].x();
                scalar centerZ = thermocouplesPos[baseI].z();
                forAll(thermocouplesPos, TCi)
                {
                    scalar tcX = thermocouplesPos[TCi].x();
                    scalar tcZ = thermocouplesPos[TCi].z();
                    double radius = Foam::sqrt((tcX - centerX) * (tcX - centerX) +
                                               (tcZ - centerZ) * (tcZ - centerZ));
                    Phi(TCi, baseI) = Foam::exp(- (shapeParameter *
                                                   shapeParameter
                                                   * radius * radius));
                }
 
                coeff(baseI) = k * (b * thermocouplesPos[baseI].x() + c);
            }
 
            return Phi.fullPivLu().solve(coeff);
        }
 
};
 
 

class IHTP01inverseLaplacian_CG: public inverseLaplacianProblem_CG
{
    public:
        explicit IHTP01inverseLaplacian_CG(int argc, char* argv[])
            :
            inverseLaplacianProblem_CG(argc, argv),
            T(_T()),
            lambda(_lambda()),
            deltaT(_deltaT()),
            mesh(_mesh()),
            runTime(_runTime())
        {
            hotSide_ind = mesh.boundaryMesh().findPatchID("hotSide");
            coldSide_ind = mesh.boundaryMesh().findPatchID("coldSide");
            interpolationPlaneDefined = 0;
            cgIter = 0;
            thermocouplesRead = 0;
        }

        volScalarField& T;

        volScalarField& lambda;

        volScalarField& deltaT;

        fvMesh& mesh;

        Time& runTime;

        double a;

        double b;

        double c;

        double d;

        List<scalar> heatFlux_gammaEx1;

        List<scalar> heatFlux_gammaEx2;

        List<scalar> heatFlux_gammaEx3;

        List<scalar> heatFlux_gammaEx4;

        label gammaEx1_ind;

        label gammaEx2_ind;

        label gammaEx3_ind;

        label gammaEx4_ind;

        PtrList<volScalarField> gField;
 
        //--------------------------------------------------------------------------
        void postProcess(word folder, word heatFluxFieldName, scalar innerField = 0.0)
        {
#include"postProcess.H"
        }
 
 
        //--------------------------------------------------------------------------
        void assignDirectBC()
        {
#include"directBC.H"
        }
 
        //--------------------------------------------------------------------------
        void solveTrue()
        {
            word outputFolder = "./ITHACAoutput/true/";
#include"solveTrue.H"
        }
};
 
#endif