xmdsmomentgroup.cc

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

	    fprintf(outfile,"%s += _i%li*_%s_dx%li;\n",dimName,i,fieldName,i);

	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _i%li*_%s_%s_ncomponents",
		      fieldName,pXMLString->c_str(),i,fieldName,pXMLString->c_str());
	      for(j=i+1;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,";\n");
	    }
	  }
	}
        else {
	  // normal sampling
	  for(j=0;j<i;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,main2PostDim(i),i);
	  fprintf(outfile,"\n");
	}

        psamplingLatticeI++;
      }
    }
    
  fprintf(outfile,"// ************ Sampling code! ******************\n");
  fprintf(outfile,"%s\n",mySamplingCode.c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  if(myRequiresIntegrations|swapped) {
    for(j=0;j<nDims;j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"unsigned long _%s_raw_index_pointer=(_%s_sample_pointer",mgFieldName,mgFieldName);
    for(i=1;i<myNDims;i++) {
      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
    }
    for(i=1;i<myNDims;i++) {
      fprintf(outfile," + _i%li",post2MainDim(i));
      for(j=i+1;j<myNDims;j++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,j);
      }
    }
    fprintf(outfile,")*_%s_raw_ncomponents;\n",mgFieldName);
    fprintf(outfile,"\n");
  }

  i=0;
  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    for(j=0;j<nDims;j++) {
	fprintf(outfile,"	");
    }
    if(myRequiresIntegrations) {
      fprintf(outfile,"_%s_raw[_%s_raw_index_pointer+%li] += %s",mgFieldName,mgFieldName,i,pXMLString->c_str());
      j=0;
      psamplingLatticeI = mySamplingLatticeList.begin();
      while(j<nDims) {
	if(*psamplingLatticeI==0) {
	  if(samplingSpace(j)) {
	    fprintf(outfile,"*_%s_dk%li",fieldName,j);
	  }
	  else {
	    fprintf(outfile,"*_%s_dx%li",fieldName,j);
	  }
	}
	psamplingLatticeI++;
	j++;
      }
      fprintf(outfile,";\n");
    }
    else
      fprintf(outfile,"_%s_raw[_%s_raw_index_pointer+%li] = %s;\n",mgFieldName,mgFieldName,i,pXMLString->c_str());
    i++;
  }

  // close nested loops

  if(!myRequiresIntegrations&!swapped) {
    if(nDims>0) {
      for(i=0;i<nDims;i++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_raw_index_pointer += _%s_raw_ncomponents;\n",mgFieldName,mgFieldName);
      fprintf(outfile,"\n");
    }
  }

  list<unsigned long>::const_reverse_iterator psamplingLatticeIr = mySamplingLatticeList.rbegin();

  if(nDims>0) {
    if(swapped) {
        for(i=nDims;i>2;i--) {

	  const char* dimName = simulation()->field()->geometry()->dimension(i-1)->name.c_str();

	  fprintf(outfile,"\n");
	  if(*psamplingLatticeIr==0) {
	    // integration
	    if(i==nDims) {
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	    }
	    if(samplingSpace(i-1)) {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li)\n",dimName,fieldName,i-1,fieldName,i-1);	
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,fieldName,i-1,fieldName,i-1);
	    }
	    else {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,fieldName,i-1);
	    }
	  }
	  else if(*psamplingLatticeIr==1) {
	    // cross-section
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li",fieldName,pXMLString->c_str(),fieldName,i-1);
	      if(!samplingSpace(i-1)) {
		fprintf(outfile,"-_i%li",i-1);
	      }
	      if(!(i==nDims)) {
		fprintf(outfile,"-1");
	      }
	      fprintf(outfile,")");
	      for(j=i;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	    }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(i-1)) {
	      // narrow k-space window
	      if(i==nDims) {
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"\n");
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li) {\n",
		      dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice%li-_%s_lattice%li)",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	}\n");
	    }
	    else {
	      // coarse x-space sampling
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li/_%s_lattice%li",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		if(!(i==nDims)) {
		  fprintf(outfile,"-1");
		}
		fprintf(outfile,")");
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,mgFieldName,main2PostDim(i-1));
	    }
	  }

	  if(!(*psamplingLatticeIr==1)) {
	    // if not cross-section
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"}\n");
	  }

	  psamplingLatticeIr++;
        }

        psamplingLatticeIr++;
        
        const char* dimName = simulation()->field()->geometry()->dimension(0)->name.c_str();

        fprintf(outfile,"\n");
        if(*psamplingLatticeIr==0) {
	  // integration
	  if(2==nDims) {
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<2;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
		      fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
	  }
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk0;\n",dimName,fieldName);
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",dimName,fieldName,fieldName);	
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",dimName,fieldName,fieldName);
	}
        else if(*psamplingLatticeIr==1) {
	  // cross-section
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    for(j=0;j<2;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice0",fieldName,pXMLString->c_str(),fieldName);
	    if(!(2==nDims)) {
	      fprintf(outfile,"-1");
	    }
	    fprintf(outfile,")");
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	}
        else {
	  // normal sampling
	  // narrow k-space window
	  if(2==nDims) {
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<2;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
		      fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
	  }
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk0;\n",dimName,fieldName);
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk0) {\n",
		  dimName,mgFieldName,main2PostDim(0),fieldName);
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    for(j=0;j<2;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice0-_%s_lattice%li)",
		    fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(0));
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk0;\n",dimName,mgFieldName,main2PostDim(0),fieldName);
	  for(j=0;j<2;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	}\n");
	}

        if(!(*psamplingLatticeIr==1)) {
	  // if not cross-section
	  for(j=0;j<2;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"}\n");
	}

        psamplingLatticeIr--;
                
        fprintf(outfile,"\n");
        if(*psamplingLatticeIr==0) {
	  // integration
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);	
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
	}
        else if(*psamplingLatticeIr==1) {
	  // cross-section
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    fprintf(outfile,"	");
	    fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice1",fieldName,pXMLString->c_str(),fieldName);
	    fprintf(outfile,"-1)");
	    fprintf(outfile,"*_%s_lattice0",fieldName);
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	}
        else {
	  // normal sampling
	  // narrow k-space window
	  fprintf(outfile,"	");
	  fprintf(outfile,"k%s +=	_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);