/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2006-7 H. Jasak All rights reserved
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "mixerGgiFvMesh_modified.H"
#include "Time.H"
#include "regionSplit.H"
#include "ggiPolyPatch.H"
#include "polyPatchID.H"
#include "addToRunTimeSelectionTable.H"
#include "mapPolyMesh.H"
#include <iostream> //neu
#include <fstream> //neu
#include <stdio.h> //neu <for writing the log files
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(mixerGgiFvMesh_modified, 0);
    addToRunTimeSelectionTable(dynamicFvMesh, mixerGgiFvMesh_modified, IOobject);
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::mixerGgiFvMesh_modified::addZonesAndModifiers()
{
    // Add zones and modifiers for motion action

    if (cellZones().size() > 0)
    {
        Info<< "void mixerGgiFvMesh_modified::addZonesAndModifiers() : "
            << "Zones and modifiers already present.  Skipping."
            << endl;

        return;
    }

    Info<< "Time = " << time().timeName() << endl
        << "Adding zones and modifiers to the mesh" << endl;

    // Add zones
    List<pointZone*> pz(0);
    List<faceZone*> fz(0);
    List<cellZone*> cz(1);

    // Copy the face zones associated with the GGI interfaces
    if (faceZones().size() > 0)
    {
        // Copy face zones
        Info << "Copying existing face zones" << endl;

        fz.setSize(faceZones().size());

        forAll (faceZones(), i)
        {
            fz[i] = faceZones()[i].clone(faceZones()).ptr();
        }
    }

    regionSplit rs(*this);

    // Get the region of the cell containing the origin.
    label originRegion = rs[findNearestCell(cs().origin())];

    labelList movingCells(nCells());
    label nMovingCells = 0;

    forAll(rs, cellI)
    {
        if (rs[cellI] == originRegion)
        {
            movingCells[nMovingCells] = cellI;
            nMovingCells++;
        }
    }

    movingCells.setSize(nMovingCells);
    Info << "Number of cells in the moving region: " << nMovingCells << endl;

    cz[0] = new cellZone
    (
        "movingCells",
        movingCells,
        0,
        cellZones()
    );

    Info << "Adding point, face and cell zones" << endl;
    removeZones();
    addZones(pz, fz, cz);

    // Write mesh
    syncUpdateMesh();
    write();
}


void Foam::mixerGgiFvMesh_modified::calcMovingMasks() const
{
    if (debug)
    {
        Info<< "void mixerGgiFvMesh_modified::calcMovingMasks() const : "
            << "Calculating point and cell masks"
            << endl;
    }

    if (movingPointsMaskPtr_)
    {
        FatalErrorIn("void mixerGgiFvMesh_modified::calcMovingMasks() const")
            << "point mask already calculated"
            << abort(FatalError);
    }

    // Set the point mask
    movingPointsMaskPtr_ = new scalarField(allPoints().size(), 0);
    scalarField& movingPointsMask = *movingPointsMaskPtr_;

    const cellList& c = cells();
    const faceList& f = allFaces();

    label movingCellsID = cellZones().findZoneID("movingCells");

    if (movingCellsID < 0)
    {
        FatalErrorIn("void mixerGgiFvMesh_modified::calcMovingMasks() const")
            << "Cannot find moving cell zone ID"
            << abort(FatalError);
    }

    const labelList& cellAddr = cellZones()[movingCellsID];

    forAll (cellAddr, cellI)
    {
        const cell& curCell = c[cellAddr[cellI]];

        forAll (curCell, faceI)
        {
            // Mark all the points as moving
            const face& curFace = f[curCell[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 1;
            }
        }
    }

    // Grab the ggi patches on the moving side
    wordList movingPatches(dict_.subDict("slider").lookup("moving"));

    forAll (movingPatches, patchI)
    {
        const label movingSliderID =
            boundaryMesh().findPatchID(movingPatches[patchI]);

        if (movingSliderID < 0)
        {
            FatalErrorIn("void mixerGgiFvMesh_modified::calcMovingMasks() const")
                << "Moving slider named " << movingPatches[patchI]
                << " not found.  Valid patch names: "
                << boundaryMesh().names() << abort(FatalError);
        }

        const ggiPolyPatch& movingGgiPatch =
            refCast<const ggiPolyPatch>(boundaryMesh()[movingSliderID]);

        const labelList& movingSliderAddr = movingGgiPatch.zone();

        forAll (movingSliderAddr, faceI)
        {
            const face& curFace = f[movingSliderAddr[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 1;
            }
        }
    }

    // Grab the ggi patches on the static side
    wordList staticPatches(dict_.subDict("slider").lookup("static"));

    forAll (staticPatches, patchI)
    {
        const label staticSliderID =
            boundaryMesh().findPatchID(staticPatches[patchI]);

        if (staticSliderID < 0)
        {
            FatalErrorIn("void mixerGgiFvMesh_modified::calcMovingMasks() const")
                << "Static slider named " << staticPatches[patchI]
                << " not found.  Valid patch names: "
                << boundaryMesh().names() << abort(FatalError);
        }

        const ggiPolyPatch& staticGgiPatch =
            refCast<const ggiPolyPatch>(boundaryMesh()[staticSliderID]);

        const labelList& staticSliderAddr = staticGgiPatch.zone();

        forAll (staticSliderAddr, faceI)
        {
            const face& curFace = f[staticSliderAddr[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 0;
            }
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

// Construct from components
Foam::mixerGgiFvMesh_modified::mixerGgiFvMesh_modified
(
    const IOobject& io
)
:
    dynamicFvMesh(io),
    dict_
    (
        IOdictionary
        (
            IOobject
            (
                "dynamicMeshDict",
                time().constant(),
                *this,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        ).subDict(typeName + "Coeffs")
    ),
    csPtr_
    (
        coordinateSystem::New
        (
            "coordinateSystem",
            dict_.subDict("coordinateSystem")
        )
    ),
    rpm_(readScalar(dict_.lookup("rpm"))),
    IXX_(readScalar(dict_.lookup("IXX"))),
    IXY_(readScalar(dict_.lookup("IXY"))),
    IXZ_(readScalar(dict_.lookup("IXZ"))),
    IYY_(readScalar(dict_.lookup("IYY"))),
    IYZ_(readScalar(dict_.lookup("IYZ"))),
    IZZ_(readScalar(dict_.lookup("IZZ"))),
    MXX_(readScalar(dict_.lookup("MXX"))),
    MYY_(readScalar(dict_.lookup("MYY"))),    
    MZZ_(readScalar(dict_.lookup("MZZ"))),
    MaxMXX(readScalar(dict_.lookup("MaxMXX"))),
    MaxMYY(readScalar(dict_.lookup("MaxMYY"))),
    MaxMZZ(readScalar(dict_.lookup("MaxMZZ"))),
    movingPointsMaskPtr_(NULL)

{
    addZonesAndModifiers();

    Info<< "Mixer mesh:" << nl
        << "    origin: " << cs().origin() << nl
        << "    axis  : " << cs().axis() << nl
        << "    rpm   : " << rpm_ << nl

	<< "moment of inertia tensor" << nl
        << "\t\t" << IXX_ << "\t\t" << IXY_ << "\t\t" << IXZ_ << nl
        << "\t\t" << IXY_ << "\t\t" << IYY_ << "\t\t" << IYZ_ << nl
        << "\t\t" << IXZ_ << "\t\t" << IYZ_ << "\t\t" << IZZ_ << endl;

//calculation of inverse of moment of inertie tensor
detI_ = IXX_*IYY_*IZZ_+IXY_*IYZ_*IXZ_+IXZ_*IXY_*IYZ_-IXX_*IYZ_*IYZ_-IXY_*IXY_*IZZ_-IXZ_*IYY_*IXZ_;
iIXX_ = (IYY_*IZZ_-IYZ_*IYZ_)/detI_;
iIXY_ = (IXZ_*IYZ_-IXY_*IZZ_)/detI_;
iIXZ_ = (IXY_*IYZ_-IXZ_*IYY_)/detI_;
iIYY_ = (IXX_*IZZ_-IXZ_*IXZ_)/detI_;
iIYZ_ = (IXZ_*IXY_-IXX_*IYZ_)/detI_;
iIZZ_ = (IXX_*IYY_-IXY_*IXY_)/detI_;

        Info<< "    detI   : " << detI_ << nl

	<< "inverse of moment of inertia tensor" << nl
        << "\t\t" << iIXX_ << "\t\t" << iIXY_ << "\t\t" << iIXZ_ << nl
        << "\t\t" << iIXY_ << "\t\t" << iIYY_ << "\t\t" << iIYZ_ << nl
        << "\t\t" << iIXZ_ << "\t\t" << iIYZ_ << "\t\t" << iIZZ_ << nl

        << "    MXX   : " << MXX_ << nl
        << "    MYY   : " << MYY_ << nl
        << "    MZZ   : " << MZZ_ << nl

        << "   MaxMXX : " << MaxMXX << nl
        << "   MaxMYY : " << MaxMYY << nl
        << "   MaxMZZ : " << MaxMZZ << endl;
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::mixerGgiFvMesh_modified::~mixerGgiFvMesh_modified()
{
    deleteDemandDrivenData(movingPointsMaskPtr_);
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

// Return moving points mask.  Moving points marked with 1
const Foam::scalarField& Foam::mixerGgiFvMesh_modified::movingPointsMask() const
{
    if (!movingPointsMaskPtr_)
    {
        calcMovingMasks();
    }

    return *movingPointsMaskPtr_;
}


bool Foam::mixerGgiFvMesh_modified::update()
{



if (time().deltaT().value()==time().value())
{

//	rpm_alpha = MZZ_/IZZ_*time().deltaT().value()*60/2/3.14159265;
//	printf("RPMalpha = %f \n", rpm_alpha);
//	printf("RPM0 = %f \n", rpm_);
//	rpm_ = rpm_ + rpm_alpha;
//	printf("RPM = %f \n", rpm_);


// rotation around X-Axis
	rpmX_alpha = (iIXX_*MXX_+iIXY_*MYY_+iIXZ_*MZZ_)*time().deltaT().value()*60/2/3.14159265;
	printf("RPMXalpha = %f \n", rpmX_alpha);
	printf("RPMX0 = %f \n", rpmX_);
	rpmX_ = rpmX_ + rpmX_alpha;
	printf("RPMX = %f \n", rpmX_);


	//angle speed is written to file RPM.txt
		std::ofstream writeRPMXfile ("RPMX.txt");
		if (writeRPMXfile.is_open())
			{
			writeRPMXfile << rpmX_;
			writeRPMXfile.close();
			}
    		else
			{
			printf("unable to open file: RPMX.txt \n");
			}

	//creating log of rotatinal speed in file RPMX.log
	FILE *RPMXlog;
	RPMXlog = fopen("RPMX.log","a+");
	fprintf(RPMXlog,"Time: %lf   RPM: %f\n", time().value(), rpmX_);
	fclose(RPMXlog);


// rotation around Y-Axis
	rpmY_alpha = (iIXY_*MXX_+iIYY_*MYY_+iIYZ_*MZZ_)*time().deltaT().value()*60/2/3.14159265;
	printf("RPMYalpha = %f \n", rpmY_alpha);
	printf("RPMY0 = %f \n", rpmY_);
	rpmY_ = rpmY_ + rpmY_alpha;
	printf("RPMY = %f \n", rpmY_);


	//angle speed is written to file RPM.txt
		std::ofstream writeRPMYfile ("RPMY.txt");
		if (writeRPMYfile.is_open())
			{
			writeRPMYfile << rpmY_;
			writeRPMYfile.close();
			}
    		else
			{
			printf("unable to open file: RPMY.txt \n");
			}

	//creating log of rotatinal speed in file RPMY.log
	FILE *RPMYlog;
	RPMYlog = fopen("RPMY.log","a+");
	fprintf(RPMYlog,"Time: %lf   RPM: %f\n", time().value(), rpmY_);
	fclose(RPMYlog);


// rotation around Z-Axis
	rpmZ_alpha = (iIXZ_*MXX_+iIYZ_*MYY_+iIZZ_*MZZ_)*time().deltaT().value()*60/2/3.14159265;
	printf("RPMZalpha = %f \n", rpmZ_alpha);
	printf("RPMZ0 = %f \n", rpmZ_);
	rpmZ_ = rpmZ_ + rpmZ_alpha;
	printf("RPMZ = %f \n", rpmZ_);


	//angle speed is written to file RPM.txt
		std::ofstream writeRPMZfile ("RPMZ.txt");
		if (writeRPMZfile.is_open())
			{
			writeRPMZfile << rpmZ_;
			writeRPMZfile.close();
			}
    		else
			{
			printf("unable to open file: RPMZ.txt \n");
			}

	//creating log of rotatinal speed in file RPMZ.log
	FILE *RPMZlog;
	RPMZlog = fopen("RPMZ.log","a+");
	fprintf(RPMZlog,"Time: %lf   RPM: %f\n", time().value(), rpmZ_);
	fclose(RPMZlog);

// rotation around x-axis
    Info<< "cs().axis():" << csPtr_().axis() << nl
    << "cs().origin():" << csPtr_().origin() << nl
    << "cs().type():" << csPtr_().type() << nl
    << "cs().direction():" << csPtr_().direction() << nl
    << "cs().axis()[0]:" << csPtr_().axis()[0] << nl
    << "cs().axis()[1]:" << csPtr_().axis()[1] << nl
    << "cs().axis()[2]:" << csPtr_().axis()[2] << nl << nl << endl;

//-------------------------------




csPtr_().origin()[0] = 5;
csPtr_().origin()[1] = 6;
csPtr_().origin()[2] = 7;

std::cout << "csPtr_().origin()[0]: " << csPtr_().origin()[0] << nl;
std::cout << "csPtr_().origin()[1]: " << csPtr_().origin()[1] << nl;
std::cout << "csPtr_().origin()[2]: " << csPtr_().origin()[2] << nl;


csPtr_().axis()[0] = 2;
csPtr_().axis()[1] = 0;
csPtr_().axis()[2] = 0;

std::cout << "csPtr_().axis()[0]: " << csPtr_().axis()[0] << nl;
std::cout << "csPtr_().axis()[1]: " << csPtr_().axis()[1] << nl;
std::cout << "csPtr_().axis()[2]: " << csPtr_().axis()[2] << nl;


//csPtr_().axis()[0] = 2;
//csPtr_().axis()[1] = 0;
//csPtr_().axis()[2] = 0;


//    Info<< "cs().axis():" << csPtr_().axis() << nl
  //  << "cs().origin():" << csPtr_().origin() << nl
    //<< "cs().type():" << csPtr_().type() << nl
//    << "cs().direction():" << csPtr_().direction() << nl
  //  << "cs().axis()[0]:" << csPtr_().axis()[0] << nl
    //<< "cs().axis()[1]:" << csPtr_().axis()[1] << nl
    //<< "cs().axis()[2]:" << csPtr_().axis()[2] << nl << nl << endl;

Foam::mixerGgiFvMesh_modified
(
    const IOobject& io
)
:
    dynamicFvMesh(io),
    dict_
    (
        IOdictionary
        (
            IOobject
            (
                "dynamicMeshDict",
                time().constant(),
                *this,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        ).subDict(typeName + "Coeffs")
    ),
    csPtr_
    (
        coordinateSystem::New
        (
            "coordinateSystem",
            dict_.subDict("coordinateSystem")
        )
    );

std::cout << "eingelsen. glaub ich... ";

//------------------------------

movePoints
    (
        csPtr_->globalPosition
        (
            csPtr_->localPosition(allPoints())
          + vector(0, rpm_*360.0*time().deltaT().value()/60.0, 0)
            *movingPointsMask()
        )
    );
}

else
{
//get moments from files
	double MXX;	
	std::ifstream mxxfile ("MXX.txt");
	if (mxxfile.is_open())
		{
		while ( mxxfile.good() )
			{
			mxxfile >> MXX;

			}
		mxxfile.close();
		}
    	else
		{
		printf("unable to open file: MXX.txt \n");
		}



	double MYY;	
	std::ifstream myyfile ("MYY.txt");
	if (myyfile.is_open())
		{
		while ( myyfile.good() )
			{
			myyfile >> MYY;

			}
		myyfile.close();
		}
    	else
		{
		printf("unable to open file: MYY.txt \n");
		}



	double MZZ;	
	std::ifstream mzzfile ("MZZ.txt");
	if (mzzfile.is_open())
		{
		while ( mzzfile.good() )
			{
			mzzfile >> MZZ;

			}
		mzzfile.close();
		}
    	else
		{
		printf("unable to open file: MZZ.txt \n");
		}

//check Values of moments
if (MXX > abs(MaxMXX))
	{
	printf("high value of MXX \n");
	printf("MXX = %f \n", MXX);
	printf("MXX is set to MaxMXX \n");
	MXX = abs(MaxMXX);
	printf("MXX = %f \n", MXX);
	}
if (MXX < -abs(MaxMXX))
	{
	printf("high value of MXX \n");
	printf("MXX = %f \n", MXX);
	printf("MXX is set to MaxMXX \n");
	MXX = -abs(MaxMXX);
	printf("MXX = %f \n", MXX);
	}
else
	{
	printf("MXX = %f \n", MXX);
	}

if (MYY > abs(MaxMYY))
	{
	printf("high value of MYY \n");
	printf("MYY = %f \n", MYY);
	printf("MYY is set to MaxMYY \n");
	MYY = abs(MaxMYY);
	printf("MYY = %f \n", MYY);
	}
if (MYY < -abs(MaxMYY))
	{
	printf("high value of MYY \n");
	printf("MYY = %f \n", MYY);
	printf("MYY is set to MaxMYY \n");
	MYY = -abs(MaxMYY);
	printf("MYY = %f \n", MYY);
	}
else
	{
	printf("MYY = %f \n", MYY);
	}

if (MZZ > abs(MaxMZZ))
	{
	printf("high value of MZZ \n");
	printf("MZZ = %f \n", MZZ);
	printf("MZZ is set to MaxMzz \n");
	MZZ = abs(MaxMZZ);
	printf("MZZ = %f \n", MZZ);
	}
if (MZZ < -abs(MaxMZZ))
	{
	printf("high value of MZZ \n");
	printf("MZZ = %f \n", MZZ);
	printf("MZZ is set to MaxMzz \n");
	MZZ = -abs(MaxMZZ);
	printf("MZZ = %f \n", MZZ);
	}
else
	{
	printf("MZZ = %f \n", MZZ);
	}






// get RPM from file RPM.log
	std::ifstream openRPMfile ("RPM.txt");
	if (openRPMfile.is_open())
		{
		while ( openRPMfile.good() )
			{
			openRPMfile >> rpm_;

			}
		openRPMfile.close();
		}
    	else
		{
		printf("unable to open file: RPM.txt \n");
		}


//calculate angular acceleration with resulting moments

rpm_alpha = MZZ/IZZ_*time().deltaT().value()*60/2/3.14159265;

printf("RPM last time step= %f rounds/minute \n", rpm_);

rpm_ = rpm_ + rpm_alpha;

	//creating log of angular speed in file
	FILE *RPMlog;
	RPMlog= fopen("RPM.log","a+");
	fprintf(RPMlog,"Time: %lf   RPM: %f\n", time().value(), rpm_);
	fclose(RPMlog);

printf("RPM actual timestep = %f rounds/minute \n", rpm_);


	//angular speed is written to file RPM.txt
	FILE *RPMtxt;
	RPMtxt= fopen("RPM.txt","w+");
	fprintf(RPMtxt,"%f\n",rpm_);
	fclose(RPMtxt);


	//creating log files for moments and RPM
	FILE *rotationlog;
	rotationlog = fopen("rotation.log","a+");
	fprintf(rotationlog,"Time:\t %lf \tMXX:\t%f \tMYY:\t%f \tMZZ:\t%f \tRPM:\t%g \n", time().value(), MXX, MYY, MZZ, rpm_);
	fclose(rotationlog);

    movePoints
    (
        csPtr_->globalPosition
        (
            csPtr_->localPosition(allPoints())
          + vector(0, rpm_*360.0*time().deltaT().value()/60.0, 0)
            *movingPointsMask()
        )
    );
}


    // The mesh is not morphing
    return false;
}



// ************************************************************************* //