/*
 Copyright (C) 2000-2004

 Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane

 This file is part of xmds.
 
 This program 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.

 This program 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 this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

/*
  $Id: xmdsintegrateex.cc,v 1.14 2004/07/13 05:29:38 paultcochrane Exp $
*/

/*! @file xmdsintegrateex.cc
  @brief Integrate element parsing classes and methods; explicit picture

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateEX public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateEXs=0;  //!< Number of xmds integrate EX objects

// ******************************************************************************
xmdsIntegrateEX::xmdsIntegrateEX(
				 const xmdsSimulation *const yourSimulation,
				 const bool& yourVerboseMode) :
  xmdsIntegrate(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsIntegrateEXs++;
    printf("xmdsIntegrateEX::xmdsIntegrateEX\n");
    printf("nxmdsIntegrateEXs=%li\n",nxmdsIntegrateEXs);
  }
};

// ******************************************************************************
xmdsIntegrateEX::~xmdsIntegrateEX() {
  if(debugFlag) {
    nxmdsIntegrateEXs--;
    printf("xmdsIntegrateEX::~xmdsIntegrateEX\n");
    printf("nxmdsIntegrateEXs=%li\n",nxmdsIntegrateEXs);
  }
};

// ******************************************************************************
void xmdsIntegrateEX::processElement(
				     const Element *const yourElement) {
  if(debugFlag) {
    printf("xmdsIntegrateEX::processElement\n");
  }

  // ************************************
  // parse code for operators

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateEX::processElement: cannot find 'main' vector");
  }

  XMLString* theCode = propagationCode();

  XMLString nextOperatorName;
  XMLString nextComponentName;
  unsigned long start=0;
  unsigned long end=0;
  char s[256];

  list<XMLString> coVectorComponentNamesList;

  // first need to get all co pairs

  while(end<theCode->length()) {
    if(findNextcoPair(nextOperatorName,nextComponentName,start,end)) {

      unsigned long nextComponentNumber;

      if(!mainVector->getComponent(nextComponentName,nextComponentNumber)) {
	sprintf(errorMessage(),"[%s] is not a component of the main vector",nextComponentName.c_str());
	throw xmdsException(yourElement,errorMessage());
      }

      unsigned long nextOperatorNumber;

      if(!getKOperator(nextOperatorName,nextOperatorNumber)) {
	sprintf(errorMessage(),"'%s' was not defined in <k_operators>",nextOperatorName.c_str());
	throw xmdsException(yourElement,errorMessage());
      }

      sprintf(s,"_segment%li_co_term_%s_%s",segmentNumber,nextOperatorName.c_str(),nextComponentName.c_str());

      const XMLString nextVectorComponentName = s;

      unsigned long coKey;

      if(!getcoKey(nextComponentNumber,nextOperatorNumber,coKey)) {

	if(verbose()) {
	  printf("adding operator-component pair: %s[%s]\n",nextOperatorName.c_str(),nextComponentName.c_str());
	}

	coKey = addcoPair(nextComponentNumber,nextOperatorNumber);

	coVectorComponentNamesList.push_back(nextVectorComponentName);
      }

      theCode->replaceData(start,end-start,nextVectorComponentName);

      start = start+nextVectorComponentName.length();
    }
    else {
      start = end;
    }
  }

  if(coVectorComponentNamesList.size()>0) {

    xmdsVector* coVector = simulation()->field()->createxmdsVector();

    sprintf(s,"segment%li_co_terms",segmentNumber);

    coVector->setName(s);

    coVector->setVectorType(COMPLEX);

    coVector->setComponents(coVectorComponentNamesList);

    list<XMLString> tempVectorNamesList;

    tempVectorNamesList.push_back("main");
    tempVectorNamesList.push_back(s);

    simulation()->field()->processVectors(tempVectorNamesList,simulation()->field()->geometry()->fullSpace());

    vectorNamesList()->push_back(s);
  }

};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateEX protected
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsIntegrateEX::writePrototypes(
				      FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateEX::writePrototypes\n");
  }

  if(usesKOperators()) {

    fprintf(outfile,"// integrate (EX) prototypes\n");
    fprintf(outfile,"\n");
    if(constantK()) {
      fprintf(outfile,"void _segment%li_calculate_k_operator_field();\n",segmentNumber);
      fprintf(outfile,"\n");
    }
    fprintf(outfile,"void _segment%li_calculate_co_terms();\n",segmentNumber);
  }
};

// ******************************************************************************
void xmdsIntegrateEX::writeRoutines(
				    FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateEX::writeRoutines\n");
  }

  if(usesKOperators()) {

    if(constantK()) {
      writeCalculatekOperatorFieldRoutine(outfile);
      writeCalculateCOTermsConstantKRoutine(outfile);
    }
    else
      writeCalculateCOTermsTimeDepKRoutine(outfile);
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateEX private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsIntegrateEX::writeCalculatekOperatorFieldRoutine(
							  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateEX::writeCalculatekOperatorFieldRoutine\n");
  }

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();

  unsigned long i;

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_k_operator_field() {\n",segmentNumber);
  fprintf(outfile,"\n");
  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      fprintf(outfile,"_segment%li_k_operator_field = new complex[total_local_size*_segment%li_nkoperators];\n",
	      segmentNumber,segmentNumber);
      fprintf(outfile,"\n");
    }

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");

  fprintf(outfile,"unsigned long _k_operator_index_pointer=0;\n");

  simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());

  char indent[64];
  for(i=0;i<nDims;i++) {
    indent[i]=0x09;
  }
  indent[nDims]=0;

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"%s_segment%li_k_operator_field[_k_operator_index_pointer + %li] = %s;\n",
	    indent,segmentNumber,i,KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
  simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());

  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsIntegrateEX::writeCalculateCOTermsConstantKRoutine(
							    FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateEX::writeCalculateCOTermsConstantKRoutine\n");
  }

  const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
  const char *const fieldName = simulation()->field()->name()->c_str();

  char coVectorName[256];
  sprintf(coVectorName,"segment%li_co_terms",segmentNumber);

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateEX::writeCalculateCOTermsConstantKRoutine: cannot find 'main' vector");
  }

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_co_terms() {\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"double _temp;\n");
  fprintf(outfile,"unsigned long _segment%li_kop_index_pointer=0;\n",segmentNumber);
  fprintf(outfile,"unsigned long _active_%s_index_pointer=0;\n",fieldName);
  fprintf(outfile,"unsigned long _%s_segment_%li_co_terms_index_pointer=0;\n",fieldName,segmentNumber);
  fprintf(outfile,"\n");

  list<XMLString> tempVectorNamesList;
  tempVectorNamesList.push_back("main");

  simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,"");

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"for(long _i0=0;_i0<total_local_size;_i0++) {\n");
  }
  else {
    fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size;_i0++) {\n",fieldName);
  }
  fprintf(outfile,"\n");

  for(unsigned long i=0; i<mainVector->nComponents(); i++) {

    const coStruct* thecoStruct;

    if(getcoStruct(i,thecoStruct)) {

      list<unsigned long>::const_iterator pULong1 = thecoStruct->operatorNumbersList.begin();
      list<unsigned long>::const_iterator pULong2 = thecoStruct->coKeysList.begin();

      while(pULong1 != thecoStruct->operatorNumbersList.end()) {

	fprintf(outfile,"\n");
	fprintf(outfile,"	_temp = _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
		segmentNumber,segmentNumber,*pULong1);
	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].re\n",fieldName,fieldName,i);
	fprintf(outfile,"	      - _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,segmentNumber,*pULong1);
	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].im;\n",fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"	_%s_segment%li_co_terms[_%s_segment_%li_co_terms_index_pointer + %li].im =\n",
		fieldName,segmentNumber,fieldName,segmentNumber,*pULong2);
	fprintf(outfile,"	     _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
		segmentNumber,segmentNumber,*pULong1);
	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].im\n",fieldName,fieldName,i);
	fprintf(outfile,"	   + _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,segmentNumber,*pULong1);
	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].re;\n",fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"	_%s_segment%li_co_terms[_%s_segment_%li_co_terms_index_pointer + %li].re=_temp;\n",
		fieldName,segmentNumber,fieldName,segmentNumber,*pULong2);

	pULong1++;
	pULong2++;
      }
    }
  }

  fprintf(outfile,"  _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
  fprintf(outfile,"  _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
  fprintf(outfile,"  _%s_segment_%li_co_terms_index_pointer+=_%s_segment%li_co_terms_ncomponents;\n",fieldName,segmentNumber,fieldName,segmentNumber);
  fprintf(outfile,"	}\n");

  fprintf(outfile,"\n");
  fprintf(outfile,"_%s_segment%li_co_terms_space=%li;\n",fieldName,segmentNumber,fullSpace);
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsIntegrateEX::writeCalculateCOTermsTimeDepKRoutine(
							   FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateIP::writeCalculateCOTermsTimeDepKRoutine\n");
  }

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
  const char *const fieldName = simulation()->field()->name()->c_str();

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateIP::writeCalculateCOTermsTimeDepKRoutine: cannot find 'main' vector");
  }

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_co_terms() {\n",segmentNumber);
  fprintf(outfile,"\n");
  for(unsigned long i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"double _temp;\n");
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");

  list<XMLString> tempVectorNamesList = *KVectorNamesList();
  char coVectorName[256];
  sprintf(coVectorName,"segment%li_co_terms",segmentNumber);
  tempVectorNamesList.push_back(coVectorName);

  simulation()->field()->openLoops(outfile,fullSpace,tempVectorNamesList);

  char indent[64];
  for(unsigned long i=0;i<nDims;i++) {
    indent[i]=0x09;
  }
  indent[nDims]=0;

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  for(unsigned long i=0; i<mainVector->nComponents(); i++) {

    const coStruct* thecoStruct;

    if(getcoStruct(i,thecoStruct)) {

      for (list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {

	fprintf(outfile,"%s_temp =  %s.re*%s.re - %s.im*%s.im;\n",
		indent,KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str(),
		KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str());
	fprintf(outfile,"%s_segment%li_co_term_%s_%s.im =  %s.re*%s.im + %s.im*%s.re;\n",
		indent,segmentNumber,KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str(),
		KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str(),
		KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str());
	fprintf(outfile,"%s_segment%li_co_term_%s_%s.re=_temp;\n",
		indent,segmentNumber,KOperator(*pULong)->c_str(),mainVector->componentName(i)->c_str());
	fprintf(outfile,"\n");
      }
    }
  }

  simulation()->field()->closeLoops(outfile,fullSpace,tempVectorNamesList);

  fprintf(outfile,"\n");
  fprintf(outfile,"_%s_segment%li_co_terms_space=%li;\n",fieldName,segmentNumber,fullSpace);
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};