nucleationParameters.h

#ifndef INCLUDE_NUCLEATIONPARAMETERS_H_
#define INCLUDE_NUCLEATIONPARAMETERS_H_
 
#include <deal.II/base/tensor.h>
 
#include <utility>
#include <vector>
 
template <int dim>
class nucleationParameters
{
public:
  unsigned int                   var_index;
  std::vector<double>            semiaxes;
  std::vector<double>            ellipsoid_rotation;
  std::vector<double>            freeze_semiaxes;
  double                         no_nucleation_border_thickness;
  double                         hold_time;
  dealii::Tensor<2, dim, double> rotation_matrix;
 
  nucleationParameters(unsigned int        _var_index,
                       std::vector<double> _semiaxes,
                       std::vector<double> _freeze_semiaxes,
                       std::vector<double> _ellipsoid_rotation,
                       double              _hold_time,
                       double              _no_nucleation_border_thickness)
    : var_index(_var_index)
    , semiaxes(std::move(_semiaxes))
    , ellipsoid_rotation(std::move(_ellipsoid_rotation))
    , freeze_semiaxes(std::move(_freeze_semiaxes))
    , no_nucleation_border_thickness(_no_nucleation_border_thickness)
    , hold_time(_hold_time)
  {
    set_rotation_matrix();
  };
 
  void
  set_rotation_matrix()
  {
    // Rotation conventions:
    // Rx refers to rotations about the x axis, Ry refers to rotations about the
    // y axis, and Rz refers to rotations about the z axis. A positive rotation
    // angle about the z axis corresponds to a clockwise rotation of the
    // particle in a 2D calculation.
 
    double degrees_to_rad = std::acos(0.0) / 90.0;
 
    dealii::Tensor<2, dim, double> Rx;
    dealii::Tensor<2, dim, double> Ry;
    dealii::Tensor<2, dim, double> Rz;
 
    Rx[0][0] = 1.0;
    Rx[1][1] = std::cos(ellipsoid_rotation.at(0) * degrees_to_rad);
 
    Ry[0][0] = std::cos(ellipsoid_rotation.at(1) * degrees_to_rad);
    Ry[1][1] = 1.0;
 
    Rz[0][0] = std::cos(ellipsoid_rotation.at(2) * degrees_to_rad);
    Rz[1][0] = std::sin(ellipsoid_rotation.at(2) * degrees_to_rad);
    Rz[0][1] = -std::sin(ellipsoid_rotation.at(2) * degrees_to_rad);
    Rz[1][1] = std::cos(ellipsoid_rotation.at(2) * degrees_to_rad);
 
    if (dim == 3)
      {
        Rx[1][2] = -std::sin(ellipsoid_rotation.at(0) * degrees_to_rad);
        Rx[2][1] = std::sin(ellipsoid_rotation.at(0) * degrees_to_rad);
        Rx[2][2] = std::cos(ellipsoid_rotation.at(0) * degrees_to_rad);
 
        Ry[0][2] = std::sin(ellipsoid_rotation.at(1) * degrees_to_rad);
        Ry[2][0] = -std::sin(ellipsoid_rotation.at(1) * degrees_to_rad);
        Ry[2][2] = std::cos(ellipsoid_rotation.at(1) * degrees_to_rad);
 
        Rz[2][2] = 1.0;
      }
 
    rotation_matrix = Rx * Ry * Rz; // Note: these are tensor multiplications
  };
};
 
#endif