/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/
 
// ------------------------------------------------------------
//
// variableFFT.cc
//
// ------------------------------------------------------------

#include <string>
#include <stdexcept>
#include "../../../program/source/mympi.h"  // von 3D UGBlocks
#include "../abbrevi2D.h"
//#include "../../../linAlgExptemp/source/linAlg.h"
//#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../grid/elements2D.h"
#include "../grid/parti2D.h"
#include "../grid/ug2D.h"
#include "../grid/blockgrid2D.h"
#include "../grid/marker2D.h"
//#include "parallel.h"
#include "extemp2D.h"
#include "variable2D.h"
#include "variableFFT.h"
#include "variable2D_cc.h"
#include "update2D.h"
#include "print_var2D_vtk_cc.h"
#include "../grid/examples_ug2D.h"
#include "../../../common_source/mathlib/fouriertransform.h"
#include "../interpol/interpolTwoD.h"

int VariableFFT::powerTwo(int t) {
    int n=1;
    for(int l=0;l<t;++l) n=n*2;
    return n;
}

void VariableFFT::operator=(VariableFFT& v) {
    Variable2D<std::complex<double> >*pv = static_cast<Variable2D<std::complex<double> >* >(&v);
    static_cast<Variable2D<std::complex<double> >* >(this)->operator=(*pv); 
}

void VariableFFT::operator=(std::complex<double> x) {
    static_cast<Variable2D<std::complex<double> >* >(this)->operator=(x); 
}

void VariableFFT::interpolate(Variable2D<std::complex<double> >& v, 
                              std::complex<double> defaultInterpolation) {
     bool createInterpolator = (assignmentBlockgrid==NULL);
     if(createInterpolator==false) {
          if(assignmentBlockgrid->getId() != v.Give_blockgrid()->getId()) {
             delete assignmentInterpolator;
             createInterpolator=true;
          }
     }
     
     if(createInterpolator) {
        assignmentBlockgrid = v.Give_blockgrid();
        assignmentInterpolator 
           = new Interpolate_on_structured_2Dgrid(Nx+1,Ny+1,
				     D2vector(ug->MinimumX(),ug->MinimumY()),
                     D2vector(ug->MaximumX(),ug->MaximumY()),                                                  
                     *assignmentBlockgrid);
     }
     assignmentInterpolator->interpolate(v,data_rectangles[0],defaultInterpolation);
     Update_back(0);       
}

 
VariableFFT::VariableFFT(int tx, int ty, Rechteck& ugRectangle) :
    Variable2D<std::complex<double> >(
          *(blockgrid = new Blockgrid2D(&ugRectangle,powerTwo(tx)-1,powerTwo(ty)-1))) {

    ownBlockGrid = true;    
    init();    

    Hx = ugRectangle.getSizeX() / Nx;
    Hy = ugRectangle.getSizeY() / Ny;
}

VariableFFT::VariableFFT(VariableFFT* other) :
    Variable2D<std::complex<double> >(*(other->Give_blockgrid())) {
    ownBlockGrid = false;

    Hx = other->getHx();
    Hy = other->getHy();

    init();
}
 
 VariableFFT::VariableFFT(const VariableFFT& other) :
    Variable2D<std::complex<double> >(*(other.Give_blockgrid())) {
    ownBlockGrid = false;

    Hx = other.getHx();
    Hy = other.getHy();

    init();
}
 
void VariableFFT::init() {
    Nx = blockgrid->Give_Nx_rectangle(0);
    Ny = blockgrid->Give_Ny_rectangle(0);       
    
    startValueY = new complex<double>[Ny+1];  

    FourierTransform intF(MAX(Nx+1,Ny+1));
    assignmentInterpolator = NULL;
    assignmentBlockgrid    = NULL;
}


VariableFFT::~VariableFFT() {
    if(ownBlockGrid) blockgrid;
    delete[] startValueY;    
}

    
void VariableFFT::FFT() {
    Update<rectangleEl>(0);  
    
    for(int j=0;j<=Ny;++j) {
        std::complex<double>* startValue = &(data_rectangles[0][( Nx+1 ) *j]);
        FourierTransform::FFT(startValue,Nx+1);
    }
    
    for(int i=0;i<=Nx;++i) {
        for(int j=0;j<=Ny;++j) {
            startValueY[j] = data_rectangles[0][i + ( Nx+1 ) *j];
        }
        FourierTransform::FFT(startValueY,Ny+1);
        for(int j=0;j<=Ny;++j) {
           data_rectangles[0][i + ( Nx+1 ) *j] = startValueY[j];
        }
    }
    
    Update_back(0);
}
    

void VariableFFT::inversFFT() {    
    Update<rectangleEl>(0);  
    
    for(int i=0;i<=Nx;++i) {
        for(int j=0;j<=Ny;++j) {
            startValueY[j] = data_rectangles[0][i + ( Nx+1 ) *j];
        }
        FourierTransform::iFFT(startValueY,Nx+1);
        for(int j=0;j<=Ny;++j) {
            data_rectangles[0][i + ( Nx+1 ) *j] = startValueY[j];
        }
    }  
    
    for(int j=0;j<=Ny;++j) {
        complex<double>* startValue = &(data_rectangles[0][( Nx+1 ) *j]);
        FourierTransform::iFFT(startValue,Nx+1);
    }    
    
    Update_back(0);    
}