Merge branch 'mesh_improvements' into 'master' (d610abf2) · Commits · XRFLab109 / Allpix Squared

tools/mesh_converter/MeshConverter.cpp

+62 −62

Original line number	Diff line number	Diff line
		@@ -32,7 +32,9 @@

		using namespace mesh_converter;
		using namespace ROOT::Math;
		using allpix::Log;
		using allpix::ThreadPool;
		using allpix::Units;

		void interrupt_handler(int);

		@@ -41,7 +43,7 @@ void interrupt_handler(int);
		*/
		void interrupt_handler(int) {
		LOG(STATUS) << "Interrupted! Aborting conversion...";
		allpix::Log::finish();
		Log::finish();
		std::exit(0);
		}

		@@ -66,7 +68,7 @@ int main(int argc, char** argv) {
		}

		// Add stream and set default logging level
		allpix::Log::addStream(std::cout);
		Log::addStream(std::cout);

		// Install abort handler (CTRL+\) and interrupt handler (CTRL+C)
		std::signal(SIGQUIT, interrupt_handler);
		@@ -84,7 +86,7 @@ int main(int argc, char** argv) {
		print_help = true;
		} else if(strcmp(argv[i], "-v") == 0 && (i + 1 < argc)) {
		try {
		log_level = allpix::Log::getLevelFromString(std::string(argv[++i]));
		log_level = Log::getLevelFromString(std::string(argv[++i]));
		} catch(std::invalid_argument& e) {
		LOG(ERROR) << "Invalid verbosity level \"" << std::string(argv[i]) << "\", ignoring overwrite";
		return_code = 1;
		@@ -136,7 +138,7 @@ int main(int argc, char** argv) {
		<< std::endl;
		std::cout << "\t -v <level> verbosity level (default reporiting level is INFO)" << std::endl;

		allpix::Log::finish();
		Log::finish();
		return return_code;
		}

		@@ -148,7 +150,7 @@ int main(int argc, char** argv) {
		auto log_level_string = config.get<std::string>("log_level", "INFO");
		std::transform(log_level_string.begin(), log_level_string.end(), log_level_string.begin(), ::toupper);
		try {
		log_level = allpix::Log::getLevelFromString(log_level_string);
		log_level = Log::getLevelFromString(log_level_string);
		} catch(std::invalid_argument& e) {
		LOG(ERROR) << "Log level \"" << log_level_string
		<< "\" specified in the configuration is invalid, defaulting to INFO instead";
		@@ -157,7 +159,7 @@ int main(int argc, char** argv) {
		}

		// Set log level:
		allpix::Log::setReportingLevel(log_level);
		Log::setReportingLevel(log_level);

		// NOTE: this stream should be available for the duration of the logging
		std::ofstream log_file;
		@@ -165,10 +167,10 @@ int main(int argc, char** argv) {
		log_file.open(log_file_name, std::ios_base::out \| std::ios_base::trunc);
		if(!log_file.good()) {
		LOG(FATAL) << "Cannot write to provided log file! Check if permissions are sufficient.";
		allpix::Log::finish();
		Log::finish();
		return 1;
		}
		allpix::Log::addStream(log_file);
		Log::addStream(log_file);
		}

		LOG(STATUS) << "Welcome to the Mesh Converter Tool of Allpix^2 " << ALLPIX_PROJECT_VERSION;
		@@ -188,28 +190,16 @@ int main(int argc, char** argv) {
		// Region, observable and binning of output field
		auto regions = config.getArray<std::string>("region");
		auto observable = config.get<std::string>("observable");

		const auto radius_step = config.get<double>("radius_step", 0.5);
		const auto max_radius = config.get<double>("max_radius", 50);
		const auto volume_cut = config.get<double>("volume_cut", 10e-9);
		const auto units = config.get<std::string>("observable_units");
		const auto vector_field = config.get<bool>("vector_field", (observable == "ElectricField"));

		XYZVectorUInt divisions;
		const auto dimension = config.get<size_t>("dimension", 3);
		if(dimension == 2) {
		auto divisions_yz = config.get<XYVectorUInt>("divisions", XYVectorUInt(100, 100));
		divisions = XYZVectorUInt(1, divisions_yz.x(), divisions_yz.y());
		} else if(dimension == 3) {
		divisions = config.get<XYZVectorUInt>("divisions", XYZVectorUInt(100, 100, 100));
		} else {
		throw allpix::InvalidValueError(config, "dimension", "only two or three dimensional fields are supported");
		}
		const auto radius_step = config.get<double>("radius_step", 0.5);
		const auto volume_cut = config.get<double>("volume_cut", 10e-9);

		// Swapping elements
		auto rot = config.getArray<std::string>("xyz", {"x", "y", "z"});
		if(rot.size() != 3) {
		throw allpix::InvalidValueError(config, "xyz", "three entries required");
		throw allpix::InvalidValueError(config, "xyz", "Three entries required");
		}

		auto start = std::chrono::system_clock::now();
		@@ -217,6 +207,22 @@ int main(int argc, char** argv) {
		std::string grid_file = file_prefix + ".grd";
		std::vector<Point> points = parser->getMesh(grid_file, regions);

		// Obtain number of mesh dimensions from mesh point:
		XYZVectorUInt divisions;
		const auto dimension = points.front().dim;
		LOG(STATUS) << "Determined dimensionality of input mesh to be " << dimension << "D";
		if(dimension == 2) {
		auto divisions_yz = config.get<XYVectorUInt>("divisions", XYVectorUInt(100, 100));
		divisions = XYZVectorUInt(1, divisions_yz.x(), divisions_yz.y());
		} else if(dimension == 3) {
		divisions = config.get<XYZVectorUInt>("divisions", XYZVectorUInt(100, 100, 100));
		}

		// If we are looking at a 2D mesh, we force x to be the coordinate out of range:
		if(dimension == 2 && rot.at(0) != "-x" && rot.at(0) != "x") {
		throw allpix::InvalidValueError(config, "xyz", "For 2D meshes the first coordinate has to remain 'x'");
		}

		std::string data_file = file_prefix + ".dat";
		std::vector<Point> field = parser->getField(data_file, observable, regions);

		@@ -267,27 +273,14 @@ int main(int argc, char** argv) {
		field = field_temp;

		// Find minimum and maximum from mesh coordinates
		double minx = DBL_MAX, miny = DBL_MAX, minz = DBL_MAX;
		double maxx = DBL_MIN, maxy = DBL_MIN, maxz = DBL_MIN;
		for(auto& point : points) {
		if(dimension == 2) {
		maxx = 1;
		minx = 0;
		maxy = std::max(maxy, point.y);
		maxz = std::max(maxz, point.z);
		miny = std::min(miny, point.y);
		minz = std::min(minz, point.z);
		}

		if(dimension == 3) {
		maxx = std::max(maxx, point.x);
		maxy = std::max(maxy, point.y);
		maxz = std::max(maxz, point.z);
		minx = std::min(minx, point.x);
		miny = std::min(miny, point.y);
		minz = std::min(minz, point.z);
		}
		}
		auto maxx =
		(dimension == 2 ? 1
		: std::max_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.x < b.x; })->x);
		auto minx = std::min_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.x < b.x; })->x;
		auto maxy = std::max_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.y < b.y; })->y;
		auto miny = std::min_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.y < b.y; })->y;
		auto maxz = std::max_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.z < b.z; })->z;
		auto minz = std::min_element(points.begin(), points.end(), [](auto& a, auto& b) { return a.z < b.z; })->z;

		// Creating a new mesh points cloud with a regular pitch
		const double xstep = (maxx - minx) / static_cast<double>(divisions.x());
		@@ -297,7 +290,9 @@ int main(int argc, char** argv) {

		// Using the minimal cell dimension as initial search radius for the point cloud:
		const auto initial_radius = config.get<double>("initial_radius", std::min({xstep, ystep, zstep}));
		LOG(INFO) << "Using initial neighbor search radius of " << initial_radius;
		const auto max_radius = config.get<double>("max_radius", std::max(initial_radius, 50.));
		LOG(INFO) << "Using initial neighbor search radius of " << initial_radius << " and maximum search radius of "
		<< max_radius;

		if(rot.at(0) != "x" \|\| rot.at(1) != "y" \|\| rot.at(2) != "z") {
		LOG(STATUS) << "TCAD mesh (x,y,z) coords. transformation into: (" << rot.at(0) << "," << rot.at(1) << ","
		@@ -333,10 +328,6 @@ int main(int argc, char** argv) {
		}
		}

		rot.at(0).erase(std::remove(rot.at(0).begin(), rot.at(0).end(), '-'), rot.at(0).end());
		rot.at(1).erase(std::remove(rot.at(1).begin(), rot.at(1).end(), '-'), rot.at(1).end());
		rot.at(2).erase(std::remove(rot.at(2).begin(), rot.at(2).end(), '-'), rot.at(2).end());

		auto end = std::chrono::system_clock::now();
		auto elapsed_seconds = std::chrono::duration_cast<std::chrono::seconds>(end - start).count();
		LOG(INFO) << "Reading the files took " << elapsed_seconds << " seconds.";
		@@ -347,7 +338,7 @@ int main(int argc, char** argv) {

		unsigned int mesh_points_done = 0;
		auto mesh_section = [&](double x, double y) {
		allpix::Log::setReportingLevel(log_level);
		Log::setReportingLevel(log_level);

		// New mesh slice
		std::vector<Point> new_mesh;
		@@ -355,19 +346,29 @@ int main(int argc, char** argv) {
		double z = minz + zstep / 2.0;
		for(unsigned int k = 0; k < divisions.z(); ++k) {
		// New mesh vertex and field
		Point q(dimension == 2 ? -1 : x, y, z), e;
		auto q = (dimension == 2 ? Point(y, z) : Point(x, y, z));
		Point e;
		bool valid = false;

		size_t prev_neighbours = 0;
		double radius = initial_radius;

		while(radius < max_radius) {
		while(radius <= max_radius) {
		LOG(DEBUG) << "Search radius: " << radius;
		// Calling octree neighbours search and sorting the results list with the closest neighbours first
		std::vector<unsigned int> results;
		octree.radiusNeighbors<unibn::L2Distance<Point>>(q, radius, results);
		LOG(DEBUG) << "Number of vertices found: " << results.size();

		if(radius == initial_radius && results.size() > 100) {
		LOG(WARNING) << "Found " << results.size() << " mesh vertices within initial search radius of "
		<< initial_radius << "um." << std::endl
		<< "This might indicate that the output mesh granularity is too low and field features "
		"might be missed."
		<< std::endl
		<< "Consider increasing output mesh granularity.";
		}

		// If after a radius step no new neighbours are found, go to the next radius step
		if(results.size() <= prev_neighbours \|\| results.empty()) {
		prev_neighbours = results.size();
		@@ -430,11 +431,11 @@ int main(int argc, char** argv) {
		std::vector<Point> e_field_new_mesh;

		// clang-format off
		auto init_function = [log_level = allpix::Log::getReportingLevel(), log_format = allpix::Log::getFormat()]() {
		auto init_function = [log_level = Log::getReportingLevel(), log_format = Log::getFormat()]() {
		// clang-format on
		// Initialize the threads to the same log level and format as the master setting
		allpix::Log::setReportingLevel(log_level);
		allpix::Log::setFormat(log_format);
		Log::setReportingLevel(log_level);
		Log::setFormat(log_format);
		};

		ThreadPool pool(num_threads, num_threads * 1024, init_function);
		@@ -465,9 +466,8 @@ int main(int argc, char** argv) {
		// Prepare header and auxiliary information:
		std::string header =
		"Allpix Squared " + std::string(ALLPIX_PROJECT_VERSION) + " TCAD Mesh Converter, observable: " + observable;
		std::array<double, 3> size{{allpix::Units::get(maxx - minx, "um"),
		allpix::Units::get(maxy - miny, "um"),
		allpix::Units::get(maxz - minz, "um")}};
		std::array<double, 3> size{
		{Units::get(maxx - minx, "um"), Units::get(maxy - miny, "um"), Units::get(maxz - minz, "um")}};
		std::array<size_t, 3> gridsize{
		{static_cast<size_t>(divisions.x()), static_cast<size_t>(divisions.y()), static_cast<size_t>(divisions.z())}};

		@@ -480,11 +480,11 @@ int main(int argc, char** argv) {
		for(unsigned int k = 0; k < divisions.z(); ++k) {
		auto& point = e_field_new_mesh[i * divisions.y() * divisions.z() + j * divisions.z() + k];
		// We need to convert to framework-internal units:
		data->push_back(allpix::Units::get(point.x, units));
		data->push_back(Units::get(point.x, units));
		// For a vector field, we push three values:
		if(quantity == FieldQuantity::VECTOR) {
		data->push_back(allpix::Units::get(point.y, units));
		data->push_back(allpix::Units::get(point.z, units));
		data->push_back(Units::get(point.y, units));
		data->push_back(Units::get(point.z, units));
		}
		}
		}
		@@ -512,6 +512,6 @@ int main(int argc, char** argv) {
		}

		// Finish the logging
		allpix::Log::finish();
		Log::finish();
		return return_code;
		}

tools/mesh_converter/MeshElement.cpp

+3 −3

Original line number	Diff line number	Diff line
		@@ -58,12 +58,12 @@ double MeshElement::get_distance(size_t index, Point& qp) const {
		return unibn::L2Distance<Point>::compute(vertices_[index], qp);
		}

		bool MeshElement::validElement(double volume_cut, Point& qp) const {
		bool MeshElement::isValid(double volume_cut, Point& qp) const {
		if(std::fabs(volume_) < MIN_VOLUME) {
		LOG(TRACE) << "Invalid tetrahedron with coplanar(3D)/colinear(2D) vertices.";
		LOG(TRACE) << "Invalid tetrahedron, all vertices are " << (dimension_ == 3 ? "coplanar" : "colinear");
		return false;
		} else if(std::fabs(volume_) <= volume_cut) {
		LOG(TRACE) << "Tetrahedron volume smaller than volume cut.";
		LOG(TRACE) << "Invalid tetrahedron with volume " << std::fabs(volume_) << " smaller than volume cut " << volume_cut;
		return false;
		}

tools/mesh_converter/MeshElement.hpp

+9 −5

Original line number	Diff line number	Diff line
		@@ -13,9 +13,11 @@ namespace mesh_converter {
		class Point {
		public:
		Point() noexcept = default;
		Point(double px, double py, double pz = 0) noexcept : x(px), y(py), z(pz) {}
		Point(double px, double py, double pz) noexcept : x(px), y(py), z(pz), dim(3){};
		Point(double py, double pz) noexcept : y(py), z(pz), dim(2){};

		double x{0}, y{0}, z{0};
		unsigned int dim{0};

		friend std::ostream& operator<<(std::ostream& out, const Point& pt) {
		out << "(" << pt.x << "," << pt.y << "," << pt.z << ")";
		@@ -53,7 +55,7 @@ namespace mesh_converter {
		* @param volume_cut Threshold for the minimum tetrahedron volume
		* @param qp Desired point for the interpolation
		*/
		bool validElement(double volume_cut, Point& qp) const;
		bool isValid(double volume_cut, Point& qp) const;

		/**
		* @brief Barycentric interpolation implementation
		@@ -129,20 +131,22 @@ namespace mesh_converter {
		// Dimensionality is number of iterator elements minus one:
		size_t dimensions = static_cast<size_t>(end - begin) - 1;
		size_t idx = 0;

		LOG(TRACE) << "Constructing " << dimensions << "D element at " << reference_ << " with mesh points:";
		for(; begin < end; begin++) {
		LOG(TRACE) << "\t\t" << (grid_)[begin];
		grid_elements[idx] = (grid_)[begin];
		field_elements[idx++] = (field_)[begin];
		}

		LOG(TRACE) << "Constructing element with dim " << dimensions << " at " << reference_;
		MeshElement element(dimensions, grid_elements, field_elements);
		valid_ = element.validElement(cut_, reference_);
		valid_ = element.isValid(cut_, reference_);
		if(valid_) {
		LOG(DEBUG) << element.print(reference_);
		result_ = element.getObservable(reference_);
		}

		return valid_; // Don't break out of the loop if element is invalid
		return valid_;
		}

		/**

tools/mesh_converter/README.md

+1 −1

Original line number	Diff line number	Diff line
		@@ -61,6 +61,7 @@ It should be noted that the Mesh Converter depends on the core utilities of the

		## Features
		- TCAD DF-ISE file format parser.
		- Automatic determination of the input mesh dimensionality (2D/3D).
		- Fast radius neighbor search for three-dimensional point clouds.
		- Barycentric interpolation between non-regular mesh points.
		- Several cuts available on the interpolation algorithm variables.
		@@ -69,7 +70,6 @@ It should be noted that the Mesh Converter depends on the core utilities of the
		### Parameters
		* `model`: Field file format to use, can be INIT or APF, defaults to APF (binary format).
		* `parser`: Parser class to interpret input data in. Currently, only DF-ISE is supported and used as default.
		* `dimension`: Specify mesh dimensionality (defaults to 3).
		* `region`: Region name or list of region names to be meshed, such as for example `bulk` or `"bulk","epi"` (No default value; required parameter).
		* `observable`: Observable to be interpolated, such as for example `ElectricField` (No default value; required parameter).
		* `observable_units`: Units in which the observable is stored in the input file (No default value; required parameter).

tools/mesh_converter/parsers/DFISEParser.cpp

+6 −12

Original line number	Diff line number	Diff line
		@@ -38,8 +38,6 @@ MeshMap DFISEParser::read_meshes(const std::string& file_name) {

		std::map<std::string, std::vector<long unsigned int>> regions_vertices;

		Point point(-1.0, -1.0, -1.0);

		std::string region;
		long unsigned int dimension = 1;
		long unsigned int data_count = 0;
		@@ -192,19 +190,15 @@ MeshMap DFISEParser::read_meshes(const std::string& file_name) {
		case DFSection::VERTICES: {
		// Read vertex points
		if(dimension == 3) {
		point.x = -1.0;
		point.y = -1.0;
		point.z = -1.0;
		while(sstr >> point.x >> point.y >> point.z) {
		vertices.push_back(point);
		double x = 0, y = 0, z = 0;
		while(sstr >> x >> y >> z) {
		vertices.emplace_back(x, y, z);
		}
		}
		if(dimension == 2) {
		point.x = -1.0;
		point.y = -1.0;
		point.z = -1.0;
		while(sstr >> point.y >> point.z) {
		vertices.push_back(point);
		double y = 0, z = 0;
		while(sstr >> y >> z) {
		vertices.emplace_back(y, z);
		}
		}
		} break;