Implement test on overlapping speheres with tolerance (b06ba108) · Commits · Giacomo Mulas / NP_TMcode

src/include/Configuration.h

+40 −2

Original line number	Original line	Diff line number	Diff line
	@@ -117,6 +117,8 @@ protected:
	double _host_ram_gb;		double _host_ram_gb;
	//! \brief GPU RAM in GB		//! \brief GPU RAM in GB
	double _gpu_ram_gb;		double _gpu_ram_gb;
			//! \brief Sphere overlap tolerance;
			double _tolerance;

	public:		public:
	//! \brief Read-only view on number of spherical components.		//! \brief Read-only view on number of spherical components.
	@@ -173,6 +175,8 @@ public:
	const double& host_ram_gb = _host_ram_gb;		const double& host_ram_gb = _host_ram_gb;
	//! \brief Read-only view on GPU RAM in GB		//! \brief Read-only view on GPU RAM in GB
	const double& gpu_ram_gb = _gpu_ram_gb;		const double& gpu_ram_gb = _gpu_ram_gb;
			//! \brief Read-only view on sphere overlap tolerance
			const double& tolerance = _tolerance;

	/*! \brief Build a scattering geometry configuration structure.		/*! \brief Build a scattering geometry configuration structure.
	*		*
	@@ -533,6 +537,12 @@ public:
	*/		*/
	double get_max_radius();		double get_max_radius();

			/*! \brief Get the minimum radius of the sphere components.
			*
			* \return radius: `double` The radius of the smallest sphere.
			*/
			double get_min_radius();

	/*! \brief Get the number of layers for a given configuration.		/*! \brief Get the number of layers for a given configuration.
	*		*
	* This is a specialized function to get the number of layers in a specific		* This is a specialized function to get the number of layers in a specific
	@@ -558,7 +568,7 @@ public:
	* \param index: `int` Index of the ID to be retrieved.		* \param index: `int` Index of the ID to be retrieved.
	* \return radius: `double` The requested sphere radius.		* \return radius: `double` The requested sphere radius.
	*/		*/
	double get_radius(int index) { return _radii_of_spheres[index]; }		double get_radius(int index);

	/*! \brief Get the value of a scale by its index.		/*! \brief Get the value of a scale by its index.
	*		*
	@@ -581,6 +591,16 @@ public:
	*/		*/
	double get_scale(int index) { return _scale_vec[index]; }		double get_scale(int index) { return _scale_vec[index]; }

			/*! \brief Get the radius of a sphere type by its index.
			*
			* This is a specialized function to get the radius of a sphere type through its
			* index. Sphere types range from 0 to NUM_CONFIGURATIONS - 1.
			*
			* \param index: `int` Index of the ID to be retrieved.
			* \return radius: `double` The requested sphere radius.
			*/
			double get_type_radius(int index) { return _radii_of_spheres[index]; }

	/*! \brief Print the contents of the configuration object to terminal.		/*! \brief Print the contents of the configuration object to terminal.
	*		*
	* In case of quick debug testing, `ScattererConfiguration.print()` allows printing		* In case of quick debug testing, `ScattererConfiguration.print()` allows printing
	@@ -630,4 +650,22 @@ public:

	};		};

	#endif		/*! \brief Check the presence of overlapping spheres in an aggregate.
			*
			* The theoretical implementation of the code is based on the assumption that the
			* spherical monomers that compose the aggregate do not overlap in space. Small
			* violations of this assumption do not critically affect the results, but they
			* need to be reported. This function performs a scan of the aggregate and returns
			* a summary of any potential overlap, if detected.
			*
			* \param sconf: `ScattererConfiguration *` Pointer to a ScattererConfiguration instance.
			* \param gconf: `GeometryConfiguration *` Pointer to a GeometryConfiguration instance.
			* \param tolerance: `const double` A tolerance value below which overlap is ignored.
			* \return overlaps: `int *` A vector containing the number of overlaps and a list of
			* overlapping sphere indices.
			*/
			int *check_overlaps(
			ScattererConfiguration sconf, GeometryConfiguration gconf, const double tolerance=0.0
			);

			#endif // INCLUDE_CONFIGURATION_H_

src/libnptm/Commons.cpp

+1 −1

Original line number	Original line	Diff line number	Diff line
	@@ -133,7 +133,7 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
	ros = new double[num_configurations]();		ros = new double[num_configurations]();
	nshl = new int[num_configurations]();		nshl = new int[num_configurations]();
	for (int ci = 0; ci < num_configurations; ci++) {		for (int ci = 0; ci < num_configurations; ci++) {
	ros[ci] = sconf->get_radius(ci);		ros[ci] = sconf->get_type_radius(ci);
	nshl[ci] = sconf->get_nshl(ci);		nshl[ci] = sconf->get_nshl(ci);
	}		}
	_npnt = gconf->npnt;		_npnt = gconf->npnt;

src/libnptm/Configuration.cpp

+84 −5

Original line number	Original line	Diff line number	Diff line
	@@ -101,6 +101,7 @@ GeometryConfiguration::GeometryConfiguration(
	_dyn_order_flag = 1;		_dyn_order_flag = 1;
	_host_ram_gb = 0.0;		_host_ram_gb = 0.0;
	_gpu_ram_gb = 0.0;		_gpu_ram_gb = 0.0;
			_tolerance = -1.0;
	}		}

	GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)		GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)
	@@ -140,6 +141,7 @@ GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)
	_dyn_order_flag = rhs._dyn_order_flag;		_dyn_order_flag = rhs._dyn_order_flag;
	_host_ram_gb = rhs._host_ram_gb;		_host_ram_gb = rhs._host_ram_gb;
	_gpu_ram_gb = rhs._gpu_ram_gb;		_gpu_ram_gb = rhs._gpu_ram_gb;
			_tolerance = rhs._tolerance;
	}		}

	#ifdef MPI_VERSION		#ifdef MPI_VERSION
	@@ -179,6 +181,7 @@ GeometryConfiguration::GeometryConfiguration(const mixMPI *mpidata) {
	MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD);		MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_host_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_host_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_gpu_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_gpu_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
			MPI_Bcast(&_tolerance, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	}		}

	void GeometryConfiguration::mpibcast(const mixMPI *mpidata) {		void GeometryConfiguration::mpibcast(const mixMPI *mpidata) {
	@@ -214,6 +217,7 @@ void GeometryConfiguration::mpibcast(const mixMPI *mpidata) {
	MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD);		MPI_Bcast(&_dyn_order_flag, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_host_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_host_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&_gpu_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_gpu_ram_gb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
			MPI_Bcast(&_tolerance, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	}		}
	#endif		#endif

	@@ -374,6 +378,12 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(const std::string& fil
	conf->_gpu_ram_gb = ram_gb;		conf->_gpu_ram_gb = ram_gb;
	}		}
	}		}
			if (str_target.size() > 10) {
			if (str_target.substr(0, 10).compare("TOLERANCE=") == 0) {
			double tolerance = (double)stod(str_target.substr(10, str_target.length()));
			conf->_tolerance = tolerance;
			}
			}
	}		}

	// Clean up memory and return configuration object.		// Clean up memory and return configuration object.
	@@ -439,7 +449,7 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs
	_xip = rhs._xip;		_xip = rhs._xip;
	_max_layers = rhs._max_layers;		_max_layers = rhs._max_layers;
	_iog_vec = new int[_number_of_spheres]();		_iog_vec = new int[_number_of_spheres]();
	_radii_of_spheres = new double[_number_of_spheres]();		_radii_of_spheres = new double[_configurations]();
	_nshl_vec = new int[_configurations]();		_nshl_vec = new int[_configurations]();
	_rcf = new double*[_configurations];		_rcf = new double*[_configurations];
	_scale_vec = new double[_number_of_scales]();		_scale_vec = new double[_number_of_scales]();
	@@ -1019,6 +1029,15 @@ double ScattererConfiguration::get_max_radius() {
	return result;		return result;
	}		}

			double ScattererConfiguration::get_min_radius() {
			double result = -1.0;
			for (int ci = 0; ci < _configurations; ci++) {
			if (_radii_of_spheres[ci] < result or result <= 0.0)
			result = _radii_of_spheres[ci];
			}
			return result;
			}

	double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {		double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {
	double result = 0.0;		double result = 0.0;
	if (_use_external_sphere) {		if (_use_external_sphere) {
	@@ -1036,14 +1055,15 @@ double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {
	avgX /= _number_of_spheres;		avgX /= _number_of_spheres;
	avgY /= _number_of_spheres;		avgY /= _number_of_spheres;
	avgZ /= _number_of_spheres;		avgZ /= _number_of_spheres;
	int configuration = 0;		// int configuration = 0;
	for (int si = 0; si < _number_of_spheres; si++) {		for (int si = 0; si < _number_of_spheres; si++) {
	int iogi = _iog_vec[si];		// int iogi = _iog_vec[si];
	if (iogi >= si + 1) configuration++;		// if (iogi >= si + 1) configuration++;
	double dX = gc->get_sph_x(si) - avgX;		double dX = gc->get_sph_x(si) - avgX;
	double dY = gc->get_sph_y(si) - avgY;		double dY = gc->get_sph_y(si) - avgY;
	double dZ = gc->get_sph_z(si) - avgZ;		double dZ = gc->get_sph_z(si) - avgZ;
	double radius = _radii_of_spheres[configuration - 1];		// double radius = _radii_of_spheres[configuration - 1];
			double radius = get_radius(si);
	dist2 = dX * dX + dY * dY + dZ * dZ + radius * radius;		dist2 = dX * dX + dY * dY + dZ * dZ + radius * radius;
	if (dist2 > max_dist2) {		if (dist2 > max_dist2) {
	max_dist2 = dist2;		max_dist2 = dist2;
	@@ -1055,6 +1075,12 @@ double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {
	return result;		return result;
	}		}

			double ScattererConfiguration::get_radius(int index) {
			int type_id = _iog_vec[index];
			double radius = _radii_of_spheres[type_id - 1];
			return radius;
			}

	void ScattererConfiguration::print() {		void ScattererConfiguration::print() {
	int ies = (_use_external_sphere)? 1 : 0;		int ies = (_use_external_sphere)? 1 : 0;
	printf("### CONFIGURATION DATA ###\n");		printf("### CONFIGURATION DATA ###\n");
	@@ -1516,3 +1542,56 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) {
	}		}
	return true;		return true;
	}		}

			int *check_overlaps(
			ScattererConfiguration sconf, GeometryConfiguration gconf,
			const double tolerance
			) {
			const int nsphs = sconf->number_of_spheres;
			const int nsphg = gconf->number_of_spheres;
			const int nn = nsphs * nsphs;
			if (nsphs != nsphg)
			throw(UnrecognizedConfigurationException("ERROR: mismatch in number of spheres!"));
			int num_overlaps = 0;
			int *index_matrix = new int[nn]();
			#pragma omp parallel for reduction(+: num_overlaps)
			for (int si = 0; si < nn; si++) {
			const int i0 = si / nsphs;
			const int i1 = si % nsphs;
			if (i0 != i1) {
			const double x0 = gconf->get_sph_x(i0);
			const double y0 = gconf->get_sph_y(i0);
			const double z0 = gconf->get_sph_z(i0);
			const double x1 = gconf->get_sph_x(i1);
			const double y1 = gconf->get_sph_y(i1);
			const double z1 = gconf->get_sph_z(i1);
			const double r0 = sconf->get_radius(i0);
			const double r1 = sconf->get_radius(i1);
			const double r_sum = (r0 + r1);
			const double dist2 = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)
			+ (z1 - z0) * (z1 - z0);
			const double dist = sqrt(dist2);
			const double min_dist = r0 + r1 - tolerance;
			if (dist < min_dist) {
			num_overlaps++;
			index_matrix[nsphs * i0 + i1] = 1;
			}
			}
			}
			int *vec_overlaps = new int[1 + num_overlaps];
			vec_overlaps[0] = num_overlaps / 2; // we do not count symmetric overlaps
			int last_index = 1;
			for (int si = 0; si < nn; si++) {
			if (index_matrix[si] != 0) {
			const int i0 = si / nsphs;
			const int i1 = si % nsphs;
			if (i0 != i1) {
			vec_overlaps[last_index++] = i0 + 1;
			vec_overlaps[last_index++] = i1 + 1;
			if (last_index > num_overlaps) break; // for loop
			}
			}
			}
			delete[] index_matrix;
			return vec_overlaps;
			}