cpu_acc_intersections.cpp

Go to the documentation of this file.

/*
 * This file is part of Vlasiator.
 * Copyright 2010-2016 Finnish Meteorological Institute
 *
 * For details of usage, see the COPYING file and read the "Rules of the Road"
 * at http://www.physics.helsinki.fi/vlasiator/
 *
 * This program 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 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */
 
#include <algorithm>
#include <cmath>
#include <utility>
 
#include <Eigen/Geometry>
#include "../common.h"
#include "cpu_acc_intersections.hpp"
#include "cpu_acc_transform.hpp"
 
using namespace std;
using namespace Eigen;

void compute_cell_intersections(
        spatial_cell::SpatialCell* spatial_cell,
        const uint popID,
        const uint map_order,
        const Real& dt,
        int intersections_id
   ) {
   vmesh::VelocityMesh* vmesh    = spatial_cell->get_velocity_mesh(popID);
   spatial_cell::Population& pop = spatial_cell->get_population(popID);
 
   // compute transform, forward in time and backward in time, performed in this acceleration
   Transform<Real,3,Affine> fwd_transform= compute_acceleration_transformation(spatial_cell,popID,dt);
   Transform<Real,3,Affine> bwd_transform= fwd_transform.inverse();
 
   phiprof::Timer intersectionsTimer {intersections_id};
   switch(map_order){
      case 0: {
         //Map order XYZ
         compute_intersections_1st(vmesh,bwd_transform, fwd_transform, 0,
                                   pop.intersection_x,pop.intersection_x_di,
                                   pop.intersection_x_dj,pop.intersection_x_dk);
         compute_intersections_2nd(vmesh,bwd_transform, fwd_transform, 1,
                                   pop.intersection_y,pop.intersection_y_di,
                                   pop.intersection_y_dj,pop.intersection_y_dk);
         compute_intersections_3rd(vmesh,bwd_transform, fwd_transform, 2,
                                   pop.intersection_z,pop.intersection_z_di,
                                   pop.intersection_z_dj,pop.intersection_z_dk);
         break;
      }
      case 1: {
         //Map order YZX
         compute_intersections_1st(vmesh, bwd_transform, fwd_transform, 1,
                                   pop.intersection_y,pop.intersection_y_di,
                                   pop.intersection_y_dj,pop.intersection_y_dk);
         compute_intersections_2nd(vmesh, bwd_transform, fwd_transform, 2,
                                   pop.intersection_z,pop.intersection_z_di,
                                   pop.intersection_z_dj,pop.intersection_z_dk);
         compute_intersections_3rd(vmesh, bwd_transform, fwd_transform, 0,
                                   pop.intersection_x,pop.intersection_x_di,
                                   pop.intersection_x_dj,pop.intersection_x_dk);
         break;
      }
      case 2: {
         //Map order Z X Y
         compute_intersections_1st(vmesh, bwd_transform, fwd_transform, 2,
                                   pop.intersection_z,pop.intersection_z_di,
                                   pop.intersection_z_dj,pop.intersection_z_dk);
         compute_intersections_2nd(vmesh, bwd_transform, fwd_transform, 0,
                                   pop.intersection_x,pop.intersection_x_di,
                                   pop.intersection_x_dj,pop.intersection_x_dk);
         compute_intersections_3rd(vmesh, bwd_transform, fwd_transform, 1,
                                   pop.intersection_y,pop.intersection_y_di,
                                   pop.intersection_y_dj,pop.intersection_y_dk);
         break;
      }
   }
   intersectionsTimer.stop();
}

Eigen::Matrix<Real,3,1> line_plane_intersection(const Eigen::Matrix<Real,3,1>& l_point,const Eigen::Matrix<Real,3,1>& l_direction,
                                                const Eigen::Matrix<Real,3,1>& p_point,const Eigen::Matrix<Real,3,1>& p_normal){
  const Real nom=p_normal.dot(p_point-l_point);
  const Real dem=p_normal.dot(l_direction);
  return l_point+(nom/dem)*l_direction;
}

void compute_intersections_1st(
        const vmesh::VelocityMesh* vmesh,
        const Transform<Real,3,Affine>& bwd_transform,const Transform<Real,3,Affine>& fwd_transform,
        const uint dimension,
        Real& intersection,Real& intersection_di,Real& intersection_dj,Real& intersection_dk) {
 
    if (dimension == 0) { //Prepare intersections for mapping along X first (mapping order X-Y-Z)
        // Normal of Lagrangian planes
        const Eigen::Matrix<Real,3,1> plane_normal = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(1.0, 0.0, 0.0);
        // Point on lowest potential Lagrangian plane
        const Eigen::Matrix<Real,3,1> plane_point
            = bwd_transform*Eigen::Matrix<Real,3,1>(vmesh->getMeshMinLimits()[0], 0.0, 0.0);
        // line along Euclidian x direction, unit vector
        const Eigen::Matrix<Real,3,1> line_direction = Eigen::Matrix<Real,3,1>(1.0, 0.0, 0.0);
        const Eigen::Matrix<Real,3,1> line_point(
            0.0,
            0.5*vmesh->getCellSize()[1]+vmesh->getMeshMinLimits()[1],
            0.5*vmesh->getCellSize()[2]+vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> lagrangian_di = bwd_transform.linear()
            * Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0],0,0.0);
        const Eigen::Matrix<Real,3,1> euclidian_dj = Eigen::Matrix<Real,3,1>(0,vmesh->getCellSize()[1],0.0);
        const Eigen::Matrix<Real,3,1> euclidian_dk = Eigen::Matrix<Real,3,1>(0.0,0.0,vmesh->getCellSize()[2]);
 
        // compute intersections, varying lines and plane in i,j,k
        const Eigen::Matrix<Real,3,1> intersection_0_0_0 = line_plane_intersection(line_point, line_direction, plane_point, plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_1_0_0 = line_plane_intersection(line_point, line_direction, plane_point + lagrangian_di, plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_0_1_0 = line_plane_intersection(line_point + euclidian_dj, line_direction, plane_point, plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_0_0_1 = line_plane_intersection(line_point + euclidian_dk, line_direction, plane_point, plane_normal);
        intersection=intersection_0_0_0[dimension];
        intersection_di=intersection_1_0_0[dimension]-intersection_0_0_0[dimension];
        intersection_dj=intersection_0_1_0[dimension]-intersection_0_0_0[dimension];
        intersection_dk=intersection_0_0_1[dimension]-intersection_0_0_0[dimension];
    }
    if (dimension == 1) { //Prepare intersections for mapping along Y first (mapping order Y-Z-X)
        // Normal of Lagrangian planes
        const Eigen::Matrix<Real,3,1> plane_normal = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(0.0, 1.0, 0.0);
        // Point on lowest possible Lagrangian plane
        const Eigen::Matrix<Real,3,1> plane_point
            = bwd_transform*Eigen::Matrix<Real,3,1>(0.0, vmesh->getMeshMinLimits()[1], 0.0);
        // line along Euclidian y direction, unit vector
        const Eigen::Matrix<Real,3,1> line_direction = Eigen::Matrix<Real,3,1>(0.0, 1.0, 0.0);
        const Eigen::Matrix<Real,3,1> line_point(
            0.5*vmesh->getCellSize()[0]+vmesh->getMeshMinLimits()[0],
            0.0,
            0.5*vmesh->getCellSize()[2]+vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> euclidian_di
            = Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0], 0.0, 0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dj
            = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(0.0 ,vmesh->getCellSize()[1], 0.0);
        const Eigen::Matrix<Real,3,1> euclidian_dk
            = Eigen::Matrix<Real,3,1>(0.0 , 0.0 ,vmesh->getCellSize()[2]);
 
        // compute intersections, varying lines and plane in i,j,k
        const Eigen::Matrix<Real,3,1> intersection_0_0_0 = line_plane_intersection(line_point,line_direction,plane_point,plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_1_0_0 = line_plane_intersection(line_point + euclidian_di, line_direction, plane_point, plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_0_1_0 = line_plane_intersection(line_point, line_direction, plane_point + lagrangian_dj, plane_normal);
        const Eigen::Matrix<Real,3,1> intersection_0_0_1 = line_plane_intersection(line_point + euclidian_dk, line_direction, plane_point, plane_normal);
 
        intersection=intersection_0_0_0[dimension];
        intersection_di=intersection_1_0_0[dimension]-intersection_0_0_0[dimension];
        intersection_dj=intersection_0_1_0[dimension]-intersection_0_0_0[dimension];
        intersection_dk=intersection_0_0_1[dimension]-intersection_0_0_0[dimension];
    }
 
    if (dimension == 2) {
        // This is the  case presented in the Slice 3D article
        // Prepare intersections for mapping along Z first (mapping order Z-X-Y)
 
        //Normal of Lagrangian planes
        const Eigen::Matrix<Real,3,1> plane_normal
            = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(0,0,1.0);
 
        // Point on lowest possible Lagrangian plane
        const Eigen::Matrix<Real,3,1> plane_point
            = bwd_transform*Eigen::Matrix<Real,3,1>(0.0,0.0,vmesh->getMeshMinLimits()[2]);
 
        // line along Euclidian z direction, unit vector
        const Eigen::Matrix<Real,3,1> line_direction = Eigen::Matrix<Real,3,1>(0,0,1.0);
        const Eigen::Matrix<Real,3,1> line_point(
            0.5*vmesh->getCellSize()[0]+vmesh->getMeshMinLimits()[0],
            0.5*vmesh->getCellSize()[1]+vmesh->getMeshMinLimits()[1],
            0.0);
        const Eigen::Matrix<Real,3,1> euclidian_di
            = Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0],0,0.0);
        const Eigen::Matrix<Real,3,1> euclidian_dj
            = Eigen::Matrix<Real,3,1>(0,vmesh->getCellSize()[1],0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dk
            = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(0.0,0.0,vmesh->getCellSize()[2]);
 
      // compute intersections, varying lines and plane in i,j,k
      const Eigen::Matrix<Real,3,1> intersection_0_0_0 = line_plane_intersection(line_point,line_direction,plane_point,plane_normal);
      const Eigen::Matrix<Real,3,1> intersection_1_0_0 = line_plane_intersection(line_point + euclidian_di, line_direction, plane_point, plane_normal);
      const Eigen::Matrix<Real,3,1> intersection_0_1_0 = line_plane_intersection(line_point + euclidian_dj, line_direction, plane_point, plane_normal);
      const Eigen::Matrix<Real,3,1> intersection_0_0_1 = line_plane_intersection(line_point, line_direction, plane_point + lagrangian_dk, plane_normal);
      intersection=intersection_0_0_0[dimension];
      intersection_di=intersection_1_0_0[dimension]-intersection_0_0_0[dimension];
      intersection_dj=intersection_0_1_0[dimension]-intersection_0_0_0[dimension];
      intersection_dk=intersection_0_0_1[dimension]-intersection_0_0_0[dimension];
   }
}

void compute_intersections_2nd(
        const vmesh::VelocityMesh* vmesh,
        const Transform<Real,3,Affine>& bwd_transform,const Transform<Real,3,Affine>& fwd_transform,
        const uint dimension,
        Real& intersection,Real& intersection_di,Real& intersection_dj,Real& intersection_dk){
 
    if (dimension == 0) { // Prepare intersections for mapping along X second (mapping order Z-X-Y)
        // This is the case presented in the Slice 3D article,
        // data along z has been moved to Lagrangian coordinates.
 
        // Normal of Euclidian y-plane
        const Eigen::Matrix<Real,3,1> plane_normal = Eigen::Matrix<Real,3,1>(0.0, 1.0, 0.0);
 
        //Point on lowest Euclidian y-plane through middle of cells
        Eigen::Matrix<Real,3,1> plane_point
            = Eigen::Matrix<Real,3,1>(0,vmesh->getMeshMinLimits()[1]+vmesh->getCellSize()[1]*0.5,0);
        const Eigen::Matrix<Real,3,1> lagrangian_di
            = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0],0,0.0);
 
        // Distance between Euclidian planes
        const Eigen::Matrix<Real,3,1> euclidian_dj
            = Eigen::Matrix<Real,3,1>(0,vmesh->getCellSize()[1],0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dk
            = bwd_transform.linear()*Eigen::Matrix<Real,3,1>(0.0,0.0,vmesh->getCellSize()[2]);
 
        // line along Lagrangian y line, unit vector. Only rotation here, not translation
        const Eigen::Matrix<Real,3,1> line_direction = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(0,1.0,0.0);
        const Eigen::Matrix<Real,3,1> line_point = bwd_transform * Eigen::Matrix<Real,3,1>(
            vmesh->getMeshMinLimits()[0],
            0.5*vmesh->getCellSize()[1]+vmesh->getMeshMinLimits()[1],
            0.5*vmesh->getCellSize()[2]+vmesh->getMeshMinLimits()[2]);
 
        // Compute two intersections between Lagrangian line (absolute position
        // does not matter so set to 0,0,0, and two Euclidian planes.
        Eigen::Matrix<Real,3,1> intersect_0_0_0 = line_plane_intersection(line_point,line_direction,plane_point,plane_normal);
        Eigen::Matrix<Real,3,1> intersect_1_0_0 = line_plane_intersection(line_point + lagrangian_di, line_direction, plane_point, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_1_0 = line_plane_intersection(line_point, line_direction, plane_point + euclidian_dj, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_0_1 = line_plane_intersection(line_point + lagrangian_dk, line_direction, plane_point, plane_normal);
 
        intersection=intersect_0_0_0[dimension];
        intersection_di = intersect_1_0_0[dimension] - intersect_0_0_0[dimension];
        intersection_dj = intersect_0_1_0[dimension] - intersect_0_0_0[dimension];
        intersection_dk = intersect_0_0_1[dimension] - intersect_0_0_0[dimension];
    }
    if (dimension == 1) { //Prepare intersections for mapping along Y second (mapping order X-Y-Z)
        // Normal of Euclidian z-plane
        const Eigen::Matrix<Real,3,1> plane_normal = Eigen::Matrix<Real,3,1>(0.0, 0.0, 1.0);
 
        // Point on lowest Euclidian z-plane through middle of cells
        Eigen::Matrix<Real,3,1> plane_point
            = Eigen::Matrix<Real,3,1>(0.0, 0.0,vmesh->getMeshMinLimits()[2]+vmesh->getCellSize()[2] * 0.5);
 
        const Eigen::Matrix<Real,3,1> lagrangian_di
            = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0], 0.0,  0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dj
            = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(0.0, vmesh->getCellSize()[1], 0.0);
        // Distance between Euclidian planes
        const Eigen::Matrix<Real,3,1> euclidian_dk
            = Eigen::Matrix<Real,3,1>(0.0, 0.0, vmesh->getCellSize()[2]);
 
        // line along Lagrangian z line, unit vector. Only rotation here, not translation
        const Eigen::Matrix<Real,3,1> line_direction
            = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(0.0, 0.0, 1.0);
        const Eigen::Matrix<Real,3,1> line_point = bwd_transform * Eigen::Matrix<Real,3,1>(
            0.5*vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
            vmesh->getMeshMinLimits()[1],
            0.5*vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
 
        // Compute two intersections between Lagrangian line (absolute position
        // does not matter so set to 0,0,0, and two Euclidian planes.
        Eigen::Matrix<Real,3,1> intersect_0_0_0 = line_plane_intersection(line_point,line_direction,plane_point,plane_normal);
        Eigen::Matrix<Real,3,1> intersect_1_0_0 = line_plane_intersection(line_point + lagrangian_di, line_direction, plane_point, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_1_0 = line_plane_intersection(line_point + lagrangian_dj, line_direction, plane_point, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_0_1 = line_plane_intersection(line_point, line_direction, plane_point + euclidian_dk, plane_normal);
 
        intersection=intersect_0_0_0[dimension];
        intersection_di = intersect_1_0_0[dimension] - intersect_0_0_0[dimension];
        intersection_dj = intersect_0_1_0[dimension] - intersect_0_0_0[dimension];
        intersection_dk = intersect_0_0_1[dimension] - intersect_0_0_0[dimension];
 
    }
    if (dimension == 2) { //Prepare intersections for mapping along Z second (mapping order Y-Z-X)
        // Normal of Euclidian x-plane
        const Eigen::Matrix<Real,3,1> plane_normal = Eigen::Matrix<Real,3,1>(1.0, 0.0, 0.0);
        //Point on lowest Euclidian x-plane through middle of cells
        Eigen::Matrix<Real,3,1> plane_point
            = Eigen::Matrix<Real,3,1>(vmesh->getMeshMinLimits()[0]+vmesh->getCellSize()[0]*0.5, 0.0, 0.0);
        // Distance between Euclidian planes
        const Eigen::Matrix<Real,3,1> euclidian_di = Eigen::Matrix<Real,3,1>(vmesh->getCellSize()[0], 0.0, 0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dj
            = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(0.0, vmesh->getCellSize()[1], 0.0);
        const Eigen::Matrix<Real,3,1> lagrangian_dk = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(0.0, 0.0, vmesh->getCellSize()[2]);
 
        // line along Lagrangian x line, unit vector. Only rotation here, not translation
        const Eigen::Matrix<Real,3,1> line_direction = bwd_transform.linear() * Eigen::Matrix<Real,3,1>(1.0, 0.0, 0.0);
        const Eigen::Matrix<Real,3,1> line_point = bwd_transform * Eigen::Matrix<Real,3,1>(
            0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
            0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
            vmesh->getMeshMinLimits()[2]);
 
        // Compute two intersections between Lagrangian line (absolute position
        // does not matter so set to 0,0,0, and two Euclidian planes.
        Eigen::Matrix<Real,3,1> intersect_0_0_0 = line_plane_intersection(line_point,line_direction,plane_point,plane_normal);
        Eigen::Matrix<Real,3,1> intersect_1_0_0 = line_plane_intersection(line_point, line_direction, plane_point + euclidian_di, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_1_0 = line_plane_intersection(line_point + lagrangian_dj, line_direction, plane_point, plane_normal);
        Eigen::Matrix<Real,3,1> intersect_0_0_1 = line_plane_intersection(line_point + lagrangian_dk, line_direction, plane_point, plane_normal);
 
        intersection=intersect_0_0_0[dimension];
        intersection_di = intersect_1_0_0[dimension] - intersect_0_0_0[dimension];
        intersection_dj = intersect_0_1_0[dimension] - intersect_0_0_0[dimension];
        intersection_dk = intersect_0_0_1[dimension] - intersect_0_0_0[dimension];
    }
}

void compute_intersections_3rd(
    const vmesh::VelocityMesh* vmesh,
    const Transform<Real,3,Affine>& bwd_transform,const Transform<Real,3,Affine>& fwd_transform,
    const uint dimension,
    Real& intersection,Real& intersection_di,Real& intersection_dj,Real& intersection_dk) {
 
    if (dimension == 0) { //Prepare intersections for mapping along X third (mapping order Y-Z-X)
        const Eigen::Matrix<Real,3,1> point_0_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.0 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_1_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(1.0 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_1_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.0 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      1.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_0_1 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.0 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      1.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        intersection = point_0_0_0[dimension];
        intersection_di = point_1_0_0[dimension]-point_0_0_0[dimension];
        intersection_dj = point_0_1_0[dimension]-point_0_0_0[dimension];
        intersection_dk = point_0_0_1[dimension]-point_0_0_0[dimension];
   }
   if (dimension == 1) { //Prepare intersections for mapping along Y third (mapping order Z-X-Y)
        // This is the case presented in the Slice 3D article,
        // data along z has been moved to Lagrangian coordinates
        const Eigen::Matrix<Real,3,1> point_0_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_1_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(1.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_1_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      1.0 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_0_1 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      vmesh->getMeshMinLimits()[1],
                                      1.5 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        intersection = point_0_0_0[dimension];
        intersection_di = point_1_0_0[dimension]-point_0_0_0[dimension];
        intersection_dj = point_0_1_0[dimension]-point_0_0_0[dimension];
        intersection_dk = point_0_0_1[dimension]-point_0_0_0[dimension];
    }
    if (dimension == 2) { //Prepare intersections for mapping along Z third (mapping order X-Y-Z)
        const Eigen::Matrix<Real,3,1> point_0_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.0 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_1_0_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(1.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.0 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_1_0 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      1.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      0.0 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        const Eigen::Matrix<Real,3,1> point_0_0_1 = bwd_transform
            * Eigen::Matrix<Real,3,1>(0.5 * vmesh->getCellSize()[0] + vmesh->getMeshMinLimits()[0],
                                      0.5 * vmesh->getCellSize()[1] + vmesh->getMeshMinLimits()[1],
                                      1.0 * vmesh->getCellSize()[2] + vmesh->getMeshMinLimits()[2]);
        intersection = point_0_0_0[dimension];
        intersection_di = point_1_0_0[dimension]-point_0_0_0[dimension];
        intersection_dj = point_0_1_0[dimension]-point_0_0_0[dimension];
        intersection_dk = point_0_0_1[dimension]-point_0_0_0[dimension];
    }
}