diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index d15f3cf7707..a346007fd20 100644 --- a/doc/sphinx/analysis.rst +++ b/doc/sphinx/analysis.rst @@ -411,10 +411,21 @@ or bin edges for the axes. Example:: density_profile.min_y, density_profile.max_y]) plt.show() +Observables based on cylindrical coordinates are also available. +They require special parameters if the cylindrical coordinate system is non-standard, e.g. if you want the origin of the cylindrical coordinates to be at a special location of the box or if you want to make use of symmetries along an axis that is not parallel to the z-axis. +For this purpose, use :class:`espressomd.math.CylindricalTransformationParameters` to create a consistent set of the parameters needed. Example:: + + import espressomd.math + + # shifted and rotated cylindrical coordinates + cyl_transform_params = espressomd.math.CylindricalTransformationParameters( + center=[5.0, 5.0, 0.0], axis=[0, 1, 0], orientation=[0, 0, 1]) + # histogram in cylindrical coordinates density_profile = espressomd.observables.CylindricalDensityProfile( - ids=[0, 1], center=[5.0, 5.0, 0.0], axis=[0, 0, 1], - n_r_bins=8, min_r=0.0, max_r=4.0, + ids=[0, 1], + transform_params = cyl_transform_params, + n_r_bins=8, min_r=1.0, max_r=4.0, n_phi_bins=16, min_phi=-np.pi, max_phi=np.pi, n_z_bins=4, min_z=4.0, max_z=8.0) obs_data = density_profile.calculate() @@ -779,4 +790,3 @@ Note that the cluster objects do not contain copies of the particles, but refer - diff --git a/doc/tutorials/charged_system/charged_system-1.ipynb b/doc/tutorials/charged_system/charged_system-1.ipynb index d743b8731db..7f223216e8e 100644 --- a/doc/tutorials/charged_system/charged_system-1.ipynb +++ b/doc/tutorials/charged_system/charged_system-1.ipynb @@ -36,7 +36,7 @@ "import espressomd\n", "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", "\n", - "from espressomd import System, interactions, electrostatics, observables, accumulators\n", + "from espressomd import System, interactions, electrostatics, observables, accumulators, math\n", "\n", "import numpy as np\n", "from scipy import optimize\n", @@ -438,20 +438,18 @@ "```python\n", "def setup_profile_calculation(system, delta_N, ion_types, r_min, n_radial_bins):\n", " radial_profile_accumulators = {}\n", + " ctp = math.CylindricalTransformationParameters(center = np.array(system.box_l) / 2.,\n", + " axis = [0, 0, 1],\n", + " orientation = [1, 0, 0])\n", " for ion_type in ion_types:\n", " ion_ids = system.part.select(type=ion_type).id\n", " radial_profile_obs = observables.CylindricalDensityProfile(\n", " ids=ion_ids,\n", - " center=np.array(system.box_l) / 2.,\n", - " axis=[0, 0, 1, ],\n", + " transform_params = ctp,\n", " n_r_bins=n_radial_bins,\n", - " n_phi_bins=1,\n", - " n_z_bins=1,\n", " min_r=r_min,\n", - " min_phi=-np.pi,\n", " min_z=-system.box_l[2] / 2.,\n", " max_r=system.box_l[0] / 2.,\n", - " max_phi=np.pi,\n", " max_z=system.box_l[2] / 2.)\n", "\n", " bin_edges = radial_profile_obs.bin_edges()\n", @@ -945,7 +943,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/samples/lb_profile.py b/samples/lb_profile.py index 7f4f6235dd6..06e6e05f516 100644 --- a/samples/lb_profile.py +++ b/samples/lb_profile.py @@ -32,6 +32,7 @@ import espressomd.shapes import espressomd.lbboundaries import espressomd.accumulators +import espressomd.math system = espressomd.System(box_l=[10.0, 10.0, 5.0]) system.time_step = 0.01 @@ -42,16 +43,14 @@ agrid=1.0, dens=1.0, visc=1.0, tau=0.01, ext_force_density=[0, 0, 0.15], kT=1.0, seed=32) system.actors.add(lb_fluid) system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23) -fluid_obs = espressomd.observables.CylindricalLBVelocityProfile( +ctp = espressomd.math.CylindricalTransformationParameters( center=[5.0, 5.0, 0.0], axis=[0, 0, 1], + orientation=[1, 0, 0]) +fluid_obs = espressomd.observables.CylindricalLBVelocityProfile( + transform_params=ctp, n_r_bins=100, - n_phi_bins=1, - n_z_bins=1, - min_r=0.0, max_r=4.0, - min_phi=-np.pi, - max_phi=np.pi, min_z=0.0, max_z=10.0, sampling_density=0.1) diff --git a/src/core/grid_based_algorithms/lb_collective_interface.cpp b/src/core/grid_based_algorithms/lb_collective_interface.cpp index f6a5b9aa69a..a2e1cb8e13c 100644 --- a/src/core/grid_based_algorithms/lb_collective_interface.cpp +++ b/src/core/grid_based_algorithms/lb_collective_interface.cpp @@ -81,6 +81,15 @@ mpi_lb_get_interpolated_velocity(Utils::Vector3d const &pos) { REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_interpolated_velocity) +boost::optional +mpi_lb_get_interpolated_density(Utils::Vector3d const &pos) { + return detail::lb_calc_for_pos(pos, [&](auto pos) { + return lb_lbinterpolation_get_interpolated_density(pos); + }); +} + +REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_interpolated_density) + auto mpi_lb_get_density(Utils::Vector3i const &index) { return detail::lb_calc_fluid_kernel( index, [&](auto const &modes, auto const &force_density) { diff --git a/src/core/grid_based_algorithms/lb_collective_interface.hpp b/src/core/grid_based_algorithms/lb_collective_interface.hpp index 1c9afa8fb31..4b6b0272ab4 100644 --- a/src/core/grid_based_algorithms/lb_collective_interface.hpp +++ b/src/core/grid_based_algorithms/lb_collective_interface.hpp @@ -27,6 +27,8 @@ /* collective getter functions */ boost::optional mpi_lb_get_interpolated_velocity(Utils::Vector3d const &pos); +boost::optional +mpi_lb_get_interpolated_density(Utils::Vector3d const &pos); boost::optional mpi_lb_get_density(Utils::Vector3i const &index); boost::optional mpi_lb_get_populations(Utils::Vector3i const &index); diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp index d6c0e58ae89..3fabd382a99 100644 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ b/src/core/grid_based_algorithms/lb_interface.cpp @@ -1273,3 +1273,23 @@ lb_lbfluid_get_interpolated_velocity(const Utils::Vector3d &pos) { } throw NoLBActive(); } + +double lb_lbfluid_get_interpolated_density(const Utils::Vector3d &pos) { + auto const folded_pos = folded_position(pos, box_geo); + auto const interpolation_order = lb_lbinterpolation_get_interpolation_order(); + if (lattice_switch == ActiveLB::GPU) { + throw std::runtime_error( + "Density interpolation is not implemented for the GPU LB."); + } + if (lattice_switch == ActiveLB::CPU) { + switch (interpolation_order) { + case (InterpolationOrder::quadratic): + throw std::runtime_error("The non-linear interpolation scheme is not " + "implemented for the CPU LB."); + case (InterpolationOrder::linear): + return mpi_call(::Communication::Result::one_rank, + mpi_lb_get_interpolated_density, folded_pos); + } + } + throw NoLBActive(); +} diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp index c72ef9b7f8b..9b831197a68 100644 --- a/src/core/grid_based_algorithms/lb_interface.hpp +++ b/src/core/grid_based_algorithms/lb_interface.hpp @@ -265,4 +265,11 @@ Utils::Vector3d lb_lbfluid_calc_fluid_momentum(); const Utils::Vector3d lb_lbfluid_get_interpolated_velocity(const Utils::Vector3d &pos); +/** + * @brief Calculates the interpolated fluid density on the master process. + * @param pos Position at which the density is to be calculated. + * @retval interpolated fluid density. + */ +double lb_lbfluid_get_interpolated_density(const Utils::Vector3d &pos); + #endif diff --git a/src/core/grid_based_algorithms/lb_interpolation.cpp b/src/core/grid_based_algorithms/lb_interpolation.cpp index ea0cd46e1dd..d0db33f33ae 100644 --- a/src/core/grid_based_algorithms/lb_interpolation.cpp +++ b/src/core/grid_based_algorithms/lb_interpolation.cpp @@ -82,6 +82,16 @@ Utils::Vector3d node_u(Lattice::index_t index) { return Utils::Vector3d{modes[1], modes[2], modes[3]} / local_density; } +double node_dens(Lattice::index_t index) { +#ifdef LB_BOUNDARIES + if (lbfields[index].boundary) { + return lbpar.density; + } +#endif // LB_BOUNDARIES + auto const modes = lb_calc_modes(index, lbfluid); + return lbpar.density + modes[0]; +} + } // namespace const Utils::Vector3d @@ -98,6 +108,19 @@ lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &pos) { return interpolated_u; } +double lb_lbinterpolation_get_interpolated_density(const Utils::Vector3d &pos) { + double interpolated_dens = 0.; + + /* Calculate fluid density at the position. + This is done by linear interpolation (eq. (11) @cite ahlrichs99a) */ + lattice_interpolation(lblattice, pos, + [&interpolated_dens](Lattice::index_t index, double w) { + interpolated_dens += w * node_dens(index); + }); + + return interpolated_dens; +} + void lb_lbinterpolation_add_force_density( const Utils::Vector3d &pos, const Utils::Vector3d &force_density) { switch (interpolation_order) { diff --git a/src/core/grid_based_algorithms/lb_interpolation.hpp b/src/core/grid_based_algorithms/lb_interpolation.hpp index 68d9c6e1be6..28544afa9ef 100644 --- a/src/core/grid_based_algorithms/lb_interpolation.hpp +++ b/src/core/grid_based_algorithms/lb_interpolation.hpp @@ -41,6 +41,13 @@ InterpolationOrder lb_lbinterpolation_get_interpolation_order(); const Utils::Vector3d lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &p); +/** + * @brief Calculates the fluid density at a given position of the + * lattice. + * @note It can lead to undefined behaviour if the + * position is not within the local lattice. */ +double lb_lbinterpolation_get_interpolated_density(const Utils::Vector3d &p); + /** * @brief Add a force density to the fluid at the given position. */ diff --git a/src/core/observables/CylindricalDensityProfile.hpp b/src/core/observables/CylindricalDensityProfile.hpp index d18ea8baeac..a7d1de6e312 100644 --- a/src/core/observables/CylindricalDensityProfile.hpp +++ b/src/core/observables/CylindricalDensityProfile.hpp @@ -41,7 +41,9 @@ class CylindricalDensityProfile : public CylindricalPidProfileObservable { for (auto p : particles) { histogram.update(Utils::transform_coordinate_cartesian_to_cylinder( - folded_position(traits.position(p), box_geo) - center, axis)); + folded_position(traits.position(p), box_geo) - + transform_params->center(), + transform_params->axis(), transform_params->orientation())); } histogram.normalize(); diff --git a/src/core/observables/CylindricalFluxDensityProfile.hpp b/src/core/observables/CylindricalFluxDensityProfile.hpp index 13ef33b9178..73b65f08c81 100644 --- a/src/core/observables/CylindricalFluxDensityProfile.hpp +++ b/src/core/observables/CylindricalFluxDensityProfile.hpp @@ -43,11 +43,13 @@ class CylindricalFluxDensityProfile : public CylindricalPidProfileObservable { // Write data to the histogram for (auto p : particles) { - auto const pos = folded_position(traits.position(p), box_geo) - center; + auto const pos = folded_position(traits.position(p), box_geo) - + transform_params->center(); histogram.update( - Utils::transform_coordinate_cartesian_to_cylinder(pos, axis), - Utils::transform_vector_cartesian_to_cylinder(traits.velocity(p), - axis, pos)); + Utils::transform_coordinate_cartesian_to_cylinder( + pos, transform_params->axis(), transform_params->orientation()), + Utils::transform_vector_cartesian_to_cylinder( + traits.velocity(p), transform_params->axis(), pos)); } histogram.normalize(); return histogram.get_histogram(); diff --git a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp index de51d4a9e10..5ef211f40db 100644 --- a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp @@ -42,13 +42,24 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate( auto const pos = folded_position(traits.position(p), box_geo); auto const v = lb_lbfluid_get_interpolated_velocity(pos) * lb_lbfluid_get_lattice_speed(); + auto const flux_dens = lb_lbfluid_get_interpolated_density(pos) * v; - histogram.update( - Utils::transform_coordinate_cartesian_to_cylinder(pos - center, axis), - Utils::transform_vector_cartesian_to_cylinder(v, axis, pos - center)); + histogram.update(Utils::transform_coordinate_cartesian_to_cylinder( + pos - transform_params->center(), + transform_params->axis(), + transform_params->orientation()), + Utils::transform_vector_cartesian_to_cylinder( + flux_dens, transform_params->axis(), + pos - transform_params->center())); } - histogram.normalize(); - return histogram.get_histogram(); + // normalize by number of hits per bin + auto hist_tmp = histogram.get_histogram(); + auto tot_count = histogram.get_tot_count(); + std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(), + hist_tmp.begin(), [](auto hi, auto ci) { + return ci > 0 ? hi / static_cast(ci) : 0.; + }); + return hist_tmp; } } // namespace Observables diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp index c1d5eea4d7e..df4e66e5936 100644 --- a/src/core/observables/CylindricalLBProfileObservable.hpp +++ b/src/core/observables/CylindricalLBProfileObservable.hpp @@ -21,6 +21,7 @@ #include "CylindricalProfileObservable.hpp" +#include #include #include #include @@ -30,15 +31,15 @@ namespace Observables { class CylindricalLBProfileObservable : public CylindricalProfileObservable { public: - CylindricalLBProfileObservable(Utils::Vector3d const ¢er, - Utils::Vector3d const &axis, int n_r_bins, - int n_phi_bins, int n_z_bins, double min_r, - double max_r, double min_phi, double max_phi, - double min_z, double max_z, - double sampling_density) - : CylindricalProfileObservable(center, axis, n_r_bins, n_phi_bins, - n_z_bins, min_r, max_r, min_phi, max_phi, - min_z, max_z), + CylindricalLBProfileObservable( + std::shared_ptr + transform_params, + int n_r_bins, int n_phi_bins, int n_z_bins, double min_r, double max_r, + double min_phi, double max_phi, double min_z, double max_z, + double sampling_density) + : CylindricalProfileObservable(std::move(transform_params), n_r_bins, + n_phi_bins, n_z_bins, min_r, max_r, + min_phi, max_phi, min_z, max_z), sampling_density(sampling_density) { calculate_sampling_positions(); } @@ -47,17 +48,16 @@ class CylindricalLBProfileObservable : public CylindricalProfileObservable { limits[0], limits[1], limits[2], n_bins[0], n_bins[1], n_bins[2], sampling_density); for (auto &p : sampling_positions) { - double theta; - Utils::Vector3d rotation_axis; - auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian( - p, Utils::Vector3d{{0.0, 0.0, 1.0}}); + auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian(p); // We have to rotate the coordinates since the utils function assumes // z-axis symmetry. - std::tie(theta, rotation_axis) = - Utils::rotation_params(Utils::Vector3d{{0.0, 0.0, 1.0}}, axis); + constexpr Utils::Vector3d z_axis{{0.0, 0.0, 1.0}}; + auto const theta = Utils::angle_between(z_axis, transform_params->axis()); + auto const rot_axis = + Utils::vector_product(z_axis, transform_params->axis()).normalize(); if (theta > std::numeric_limits::epsilon()) - p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart); - p = p_cart + center; + p_cart = Utils::vec_rotate(rot_axis, theta, p_cart); + p = p_cart + transform_params->center(); } } std::vector sampling_positions; diff --git a/src/core/observables/CylindricalLBVelocityProfile.cpp b/src/core/observables/CylindricalLBVelocityProfile.cpp index 28148f5d493..791417859b8 100644 --- a/src/core/observables/CylindricalLBVelocityProfile.cpp +++ b/src/core/observables/CylindricalLBVelocityProfile.cpp @@ -35,11 +35,12 @@ std::vector CylindricalLBVelocityProfile::operator()() const { for (auto const &p : sampling_positions) { auto const velocity = lb_lbfluid_get_interpolated_velocity(p) * lb_lbfluid_get_lattice_speed(); - auto const pos_shifted = p - center; - auto const pos_cyl = - Utils::transform_coordinate_cartesian_to_cylinder(pos_shifted, axis); - histogram.update(pos_cyl, Utils::transform_vector_cartesian_to_cylinder( - velocity, axis, pos_shifted)); + auto const pos_shifted = p - transform_params->center(); + auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( + pos_shifted, transform_params->axis(), transform_params->orientation()); + histogram.update(pos_cyl, + Utils::transform_vector_cartesian_to_cylinder( + velocity, transform_params->axis(), pos_shifted)); } auto hist_data = histogram.get_histogram(); auto const tot_count = histogram.get_tot_count(); diff --git a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp index e1106574cfa..9650c7dda6f 100644 --- a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp @@ -41,17 +41,20 @@ std::vector CylindricalLBVelocityProfileAtParticlePositions::evaluate( lb_lbfluid_get_lattice_speed(); histogram.update( - Utils::transform_coordinate_cartesian_to_cylinder(pos - center, axis), - Utils::transform_vector_cartesian_to_cylinder(v, axis, pos - center)); + Utils::transform_coordinate_cartesian_to_cylinder( + pos - transform_params->center(), transform_params->axis(), + transform_params->orientation()), + Utils::transform_vector_cartesian_to_cylinder( + v, transform_params->axis(), pos - transform_params->center())); } + // normalize by number of hits per bin auto hist_tmp = histogram.get_histogram(); auto tot_count = histogram.get_tot_count(); - for (size_t ind = 0; ind < hist_tmp.size(); ++ind) { - if (tot_count[ind] > 0) { - hist_tmp[ind] /= static_cast(tot_count[ind]); - } - } + std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(), + hist_tmp.begin(), [](auto hi, auto ci) { + return ci > 0 ? hi / static_cast(ci) : 0.; + }); return hist_tmp; } diff --git a/src/core/observables/CylindricalPidProfileObservable.hpp b/src/core/observables/CylindricalPidProfileObservable.hpp index bbece2be6e5..8a3ae7886e1 100644 --- a/src/core/observables/CylindricalPidProfileObservable.hpp +++ b/src/core/observables/CylindricalPidProfileObservable.hpp @@ -19,6 +19,8 @@ #ifndef OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP +#include + #include "CylindricalProfileObservable.hpp" #include "PidObservable.hpp" @@ -27,16 +29,16 @@ namespace Observables { class CylindricalPidProfileObservable : public PidObservable, public CylindricalProfileObservable { public: - CylindricalPidProfileObservable(std::vector const &ids, - Utils::Vector3d const ¢er, - Utils::Vector3d const &axis, int n_r_bins, - int n_phi_bins, int n_z_bins, double min_r, - double max_r, double min_phi, double max_phi, - double min_z, double max_z) + CylindricalPidProfileObservable( + std::vector const &ids, + std::shared_ptr + transform_params, + int n_r_bins, int n_phi_bins, int n_z_bins, double min_r, double max_r, + double min_phi, double max_phi, double min_z, double max_z) : PidObservable(ids), - CylindricalProfileObservable(center, axis, n_r_bins, n_phi_bins, - n_z_bins, min_r, max_r, min_phi, max_phi, - min_z, max_z) {} + CylindricalProfileObservable(std::move(transform_params), n_r_bins, + n_phi_bins, n_z_bins, min_r, max_r, + min_phi, max_phi, min_z, max_z) {} }; } // Namespace Observables diff --git a/src/core/observables/CylindricalProfileObservable.hpp b/src/core/observables/CylindricalProfileObservable.hpp index c28669a4e62..224856f0cfd 100644 --- a/src/core/observables/CylindricalProfileObservable.hpp +++ b/src/core/observables/CylindricalProfileObservable.hpp @@ -22,12 +22,17 @@ #include "ProfileObservable.hpp" #include +#include +#include #include #include #include #include +#include +#include +#include #include namespace Observables { @@ -35,16 +40,16 @@ namespace Observables { /** Cylindrical profile observable */ class CylindricalProfileObservable : public ProfileObservable { public: - CylindricalProfileObservable(Utils::Vector3d const ¢er, - Utils::Vector3d const &axis, int n_r_bins, - int n_phi_bins, int n_z_bins, double min_r, - double max_r, double min_phi, double max_phi, - double min_z, double max_z) + CylindricalProfileObservable( + std::shared_ptr + transform_params, + int n_r_bins, int n_phi_bins, int n_z_bins, double min_r, double max_r, + double min_phi, double max_phi, double min_z, double max_z) : ProfileObservable(n_r_bins, n_phi_bins, n_z_bins, min_r, max_r, min_phi, max_phi, min_z, max_z), - center(center), axis(axis) {} - Utils::Vector3d center; - Utils::Vector3d axis; + transform_params(std::move(transform_params)) {} + + std::shared_ptr transform_params; }; } // Namespace Observables diff --git a/src/core/observables/CylindricalVelocityProfile.hpp b/src/core/observables/CylindricalVelocityProfile.hpp index 58dfa0bc5fc..6b70dd2feaf 100644 --- a/src/core/observables/CylindricalVelocityProfile.hpp +++ b/src/core/observables/CylindricalVelocityProfile.hpp @@ -43,11 +43,13 @@ class CylindricalVelocityProfile : public CylindricalPidProfileObservable { Utils::CylindricalHistogram histogram(n_bins, 3, limits); for (auto p : particles) { - auto const pos = folded_position(traits.position(p), box_geo) - center; + auto const pos = folded_position(traits.position(p), box_geo) - + transform_params->center(); histogram.update( - Utils::transform_coordinate_cartesian_to_cylinder(pos, axis), - Utils::transform_vector_cartesian_to_cylinder(traits.velocity(p), - axis, pos)); + Utils::transform_coordinate_cartesian_to_cylinder( + pos, transform_params->axis(), transform_params->orientation()), + Utils::transform_vector_cartesian_to_cylinder( + traits.velocity(p), transform_params->axis(), pos)); } auto hist_tmp = histogram.get_histogram(); diff --git a/src/python/espressomd/math.py b/src/python/espressomd/math.py new file mode 100644 index 00000000000..5b7ff893eb4 --- /dev/null +++ b/src/python/espressomd/math.py @@ -0,0 +1,36 @@ +# Copyright (C) 2010-2019 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from .script_interface import ScriptInterfaceHelper, script_interface_register + + +@script_interface_register +class CylindricalTransformationParameters(ScriptInterfaceHelper): + """ + Class to hold and validate the parameters needed for a cylindrical transformation. + The three parameters are available as attributes but are read-only. + + Parameters + ---------- + center : (3,) array_like of :obj:`float`, default = [0, 0, 0] + Position of the origin of the cylindrical coordinate system. + axis : (3,) array_like of :obj:`float`, default = [0, 0, 1] + Orientation vector of the ``z``-axis of the cylindrical coordinate system. + orientation: (3,) array_like of :obj:`float`, default = [1, 0, 0] + The axis on which ``phi = 0``. + """ + _so_name = "CylindricalTransformationParameters" diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index 2169bee15e9..775ea09fbce 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -17,6 +17,7 @@ import itertools import numpy as np from .script_interface import ScriptInterfaceHelper, script_interface_register +from .math import CylindricalTransformationParameters @script_interface_register @@ -69,6 +70,18 @@ def bin_centers(self): return np.array(list(itertools.product(*edges))).reshape(shape) +class CylindricalProfileObservable(ProfileObservable): + """ + Base class for observables that work with cylinder coordinates + """ + + def __init__( + self, transform_params=CylindricalTransformationParameters(), **kwargs): + # Provide default transformation parameters if not user-provided + kwargs['transform_params'] = transform_params + super().__init__(**kwargs) + + @script_interface_register class ComPosition(Observable): @@ -636,7 +649,7 @@ class DPDStress(Observable): @script_interface_register -class CylindricalDensityProfile(ProfileObservable): +class CylindricalDensityProfile(CylindricalProfileObservable): """Calculates the particle density in cylindrical coordinates. @@ -644,26 +657,24 @@ class CylindricalDensityProfile(ProfileObservable): ---------- ids : array_like of :obj:`int` The ids of (existing) particles to take into account. - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. @@ -676,7 +687,7 @@ class CylindricalDensityProfile(ProfileObservable): @script_interface_register -class CylindricalFluxDensityProfile(ProfileObservable): +class CylindricalFluxDensityProfile(CylindricalProfileObservable): """Calculates the particle flux density in cylindrical coordinates. @@ -684,26 +695,24 @@ class CylindricalFluxDensityProfile(ProfileObservable): ---------- ids : array_like of :obj:`int` The ids of (existing) particles to take into account. - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. @@ -718,7 +727,8 @@ class CylindricalFluxDensityProfile(ProfileObservable): @script_interface_register -class CylindricalLBFluxDensityProfileAtParticlePositions(ProfileObservable): +class CylindricalLBFluxDensityProfileAtParticlePositions( + CylindricalProfileObservable): """Calculates the LB fluid flux density at the particle positions in cylindrical coordinates. @@ -727,26 +737,24 @@ class CylindricalLBFluxDensityProfileAtParticlePositions(ProfileObservable): ---------- ids : array_like of :obj:`int` The ids of (existing) particles to take into account. - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. @@ -761,7 +769,8 @@ class CylindricalLBFluxDensityProfileAtParticlePositions(ProfileObservable): @script_interface_register -class CylindricalLBVelocityProfileAtParticlePositions(ProfileObservable): +class CylindricalLBVelocityProfileAtParticlePositions( + CylindricalProfileObservable): """Calculates the LB fluid velocity at the particle positions in cylindrical coordinates. @@ -770,26 +779,24 @@ class CylindricalLBVelocityProfileAtParticlePositions(ProfileObservable): ---------- ids : array_like of :obj:`int` The ids of (existing) particles to take into account. - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. @@ -804,7 +811,7 @@ class CylindricalLBVelocityProfileAtParticlePositions(ProfileObservable): @script_interface_register -class CylindricalVelocityProfile(ProfileObservable): +class CylindricalVelocityProfile(CylindricalProfileObservable): """Calculates the particle velocity profile in cylindrical coordinates. @@ -812,26 +819,24 @@ class CylindricalVelocityProfile(ProfileObservable): ---------- ids : array_like of :obj:`int` The ids of (existing) particles to take into account. - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. @@ -846,7 +851,7 @@ class CylindricalVelocityProfile(ProfileObservable): @script_interface_register -class CylindricalLBVelocityProfile(ProfileObservable): +class CylindricalLBVelocityProfile(CylindricalProfileObservable): """Calculates the LB fluid velocity profile in cylindrical coordinates. @@ -856,26 +861,24 @@ class CylindricalLBVelocityProfile(ProfileObservable): Parameters ---------- - center : (3,) array_like of :obj:`float` - Position of the center of the cylindrical coordinate system for the histogram. - axis : (3,) array_like of :obj:`float` - Orientation vector of the ``z``-axis of the cylindrical coordinate system for the histogram. - n_r_bins : :obj:`int` + transform_params : :class:`espressomd.math.CylindricalTransformationParameters`, optional + Parameters of the cylinder transformation. Defaults to the default of :class:`espressomd.math.CylindricalTransformationParameters` + n_r_bins : :obj:`int`, default = 1 Number of bins in radial direction. - n_phi_bins : :obj:`int` + n_phi_bins : :obj:`int`, default = 1 Number of bins for the azimuthal direction. - n_z_bins : :obj:`int` + n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. - min_r : :obj:`float` + min_r : :obj:`float`, default = 0 Minimum ``r`` to consider. - min_phi : :obj:`float` - Minimum ``phi`` to consider. + min_phi : :obj:`float`, default = -pi + Minimum ``phi`` to consider. Must be in [-pi,pi). min_z : :obj:`float` Minimum ``z`` to consider. max_r : :obj:`float` Maximum ``r`` to consider. - max_phi : :obj:`float` - Maximum ``phi`` to consider. + max_phi : :obj:`float`, default = pi + Maximum ``phi`` to consider. Must be in (-pi,pi]. max_z : :obj:`float` Maximum ``z`` to consider. sampling_density : :obj:`float` diff --git a/src/script_interface/CylindricalTransformationParameters.hpp b/src/script_interface/CylindricalTransformationParameters.hpp new file mode 100644 index 00000000000..89ed4b840d0 --- /dev/null +++ b/src/script_interface/CylindricalTransformationParameters.hpp @@ -0,0 +1,62 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef SCRIPT_INTERFACE_CYL_TRANSFORM_PARAMS_HPP +#define SCRIPT_INTERFACE_CYL_TRANSFORM_PARAMS_HPP + +#include "script_interface/ScriptInterface.hpp" + +#include "utils/math/cylindrical_transformation_parameters.hpp" + +namespace ScriptInterface { + +class CylindricalTransformationParameters + : public AutoParameters { +public: + CylindricalTransformationParameters() { + add_parameters({{"center", AutoParameter::read_only, + [this]() { return m_transform_params->center(); }}, + {"axis", AutoParameter::read_only, + [this]() { return m_transform_params->axis(); }}, + {"orientation", AutoParameter::read_only, + [this]() { return m_transform_params->orientation(); }}}); + } + std::shared_ptr<::Utils::CylindricalTransformationParameters> + cyl_transform_params() { + return m_transform_params; + } + void do_construct(VariantMap const ¶ms) override { + m_transform_params = + std::make_shared( + get_value_or(params, "center", + Utils::Vector3d{{0, 0, 0}}), + get_value_or(params, "axis", + Utils::Vector3d{{0, 0, 1}}), + get_value_or(params, "orientation", + Utils::Vector3d{{1, 0, 0}})); + } + +private: + std::shared_ptr + m_transform_params; +}; +} // namespace ScriptInterface +#endif diff --git a/src/script_interface/initialize.cpp b/src/script_interface/initialize.cpp index 88a7048bc4b..f971b29ae12 100644 --- a/src/script_interface/initialize.cpp +++ b/src/script_interface/initialize.cpp @@ -29,6 +29,7 @@ #include "h5md/initialize.hpp" #endif #include "ComFixed.hpp" +#include "CylindricalTransformationParameters.hpp" #include "accumulators/initialize.hpp" #include "collision_detection/initialize.hpp" #include "lbboundaries/initialize.hpp" @@ -53,6 +54,8 @@ void initialize(Utils::Factory *f) { CollisionDetection::initialize(f); f->register_new("ComFixed"); + f->register_new( + "CylindricalTransformationParameters"); } } /* namespace ScriptInterface */ diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp index f0f12e3e639..b0eafa3942a 100644 --- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp @@ -28,6 +28,8 @@ #include "core/observables/CylindricalLBProfileObservable.hpp" #include "script_interface/get_value.hpp" +#include "script_interface/CylindricalTransformationParameters.hpp" + #include #include @@ -53,18 +55,7 @@ class CylindricalLBProfileObservable using Base::Base; CylindricalLBProfileObservable() { this->add_parameters({ - {"center", - [this](const Variant &v) { - cylindrical_profile_observable()->center = - get_value<::Utils::Vector3d>(v); - }, - [this]() { return cylindrical_profile_observable()->center; }}, - {"axis", - [this](const Variant &v) { - cylindrical_profile_observable()->axis = - get_value(v); - }, - [this]() { return cylindrical_profile_observable()->axis; }}, + {"transform_params", m_transform_params}, {"n_r_bins", [this](const Variant &v) { cylindrical_profile_observable()->n_bins[0] = @@ -149,13 +140,21 @@ class CylindricalLBProfileObservable } void do_construct(VariantMap const ¶ms) override { - m_observable = - make_shared_from_args( - params, "center", "axis", "n_r_bins", "n_phi_bins", "n_z_bins", - "min_r", "max_r", "min_phi", "max_phi", "min_z", "max_z", - "sampling_density"); + set_from_args(m_transform_params, params, "transform_params"); + + if (m_transform_params) + m_observable = std::make_shared( + m_transform_params->cyl_transform_params(), + get_value_or(params, "n_r_bins", 1), + get_value_or(params, "n_phi_bins", 1), + get_value_or(params, "n_z_bins", 1), + get_value_or(params, "min_r", 0.), + get_value(params, "max_r"), + get_value_or(params, "min_phi", -Utils::pi()), + get_value_or(params, "max_phi", Utils::pi()), + get_value(params, "min_z"), + get_value(params, "max_z"), + get_value(params, "sampling_density")); } Variant do_call_method(std::string const &method, @@ -180,6 +179,7 @@ class CylindricalLBProfileObservable private: std::shared_ptr m_observable; + std::shared_ptr m_transform_params; }; } /* namespace Observables */ diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp index 9d22325af17..d541880b367 100644 --- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp @@ -27,11 +27,14 @@ #include "Observable.hpp" #include "core/observables/CylindricalPidProfileObservable.hpp" -#include +#include +#include +#include #include #include #include + #include #include @@ -58,18 +61,7 @@ class CylindricalPidProfileObservable get_value>(v); }, [this]() { return cylindrical_pid_profile_observable()->ids(); }}, - {"center", - [this](const Variant &v) { - cylindrical_pid_profile_observable()->center = - get_value<::Utils::Vector3d>(v); - }, - [this]() { return cylindrical_pid_profile_observable()->center; }}, - {"axis", - [this](const Variant &v) { - cylindrical_pid_profile_observable()->axis = - get_value(v); - }, - [this]() { return cylindrical_pid_profile_observable()->axis; }}, + {"transform_params", m_transform_params}, {"n_r_bins", [this](const Variant &v) { cylindrical_pid_profile_observable()->n_bins[0] = @@ -149,13 +141,21 @@ class CylindricalPidProfileObservable }; void do_construct(VariantMap const ¶ms) override { - m_observable = - make_shared_from_args, Utils::Vector3d, - Utils::Vector3d, int, int, int, double, double, - double, double, double, double>( - params, "ids", "center", "axis", "n_r_bins", "n_phi_bins", - "n_z_bins", "min_r", "max_r", "min_phi", "max_phi", "min_z", - "max_z"); + set_from_args(m_transform_params, params, "transform_params"); + + if (m_transform_params) + m_observable = std::make_shared( + get_value>(params, "ids"), + m_transform_params->cyl_transform_params(), + get_value_or(params, "n_r_bins", 1), + get_value_or(params, "n_phi_bins", 1), + get_value_or(params, "n_z_bins", 1), + get_value_or(params, "min_r", 0.), + get_value(params, "max_r"), + get_value_or(params, "min_phi", -Utils::pi()), + get_value_or(params, "max_phi", Utils::pi()), + get_value(params, "min_z"), + get_value(params, "max_z")); } Variant do_call_method(std::string const &method, @@ -180,6 +180,7 @@ class CylindricalPidProfileObservable private: std::shared_ptr m_observable; + std::shared_ptr m_transform_params; }; } /* namespace Observables */ diff --git a/src/shapes/include/shapes/HollowConicalFrustum.hpp b/src/shapes/include/shapes/HollowConicalFrustum.hpp index ac846e0e545..1deb0d98f2e 100644 --- a/src/shapes/include/shapes/HollowConicalFrustum.hpp +++ b/src/shapes/include/shapes/HollowConicalFrustum.hpp @@ -22,6 +22,9 @@ #include "Shape.hpp" #include +#include + +#include namespace Shapes { @@ -48,15 +51,20 @@ class HollowConicalFrustum : public Shape { HollowConicalFrustum() : m_r1(0.0), m_r2(0.0), m_length(0.0), m_thickness(0.0), m_direction(1), m_center{Utils::Vector3d{}}, m_axis{Utils::Vector3d{ - 0, 0, 1}} {} - - void set_r1(double radius) { m_r1 = radius; } - void set_r2(double radius) { m_r2 = radius; } - void set_length(double length) { m_length = length; } - void set_thickness(double thickness) { m_thickness = thickness; } - void set_direction(int dir) { m_direction = dir; } - void set_axis(Utils::Vector3d const &axis) { m_axis = axis; } + 0., 0., 1.}}, + m_orientation{Utils::Vector3d{1., 0., 0.}} {} + void set_r1(double const radius) { m_r1 = radius; } + void set_r2(double const radius) { m_r2 = radius; } + void set_length(double const length) { m_length = length; } + void set_thickness(double const thickness) { m_thickness = thickness; } + void set_direction(int const dir) { m_direction = dir; } + void set_axis(Utils::Vector3d const &axis) { + m_axis = axis; + // Even though the HCF is cylinder-symmetric, it needs a well defined phi=0 + // orientation for the coordinate transformation. + m_orientation = Utils::calc_orthonormal_vector(axis); + } void set_center(Utils::Vector3d const ¢er) { m_center = center; } /// Get radius 1 perpendicular to axis. @@ -92,6 +100,7 @@ class HollowConicalFrustum : public Shape { int m_direction; Utils::Vector3d m_center; Utils::Vector3d m_axis; + Utils::Vector3d m_orientation; }; } // namespace Shapes diff --git a/src/shapes/src/HollowConicalFrustum.cpp b/src/shapes/src/HollowConicalFrustum.cpp index 3db5fb59fc0..ea384252cc4 100644 --- a/src/shapes/src/HollowConicalFrustum.cpp +++ b/src/shapes/src/HollowConicalFrustum.cpp @@ -33,8 +33,9 @@ void HollowConicalFrustum::calculate_dist(const Utils::Vector3d &pos, Utils::Vector3d &vec) const { // transform given position to cylindrical coordinates in the reference frame // of the cone - auto const pos_cyl = - Utils::transform_coordinate_cartesian_to_cylinder(pos - m_center, m_axis); + auto const v = pos - m_center; + auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( + v, m_axis, m_orientation); // clang-format off /* * the following implementation is based on: @@ -61,7 +62,7 @@ void HollowConicalFrustum::calculate_dist(const Utils::Vector3d &pos, // Transform back to cartesian coordinates. auto const pos_intersection = Utils::transform_coordinate_cylinder_to_cartesian( - {r_intersection, pos_cyl[1], z_intersection}, m_axis) + + {r_intersection, pos_cyl[1], z_intersection}, m_axis, m_orientation) + m_center; auto const u = (pos - pos_intersection).normalize(); diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index db31193b062..07f84265f8b 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -19,64 +19,141 @@ #ifndef UTILS_COORDINATE_TRANSFORMATION_HPP #define UTILS_COORDINATE_TRANSFORMATION_HPP +/** + * @file + * Convert coordinates from the Cartesian system to the cylindrical system. + * The transformation functions are provided with three overloads: + * - one function for the trivial Cartesian <-> cylindrical transformation + * - one function to transform from/to a cylindrical system with custom axis + * (extra @p axis argument, keep in mind the angle phi is under-defined) + * - one function to transform from/to an oriented cylindrical system with + * custom axis (extra @p orientation argument, the angle phi is well-defined) + */ + #include "utils/Vector.hpp" #include "utils/constants.hpp" #include "utils/math/vec_rotate.hpp" +#include "utils/matrix.hpp" #include "utils/quaternion.hpp" +#include +#include + namespace Utils { -/** \brief Transform the given 3D position to cylinder coordinates with - * longitudinal axis aligned with axis parameter. +/** + * @brief Basis change. */ -inline Vector3d -transform_coordinate_cartesian_to_cylinder(const Vector3d &pos, - const Vector3d &axis) { - static auto const z_axis = Vector3d{{0, 0, 1}}; - double theta; - Vector3d rotation_axis; - auto r = [](auto const &pos) { - return std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]); - }; - auto phi = [](auto const &pos) { return std::atan2(pos[1], pos[0]); }; - if (axis != z_axis) { - std::tie(theta, rotation_axis) = rotation_params(axis, z_axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, pos); - return {r(rotated_pos), phi(rotated_pos), rotated_pos[2]}; +inline Vector3d basis_change(Vector3d const &b1, Vector3d const &b2, + Vector3d const &b3, Vector3d const &v, + bool reverse = false) { + auto const e_x = b1.normalized(); + auto const e_y = b2.normalized(); + auto const e_z = b3.normalized(); + auto const M = Matrix{ + {e_x[0], e_x[1], e_x[2]}, + {e_y[0], e_y[1], e_y[2]}, + {e_z[0], e_z[1], + e_z[2]}}.transposed(); + if (reverse) { + return M * v; } - return {r(pos), phi(pos), pos[2]}; + return M.inversed() * v; } /** - * @brief Coordinate transformation from cylinder to cartesian coordinates. + * @brief Coordinate transformation from Cartesian to cylindrical coordinates. + * The origins and z-axis of the coordinate systems co-incide. + * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the + * original coordinate system. + * @param pos %Vector to transform */ inline Vector3d -transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, - Vector3d const &axis) { - Vector3d const transformed{ - {pos[0] * std::cos(pos[1]), pos[0] * std::sin(pos[1]), pos[2]}}; - static auto const z_axis = Vector3d{{0, 0, 1}}; - if (axis == z_axis) - return transformed; - double theta; - Vector3d rotation_axis; - std::tie(theta, rotation_axis) = rotation_params(z_axis, axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, transformed); - return rotated_pos; +transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) { + auto const r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]); + auto const phi = std::atan2(pos[1], pos[0]); + return {r, phi, pos[2]}; +} + +/** + * @brief Coordinate transformation from Cartesian to cylindrical coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * + * If the parameter @p axis is not equal to [0, 0, 1], the value + * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined. + * To fully define it, it is necessary to provide an orientation vector + * in Cartesian coordinates that will be used as the reference point + * (i.e. such that @f$ \phi = 0 @f$), by default it is the x-axis. + * + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param orientation Reference point (in untransformed coordinates) for + * which @f$ \phi = 0 @f$ + */ +inline Vector3d transform_coordinate_cartesian_to_cylinder( + Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) { + // check that axis and orientation are orthogonal + assert(std::abs(axis * orientation) < + 5 * std::numeric_limits::epsilon()); + auto const rotation_axis = vector_product(axis, orientation); + auto const pos_t = basis_change(orientation, rotation_axis, axis, pos); + return transform_coordinate_cartesian_to_cylinder(pos_t); } -/** \brief Transform the given 3D vector to cylinder coordinates with - * symmetry axis aligned with axis parameter. +/** + * @brief Coordinate transformation from cylindrical to Cartesian coordinates. + * The origins and z-axis of the coordinate systems co-incide. + * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the + * transformed coordinate system. + * @param pos %Vector to transform + */ +inline Vector3d +transform_coordinate_cylinder_to_cartesian(Vector3d const &pos) { + auto const &rho = pos[0]; + auto const &phi = pos[1]; + auto const &z = pos[2]; + return {rho * std::cos(phi), rho * std::sin(phi), z}; +} + +/** + * @brief Coordinate transformation from cylindrical to Cartesian coordinates + * with change of basis. The origins of the coordinate systems co-incide. + * + * If the parameter @p axis is not equal to [0, 0, 1], the value + * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined. + * To fully define it, it is necessary to provide an orientation vector + * in Cartesian coordinates that will be used as the reference point + * (i.e. such that @f$ \phi = 0 @f$). + * + * @param pos %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param orientation Reference point (in Cartesian coordinates) for + * which @f$ \phi = 0 @f$ + */ +inline Vector3d transform_coordinate_cylinder_to_cartesian( + Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) { + // check that axis and orientation are orthogonal + assert(std::abs(axis * orientation) < + 5 * std::numeric_limits::epsilon()); + auto const rotation_axis = vector_product(axis, orientation); + auto const pos_t = transform_coordinate_cylinder_to_cartesian(pos); + return basis_change(orientation, rotation_axis, axis, pos_t, true); +} + +/** + * @brief Vector transformation from Cartesian to cylindrical coordinates. + * @param vec %Vector to transform + * @param axis Longitudinal axis of the cylindrical coordinates + * @param pos Origin of the vector */ inline Vector3d transform_vector_cartesian_to_cylinder(Vector3d const &vec, Vector3d const &axis, Vector3d const &pos) { static auto const z_axis = Vector3d{{0, 0, 1}}; - double theta; - Vector3d rotation_axis; - std::tie(theta, rotation_axis) = rotation_params(axis, z_axis); - auto const rotated_pos = vec_rotate(rotation_axis, theta, pos); - auto const rotated_vec = vec_rotate(rotation_axis, theta, vec); + auto const angle = angle_between(axis, z_axis); + auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize(); + auto const rotated_pos = vec_rotate(rotation_axis, angle, pos); + auto const rotated_vec = vec_rotate(rotation_axis, angle, vec); auto const r = std::sqrt(rotated_pos[0] * rotated_pos[0] + rotated_pos[1] * rotated_pos[1]); // v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2) diff --git a/src/utils/include/utils/math/cylindrical_transformation_parameters.hpp b/src/utils/include/utils/math/cylindrical_transformation_parameters.hpp new file mode 100644 index 00000000000..20f18d78a77 --- /dev/null +++ b/src/utils/include/utils/math/cylindrical_transformation_parameters.hpp @@ -0,0 +1,80 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef ESPRESSO_CYLINDER_TRANSFORMATION_PARAMETERS_HPP +#define ESPRESSO_CYLINDER_TRANSFORMATION_PARAMETERS_HPP + +#include +#include + +#include + +namespace Utils { + +/** + * @brief A class to hold and validate parameters for a cylindrical coordinate + * transformations. + * + * @param center The origin of the cylindrical coordinates. + * @param axis The "z"-axis. Must be normalized. + * @param orientation The axis along which phi = 0. Must be normalized and + * orthogonal to axis. + */ +class CylindricalTransformationParameters { +public: + CylindricalTransformationParameters() = default; + CylindricalTransformationParameters(Utils::Vector3d const ¢er, + Utils::Vector3d const &axis, + Utils::Vector3d const &orientation) + : m_center(center), m_axis(axis), m_orientation(orientation) { + validate(); + } + + Utils::Vector3d center() const { return m_center; } + Utils::Vector3d axis() const { return m_axis; } + Utils::Vector3d orientation() const { return m_orientation; } + +private: + void validate() const { + auto constexpr eps = 10 * std::numeric_limits::epsilon(); + if (Utils::abs(m_orientation * m_axis) > eps) { + throw std::runtime_error( + "CylindricalTransformationParameters: Axis and orientation must be " + "orthogonal. Scalar product is " + + std::to_string(m_orientation * m_axis)); + } + if (Utils::abs(m_axis.norm() - 1) > eps) { + throw std::runtime_error("CylindricalTransformationParameters: Axis must " + "be normalized. Norm is " + + std::to_string(m_axis.norm())); + } + if (Utils::abs(m_orientation.norm() - 1) > eps) { + throw std::runtime_error("CylindricalTransformationParameters: " + "orientation must be normalized. Norm is " + + std::to_string(m_orientation.norm())); + } + } + + const Utils::Vector3d m_center{}; + const Utils::Vector3d m_axis{0, 0, 1}; + const Utils::Vector3d m_orientation{1, 0, 0}; +}; + +} // namespace Utils + +#endif // ESPRESSO_CYLINDER_TRANSFORMATION_PARAMETERS_HPP diff --git a/src/utils/include/utils/math/orthonormal_vec.hpp b/src/utils/include/utils/math/orthonormal_vec.hpp new file mode 100644 index 00000000000..57f2637fab4 --- /dev/null +++ b/src/utils/include/utils/math/orthonormal_vec.hpp @@ -0,0 +1,54 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef ESPRESSO_ORTHONORMAL_VEC_HPP +#define ESPRESSO_ORTHONORMAL_VEC_HPP + +#include "utils/Vector.hpp" +#include "utils/constants.hpp" + +namespace Utils { +/** + * @brief Return a vector that is orthonormal to vec + */ +template +Vector calc_orthonormal_vector(Vector const &vec) { + /* Calculate orthonormal vector using Gram-Schmidt orthogonalization of a + trial vector. Only works if the trial vector is not parallel, so we have to + try a second one in that case + */ + Vector, 2> try_vectors = {Vector::broadcast(0), + Vector::broadcast(0)}; + try_vectors[0][0] = 1; + try_vectors[1][1] = 1; + + Vector ret; + for (auto v : try_vectors) { + auto orth_component = v - (v * vec) / vec.norm2() * vec; + auto norm = orth_component.norm(); + if (norm >= 1. / Utils::sqrt_2()) { + ret = orth_component / norm; + break; + } + } + return ret; +} + +} // namespace Utils + +#endif // ESPRESSO_ORTHONORMAL_VEC_HPP \ No newline at end of file diff --git a/src/utils/include/utils/math/vec_rotate.hpp b/src/utils/include/utils/math/vec_rotate.hpp index 32c1cc23d8d..029852ddd56 100644 --- a/src/utils/include/utils/math/vec_rotate.hpp +++ b/src/utils/include/utils/math/vec_rotate.hpp @@ -47,21 +47,10 @@ inline Vector3d vec_rotate(const Vector3d &axis, double angle, } /** - * @brief Determine rotation angle and axis for rotating vec onto target_vec. - * @param vec Vector to be rotated - * @param target_vec Target vector - * @return rotation angle and rotation axis + * @brief Determine the angle between two vectors. */ -inline std::tuple -rotation_params(Vector3d const &vec, Vector3d const &target_vec) { - if (vec.normalized() != target_vec.normalized()) { - auto const theta = - std::acos(vec * target_vec / (vec.norm() * target_vec.norm())); - auto const rotation_axis = - Utils::vector_product(vec, target_vec).normalize(); - return std::make_tuple(theta, rotation_axis); - } - return std::make_tuple(0.0, Vector3d{}); +inline double angle_between(Vector3d const &v1, Vector3d const &v2) { + return std::acos(v1 * v2 / std::sqrt(v1.norm2() * v2.norm2())); } } // namespace Utils diff --git a/src/utils/tests/CMakeLists.txt b/src/utils/tests/CMakeLists.txt index 90508ce2539..a230d4b84ba 100644 --- a/src/utils/tests/CMakeLists.txt +++ b/src/utils/tests/CMakeLists.txt @@ -78,3 +78,5 @@ unit_test(NAME sendrecv_test SRC sendrecv_test.cpp DEPENDS EspressoUtils Boost::mpi MPI::MPI_CXX EspressoUtils NUM_PROC 3) unit_test(NAME matrix_test SRC matrix_test.cpp DEPENDS EspressoUtils Boost::serialization NUM_PROC 1) +unit_test(NAME orthonormal_vec_test SRC orthonormal_vec_test.cpp DEPENDS + EspressoUtils Boost::serialization NUM_PROC 1) diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index b5c5d19e2ee..e70709718b6 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -25,64 +25,216 @@ #include #include +#include using Utils::Vector3d; BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) { - Vector3d const cart_coord{{1.0, 3.3, 2.0}}; - auto const transformed_x = transform_coordinate_cartesian_to_cylinder( - cart_coord, Vector3d{{1, 0, 0}}); - auto const transformed_y = transform_coordinate_cartesian_to_cylinder( - cart_coord, Vector3d{{0, 1, 0}}); - auto const transformed_z = transform_coordinate_cartesian_to_cylinder( - cart_coord, Vector3d{{0, 0, 1}}); - // For x as the symmetry axis we rotate the cartesian coordinates around the - // y-axis by -pi/2. - auto const expected_x = transform_coordinate_cartesian_to_cylinder( - vec_rotate(Vector3d{{0.0, 1.0, 0.0}}, -Utils::pi() / 2.0, cart_coord), - Vector3d{{0, 0, 1}}); - // For y as the symmetry axis we rotate the cartesian coordinates around the - // x-axis by pi/2. - auto const expected_y = transform_coordinate_cartesian_to_cylinder( - vec_rotate(Vector3d{{1.0, 0.0, 0.0}}, Utils::pi() / 2.0, cart_coord), - Vector3d{{0, 0, 1}}); - auto const expected_z = Vector3d{ - {std::sqrt(cart_coord[0] * cart_coord[0] + cart_coord[1] * cart_coord[1]), - std::atan2(cart_coord[1], cart_coord[0]), cart_coord[2]}}; + constexpr auto eps = 1e-14; + auto const pos = Vector3d{{1.0, 3.3, 2.0}}; + auto const cyl = transform_coordinate_cartesian_to_cylinder(pos); + BOOST_CHECK_SMALL(cyl[0] - std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]), eps); + BOOST_CHECK_SMALL(cyl[1] - std::atan2(pos[1], pos[0]), eps); + BOOST_CHECK_SMALL(cyl[2] - pos[2], eps); +} + +BOOST_AUTO_TEST_CASE(basis_transform_test) { + constexpr auto eps = 1e-14; + Vector3d const b_x{{1, 0, 0}}; + Vector3d const b_y{{0, 1, 0}}; + Vector3d const b_z{{0, 0, 1}}; + // identity transform + Vector3d const v{{1, 2, 3}}; + Vector3d const v_identity_transform = Utils::basis_change(b_x, b_y, b_z, v); + // identity transform (swap both the vector and coordinate system) + Vector3d const v_swap_coord_transform = + Utils::basis_change(b_z, b_y, b_x, {{v[2], v[1], v[0]}}); + // non-trivial transform + Vector3d const v1 = Vector3d{{2, 2, 2}}.normalized(); + Vector3d const v2 = Vector3d{{3, 3, -6}}.normalized(); + Vector3d const v3 = Utils::vector_product(v1, v2).normalized(); + Vector3d const v4 = basis_change(v1, v2, v3, 0.1 * v1 + 0.2 * v2 - 0.3 * v3); + Vector3d const v4_expected = Vector3d{{0.1, 0.2, -0.3}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(v_identity_transform[i] - v[i], eps); + BOOST_CHECK_SMALL(v_swap_coord_transform[i] - v[i], eps); + BOOST_CHECK_SMALL(v4[i] - v4_expected[i], eps); + } +} +BOOST_AUTO_TEST_CASE( + transform_coordinate_cartesian_to_cylinder_base_change_test) { + constexpr auto eps = 1e-14; + Vector3d const v1{{1, 3, 4}}; + Vector3d const v2{{-3.0, 7, 2}}; + Vector3d const axis = Utils::vector_product(v1, v2).normalized(); + Vector3d const v3 = + Utils::transform_coordinate_cartesian_to_cylinder(v1, axis, v1); + Vector3d const v3_ref{{v1.norm(), 0, 0}}; + auto const angle_v1_v2 = Utils::angle_between(v1, v2); + auto const v4 = Utils::transform_coordinate_cartesian_to_cylinder( + v2 + 2 * axis, axis, v1); + Vector3d v4_ref{{v2.norm(), angle_v1_v2, 2}}; + auto const v5 = Utils::transform_coordinate_cartesian_to_cylinder( + v1 + 2 * axis, axis, v2); + Vector3d v5_ref{{v1.norm(), -angle_v1_v2, 2}}; for (int i = 0; i < 3; ++i) { - BOOST_CHECK(transformed_x[i] == expected_x[i]); - BOOST_CHECK(transformed_y[i] == expected_y[i]); - BOOST_CHECK(transformed_z[i] == expected_z[i]); + BOOST_CHECK_SMALL(v3[i] - v3_ref[i], eps); + BOOST_CHECK_SMALL(v4[i] - v4_ref[i], eps); + BOOST_CHECK_SMALL(v5[i] - v5_ref[i], eps); + } +} + +BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_and_orientation_test) { + constexpr auto eps = 1e-14; + // tilted orthogonal basis + auto const y = (Vector3d{{0, 1, -1}}).normalize(); + auto const z = (Vector3d{{1, 1, 1}}).normalize(); + auto const x = Utils::vector_product(y, z); + + // check transformation with orientation (phi is random for r=0) + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y); + auto const x_ref = Vector3d{{1.0, -Utils::pi() / 2.0, 0.0}}; + auto const y_ref = Vector3d{{1.0, 0.0, 0.0}}; + auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(x_cyl[i] - x_ref[i], eps); + BOOST_CHECK_SMALL(y_cyl[i] - y_ref[i], eps); + BOOST_CHECK_SMALL(z_cyl[i] - z_ref[i], eps); + } + } + // check transformation with orientation for another angle + { + auto const u = vec_rotate(z, Utils::pi() / 3.0, x); + auto const v = vec_rotate(z, Utils::pi() / 3.0, y); + auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y); + auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y); + auto const u_ref = Vector3d{{1.0, Utils::pi() * (1. / 3. - 1. / 2.), 0.0}}; + auto const v_ref = Vector3d{{1.0, Utils::pi() / 3.0, 0.0}}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(u_cyl[i] - u_ref[i], eps); + BOOST_CHECK_SMALL(v_cyl[i] - v_ref[i], eps); + } + } + // check transformation of random vectors + { + std::subtract_with_carry_engine rng(2); + auto const r_uniform = [&rng]() { + return static_cast(rng() - rng.min()) / (rng.max() - rng.min()); + }; + for (int trial = 0; trial < 100; ++trial) { + Vector3d const v1{r_uniform(), r_uniform(), r_uniform()}; + Vector3d const v2{r_uniform(), r_uniform(), r_uniform()}; + auto const a = Utils::vector_product(v1, v2) / v1.norm() / v2.norm(); + auto const v1_v1 = transform_coordinate_cartesian_to_cylinder(v1, a, v1); + auto const v2_v1 = transform_coordinate_cartesian_to_cylinder(v2, a, v1); + auto const v1_v2 = transform_coordinate_cartesian_to_cylinder(v1, a, v2); + Vector3d const v1_v1_ref{v1.norm(), 0.0, 0.0}; + Vector3d const v2_v1_ref{v2.norm(), Utils::angle_between(v1, v2), 0.0}; + Vector3d const v1_v2_ref{v1.norm(), -Utils::angle_between(v1, v2), 0.0}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(v1_v1[i] - v1_v1_ref[i], eps); + BOOST_CHECK_SMALL(v2_v1[i] - v2_v1_ref[i], eps); + BOOST_CHECK_SMALL(v1_v2[i] - v1_v2_ref[i], eps); + } + } } } BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) { + constexpr auto eps = 1e-14; + auto const cyl = Vector3d{{1.0, Utils::pi() / 4, 2.0}}; + auto const pos = transform_coordinate_cylinder_to_cartesian(cyl); + BOOST_CHECK_SMALL(pos[0] - std::sqrt(2) / 2, eps); + BOOST_CHECK_SMALL(pos[1] - std::sqrt(2) / 2, eps); + BOOST_CHECK_SMALL(pos[2] - cyl[2], eps); +} + +BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_and_orientation_test) { + constexpr auto eps = 2e-14; Vector3d const cylinder_coord{{1.2, 3.123, 42.0}}; - auto const transformed_x = transform_coordinate_cylinder_to_cartesian( - cylinder_coord, Vector3d{{1, 0, 0}}); - auto const transformed_y = transform_coordinate_cylinder_to_cartesian( - cylinder_coord, Vector3d{{0, 1, 0}}); - auto const transformed_z = transform_coordinate_cylinder_to_cartesian( - cylinder_coord, Vector3d{{0, 0, 1}}); + auto const e_x = Vector3d{{1., 0., 0.}}; + auto const e_y = Vector3d{{0., 1., 0.}}; + auto const e_z = Vector3d{{0., 0., 1.}}; + + auto const transformed_x = + transform_coordinate_cylinder_to_cartesian(cylinder_coord, e_x, -e_z); + auto const transformed_y = + transform_coordinate_cylinder_to_cartesian(cylinder_coord, e_y, e_x); + auto const transformed_z = + transform_coordinate_cylinder_to_cartesian(cylinder_coord, e_z, e_x); // We transform from cylinder zu cartesian and have to rotate back. See test // cartesian_to_cylinder_test. - auto const expected_x = - vec_rotate(Vector3d{{0.0, 1.0, 0.0}}, Utils::pi() / 2.0, - transform_coordinate_cylinder_to_cartesian( - cylinder_coord, Vector3d{{0, 0, 1}})); - auto const expected_y = - vec_rotate(Vector3d{{1.0, 0.0, 0.0}}, -Utils::pi() / 2.0, - transform_coordinate_cylinder_to_cartesian( - cylinder_coord, Vector3d{{0, 0, 1}})); + auto const expected_x = vec_rotate( + e_y, Utils::pi() / 2.0, + transform_coordinate_cylinder_to_cartesian(cylinder_coord, e_z, e_x)); + auto const expected_y = vec_rotate( + e_x, -Utils::pi() / 2.0, + transform_coordinate_cylinder_to_cartesian(cylinder_coord, e_z, e_x)); // x = r * cos(phi); y = r * sin(phi); z = z auto const expected_z = Vector3d{ {cylinder_coord[0] * std::cos(cylinder_coord[1]), cylinder_coord[0] * std::sin(cylinder_coord[1]), cylinder_coord[2]}}; for (int i = 0; i < 3; ++i) { - BOOST_CHECK(transformed_x[i] == expected_x[i]); - BOOST_CHECK(transformed_y[i] == expected_y[i]); - BOOST_CHECK(transformed_z[i] == expected_z[i]); + BOOST_CHECK_SMALL(transformed_x[i] - expected_x[i], eps); + BOOST_CHECK_SMALL(transformed_y[i] - expected_y[i], eps); + BOOST_CHECK_SMALL(transformed_z[i] - expected_z[i], eps); + } +} + +BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_with_phi_2_test) { + constexpr auto eps = 1e-14; + // tilted orthogonal basis + auto const y = (Vector3d{{0, 1, -1}}).normalize(); + auto const z = (Vector3d{{1, 1, 1}}).normalize(); + auto const x = Utils::vector_product(y, z); + + // check transformation with orientation + { + auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y); + auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y); + auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y); + auto const x_cart = transform_coordinate_cylinder_to_cartesian(x_cyl, z, y); + auto const y_cart = transform_coordinate_cylinder_to_cartesian(y_cyl, z, y); + auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z, y); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(x_cart[i] - x[i], eps); + BOOST_CHECK_SMALL(y_cart[i] - y[i], eps); + BOOST_CHECK_SMALL(z_cart[i] - z[i], eps); + } + } + // check transformation with orientation for another angle + { + auto const u = vec_rotate(z, Utils::pi() / 3.0, x); + auto const v = vec_rotate(z, Utils::pi() / 3.0, y); + auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y); + auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y); + auto const u_cart = transform_coordinate_cylinder_to_cartesian(u_cyl, z, y); + auto const v_cart = transform_coordinate_cylinder_to_cartesian(v_cyl, z, y); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(u_cart[i] - u[i], eps); + BOOST_CHECK_SMALL(v_cart[i] - v[i], eps); + } + } + // check transformation of random vectors + { + std::subtract_with_carry_engine rng(2); + auto const r_uniform = [&rng]() { + return static_cast(rng() - rng.min()) / (rng.max() - rng.min()); + }; + for (int trial = 0; trial < 100; ++trial) { + Vector3d const v1{r_uniform(), r_uniform(), r_uniform()}; + Vector3d const v2{r_uniform(), r_uniform(), r_uniform()}; + auto const a = Utils::vector_product(v1, v2) / v1.norm() / v2.norm(); + auto const v3 = transform_coordinate_cartesian_to_cylinder(v2, a, v1); + auto const v4 = transform_coordinate_cylinder_to_cartesian(v3, a, v1); + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(v4[i] - v2[i], eps); + } + } } } diff --git a/src/utils/tests/matrix_test.cpp b/src/utils/tests/matrix_test.cpp index 8030490d2ba..d3a9533cd2f 100644 --- a/src/utils/tests/matrix_test.cpp +++ b/src/utils/tests/matrix_test.cpp @@ -1,6 +1,8 @@ /* * Copyright (C) 2018-2019 The ESPResSo project * + * This file is part of ESPResSo. + * * ESPResSo is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or diff --git a/src/utils/tests/orthonormal_vec_test.cpp b/src/utils/tests/orthonormal_vec_test.cpp new file mode 100644 index 00000000000..51f9679dfe2 --- /dev/null +++ b/src/utils/tests/orthonormal_vec_test.cpp @@ -0,0 +1,39 @@ +/* + * Copyright (C) 2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#define BOOST_TEST_MODULE Utils::orthonormal_vec test +#define BOOST_TEST_DYN_LINK +#include + +#include +#include + +BOOST_AUTO_TEST_CASE(orthonormal_vec_test) { + constexpr auto eps = 1e-14; + + auto const v0 = Utils::Vector3d{{1.1, -2.2, 3.3}}; + auto v0_orth = Utils::calc_orthonormal_vector(v0); + BOOST_CHECK_SMALL(v0 * v0_orth, eps); + BOOST_CHECK_SMALL(1 - v0_orth.norm(), eps); + + auto const v1 = Utils::VectorXd<2>{{1., 0.}}; + auto v1_orth = Utils::calc_orthonormal_vector(v1); + BOOST_CHECK_SMALL(v1 * v1_orth, eps); + BOOST_CHECK_SMALL(1 - v1_orth.norm(), eps); +} \ No newline at end of file diff --git a/src/utils/tests/vec_rotate_test.cpp b/src/utils/tests/vec_rotate_test.cpp index 18b9b14ef45..f1d5727a2cd 100644 --- a/src/utils/tests/vec_rotate_test.cpp +++ b/src/utils/tests/vec_rotate_test.cpp @@ -22,11 +22,9 @@ #include #include #include -using Utils::vec_rotate; #include #include -#include BOOST_AUTO_TEST_CASE(rotation) { using std::cos; @@ -43,20 +41,16 @@ BOOST_AUTO_TEST_CASE(rotation) { auto const expected = cos(t) * v + sin(t) * vector_product(k, v) + (1. - cos(t)) * (k * v) * k; - auto const is = vec_rotate(k, t, v); + auto const is = Utils::vec_rotate(k, t, v); auto const rel_diff = (expected - is).norm() / expected.norm(); BOOST_CHECK(rel_diff < std::numeric_limits::epsilon()); } -BOOST_AUTO_TEST_CASE(rotation_params) { - Utils::Vector3d v1 = {1.0, 0.0, 0.0}; - Utils::Vector3d v2 = {1.0, 1.0, 0.0}; +BOOST_AUTO_TEST_CASE(angle_between) { + Utils::Vector3d const v1 = {1.0, 0.0, 0.0}; + Utils::Vector3d const v2 = {1.0, 1.0, 0.0}; - double angle; - Utils::Vector3d rotation_axis; - std::tie(angle, rotation_axis) = Utils::rotation_params(v1, v2); + auto const angle = Utils::angle_between(v1, v2); BOOST_CHECK_CLOSE(angle, Utils::pi() / 4.0, 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v1), 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v2), 1e-7); } diff --git a/testsuite/python/lb_poiseuille_cylinder.py b/testsuite/python/lb_poiseuille_cylinder.py index 981c0c23e28..0607a60641f 100644 --- a/testsuite/python/lb_poiseuille_cylinder.py +++ b/testsuite/python/lb_poiseuille_cylinder.py @@ -18,6 +18,7 @@ import unittest_decorators as utx import numpy as np +import espressomd.math import espressomd.lb import espressomd.lbboundaries import espressomd.observables @@ -81,7 +82,8 @@ class LBPoiseuilleCommon: system = espressomd.System(box_l=[BOX_L] * 3) system.time_step = TIME_STEP system.cell_system.skin = 0.4 * AGRID - params = {'axis': [0, 0, 1]} + params = {'axis': [0, 0, 1], + 'orientation': [1, 0, 0]} def prepare(self): """ @@ -150,8 +152,10 @@ def prepare_obs(self): else: obs_center = [BOX_L / 2.0, BOX_L / 2.0, 0.0] local_obs_params = OBS_PARAMS.copy() - local_obs_params['center'] = obs_center - local_obs_params['axis'] = self.params['axis'] + ctp = espressomd.math.CylindricalTransformationParameters(center=obs_center, + axis=self.params['axis'], + orientation=self.params['orientation']) + local_obs_params['transform_params'] = ctp obs = espressomd.observables.CylindricalLBVelocityProfile( **local_obs_params) self.accumulator = espressomd.accumulators.MeanVarianceCalculator( @@ -178,16 +182,19 @@ def check_observable(self): def test_x(self): self.params['axis'] = [1, 0, 0] + self.params['orientation'] = [0, 0, -1] self.compare_to_analytical() self.check_observable() def test_y(self): self.params['axis'] = [0, 1, 0] + self.params['orientation'] = [1, 0, 0] self.compare_to_analytical() self.check_observable() def test_z(self): self.params['axis'] = [0, 0, 1] + self.params['orientation'] = [1, 0, 0] self.compare_to_analytical() self.check_observable() diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py index 2e6e7a79ec2..3d4cefbf0ce 100644 --- a/testsuite/python/observable_cylindrical.py +++ b/testsuite/python/observable_cylindrical.py @@ -18,6 +18,7 @@ import unittest as ut import espressomd import espressomd.observables +import espressomd.math import tests_common @@ -31,13 +32,15 @@ class TestCylindricalObservable(ut.TestCase): system.time_step = 0.01 system.cell_system.skin = 0.4 + cyl_transform_params = espressomd.math.CylindricalTransformationParameters( + center=3 * [7.5], axis=[1 / np.sqrt(2), 1 / np.sqrt(2), 0], orientation=[0, 0, 1]) + params = { - 'ids': list(range(100)), - 'center': [7.5, 7.5, 7.5], # center of the histogram - 'axis': 'y', - 'n_r_bins': 4, # number of bins in r - 'n_phi_bins': 4, # -*- in phi - 'n_z_bins': 4, # -*- in z + 'ids': None, + 'transform_params': cyl_transform_params, + 'n_r_bins': 4, + 'n_phi_bins': 3, + 'n_z_bins': 4, 'min_r': 0.0, 'min_phi': -np.pi, 'min_z': -5.0, @@ -46,180 +49,162 @@ class TestCylindricalObservable(ut.TestCase): 'max_z': 5.0, } + v_r = 0.6 + v_phi = 0.7 + v_z = 0.8 + def tearDown(self): self.system.part.clear() - def swap_axis(self, arr, axis): - if axis == 'x': - arr = np.dot( - tests_common.rotation_matrix([0, 1, 0], np.pi / 2.0), arr) - elif axis == 'y': - arr = np.dot( - tests_common.rotation_matrix([1, 0, 0], -np.pi / 2.0), arr) - return arr - - def swap_axis_inverse(self, arr, axis): - if axis == 'x': - arr = np.dot( - tests_common.rotation_matrix([0, 1, 0], -np.pi / 2.0), arr) - elif axis == 'y': - arr = np.dot( - tests_common.rotation_matrix([1, 0, 0], np.pi / 2.0), arr) - return arr - - def pol_coords(self): - positions = np.zeros((len(self.params['ids']), 3)) - for i, p in enumerate(self.system.part): - tmp = p.pos - np.array(self.params['center']) - tmp = self.swap_axis_inverse(tmp, self.params['axis']) - positions[ - i, :] = tests_common.transform_pos_from_cartesian_to_polar_coordinates(tmp) - return positions - - def set_particles(self): - self.system.part.clear() - # Parameters for an ellipse. - a = 1.0 # semi minor-axis length - b = 2.0 # semi major-axis length - # Choose the cartesian velocities such that each particle gets the same - # v_r, v_phi and v_z, respectively. - self.v_r = .75 - self.v_phi = 2.5 - self.v_z = 1.5 - for i in self.params['ids']: + def calc_ellipsis_pos_vel( + self, n_part, z_min, z_max, semi_x=1., semi_y=1.): + """ + Calculate positions on an elliptical corkscrew line. + Calculate cartesian velocities that lead to a + constant velocity in cylindrical coordinates + """ + + zs = np.linspace(z_min, z_max, num=n_part) + angles = np.linspace(-0.99 * np.pi, 0.999 * np.pi, num=n_part) + + positions = [] + velocities = [] + + for angle, z in zip(angles, zs): position = np.array( - [a * np.cos(i * 2.0 * np.pi / (len(self.params['ids']) + 1)), - b * np.sin(i * 2.0 * np.pi / (len(self.params['ids']) + 1)), - i * (self.params['max_z'] - self.params['min_z']) / - (len(self.params['ids']) + 1) - self.params['center'][2]]) - - e_z = np.array([0, 0, 1]) - e_r = position - (position * e_z) * e_z - e_r /= np.linalg.norm(e_r) - e_phi = np.cross(e_z, e_r) - velocity = e_r * self.v_r + e_phi * self.v_phi + e_z * self.v_z - - velocity = self.swap_axis(velocity, self.params['axis']) - position = self.swap_axis(position, self.params['axis']) - position += np.array(self.params['center']) - self.system.part.add(id=i, pos=position, v=velocity) - - def calculate_numpy_histogram(self): - pol_positions = self.pol_coords() + [semi_x * np.cos(angle), + semi_y * np.sin(angle), + z]) + + e_r, e_phi, e_z = tests_common.get_cylindrical_basis_vectors( + position) + velocity = self.v_r * e_r + self.v_phi * e_phi + self.v_z * e_z + + positions.append(position) + velocities.append(velocity) + + return np.array(positions), np.array(velocities) + + def align_with_observable_frame(self, vec): + """ + Rotate vectors from the original box frame to the frame of the observables. + """ + + # align original z to observable z + vec = tests_common.rodrigues_rot(vec, [1, -1, 0], -np.pi / 2.) + # original x now points along [sqrt(3),-sqrt(3),-sqrt(3)] + + # align original x to observable orientation + vec = tests_common.rodrigues_rot(vec, [1, 1, 0], -3. / 4. * np.pi) + return vec + + def setup_system_get_np_hist(self): + """ + Pick positions and velocities in the original box frame + and calculate the np histogram. + Then rotate and move the positions and velocities + to the frame of the observables. + After calculating the core observables, the result should be + the same as the np histogram obtained from the original box frame. + """ + + positions, velocities = self.calc_ellipsis_pos_vel(100, 0.99 * + self.params['min_z'], 0.9 * + self.params['max_z'], semi_x=0.9 * + self.params['max_r'], semi_y=0.2 * + self.params['max_r']) + + # first, get the numpy histogram of the cylinder coordinates + pos_cyl = [] + for pos in positions: + pos_cyl.append( + tests_common.transform_pos_from_cartesian_to_polar_coordinates(pos)) np_hist, np_edges = tests_common.get_histogram( - pol_positions, self.params, 'cylindrical') - return np_hist, np_edges - - def normalize_with_bin_volume(self, histogram): - bin_volume = tests_common.get_cylindrical_bin_volume( - self.params['n_r_bins'], - self.params['n_phi_bins'], - self.params['n_z_bins'], - self.params['min_r'], - self.params['max_r'], - self.params['min_phi'], - self.params['max_phi'], - self.params['min_z'], - self.params['max_z']) - for i in range(self.params['n_r_bins']): - histogram[i, :, :] /= bin_volume[i] - return histogram - - def density_profile_test(self): - self.set_particles() - # Set up the Observable. - local_params = self.params.copy() - if self.params['axis'] == 'x': - local_params['axis'] = [1.0, 0.0, 0.0] - elif self.params['axis'] == 'y': - local_params['axis'] = [0.0, 1.0, 0.0] - else: - local_params['axis'] = [0.0, 0.0, 1.0] - obs = espressomd.observables.CylindricalDensityProfile(**local_params) - core_hist = obs.calculate() - core_edges = obs.call_method("edges") - np_hist, np_edges = self.calculate_numpy_histogram() - np_hist = self.normalize_with_bin_volume(np_hist) - np.testing.assert_array_almost_equal(np_hist, core_hist) - for i in range(3): - np.testing.assert_array_almost_equal(np_edges[i], core_edges[i]) - self.assertEqual(np.prod(obs.shape()), len(np_hist.flatten())) - - def velocity_profile_test(self): - self.set_particles() - # Set up the Observable. - local_params = self.params.copy() - if self.params['axis'] == 'x': - local_params['axis'] = [1.0, 0.0, 0.0] - elif self.params['axis'] == 'y': - local_params['axis'] = [0.0, 1.0, 0.0] - else: - local_params['axis'] = [0.0, 0.0, 1.0] - obs = espressomd.observables.CylindricalVelocityProfile(**local_params) - core_hist = obs.calculate() + np.array(pos_cyl), self.params, 'cylindrical') + np_dens = tests_common.normalize_cylindrical_hist( + np_hist.copy(), self.params) + + # now align the positions and velocities with the frame of reference + # used in the observables + pos_aligned = [] + vel_aligned = [] + for pos, vel in zip(positions, velocities): + pos_aligned.append( + self.align_with_observable_frame(pos) + + self.cyl_transform_params.center) + vel_aligned.append(self.align_with_observable_frame(vel)) + self.system.part.add(pos=pos_aligned, v=vel_aligned) + self.params['ids'] = self.system.part[:].id + + return np_dens, np_edges + + def check_edges(self, observable, np_edges): + core_edges = observable.call_method("edges") + for core_edge, np_edge in zip(core_edges, np_edges): + np.testing.assert_array_almost_equal(core_edge, np_edge) + + def test_density_profile(self): + """ + Check that the result from the observable (in its own frame) + matches the np result from the box frame + """ + np_dens, np_edges = self.setup_system_get_np_hist() + + cyl_dens_prof = espressomd.observables.CylindricalDensityProfile( + **self.params) + core_hist = cyl_dens_prof.calculate() + np.testing.assert_array_almost_equal(np_dens, core_hist) + self.check_edges(cyl_dens_prof, np_edges) + + def test_vel_profile(self): + """ + Check that the result from the observable (in its own frame) + matches the np result from the box frame + """ + np_dens, np_edges = self.setup_system_get_np_hist() + cyl_vel_prof = espressomd.observables.CylindricalVelocityProfile( + **self.params) + core_hist = cyl_vel_prof.calculate() core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] - np_hist, _ = self.calculate_numpy_histogram() - for x in np.nditer(np_hist, op_flags=['readwrite']): - if x[...] > 0.0: - x[...] /= x[...] - np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) + np_hist_binary = np_dens + np_hist_binary[np.nonzero(np_hist_binary)] = 1 + np.testing.assert_array_almost_equal( + np_hist_binary * self.v_r, core_hist_v_r) np.testing.assert_array_almost_equal( - np_hist * self.v_phi, core_hist_v_phi) - np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z) - self.assertEqual(np.prod(obs.shape()), len(np_hist.flatten()) * 3) - - def flux_density_profile_test(self): - self.set_particles() - # Set up the Observable. - local_params = self.params.copy() - if self.params['axis'] == 'x': - local_params['axis'] = [1.0, 0.0, 0.0] - elif self.params['axis'] == 'y': - local_params['axis'] = [0.0, 1.0, 0.0] - else: - local_params['axis'] = [0.0, 0.0, 1.0] - obs = espressomd.observables.CylindricalFluxDensityProfile( - **local_params) - core_hist = obs.calculate() + np_hist_binary * self.v_phi, core_hist_v_phi) + np.testing.assert_array_almost_equal( + np_hist_binary * self.v_z, core_hist_v_z) + self.check_edges(cyl_vel_prof, np_edges) + + def test_flux_density_profile(self): + """ + Check that the result from the observable (in its own frame) + matches the np result from the box frame + """ + np_dens, np_edges = self.setup_system_get_np_hist() + cyl_flux_dens = espressomd.observables.CylindricalFluxDensityProfile( + **self.params) + core_hist = cyl_flux_dens.calculate() core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] - np_hist, _ = self.calculate_numpy_histogram() - np_hist = self.normalize_with_bin_volume(np_hist) - np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) + np.testing.assert_array_almost_equal(np_dens * self.v_r, core_hist_v_r) np.testing.assert_array_almost_equal( - np_hist * self.v_phi, core_hist_v_phi) - np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z) - self.assertEqual(np.prod(obs.shape()), len(np_hist.flatten()) * 3) - - def test_hist_x(self): - self.params['axis'] = 'x' - self.velocity_profile_test() - self.flux_density_profile_test() - self.density_profile_test() - - def test_hist_y(self): - self.params['axis'] = 'y' - self.velocity_profile_test() - self.flux_density_profile_test() - self.density_profile_test() - - def test_hist_z(self): - self.params['axis'] = 'z' - self.velocity_profile_test() - self.flux_density_profile_test() - self.density_profile_test() + np_dens * self.v_phi, core_hist_v_phi) + np.testing.assert_array_almost_equal(np_dens * self.v_z, core_hist_v_z) + self.check_edges(cyl_flux_dens, np_edges) def test_cylindrical_pid_profile_interface(self): - # test setters and getters + """ + Test setters and getters of the script interface + """ params = self.params.copy() params['n_r_bins'] = 4 params['n_phi_bins'] = 6 params['n_z_bins'] = 8 params['ids'] = [0, 1] - params['axis'] = [0.0, 1.0, 0.0] self.system.part.add(id=0, pos=[0, 0, 0], type=0) self.system.part.add(id=1, pos=[0, 0, 0], type=1) observable = espressomd.observables.CylindricalDensityProfile(**params) @@ -266,15 +251,14 @@ def test_cylindrical_pid_profile_interface(self): self.assertEqual(observable.max_z, 9) obs_bin_edges = observable.bin_edges() np.testing.assert_array_equal(obs_bin_edges[-1, -1, -1], [7, 8, 9]) - # check center - np.testing.assert_array_equal( - np.copy(observable.center), params['center']) - observable.center = [3, 2, 1] - np.testing.assert_array_equal(np.copy(observable.center), [3, 2, 1]) - # check axis - np.testing.assert_array_equal(np.copy(observable.axis), params['axis']) - observable.axis = [6, 5, 4] - np.testing.assert_array_equal(np.copy(observable.axis), [6, 5, 4]) + # check center, axis, orientation + ctp = espressomd.math.CylindricalTransformationParameters( + center=[1, 2, 3], axis=[0, 1, 0], orientation=[0, 0, 1]) + observable.transform_params = ctp + + for attr_name in ['center', 'axis', 'orientation']: + np.testing.assert_array_almost_equal(np.copy(ctp.__getattr__(attr_name)), + np.copy(observable.transform_params.__getattr__(attr_name))) if __name__ == "__main__": diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index e2ba5a665c6..45001c928e6 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -18,210 +18,158 @@ import unittest as ut import unittest_decorators as utx import espressomd +import espressomd.math import espressomd.observables import espressomd.lb import tests_common -AGRID = 1.0 -VISC = 2.7 -DENS = 1.7 -TIME_STEP = 0.1 -LB_PARAMS = {'agrid': AGRID, - 'dens': DENS, - 'visc': VISC, - 'tau': TIME_STEP, - } - - class CylindricalLBObservableCommon: """ - Testcase for the CylindricalLBObservable. + Testcase for the CylindricalLBObservables. """ lbf = None - system = espressomd.System(box_l=(10, 10, 10)) + system = espressomd.System(box_l=3 * [15]) system.time_step = 0.01 system.cell_system.skin = 0.4 positions = [] + lb_params = {'agrid': 1., + 'dens': 1.2, + 'visc': 2.7, + 'tau': 0.1, + } + cyl_transform_params = espressomd.math.CylindricalTransformationParameters( + center=3 * [7], axis=[1, 0, 0], orientation=[0, 0, 1]) + params = { - 'ids': list(range(10)), - 'center': [5.0, 5.0, 5.0], # center of the histogram - 'axis': 'y', - 'n_r_bins': 10, # number of bins in r - 'n_phi_bins': 2, # -*- in phi - 'n_z_bins': 2, # -*- in z + 'ids': None, + 'transform_params': cyl_transform_params, + 'n_r_bins': 4, + 'n_phi_bins': 3, + 'n_z_bins': 5, 'min_r': 0.0, 'min_phi': -np.pi, - 'min_z': -5.0, - 'max_r': 5.0, + 'min_z': -6.0, + 'max_r': 6.0, 'max_phi': np.pi, - 'max_z': 5.0, + 'max_z': 6.0, } - def tearDown(self): - self.system.part.clear() - - def swap_axis(self, arr, axis): - if axis == 'x': - arr = np.dot(tests_common.rotation_matrix( - [0, 1, 0], np.pi / 2.0), arr) - elif axis == 'y': - arr = np.dot(tests_common.rotation_matrix( - [1, 0, 0], -np.pi / 2.0), arr) - return arr - - def swap_axis_inverse(self, arr, axis): - if axis == 'x': - arr = np.dot(tests_common.rotation_matrix( - [0, 1, 0], -np.pi / 2.0), arr) - elif axis == 'y': - arr = np.dot(tests_common.rotation_matrix( - [1, 0, 0], np.pi / 2.0), arr) - return arr - - def pol_coords(self): - positions = np.zeros((len(self.positions), 3)) - for i, p in enumerate(self.positions): - tmp = p - np.array(self.params['center']) - tmp = self.swap_axis_inverse(tmp, self.params['axis']) - positions[i, :] = tests_common.transform_pos_from_cartesian_to_polar_coordinates( - tmp) - return positions - - def set_particles(self): - self.system.part.clear() - self.system.part.add(pos=self.positions) - - def set_fluid_velocity(self): - del self.positions[:] - # Choose the cartesian velocities such that each particle gets the same - # v_r, v_phi and v_z, respectively. - self.v_r = .75 - self.v_phi = 2.5 - self.v_z = 1.5 - node_positions = np.arange(-4.5, 5.0, 1.0) - for i, _ in enumerate(node_positions): - position = np.array( - [node_positions[i], node_positions[i], node_positions[i]]) - - e_z = np.array([0, 0, 1]) - e_r = position - (position * e_z) * e_z - e_r /= np.linalg.norm(e_r) - e_phi = np.cross(e_z, e_r) - - velocity = e_r * self.v_r + e_phi * self.v_phi + e_z * self.v_z - - velocity = self.swap_axis(velocity, self.params['axis']) - position = self.swap_axis(position, self.params['axis']) - position += np.array(self.params['center']) - self.positions.append(position) - self.lbf[np.array(position, dtype=int)].velocity = velocity - - def normalize_with_bin_volume(self, histogram): - bin_volume = tests_common.get_cylindrical_bin_volume( - self.params['n_r_bins'], - self.params['n_phi_bins'], - self.params['n_z_bins'], - self.params['min_r'], - self.params['max_r'], - self.params['min_phi'], - self.params['max_phi'], - self.params['min_z'], - self.params['max_z']) - # Normalization - for i in range(self.params['n_r_bins']): - histogram[i, :, :] /= bin_volume[i] - return histogram - - def LB_fluxdensity_profile_test(self): - self.set_fluid_velocity() - self.set_particles() - # Set up the Observable. - local_params = self.params.copy() - if self.params['axis'] == 'x': - local_params['axis'] = [1.0, 0.0, 0.0] - elif self.params['axis'] == 'y': - local_params['axis'] = [0.0, 1.0, 0.0] - else: - local_params['axis'] = [0.0, 0.0, 1.0] - p = espressomd.observables.CylindricalLBFluxDensityProfileAtParticlePositions( - **local_params) - core_hist = p.calculate() - core_hist_v_r = core_hist[:, :, :, 0] - core_hist_v_phi = core_hist[:, :, :, 1] - core_hist_v_z = core_hist[:, :, :, 2] - core_edges = p.call_method("edges") - self.pol_positions = self.pol_coords() + v_r = 0.02 + v_phi = 0.04 + v_z = 0.03 + + def calc_vel_at_pos(self, positions): + """ + In cylindrical coordinates, all velocities are the same. + In cartesian they depend on the position. + The cartesian velocities are calculated here. + """ + + vels = [] + for pos in positions: + e_r, e_phi, e_z = tests_common.get_cylindrical_basis_vectors(pos) + velocity = self.v_r * e_r + self.v_phi * e_phi + self.v_z * e_z + vels.append(velocity) + return vels + + def align_with_observable_frame(self, vec): + """ + Rotate vectors from the original box frame to + the frame of the observables. + """ + + # align original z to observable z + vec = tests_common.rodrigues_rot(vec, [0, 1, 0], np.pi / 2.) + # original x now points along [0,0,-1] + + # align original x to observable orientation + vec = tests_common.rodrigues_rot(vec, [1, 0, 0], np.pi) + return vec + + def setup_system_get_np_hist(self): + """ + Pick positions and velocities in the original box frame and + calculate the np histogram. Then rotate and move the positions + and velocities to the frame of the observables. + After calculating the core observables, the result should be + the same as the np histogram obtained from the original box frame. + """ + + nodes = np.array(np.meshgrid([1, 2], [1, 2], [ + 1, 1, 1, 1, 2])).T.reshape(-1, 3) + positions = nodes + 3 * [0.5] + velocities = self.calc_vel_at_pos(positions) + + # get the histogram from numpy + pos_cyl = [] + for pos in positions: + pos_cyl.append( + tests_common.transform_pos_from_cartesian_to_polar_coordinates(pos)) np_hist, np_edges = tests_common.get_histogram( - self.pol_positions, self.params, 'cylindrical') - np_hist = self.normalize_with_bin_volume(np_hist) - np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) + np.array(pos_cyl), self.params, 'cylindrical') + + # the particles only determine the evaluation points, not the values of + # the observables + np_hist[np.nonzero(np_hist)] = 1 + + # now align the positions and velocities with the frame of reference + # used in the observables + pos_aligned = [] + vel_aligned = [] + for pos, vel in zip(positions, velocities): + pos_aligned.append( + self.align_with_observable_frame(pos) + + self.cyl_transform_params.center) + vel_aligned.append(self.align_with_observable_frame(vel)) + node_aligned = np.array( + np.rint( + np.array(pos_aligned) - + 3 * + [0.5]), + dtype=int) + self.system.part.add(pos=pos_aligned, v=vel_aligned) + self.params['ids'] = self.system.part[:].id + + for node, vel in zip(node_aligned, vel_aligned): + self.lbf[node].velocity = vel + + return np_hist, np_edges + + def check_edges(self, observable, np_edges): + core_edges = observable.call_method("edges") + for core_edge, np_edge in zip(core_edges, np_edges): + np.testing.assert_array_almost_equal(core_edge, np_edge) + + def test_cylindrical_lb_vel_profile_obs(self): + """ + Check that the result from the observable (in its own frame) + matches the np result from the box frame + """ + + np_hist_binary, np_edges = self.setup_system_get_np_hist() + vel_obs = espressomd.observables.CylindricalLBVelocityProfileAtParticlePositions( + **self.params) + core_hist_v = vel_obs.calculate() + core_hist_v_r = core_hist_v[:, :, :, 0] + core_hist_v_phi = core_hist_v[:, :, :, 1] + core_hist_v_z = core_hist_v[:, :, :, 2] np.testing.assert_array_almost_equal( - np_hist * self.v_phi, core_hist_v_phi) - np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z) - for i in range(3): - np.testing.assert_array_almost_equal(np_edges[i], core_edges[i]) - self.assertEqual(np.prod(p.shape()), len(np_hist.flatten()) * 3) - - def LB_velocity_profile_at_particle_positions_test(self): - self.set_fluid_velocity() - self.set_particles() - # Set up the Observable. - local_params = self.params.copy() - if self.params['axis'] == 'x': - local_params['axis'] = [1.0, 0.0, 0.0] - elif self.params['axis'] == 'y': - local_params['axis'] = [0.0, 1.0, 0.0] - else: - local_params['axis'] = [0.0, 0.0, 1.0] - p = espressomd.observables.CylindricalLBVelocityProfileAtParticlePositions( - **local_params) - core_hist = p.calculate() - core_hist_v_r = core_hist[:, :, :, 0] - core_hist_v_phi = core_hist[:, :, :, 1] - core_hist_v_z = core_hist[:, :, :, 2] - self.pol_positions = self.pol_coords() - np_hist, _ = np.histogramdd( - self.pol_positions, - bins=(self.params['n_r_bins'], - self.params['n_phi_bins'], - self.params['n_z_bins']), - range=[(self.params['min_r'], - self.params['max_r']), - (self.params['min_phi'], - self.params['max_phi']), - (self.params['min_z'], - self.params['max_z'])]) - for x in np.nditer(np_hist, op_flags=['readwrite']): - if x[...] > 0.0: - x[...] /= x[...] - np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) + np_hist_binary * self.v_r, core_hist_v_r) np.testing.assert_array_almost_equal( - np_hist * self.v_phi, core_hist_v_phi) - np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z) - self.assertEqual(np.prod(p.shape()), len(np_hist.flatten()) * 3) - - def perform_tests(self): - self.LB_fluxdensity_profile_test() - self.LB_velocity_profile_at_particle_positions_test() - - def test_x_axis(self): - self.params['axis'] = 'x' - self.perform_tests() - - def test_y_axis(self): - self.params['axis'] = 'y' - self.perform_tests() - - def test_z_axis(self): - self.params['axis'] = 'z' - self.perform_tests() + np_hist_binary * self.v_phi, core_hist_v_phi) + np.testing.assert_array_almost_equal( + np_hist_binary * self.v_z, core_hist_v_z) + self.check_edges(vel_obs, np_edges) def test_cylindrical_lb_profile_interface(self): - # test setters and getters + """ + Test setters and getters of the script interface + """ + params = self.params.copy() params['n_r_bins'] = 4 params['n_phi_bins'] = 6 @@ -269,25 +217,20 @@ def test_cylindrical_lb_profile_interface(self): self.assertEqual(observable.max_z, 9) obs_bin_edges = observable.bin_edges() np.testing.assert_array_equal(obs_bin_edges[-1, -1, -1], [7, 8, 9]) - # check center - np.testing.assert_array_equal( - np.copy(observable.center), params['center']) - observable.center = [3, 2, 1] - np.testing.assert_array_equal(np.copy(observable.center), [3, 2, 1]) - # check axis - np.testing.assert_array_equal(np.copy(observable.axis), params['axis']) - observable.axis = [6, 5, 4] - np.testing.assert_array_equal(np.copy(observable.axis), [6, 5, 4]) - # check sampling_density - self.assertEqual(observable.sampling_density, 2) - observable.sampling_density = 3 - self.assertEqual(observable.sampling_density, 3) + # check center, axis, orientation + ctp = espressomd.math.CylindricalTransformationParameters( + center=[1, 2, 3], axis=[0, 1, 0], orientation=[0, 0, 1]) + observable.transform_params = ctp + + for attr_name in ['center', 'axis', 'orientation']: + np.testing.assert_array_almost_equal(np.copy(ctp.__getattr__(attr_name)), + np.copy(observable.transform_params.__getattr__(attr_name))) class CylindricalLBObservableCPU(ut.TestCase, CylindricalLBObservableCommon): def setUp(self): - self.lbf = espressomd.lb.LBFluid(**LB_PARAMS) + self.lbf = espressomd.lb.LBFluid(**self.lb_params) self.system.actors.add(self.lbf) def tearDown(self): @@ -295,16 +238,47 @@ def tearDown(self): self.system.actors.remove(self.lbf) self.system.part.clear() + def test_cylindrical_lb_flux_density_obs(self): + """ + Check that the result from the observable (in its own frame) + matches the np result from the box frame. + Only for CPU because density interpolation is not implemented for GPU LB. + """ + np_hist_binary, np_edges = self.setup_system_get_np_hist() + + flux_obs = espressomd.observables.CylindricalLBFluxDensityProfileAtParticlePositions( + **self.params) + core_hist_fl = flux_obs.calculate() + core_hist_fl_r = core_hist_fl[:, :, :, 0] + core_hist_fl_phi = core_hist_fl[:, :, :, 1] + core_hist_fl_z = core_hist_fl[:, :, :, 2] + + np.testing.assert_array_almost_equal( + np_hist_binary * + self.lb_params['dens'] * + self.v_r, + core_hist_fl_r) + np.testing.assert_array_almost_equal( + np_hist_binary * + self.lb_params['dens'] * + self.v_phi, + core_hist_fl_phi) + np.testing.assert_array_almost_equal( + np_hist_binary * + self.lb_params['dens'] * + self.v_z, + core_hist_fl_z) + self.check_edges(flux_obs, np_edges) + @utx.skipIfMissingGPU() class CylindricalLBObservableGPU(ut.TestCase, CylindricalLBObservableCommon): def setUp(self): - self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) + self.lbf = espressomd.lb.LBFluidGPU(**self.lb_params) self.system.actors.add(self.lbf) def tearDown(self): - del self.positions[:] self.system.actors.remove(self.lbf) self.system.part.clear() diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index 92d169a1620..5c259cd49b2 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -134,7 +134,7 @@ def abspath(path): def transform_pos_from_cartesian_to_polar_coordinates(pos): - """Transform the given cartesian coordinates to polar coordinates. + """Transform the given cartesian coordinates to cylindrical coordinates. Parameters ---------- @@ -167,11 +167,28 @@ def transform_vel_from_cartesian_to_polar_coordinates(pos, vel): (pos[0] * vel[1] - pos[1] * vel[0]) / np.sqrt(pos[0]**2 + pos[1]**2), vel[2]]) +def get_cylindrical_basis_vectors(pos): + phi = transform_pos_from_cartesian_to_polar_coordinates(pos)[1] + e_r = np.array([np.cos(phi), np.sin(phi), 0.]) + e_phi = np.array([-np.sin(phi), np.cos(phi), 0.]) + e_z = np.array([0., 0., 1.]) + return e_r, e_phi, e_z + + def convert_vec_body_to_space(system, part, vec): A = rotation_matrix_quat(system, part) return np.dot(A.transpose(), vec) +def rodrigues_rot(vec, axis, angle): + """ + https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Statement + """ + axis /= np.linalg.norm(axis) + return np.cos(angle) * vec + np.sin(angle) * np.cross(axis, vec) + \ + (1 - np.cos(angle)) * np.dot(axis, vec) * axis + + def rotation_matrix(axis, theta): """ Return the rotation matrix associated with counterclockwise rotation about @@ -225,55 +242,41 @@ def rotation_matrix_quat(system, part): return A -def get_cylindrical_bin_volume( - n_r_bins, - n_phi_bins, - n_z_bins, - min_r, - max_r, - min_phi, - max_phi, - min_z, - max_z): +def normalize_cylindrical_hist(histogram, cyl_obs_params): """ - Return the bin volumes for a cylindrical histogram. + normalize a histogram in cylindrical coordinates. Helper to test the output + of cylindrical histogram observables Parameters ---------- - n_r_bins : :obj:`float` - Number of bins in ``r`` direction. - n_phi_bins : :obj:`float` - Number of bins in ``phi`` direction. - n_z_bins : :obj:`float` - Number of bins in ``z`` direction. - min_r : :obj:`float` - Minimum considered value in ``r`` direction. - max_r : :obj:`float` - Maximum considered value in ``r`` direction. - min_phi : :obj:`float` - Minimum considered value in ``phi`` direction. - max_phi : :obj:`float` - Maximum considered value in ``phi`` direction. - min_z : :obj:`float` - Minimum considered value in ``z`` direction. - max_z : :obj:`float` - Maximum considered value in ``z`` direction. + histogram : (N,3) array_like of :obj:`float` + The histogram that needs to be normalized + cyl_obs_params : :obj:`dict` + A dictionary containing the common parameters of the cylindrical histogram observables. + Needs to contain the information about number and range of bins. + """ - Returns - ------- - array_like - Bin volumes. + n_r_bins = cyl_obs_params['n_r_bins'] + n_phi_bins = cyl_obs_params['n_phi_bins'] + n_z_bins = cyl_obs_params['n_z_bins'] + min_r = cyl_obs_params['min_r'] + max_r = cyl_obs_params['max_r'] + min_phi = cyl_obs_params['min_phi'] + max_phi = cyl_obs_params['max_phi'] + min_z = cyl_obs_params['min_z'] + max_z = cyl_obs_params['max_z'] - """ bin_volume = np.zeros(n_r_bins) r_bin_size = (max_r - min_r) / n_r_bins phi_bin_size = (max_phi - min_phi) / n_phi_bins z_bin_size = (max_z - min_z) / n_z_bins for i in range(n_r_bins): - bin_volume[i] = np.pi * ((min_r + r_bin_size * (i + 1))**2.0 - - (min_r + r_bin_size * i)**2.0) * \ + bin_volume = np.pi * ((min_r + r_bin_size * (i + 1))**2.0 - + (min_r + r_bin_size * i)**2.0) * \ phi_bin_size / (2.0 * np.pi) * z_bin_size - return bin_volume + histogram[i, :, :] /= bin_volume + + return histogram def get_histogram(pos, obs_params, coord_system, **kwargs):