xmdsmomentgroup.cc

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

	    fprintf(outfile,"\n");
	  }
	  else if(firstlattice==1) {
	    // cross-section in this dimension
	    fprintf(outfile,"	");
	    fprintf(outfile,"unsigned long _i0 = 0;\n");
	  }
	  else {
	    // normal sampling
	    fprintf(outfile,"	");
	    fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_lattice%li;_i0++) {\n",mgFieldName,main2PostDim(0));
	    fprintf(outfile,"\n");
	  }
                

	  for(i=2;i<nDims;i++) {
    
            const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();
    
            for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
            if(samplingSpace(i)) {
	      fprintf(outfile,"double k%s = 0;\n",dimName);
	    }
            else {
	      fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	    }
            fprintf(outfile,"\n");
    
            if(*psamplingLatticeI==0) {
	      // integration over this dimension
	      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,fieldName,i,i);
	      fprintf(outfile,"\n");
	    }
            else if(*psamplingLatticeI==1) {
	      // cross-section in this dimension
	      if(samplingSpace(i)) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	      }
	      else {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		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++;
	  }
        }
      else
        {
	  // Make the first one smaller.  Then loop.
	  // Add the k-space possibility to this first loop everywhere.

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

	  if(*psamplingLatticeI==0) {
	    // integration over this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
                
	    if(samplingSpace(0)) {
		fprintf(outfile,"double k%s = local_x_start*_%s_dk0;\n",dimName,fieldName);
		fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",dimName,fieldName,fieldName);	
		fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",dimName,fieldName,fieldName);
	      }
	    else {
	      fprintf(outfile,"double %s = _%s_xmin0 + local_x_start*_%s_dx0;\n",dimName,fieldName,fieldName);
	    }
                                    
	    fprintf(outfile,"\n");

	    fprintf(outfile,"for(long _i0=0; _i0<local_nx; _i0++) {\n");
	    fprintf(outfile,"\n");
	  }
	  else if(*psamplingLatticeI==1) {
	    // cross-section in this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    if(!myRequiresIntegrations) {
		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);
		}
		fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);
		fprintf(outfile,"\n");
	      }

	    fprintf(outfile,"\n");

	    if(samplingSpace(0)) {
		fprintf(outfile,"double k%s = 0;\n",dimName);
		fprintf(outfile,"unsigned long _i0 = 0;\n");                    
		fprintf(outfile,"if(local_x_start==0) {\n");                    
	      }
	    else
	      {
		fprintf(outfile,"double %s = _%s_xmin0;\n",dimName,fieldName);
		fprintf(outfile,"\n");

		fprintf(outfile,"long _i0=_%s_lattice0/2-local_x_start;\n",fieldName);
    
		fprintf(outfile,"if((_i0+1>0)&(_i0<local_nx)) {\n");
    
		fprintf(outfile,"  %s += (_i0+local_x_start)*_%s_dx0;\n",dimName,fieldName);
    
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"  _%s_%s_index_pointer += _i0*_%s_%s_ncomponents",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		  for(j=1;j<nDims;j++) {
		    fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		  }
		  fprintf(outfile,";\n\n");
		}
	      }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(0)) {
		fprintf(outfile,"long first_x_pointer, last_x_pointer;\n");
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"long _%s_%s_index_pointer;\n",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"double k%s;\n\n",dimName);
                    
		fprintf(outfile,"if(local_x_start<(_%s_lattice%li-1)/2+1) {\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   first_x_pointer=local_x_start;\n");
		fprintf(outfile,"   k%s=local_x_start*_%s_dk0;\n",dimName,fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"   }\n");
		fprintf(outfile,"else if(local_x_start>(_%s_lattice%li-1)/2+_%s_lattice0-_%s_lattice%li) {\n"
			,mgFieldName,main2PostDim(0),fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   first_x_pointer=local_x_start-_%s_lattice0+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   k%s=(local_x_start-_%s_lattice0)*_%s_dk0;\n",dimName,fieldName,fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"   }\n");
		fprintf(outfile,"else {\n");
		fprintf(outfile,"   first_x_pointer=(_%s_lattice%li-1)/2+1;\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   k%s=((_%s_lattice%li-1)/2-_%s_lattice%li+1)*_%s_dk0;\n",dimName,mgFieldName,main2PostDim(0),mgFieldName,main2PostDim(0),fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		    fprintf(outfile,"   _%s_%s_index_pointer=(first_x_pointer+_%s_lattice0-_%s_lattice%li-local_x_start)",fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(0));
		    for(i=1;i<nDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		    }
		    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
		  }
		fprintf(outfile,"   }\n\n");
                    
		fprintf(outfile,"if(local_x_start+local_nx-1<(_%s_lattice%li-1)/2+1)\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   last_x_pointer=local_x_start+local_nx-1;\n");
		fprintf(outfile,"else if(local_x_start+local_nx-1>(_%s_lattice%li-1)/2+_%s_lattice0-_%s_lattice%li)\n"
			,mgFieldName,main2PostDim(0),fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   last_x_pointer=local_x_start+local_nx-1-_%s_lattice0+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"else \n");
		fprintf(outfile,"   last_x_pointer=(_%s_lattice%li-1)/2;\n\n",mgFieldName,main2PostDim(0));
 
		if(!myRequiresIntegrations) {
		    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);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents + first_x_pointer",mgFieldName);
		    for(i=2;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);            
                
		    fprintf(outfile,"\n");
		  }
                       
		fprintf(outfile,"for(long _i0=first_x_pointer;_i0<last_x_pointer+1;_i0++) {\n");
		fprintf(outfile,"\n");
	      }
	    else
	      {
		fprintf(outfile,"long first_x_pointer=_%s_lattice1-(_%s_lattice1*(_%s_lattice0-local_x_start))/_%s_lattice0;\n"
			,mgFieldName,mgFieldName,fieldName,fieldName);
		fprintf(outfile,"long last_x_pointer=(_%s_lattice1*(local_x_start+local_nx-1))/_%s_lattice0;\n\n",mgFieldName,fieldName);
 
		if(!myRequiresIntegrations) {
		    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);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents + first_x_pointer",mgFieldName);
		    for(i=2;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);            
                
		    fprintf(outfile,"\n");
		  }
                       
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		    fprintf(outfile,"unsigned long _%s_%s_index_pointer=(first_x_pointer*(_%s_lattice0/_%s_lattice1)-local_x_start)"
                            ,fieldName,pXMLString->c_str(),fieldName,mgFieldName);
		    for(i=1;i<nDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		    }
		    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
		  }
		fprintf(outfile,"\n");
    
		fprintf(outfile,"double %s = _%s_xmin0 + first_x_pointer*_%s_dx1;\n",dimName,fieldName,mgFieldName);
		fprintf(outfile,"\n");
            
		fprintf(outfile,"for(long _i0=first_x_pointer;_i0<last_x_pointer+1;_i0++) {\n");
		fprintf(outfile,"\n");
	      }
        
	  }

	  psamplingLatticeI++;

	  for(i=1;i<nDims;i++) {
    
            const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();
    
            for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
            if(samplingSpace(i)) {
	      fprintf(outfile,"double k%s = 0;\n",dimName);
	    }
            else {
	      fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	    }
            fprintf(outfile,"\n");
    
            if(*psamplingLatticeI==0) {
	      // integration over this dimension
	      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,fieldName,i,i);
	      fprintf(outfile,"\n");
	    }
            else if(*psamplingLatticeI==1) {
	      // cross-section in this dimension
	      if(samplingSpace(i)) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	      }
	      else {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		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++;
	  }
        }
    } 
  else 
    {
      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
      }
      fprintf(outfile,"\n");
    
      for(i=0;i<nDims;i++) {
    
        const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();

        for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
        if(samplingSpace(i)) {
	  fprintf(outfile,"double k%s = 0;\n",dimName);
	}
        else {
	  fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	}
        fprintf(outfile,"\n");

        if(*psamplingLatticeI==0) {
	  // integration over this dimension
	  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,fieldName,i,i);
	  fprintf(outfile,"\n");
	}
        else if(*psamplingLatticeI==1) {
	  // cross-section in this dimension
	  if(samplingSpace(i)) {
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	  }
	  else {
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }