C++ Molecule::componentid方法代码示例

void StatisticsWriter::initOutput(ParticleContainer* particleContainer,
		DomainDecompBase* domainDecomp, Domain* domain) {

#ifndef NDEBUG
	if (static_cast<const void*>(&_container) != static_cast<const void*>(particleContainer)) {
		global_log->error() << "VTKGridWriter works only with PlottableLinkCells!" << std::endl;
		exit(1);
	}
#endif

	std::stringstream fileNameStream;
	fileNameStream << _outputPrefix << "_LC" << ".stat";
	ofstream outfile(fileNameStream.str().c_str());

	outfile << "Number of Cells (including halo): [" << _container._cellsPerDimension[0]
	        << "," << _container._cellsPerDimension[1] << "," << _container._cellsPerDimension[2] << "]" << endl;

	outfile << "HaloWidth in num cells: [" << _container._haloWidthInNumCells[0]
		        << "," << _container._haloWidthInNumCells[1] << "," << _container._haloWidthInNumCells[2] << "]" << endl;

	outfile << "Cell Width: [" << _container._cellLength[0] << "," << _container._cellLength[1] << "," << _container._cellLength[2] << "]" << endl;

	outfile << endl;
	outfile << "SizeOf(Molecule) = " << sizeof(Molecule) << endl;

	std::vector<bool> count;

	Molecule* tmpMolecule = particleContainer->begin();
	while (tmpMolecule != particleContainer->end()) {
		if (tmpMolecule->componentid() >= count.size()) {
			count.resize(tmpMolecule->componentid()+1, false);
		}

		if (count[tmpMolecule->componentid()] == false) {
			outfile << "Molecule(cid=" << tmpMolecule->componentid() << ") " << tmpMolecule->totalMemsize() << endl;
			count[tmpMolecule->componentid()] = true;
		}
		tmpMolecule = particleContainer->next();
	}

	outfile.close();

}

unsigned long DomainDecompDummy::countMolecules(ParticleContainer* moleculeContainer, vector<unsigned long> &compCount) {
	for (unsigned i = 0; i < compCount.size(); i++) {
		compCount[i] = 0;
	}
	Molecule* tempMolecule;
	for (tempMolecule = moleculeContainer->begin(); tempMolecule != moleculeContainer->end(); tempMolecule = moleculeContainer->next()) {
		compCount[tempMolecule->componentid()] += 1;
	}
	int numMolecules = 0;
	for (unsigned i = 0; i < compCount.size(); i++) {
		numMolecules += compCount[i];
	}
	return numMolecules;
}

Molecule ChemicalPotential::loadMolecule()
{
			assert(!this->reservoir.empty());
			Molecule tmp = this->reservoir.front();
			this->reservoir.pop_front();
			if(this->reservoir.empty())
			{
	 tmp.scale_v(-1.0);
	 this->reservoir.push_back( tmp );
			}
			assert(tmp.componentid() == componentid);
#ifndef NDEBUG
			tmp.check(tmp.id());
#endif
			return tmp;
}

bool ChemicalPotential::getDeletion(TMoleculeContainer* cell, double* minco, double* maxco)
{
	 if(this->remainingDeletions.empty()) return false;
	 if(cell->getNumberOfParticles() == 0) return false;
	 
	 unsigned idx = *this->remainingDeletions.begin();
	 this->remainingDeletions.erase(this->remainingDeletions.begin());
	 double tminco[3];
	 double tmaxco[3];
	 if(restrictedControlVolume) for(int d=0; d < 3; d++)
	 {
	    tminco[d] = (minco[d] > control_bottom[d])? minco[d]
	                                              : control_bottom[d];
	    tmaxco[d] = (maxco[d] < control_top[d])? maxco[d]
	                                           : control_top[d];
	 }
	 else for(int d=0; d < 3; d++)
	 {
	    tminco[d] = minco[d];
	    tmaxco[d] = maxco[d];
	 }
	 
	 Molecule* m = cell->begin();
	 int j=0;
	 for(unsigned i=0; (i < idx); i++)
	 {
	    while(( (m->r(0) > tmaxco[0]) || (m->r(1) > tmaxco[1]) ||
	      (m->r(2) > tmaxco[2]) || (m->r(0) < tminco[0]) ||
	      (m->r(1) < tminco[1]) || (m->r(2) < tminco[2]) ||
							(m->componentid() != this->componentid) )
					 && (m != cell->end()))
			{
				 m = cell->next();
				 if(m == cell->end())
	 {
	    if(j == 0) return false;
	    m = cell->begin();
	    j = 0;
	 }
			}
			m = cell->next();
			j++;
			if(m == cell->end())
			{
	 if(j == 0) return false;
	 m = cell->begin();
	 j = 0;
			}
	 }
	 while( (m->r(0) > tmaxco[0]) || (m->r(1) > tmaxco[1]) ||
	        (m->r(2) > tmaxco[2]) || (m->r(0) < tminco[0]) ||
	  (m->r(1) < tminco[1]) || (m->r(2) < tminco[2]) ||
					(m->componentid() != this->componentid) )
	 {
	    m = cell->next();
	    if(m == cell->end())
	    {
	 if(j == 0) return false;
	 m = cell->begin();
			}
	 }
#ifndef NDEBUG
	 global_log->debug() << "ID " << m->id() << " selected for deletion (index " << idx << ")." << std::endl;
#endif
	 assert(m->id() < nextid);
	 return true;
}

/*
 * diese Version beschleunigt in alle Raumrichtungen
 */
void PressureGradient::determineAdditionalAcceleration
(
 DomainDecompBase* domainDecomp,
 ParticleContainer* molCont, double dtConstantAcc )
{
	for( map<unsigned int, double>::iterator uAAit = _universalAdditionalAcceleration[0].begin();
			uAAit != _universalAdditionalAcceleration[0].end();
			uAAit++ )
	{
		this->_localN[uAAit->first] = 0;
		for(unsigned short int d = 0; d < 3; d++)
			this->_localVelocitySum[d][uAAit->first] = 0.0;
	}
	for(Molecule* thismol = molCont->begin(); thismol != molCont->end(); thismol = molCont->next())
	{
		unsigned int cid = thismol->componentid();
		map<unsigned int, unsigned int>::iterator uCSIDit = this->_universalComponentSetID.find(cid);
		if(uCSIDit == _universalComponentSetID.end()) continue;
		unsigned cosetid = uCSIDit->second;
		this->_localN[cosetid]++;
		for(unsigned short int d = 0; d < 3; d++)
			this->_localVelocitySum[d][cosetid] += thismol->v(d);
	}

	// domainDecomp->collectCosetVelocity(&_localN, _localVelocitySum, &_globalN, _globalVelocitySum);
	//
	domainDecomp->collCommInit( 4 * _localN.size() );
	for( map<unsigned int, unsigned int long>::iterator lNit = _localN.begin(); lNit != _localN.end(); lNit++ )
	{
		domainDecomp->collCommAppendUnsLong(lNit->second);
		for(int d = 0; d < 3; d++)
			domainDecomp->collCommAppendDouble( this->_localVelocitySum[d][lNit->first]);
	}
	domainDecomp->collCommAllreduceSum();
	for( map<unsigned int, unsigned long>::iterator lNit = _localN.begin(); lNit != _localN.end(); lNit++ )
	{
		_globalN[lNit->first] = domainDecomp->collCommGetUnsLong();
		for(int d = 0; d < 3; d++)
			_globalVelocitySum[d][lNit->first] = domainDecomp->collCommGetDouble();
	}
	domainDecomp->collCommFinalize();

	map<unsigned int, long double>::iterator gVSit;
	if(!this->_localRank)
	{
		for(gVSit = _globalVelocitySum[0].begin(); gVSit != _globalVelocitySum[0].end(); gVSit++)
		{
#ifndef NDEBUG
			global_log->debug() << "required entries in velocity queue: " << _globalVelocityQueuelength[gVSit->first] << endl;
			global_log->debug() << "entries in velocity queue: " << _globalPriorVelocitySums[0][gVSit->first].size() << endl;
#endif
			for(unsigned short int d = 0; d < 3; d++)
			{
				while(_globalPriorVelocitySums[d][gVSit->first].size() < _globalVelocityQueuelength[gVSit->first])
					_globalPriorVelocitySums[d][gVSit->first].push_back(_globalVelocitySum[d][gVSit->first]);
			}
		}
	}

	if(!this->_localRank)
	{
		for(gVSit = _globalVelocitySum[0].begin(); gVSit != _globalVelocitySum[0].end(); gVSit++)
		{
			double invgN = 1.0 / this->_globalN[gVSit->first];
			double invgtau = 1.0 / this->_universalTau[gVSit->first];
			double invgtau2 = invgtau * invgtau;
			double previousVelocity[3];
			for(unsigned short int d = 0; d < 3; d++)
			{
				previousVelocity[d] = invgN * _globalPriorVelocitySums[d][gVSit->first].front();
				this->_globalPriorVelocitySums[d][gVSit->first].pop_front();
				this->_universalAdditionalAcceleration[d][gVSit->first]
					+= dtConstantAcc * invgtau2 *
					( this->_globalTargetVelocity[d][gVSit->first]
					  - 2.0*this->_globalVelocitySum[d][gVSit->first]*invgN
					  + previousVelocity[d] );
			}
#ifndef NDEBUG
			global_log->debug() << "accelerator no. " << gVSit->first 
				<< "previous vz: " << previousVelocity[2] 
				<< "current vz: " << _globalVelocitySum[2][gVSit->first]*invgN << endl;
#endif
		}
	}

	domainDecomp->collCommInit(7*this->_localN.size());
	for( map<unsigned int, unsigned long>::iterator lNit = _localN.begin();
			lNit != this->_localN.end();
			lNit++ )
	{
		domainDecomp->collCommAppendUnsLong(_globalN[lNit->first]);
		for(int d=0; d < 3; d++)
		{
			domainDecomp->collCommAppendDouble(
					this->_universalAdditionalAcceleration[d][lNit->first]
					);
			domainDecomp->collCommAppendDouble(
//.........这里部分代码省略.........

void PovWriter::doOutput(ParticleContainer* particleContainer,
                         DomainDecompBase* domainDecomp, Domain* domain,
                         unsigned long simstep, list<ChemicalPotential>* lmu) {
	if (simstep % _writeFrequency == 0) {

		stringstream filenamestream;
		if (_filenameisdate) {
			filenamestream << "mardyn" << gettimestring();
		}
		else {
			filenamestream << _filename;
		}

		if (_incremental) {
			filenamestream << "-";
			/* align file numbers with preceding '0's in the required range from 0 to _numberOfTimesteps. */
			int num_digits = ceil(log(double(_numberOfTimesteps / _writeFrequency)) / log(10.));
			filenamestream << aligned_number(simstep / _writeFrequency, num_digits, '0');
		}
		filenamestream << ".pov";

		ofstream ostrm(filenamestream.str().c_str());

		ostrm << "// " << filenamestream.str() << endl;
		ostrm << "// moldy" << endl;
		time_t now;
		now = time(NULL);
		ostrm << "// " << ctime(&now) << endl;

		ostrm << "// bb: [0," << domain->getGlobalLength(0) << "]^3" << endl;
		ostrm << "//*PMRawBegin" << endl;
		ostrm << "background {rgb <1,1,1>}" << endl;
		ostrm << "//*PMRawEnd" << endl;
		vector<Component> dcomponents = domain->getComponents();
		for (unsigned int i = 0; i < dcomponents.size(); ++i) {
			ostringstream osstrm;
			osstrm.clear();
			osstrm.str("");
			osstrm << " pigment {color rgb <" << (i + 1) % 2 << "," << (i + 1) / 2 % 2 << "," << (i + 1) / 4 % 2 << ">}";
			osstrm << " finish{ambient 0.5 diffuse 0.4 phong 0.3 phong_size 3}";
			ostrm << "#declare T" << i << " = ";
			//ostrm << "sphere {<0,0,0>,0.5 pigment {color rgb<1,0,0>} finish{ambient 0.5 diffuse 0.4 phong 0.3 phong_size 3} scale 1.}";
			dcomponents.at(i).writePOVobjs(ostrm, osstrm.str());
			ostrm << endl;
		}
		ostrm << endl;
		ostrm << "camera { perspective" << endl;
		float xloc = -.1 * domain->getGlobalLength(0);
		float yloc = 1.1 * domain->getGlobalLength(1);
		float zloc = -1.5 * domain->getGlobalLength(2);
		ostrm << " location <" << xloc << ", " << yloc << ", " << zloc << ">" << endl;
		//ostrm << " direction <0, 0, 1>" << endl;
		//ostrm << " right <1.33333, 0, 0>" << endl;
		//ostrm << " up <0, 1, 0>" << endl;
		//ostrm << " sky <0, 1, 0>" << endl;
		ostrm << " look_at <" << .5 * domain->getGlobalLength(0) << ", " << .5 * domain->getGlobalLength(1) << ", " << .5 * domain->getGlobalLength(2) << ">" << endl;
		ostrm << "}" << endl;
		ostrm << endl;
		ostrm << "light_source { <" << xloc << ", " << yloc << ", " << zloc << ">, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <0,0,0>, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <0,0," << domain->getGlobalLength(2) << ">, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <0," << domain->getGlobalLength(1) << ",0>, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <0," << domain->getGlobalLength(1) << "," << domain->getGlobalLength(2) << ">, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <" << domain->getGlobalLength(0) << ",0,0>, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <" << domain->getGlobalLength(0) << ",0," << domain->getGlobalLength(2) << ">, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <" << domain->getGlobalLength(0) << "," << domain->getGlobalLength(1) << ",0>, color rgb <1,1,1> }" << endl;
		ostrm << "light_source { <" << domain->getGlobalLength(0) << "," << domain->getGlobalLength(1) << "," << domain->getGlobalLength(2) << ">, color rgb <1,1,1> }" << endl;
		ostrm << endl;
		ostrm << "// " << dcomponents.size() << " objects for the atoms following..." << endl;
		double mrot[3][3];
		for (Molecule* pos = particleContainer->begin(); pos != particleContainer->end(); pos = particleContainer->next()) {
			(pos->q()).getRotinvMatrix(mrot);
			//cout << "object { T0 rotate <0,0,0> translate <0,0,0>}" << endl;
			ostrm << "object { T" << pos->componentid();
			ostrm << " matrix <"
			      << mrot[0][0] << "," << mrot[0][1] << "," << mrot[0][2] << ","
			      << mrot[1][0] << "," << mrot[1][1] << "," << mrot[1][2] << ","
			      << mrot[2][0] << "," << mrot[2][1] << "," << mrot[2][2] << ","
			      << pos->r(0) << "," << pos->r(1) << "," << pos->r(2)
			      << ">";
			ostrm << "}" << endl;
		}

		// RK
		/* map cluster ID to color */
		//if ((pos->clusterid() != -1) && (m_clusters.find(pos->clusterid())->second > MINCLUSTERSIZE)) {
		//	ostrm << "  pigment {color rgb<" <<
		//			(0.25*(pos->clusterid()%5)) << "," <<
		//			(0.25*((pos->clusterid()/5)%5)) << "," <<
		//			(0.25*((pos->clusterid()/125)%5)) << "," <<
		//			">}";
		//}
		//else {
		//	ostrm << "  pigment {color rgb<0.9,0.9,0.9>}";
		//}
		// /RK

		ostrm.close();
	}
}

void DomainDecompDummy::exchangeMolecules(ParticleContainer* moleculeContainer, const vector<Component>& components, Domain* domain) {

	double rmin[3]; // lower corner of the process-specific domain //ENABLE_MPI
	double rmax[3];
	double halo_L[3]; // width of the halo strip //ENABLE_MPI
	for (int i = 0; i < 3; i++) {
		rmin[i] = moleculeContainer->getBoundingBoxMin(i);
		rmax[i] = moleculeContainer->getBoundingBoxMax(i);
		halo_L[i] = moleculeContainer->get_halo_L(i);
	}

	Molecule* currentMolecule;
	// molecules that have to be copied (because of halo), get a new position
	double new_position[3];

	double phaseSpaceSize[3];

	double low_limit; // particles below this limit have to be copied or moved to the lower process
	double high_limit; // particles above(or equal) this limit have to be copied or moved to the higher process

	for (unsigned short d = 0; d < 3; ++d) {
		phaseSpaceSize[d] = rmax[d] - rmin[d];
		// set limits (outside "inner" region)
		low_limit = rmin[d] + halo_L[d];
		high_limit = rmax[d] - halo_L[d];
		currentMolecule = moleculeContainer->begin();

		//cout << "low_limit: " << low_limit << " / high_limit: " << high_limit << endl;
		//cout << "halo_L: " << halo_L[0] << " / " << halo_L[1] << " / " << halo_L[2] << endl;
		//cout << "proc_domain_L: " << proc_domain_L[0] << " / " << proc_domain_L[1] << " / " << proc_domain_L[2] << endl;
		while (currentMolecule != moleculeContainer->end()) {
			const double& rd = currentMolecule->r(d);
			if (rd < low_limit) {
#ifdef MADPAC_LBMD_COUPLING
				if(madpac::PeanoCouplingService::getInstance().bEnabled() && d < madpac::PeanoCouplingService::getInstance().getDim())
				{
					madpac::PeanoCouplingService::getInstance().getDataContainer()->incNumInsertMolecules();
				}
				else
				{
#endif
				// determine the position for the copy of the molecule
				for (unsigned short d2 = 0; d2 < 3; d2++) {
					// when moving parallel to the coordinate d2 to another process, the
					// local coordinates in d2 change
					if (d2 == d)
						new_position[d2] = rd + phaseSpaceSize[d2];
					else
						new_position[d2] = currentMolecule->r(d2);
				}
				Molecule m1 = Molecule(currentMolecule->id(),currentMolecule->componentid(),
				                       new_position[0], new_position[1], new_position[2],
				                       currentMolecule->v(0),currentMolecule->v(1),currentMolecule->v(2),
				                       currentMolecule->q().qw(),currentMolecule->q().qx(),currentMolecule->q().qy(),currentMolecule->q().qz(),
				                       currentMolecule->D(0),currentMolecule->D(1),currentMolecule->D(2), &components);
				moleculeContainer->addParticle(m1);
				currentMolecule = moleculeContainer->next();
#ifdef MADPAC_LBMD_COUPLING
				}
#endif
			}
			else if (rd >= high_limit) {
#ifdef MADPAC_LBMD_COUPLING
				if(madpac::PeanoCouplingService::getInstance().bEnabled() && d < madpac::PeanoCouplingService::getInstance().getDim())
				{
					madpac::PeanoCouplingService::getInstance().getDataContainer()->incNumInsertMolecules();
				}
				else
				{
#endif
				// determine the position for the copy of the molecule
				for (unsigned short d2 = 0; d2 < 3; d2++) {
					// when moving parallel to the coordinate d2 to another process, the
					// local coordinates in d2 change
					if (d2 == d)
						new_position[d2] = rd - phaseSpaceSize[d2];
					else
						new_position[d2] = currentMolecule->r(d2);
				}
				Molecule m1 = Molecule(currentMolecule->id(),currentMolecule->componentid(),
				                       new_position[0], new_position[1], new_position[2],
				                       currentMolecule->v(0),currentMolecule->v(1),currentMolecule->v(2),
				                       currentMolecule->q().qw(),currentMolecule->q().qx(),currentMolecule->q().qy(),currentMolecule->q().qz(),
				                       currentMolecule->D(0),currentMolecule->D(1),currentMolecule->D(2), &components);
				moleculeContainer->addParticle(m1);
				currentMolecule = moleculeContainer->next();
#ifdef MADPAC_LBMD_COUPLING
				}
#endif
			}
			else
				currentMolecule = moleculeContainer->next();
		}
	}
}