Merge branch 'p/meshcoplanarity' into 'master' (674daf46) · Commits · XRFLab109 / Allpix Squared

tools/mesh_converter/MeshConverter.cpp

+25 −4

Original line number	Diff line number	Diff line
		@@ -203,6 +203,8 @@ int main(int argc, char** argv) {
		const auto units = config.get<std::string>("observable_units");
		const auto vector_field = config.get<bool>("vector_field", (observable == "ElectricField"));

		const auto interpolate = config.get<bool>("interpolate", true);
		const auto allow_decay = config.get<bool>("allow_coplanar_interpolation", false);
		const auto radius_step = config.get<double>("radius_step", 0.5);
		const auto volume_cut = config.get<double>("volume_cut", 10e-9);

		@@ -358,8 +360,17 @@ int main(int argc, char** argv) {
		// New mesh vertex and field
		auto q = (dimension == 2 ? Point(y, z) : Point(x, y, z));
		Point e;
		bool valid = false;

		// No interpolation requested, return nearest neighbor:
		if(!interpolate) {
		auto idx = static_cast<size_t>(octree.findNeighbor<unibn::L2Distance<Point>>(q));
		new_mesh.push_back(field.at(idx));
		z += zstep;
		continue;
		}

		bool valid = false;
		bool allow_zero_volume = false;
		size_t prev_neighbours = 0;
		double radius = initial_radius;

		@@ -394,6 +405,15 @@ int main(int argc, char** argv) {
		continue;
		}

		// If we have too many neighbors, we could decay to using lower-dimension interpolation:
		if(allow_decay && radius > initial_radius && results.size() > 100) {
		LOG_ONCE(WARNING) << "Large number of neighbors found, this hints to a quasi-co"
		<< (dimension == 3 ? "planar" : "linear") << " situation" << std::endl
		<< "Decaying to interpolation in " << (dimension == 3 ? "planar" : "linear")
		<< " space for affected points";
		allow_zero_volume = true;
		}

		// Sort by lowest distance first, this drastically reduces the number of permutations required to find a
		// valid mesh element and also ensures that this is the one with the smallest volume.
		std::sort(results.begin(), results.end(), [&](unsigned int a, unsigned int b) {
		@@ -406,7 +426,7 @@ int main(int argc, char** argv) {
		auto res = for_each_combination(results.begin(),
		results.begin() + (dimension == 3 ? 4 : 3),
		results.end(),
		Combination(&points, &field, q, volume_cut));
		Combination(&points, &field, q, allow_zero_volume ? 0 : volume_cut));
		valid = res.valid();
		if(valid) {
		e = res.result();
		@@ -434,8 +454,9 @@ int main(int argc, char** argv) {
		}

		mesh_points_done += divisions.z();
		LOG_PROGRESS(STATUS, "m") << "Interpolating new mesh: " << mesh_points_done << " of " << mesh_points_total
		<< ", " << (mesh_points_done / (mesh_points_total / 100)) << "%";
		LOG_PROGRESS(STATUS, "m") << (interpolate ? "Interpolating" : "Generating") << " new mesh: " << mesh_points_done
		<< " of " << mesh_points_total << ", "
		<< (mesh_points_done / (mesh_points_total / 100)) << "%";

		return new_mesh;
		};

tools/mesh_converter/MeshElement.cpp

+10 −6

Original line number	Diff line number	Diff line
		@@ -69,13 +69,17 @@ double MeshElement::get_distance(size_t index, Point& qp) const {
		}

		bool MeshElement::isValid(double volume_cut, Point& qp) const {
		// Check if we should apply coplanar/colinear criterions:
		if(volume_cut > 0) {
		if(std::fabs(volume_) < MIN_VOLUME) {
		LOG(TRACE) << "Invalid tetrahedron, all vertices are " << (dimension_ == 3 ? "coplanar" : "colinear");
		return false;
		} else if(std::fabs(volume_) <= volume_cut) {
		LOG(TRACE) << "Invalid tetrahedron with volume " << std::fabs(volume_) << " smaller than volume cut " << volume_cut;
		LOG(TRACE) << "Invalid tetrahedron with volume " << std::fabs(volume_) << " smaller than volume cut "
		<< volume_cut;
		return false;
		}
		}

		for(size_t i = 0; i < dimension_ + 1; i++) {
		if(volume_ * get_sub_volume(i, qp) < 0) {

tools/mesh_converter/MeshElement.hpp

+8 −1

Original line number	Diff line number	Diff line
		@@ -13,6 +13,7 @@

		#include <Eigen/Eigen>
		#include <array>
		#include <cmath>
		#include <utility>

		#include "core/utils/log.h"
		@@ -26,6 +27,8 @@ namespace mesh_converter {
		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){};

		bool isFinite() const { return std::isfinite(x) && std::isfinite(y) && std::isfinite(z); };

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

		@@ -62,7 +65,7 @@ namespace mesh_converter {

		/**
		* @brief Checks if the tetrahedron is valid for the interpolation
		* @param volume_cut Threshold for the minimum tetrahedron volume
		* @param volume_cut Threshold for the minimum tetrahedron volume. Values <= 0 disable coplanarity checks
		* @param qp Desired point for the interpolation
		*/
		bool isValid(double volume_cut, Point& qp) const;
		@@ -154,6 +157,10 @@ namespace mesh_converter {
		if(valid_) {
		LOG(DEBUG) << element.print(reference_);
		result_ = element.getObservable(reference_);
		if(!result_.isFinite()) {
		LOG(WARNING) << "Interpolated result not a finite number at " << reference_;
		return false;
		}
		}

		return valid_;

tools/mesh_converter/README.md

+24 −8

Original line number	Diff line number	Diff line
		@@ -4,11 +4,25 @@
		title: "Mesh Converter"
		---

		This code takes adaptive meshes from finite-element simulations and transforms them into a regularly spaced grid for faster field value lookup as reuqired by Monte Carlo simulations tools such as Allpix Squared.
		The input consists of two files, one containing the vertex coordinates of each input mesh node, the other providing the relevant field values associated to each of these vertices.
		One output file containing the regular interpolated mesh is produced.
		This code takes adaptive meshes from finite-element simulations and transforms them into a regularly spaced grid for faster
		field value lookup as required by Monte Carlo simulations tools such as Allpix Squared. The input consists of two files, one
		containing the vertex coordinates of each input mesh node, the other providing the relevant field values associated to each
		of these vertices. One output file containing the regular mesh is produced. This tool can work in two different modes, the
		closest-neighbor mode and interpolation mode:

		A new regular mesh is created by scanning the model volume in regular X, Y and Z steps (not necessarily coinciding with original mesh nodes) and using a barycentric interpolation method to calculate the respective electric field vector on the new point. The interpolation uses the four closest, no-coplanar, neighbor vertex nodes such, that the respective tetrahedron encloses the query point. For the neighbors search, the software uses the Octree implementation \[[@octree]\].
		### Simple Closest-Neighbor Search

		In this mode, selected by setting the parameter `interpolate = false`, no interpolation of field values is performed, but for
		every output mesh point, the values from the closest neighbor of the input mesh are taken. In most cases this approach should
		produces reasonably precise results with a granularity similar to the respective adaptive mesh granularity in the respective
		region. The tool uses the Octree `findNeighbor` algorithm \[[@octree]\] to find the closest neighbor to the query point.

		### Barycentric Interpolation Method

		In this mode, the regular mesh is created by scanning the model volume in regular X, Y and Z steps and using a barycentric
		interpolation method to calculate the respective electric field vector on the new point. The interpolation uses the four
		closest, no-coplanar, neighbor vertex nodes such, that the respective tetrahedron encloses the query point. For the neighbors
		search, the tool uses the Octree `radiusNeighbors` neighbor search algorithm \[[@octree]\].

		## File Formats

		@@ -77,11 +91,13 @@ It should be noted that the Mesh Converter depends on the core utilities of the
		* `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).
		* `interpolate`: Boolean switch to select either the barycentric interpolation method or the closest-neighbor method. Defaults to `true`, i.e. using the interpolation method.
		* `initial_radius`: Initial node neighbors search radius in micro meters. Defaults to the minimal cell dimension of the final interpolated mesh.
		* `radius_step`: Radius step if no neighbor is found (defaults to `0.5um`).
		* `max_radius`: Maximum search radius (default is `50um`).
		* `allow_failure`: Allow the interpolation of a single mesh point to fail, i.e. when no neighbors could be found. If set to `true`, the respective mesh element will be set to zero and the interpolation will continue, if `false` the interpolation will be aborted. Defaults to `false`.
		* `volume_cut`: Minimum volume for tetrahedron for non-coplanar vertices (defaults to minimum double value).
		* `radius_step`: Radius step if no neighbor is found (defaults to `0.5um`). Only used for barycentric interpolation.
		* `max_radius`: Maximum search radius (default is `50um`). Only used for barycentric interpolation.
		* `allow_coplanar_interpolation`: Allow the interpolation to use coplanar/colinear vertices if no full interpolation volume can be found after increasing the search radius and if more than 100 neighbors are found. Defaults to `false`. It should be noted that this feature is experimental and that it can produce `NaN` results for the interpolated field.
		* `allow_failure`: Allow the interpolation of a single mesh point to fail, i.e. when no neighbors could be found. If set to `true`, the respective mesh element will be set to zero and the interpolation will continue, if `false` the interpolation will be aborted. Defaults to `false`. Only used for barycentric interpolation.
		* `volume_cut`: Minimum volume for tetrahedron for non-coplanar vertices (defaults to minimum double value). Only used for barycentric interpolation.
		* `divisions`: Number of divisions of the new regular mesh for each dimension, 2D or 3D vector depending on the `dimension` setting. Defaults to 100 bins in each dimension.
		* `xyz`: Array to replace the system coordinates of the mesh. A detailed description of how to use this parameter is given below.
		* `workers`: Number of worker threads to be used for the interpolation. Defaults to the available number of cores on the machine (hardware concurrency).