xmdsmomentgroup.cc

XMDS is a code generator that integrates equations. You write them down in human readable form in a

/*
 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: xmdsmomentgroup.cc,v 1.39 2004/10/21 09:48:12 paultcochrane Exp $
*/

/*! @file xmdsmomentgroup.cc
  @brief Moment group handling classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//				xmdsMomentGroup public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsMomentGroups=0;  //!< Number of xmds moment group objects

// ******************************************************************************
xmdsMomentGroup::xmdsMomentGroup(
				 const xmdsSimulation *const yourSimulation,
				 const bool& yourVerboseMode,
				 const unsigned long& yourGroupNumber) :
  xmdsField(yourSimulation,yourVerboseMode),
  myGroupNumber(yourGroupNumber) {
  if(debugFlag) {
    nxmdsMomentGroups++;
    printf("xmdsMomentGroup::xmdsMomentGroup\n");
    printf("nxmdsMomentGroups=%li\n",nxmdsMomentGroups);
  }
  nSamples=0;
};

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

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

  if(verbose()) {
    printf("Processing moment group %li ...\n",myGroupNumber);
  }

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

  unsigned long i;
  const NodeList* candidateElements;

  myPost2MainDimList.clear();
  myMain2PostDimList.clear();

  myRequiresIntegrations=0;

  myPost2MainDimList.push_back(0);

  list<bool> rawSpaceList;
  rawSpaceList.push_back(0);

  dimensionStruct nextDimension;
  tempGeometry.addDimension(nextDimension);

  // ************************************
  // find sampling element

  candidateElements = yourElement->getElementsByTagName("sampling",0);

  if(candidateElements->length()==0) {
    throw xmdsException(yourElement,"<sampling> element required");
  }

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <sampling> elements defined");
  }

  const Element* samplingElement = dynamic_cast<const Element*> (candidateElements->item(0));

  if(verbose()) {
    printf("Processing sampling element ...\n");
  }

  // ************************************
  // get sampling field geometery

  if(nDims>0) {

    // ************************************
    // find sampling space

    list<bool> samplingSpaceList;
    getAssignmentBools(samplingElement,"fourier_space",0,simulation()->field()->geometry()->nDims(),samplingSpaceList);

    if(samplingSpaceList.size() == 0) {
      printf("Sampling space for moment group #%li defaulting to x-space.\n",myGroupNumber+1);
      mySamplingSpace = 0;
    }
    else {
      mySamplingSpace = spaceList2ULong(samplingSpaceList);

      if(verbose()) {
	for(i=0;i<nDims;i++) {
	  if(samplingSpace(i)) {
	    printf("sampling will be performed with dimension #%li in fourier space\n",i+1);
	  }
	  else {
	    printf("sampling will be performed with dimension #%li in normal space\n",i+1);
	  }
	}
      }
    }

    // ************************************
    // find sampling lattice

    getAssignmentULongs(samplingElement,"lattice",1,nDims,mySamplingLatticeList);

    list<unsigned long>::const_iterator pULong = mySamplingLatticeList.begin();
    for(i=0;i<nDims;i++) {

      const unsigned long fieldLatticeI = simulation()->field()->geometry()->dimension(i)->lattice;
      const unsigned long latticeI=*pULong;

      if(verbose()) {
	printf("sampling lattice dimension #%li has N points = %li\n",i+1,*pULong);
      }

      if(latticeI>1) {

	nextDimension.name = simulation()->field()->geometry()->dimension(i)->name;
	nextDimension.lattice = latticeI;
	nextDimension.domain = simulation()->field()->geometry()->dimension(i)->domain;

	tempGeometry.addDimension(nextDimension);

	myPost2MainDimList.push_back(i);
	rawSpaceList.push_back(samplingSpace(i));

	if(fieldLatticeI - latticeI*(fieldLatticeI/latticeI)) {
	  sprintf(errorMessage(),"moments group sampling lattice dimension[%li] does not divide evenly into same field lattice",i+1);
	  throw xmdsException(samplingElement,errorMessage());
	}

	if(verbose()) {
	  printf("sampling lattice dimension #%li divides same field lattice by %li\n",i+1,fieldLatticeI/latticeI);
	}
      }
      else if(latticeI==1) {

	// cross section
	if(verbose()) {
	  printf("transverse dimension #%li will be cross-sectioned at ",i+1);
	  if(samplingSpace(i)) {
	    printf("k = 0\n");
	  }
	  else {
	    printf("x = domain centre\n");
	  }
	}
      }

      else {
	// integrate
	myRequiresIntegrations = 1;

	if(verbose()) {
	  printf("transverse dimension #%li will be integrated in ",i+1);
	  if(samplingSpace(i)) {
	    printf(" fourier space\n");
	  }
	  else {
	    printf(" normal space\n");
	  }
	}
      }

      myMain2PostDimList.push_back(myPost2MainDimList.size()-1);
      pULong++;
    }
  }

  myRawSpace=spaceList2ULong(rawSpaceList);

  // ************************************
  // find vectors

  getAssignmentStrings(samplingElement,"vectors",0,0,myVectorNamesList);

  if(myVectorNamesList.size()==0) {
    // no vectors specified, therefore assume using only main vector
    myVectorNamesList.push_back("main");
  }

  simulation()->field()->processVectors(myVectorNamesList,mySamplingSpace);

  // ************************************
  // find sampling moments

  getAssignmentStrings(samplingElement,"moments",1,0,mySamplingMomentsList);

  if(mySamplingMomentsList.size()==0) {
    throw xmdsException(samplingElement,"No sampling moments defined!");
  }
  if(verbose()) {
    for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
      printf("adding sampling moment '%s'\n",pXMLString->c_str());
    }
  }

  // ************************************
  // find sampling code

  mySamplingCode=*samplingElement->textContent(0);

  // check for non-white space charaters in code:

  if(mySamplingCode.isAllWhiteSpace()) {
    throw xmdsException(samplingElement,"No sampling code defined!");
  }

  if(verbose()) {
    printf("moment group sampling code loaded\n");
  }

  // ************************************
  // find post_propagation element

  candidateElements = yourElement->getElementsByTagName("post_propagation",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <post_propagation> elements defined");
  }

  else if(candidateElements->length()==1) {

    // there is a post_propagation element

    const Element* postPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

    if(verbose()) {
      printf("Processing post_propagation element ...\n");
    }

    // ************************************
    // find post space

    list<bool> aSpaceList;

    getAssignmentBools(postPropagationElement,"fourier_space",0,myPost2MainDimList.size(),aSpaceList);

    if(aSpaceList.size() == 0) {
      printf("In moment group #%li, using sampling space for space of remaining post propagation dimensions.\n",myGroupNumber+1);
      for(i=0;i<myPost2MainDimList.size();i++) {
	aSpaceList.push_back(samplingSpace(post2MainDim(i)));
      }
    }

    myPostSpace = spaceList2ULong(aSpaceList);

    if(verbose()) {
      for(i=0;i<myPost2MainDimList.size();i++) {
	if((myPostSpace >> i)&1) {
	  printf("post propagation will be performed with remaining dimension #%li in fourier space\n",i+1);
	}
	else {
	  printf("post propagation will be performed with remaining dimension #%li in normal space\n",i+1);
	}
      }
    }

    // ************************************
    // find post moments

    getAssignmentStrings(postPropagationElement,"moments",1,0,myPostMomentsList);

    if(myPostMomentsList.size()==0) {
      throw xmdsException(postPropagationElement,"No post propagation moments defined!");
    }
    if(verbose()) {
      for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
	printf("post propagation moment '%s' added\n",pXMLString->c_str());
      }
    }

    // ************************************
    // find post code

    myPostCode=*postPropagationElement->textContent(0);

    // check for non-white space charaters in code:

    if(myPostCode.isAllWhiteSpace()) {
      throw xmdsException(postPropagationElement,"No post propagation code defined!");
    }

    if(verbose()) {
      printf("moment group post propagation code loaded\n");
    }
  }
  else {
    // no post_propagation element
    // need to assign across raw space and moments

    myPostSpace = myRawSpace;
    myPostMomentsList=mySamplingMomentsList;
  }

  char tempName[256];
  sprintf(tempName,"mg%li",myGroupNumber);

  setName((XMLString) tempName);

  // add moment group vectors

  // determine if raw sample vector needs to be complex

  complexRawVector = 0;

  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {

    const xmdsVector* nextVector;

    simulation()->field()->getVector(*pXMLString,nextVector);

    complexRawVector |= (nextVector->vectorType()==COMPLEX);
  }

  complexRawVector |= (myPostSpace != myRawSpace);