[cig-commits] r1285 - in trunk/aspect: include/aspect/velocity_boundary_conditions source/velocity_boundary_conditions
gassmoeller at dealii.org
gassmoeller at dealii.org
Thu Oct 18 12:58:14 PDT 2012
Author: gassmoeller
Date: 2012-10-18 13:58:14 -0600 (Thu, 18 Oct 2012)
New Revision: 1285
Modified:
trunk/aspect/include/aspect/velocity_boundary_conditions/gplates.h
trunk/aspect/source/velocity_boundary_conditions/gplates.cc
Log:
Indent
Modified: trunk/aspect/include/aspect/velocity_boundary_conditions/gplates.h
===================================================================
--- trunk/aspect/include/aspect/velocity_boundary_conditions/gplates.h 2012-10-18 19:13:21 UTC (rev 1284)
+++ trunk/aspect/include/aspect/velocity_boundary_conditions/gplates.h 2012-10-18 19:58:14 UTC (rev 1285)
@@ -35,17 +35,17 @@
namespace internal
{
- /**
- * GPlatesLookup handles all kinds of tasks around looking up a certain velocity boundary
- * condition from a gplates .gpml file. This class keeps around the contents of two sets
- * of files, corresponding to two instances in time where GPlates provides us with data;
- * the boundary values at one particular time are interpolated between the two currently
- * loaded data sets.
- */
+ /**
+ * GPlatesLookup handles all kinds of tasks around looking up a certain velocity boundary
+ * condition from a gplates .gpml file. This class keeps around the contents of two sets
+ * of files, corresponding to two instances in time where GPlates provides us with data;
+ * the boundary values at one particular time are interpolated between the two currently
+ * loaded data sets.
+ */
class GPlatesLookup
{
- public:
+ public:
/**
* Initialize all members and the two pointers referring to the actual velocities.
@@ -74,80 +74,80 @@
template <int dim>
Tensor<1,dim> surface_velocity(const Point<dim> &position_, const double time_weight) const;
- private:
+ private:
- /**
- * Tables which contain the velocities
- */
- dealii::Table<2,Tensor<1,2> > velocity_vals;
- dealii::Table<2,Tensor<1,2> > old_velocity_vals;
+ /**
+ * Tables which contain the velocities
+ */
+ dealii::Table<2,Tensor<1,2> > velocity_vals;
+ dealii::Table<2,Tensor<1,2> > old_velocity_vals;
- /**
- * Pointers to the actual tables.
- * Used to avoid unnecessary copying
- * of values.
- */
- dealii::Table<2,Tensor<1,2> > * velocity_values;
- dealii::Table<2,Tensor<1,2> > * old_velocity_values;
+ /**
+ * Pointers to the actual tables.
+ * Used to avoid unnecessary copying
+ * of values.
+ */
+ dealii::Table<2,Tensor<1,2> > * velocity_values;
+ dealii::Table<2,Tensor<1,2> > * old_velocity_values;
- /**
- * Distances between adjacent point in the Lat/Lon grid
- */
- double delta_phi,delta_theta;
+ /**
+ * Distances between adjacent point in the Lat/Lon grid
+ */
+ double delta_phi,delta_theta;
- /**
- * The rotation axis and angle around which a 2D model needs to be rotated to
- * be transformed to a plane that contains the origin and the two
- * user prescribed points. Is not used for 3D.
- */
- Tensor<1,3> rotation_axis;
- double rotation_angle;
+ /**
+ * The rotation axis and angle around which a 2D model needs to be rotated to
+ * be transformed to a plane that contains the origin and the two
+ * user prescribed points. Is not used for 3D.
+ */
+ Tensor<1,3> rotation_axis;
+ double rotation_angle;
- /**
- * A function that returns the rotated vector r' that results out of a
- * rotation from vector r around a specified rotation_axis by an defined angle
- */
- Tensor<1,3>
- rotate (const Tensor<1,3> &position,const Tensor<1,3> &rotation_axis, const double angle) const;
+ /**
+ * A function that returns the rotated vector r' that results out of a
+ * rotation from vector r around a specified rotation_axis by an defined angle
+ */
+ Tensor<1,3>
+ rotate (const Tensor<1,3> &position,const Tensor<1,3> &rotation_axis, const double angle) const;
- /**
- * Convert a tensor of rank 1 and dimension in to rank 1 and dimension out.
- * If $out < in$ the last elements will be discarded, if $out > in$ zeroes will
- * be appended to fill the tensor.
- */
- template <int in, int out>
- Tensor<1,out> convert_tensor (Tensor<1,in> old_tensor) const;
+ /**
+ * Convert a tensor of rank 1 and dimension in to rank 1 and dimension out.
+ * If $out < in$ the last elements will be discarded, if $out > in$ zeroes will
+ * be appended to fill the tensor.
+ */
+ template <int in, int out>
+ Tensor<1,out> convert_tensor (Tensor<1,in> old_tensor) const;
- /**
- * Returns spherical coordinates of a cartesian position.
- */
- Tensor<1,3> spherical_surface_coordinates(const Tensor<1,3> &position) const;
+ /**
+ * Returns spherical coordinates of a cartesian position.
+ */
+ Tensor<1,3> spherical_surface_coordinates(const Tensor<1,3> &position) const;
- /**
- * Return the cartesian coordinates of a spherical surface position
- * defined by theta (polar angle. not geographical latitude) and phi.
- */
- Tensor<1,3>
- cartesian_surface_coordinates(const Tensor<1,2> &sposition) const;
+ /**
+ * Return the cartesian coordinates of a spherical surface position
+ * defined by theta (polar angle. not geographical latitude) and phi.
+ */
+ Tensor<1,3>
+ cartesian_surface_coordinates(const Tensor<1,2> &sposition) const;
- /**
- * Returns cartesian velocities calculated from surface velocities and position in spherical coordinates
- *
- * @param s_velocities Surface velocities in spherical coordinates (theta, phi) @param s_position Position
- * in spherical coordinates (theta,phi,radius)
- */
- Tensor<1,3> sphere_to_cart_velocity(const Tensor<1,2> &s_velocities, const Tensor<1,3> &s_position) const;
+ /**
+ * Returns cartesian velocities calculated from surface velocities and position in spherical coordinates
+ *
+ * @param s_velocities Surface velocities in spherical coordinates (theta, phi) @param s_position Position
+ * in spherical coordinates (theta,phi,radius)
+ */
+ Tensor<1,3> sphere_to_cart_velocity(const Tensor<1,2> &s_velocities, const Tensor<1,3> &s_position) const;
- /**
- * calculates the phi-index given a certain phi
- */
- double get_idphi(const double phi_) const;
+ /**
+ * calculates the phi-index given a certain phi
+ */
+ double get_idphi(const double phi_) const;
- /**
- * calculates the theta-index given a certain polar angle
- */
- double get_idtheta(const double theta_) const;
+ /**
+ * calculates the theta-index given a certain polar angle
+ */
+ double get_idtheta(const double theta_) const;
};
}
Modified: trunk/aspect/source/velocity_boundary_conditions/gplates.cc
===================================================================
--- trunk/aspect/source/velocity_boundary_conditions/gplates.cc 2012-10-18 19:13:21 UTC (rev 1284)
+++ trunk/aspect/source/velocity_boundary_conditions/gplates.cc 2012-10-18 19:58:14 UTC (rev 1285)
@@ -39,287 +39,287 @@
{
namespace internal
{
- GPlatesLookup::GPlatesLookup(const Tensor<1,2> &surface_point_one, const Tensor<1,2> &surface_point_two)
- {
- velocity_vals.reinit(0,0);
- old_velocity_vals.reinit(0,0);
+ GPlatesLookup::GPlatesLookup(const Tensor<1,2> &surface_point_one, const Tensor<1,2> &surface_point_two)
+ {
+ velocity_vals.reinit(0,0);
+ old_velocity_vals.reinit(0,0);
- velocity_values = &velocity_vals;
- old_velocity_values = &old_velocity_vals;
+ velocity_values = &velocity_vals;
+ old_velocity_values = &old_velocity_vals;
- delta_phi = 0.0;
- delta_theta = 0.0;
+ delta_phi = 0.0;
+ delta_theta = 0.0;
- // get the cartesian coordinates of the points the 2D model shall lie in
- // this computation is done also for 3D since it is not expensive and the
- // template dim is currently not used here. Could be changed.
- const Tensor<1,3> point_one = cartesian_surface_coordinates(surface_point_one);
- const Tensor<1,3> point_two = cartesian_surface_coordinates(surface_point_two);
+ // get the cartesian coordinates of the points the 2D model shall lie in
+ // this computation is done also for 3D since it is not expensive and the
+ // template dim is currently not used here. Could be changed.
+ const Tensor<1,3> point_one = cartesian_surface_coordinates(surface_point_one);
+ const Tensor<1,3> point_two = cartesian_surface_coordinates(surface_point_two);
- //Set up the normal vector of an unrotated 2D spherical shell
- // that by default lies in the x-y plane.
- const double normal[3] = {0.0,0.0,1.0};
- const Tensor<1,3> unrotated_normal_vector (normal);
+ //Set up the normal vector of an unrotated 2D spherical shell
+ // that by default lies in the x-y plane.
+ const double normal[3] = {0.0,0.0,1.0};
+ const Tensor<1,3> unrotated_normal_vector (normal);
- // Compute the normal vector of the plane that contains
- // the origin and the two user-specified points
- Tensor<1,3> rotated_normal_vector;
- cross_product(rotated_normal_vector,point_one,point_two);
- rotated_normal_vector /= rotated_normal_vector.norm();
+ // Compute the normal vector of the plane that contains
+ // the origin and the two user-specified points
+ Tensor<1,3> rotated_normal_vector;
+ cross_product(rotated_normal_vector,point_one,point_two);
+ rotated_normal_vector /= rotated_normal_vector.norm();
- // Calculate the crossing line of the two normals,
- // which will be the rotation axis to transform the
- // into each other
- cross_product(rotation_axis,unrotated_normal_vector,rotated_normal_vector);
- rotation_axis /= rotation_axis.norm();
+ // Calculate the crossing line of the two normals,
+ // which will be the rotation axis to transform the
+ // into each other
+ cross_product(rotation_axis,unrotated_normal_vector,rotated_normal_vector);
+ rotation_axis /= rotation_axis.norm();
- // Calculate the rotation angle from the inner product rule
- rotation_angle = std::acos(rotated_normal_vector*unrotated_normal_vector);
+ // Calculate the rotation angle from the inner product rule
+ rotation_angle = std::acos(rotated_normal_vector*unrotated_normal_vector);
- }
+ }
- bool
- GPlatesLookup::fexists(const std::string &filename) const
- {
- std::ifstream ifile(filename.c_str());
- return ifile;
- }
+ bool
+ GPlatesLookup::fexists(const std::string &filename) const
+ {
+ std::ifstream ifile(filename.c_str());
+ return ifile;
+ }
- void
- GPlatesLookup::load_file(const std::string &filename)
- {
- using boost::property_tree::ptree;
- ptree pt;
+ void
+ GPlatesLookup::load_file(const std::string &filename)
+ {
+ using boost::property_tree::ptree;
+ ptree pt;
- // Check whether file exists, we do not want to throw
- // an exception in case it does not, because it could be by purpose
- // (i.e. the end of the boundary condition is reached)
- AssertThrow (fexists(filename),
- ExcMessage (std::string("GPlates file <")
- +
- filename
- +
- "> not found!"));
+ // Check whether file exists, we do not want to throw
+ // an exception in case it does not, because it could be by purpose
+ // (i.e. the end of the boundary condition is reached)
+ AssertThrow (fexists(filename),
+ ExcMessage (std::string("GPlates file <")
+ +
+ filename
+ +
+ "> not found!"));
- // populate tree structure pt
- read_xml(filename, pt);
+ // populate tree structure pt
+ read_xml(filename, pt);
- const int n_points = pt.get_child("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:domainSet.gml:MultiPoint").size();
- const int n_phi = static_cast<int>(std::sqrt(2.*n_points));
- const int n_theta = n_phi /2;
+ const int n_points = pt.get_child("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:domainSet.gml:MultiPoint").size();
+ const int n_phi = static_cast<int>(std::sqrt(2.*n_points));
+ const int n_theta = n_phi /2;
- if ((delta_theta != 0.0) || (delta_phi != 0.0))
- {
- Assert((delta_theta - numbers::PI / (n_theta-1)) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
- Assert((delta_phi - 2*numbers::PI / n_phi) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
- }
+ if ((delta_theta != 0.0) || (delta_phi != 0.0))
+ {
+ Assert((delta_theta - numbers::PI / (n_theta-1)) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
+ Assert((delta_phi - 2*numbers::PI / n_phi) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
+ }
- delta_theta = numbers::PI / (n_theta-1);
- delta_phi = 2*numbers::PI / n_phi;
+ delta_theta = numbers::PI / (n_theta-1);
+ delta_phi = 2*numbers::PI / n_phi;
- // swap pointers to old and new values, we overwrite the old ones
- // and the new ones become the old ones
- std::swap( velocity_values, old_velocity_values);
+ // swap pointers to old and new values, we overwrite the old ones
+ // and the new ones become the old ones
+ std::swap( velocity_values, old_velocity_values);
- (*velocity_values).reinit(n_theta,n_phi);
+ (*velocity_values).reinit(n_theta,n_phi);
- std::string velos = pt.get<std::string>("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:rangeSet.gml:DataBlock.gml:tupleList");
- std::stringstream in(velos, std::ios::in);
- AssertThrow (in,
- ExcMessage (std::string("Couldn't find velocities. Is file native gpml format for velocities?")));
+ std::string velos = pt.get<std::string>("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:rangeSet.gml:DataBlock.gml:tupleList");
+ std::stringstream in(velos, std::ios::in);
+ AssertThrow (in,
+ ExcMessage (std::string("Couldn't find velocities. Is file native gpml format for velocities?")));
- unsigned int i = 0;
- char sep;
- while (!in.eof())
- {
- double vtheta, vphi;
- const double cmyr_si = 3.1557e9;
+ unsigned int i = 0;
+ char sep;
+ while (!in.eof())
+ {
+ double vtheta, vphi;
+ const double cmyr_si = 3.1557e9;
- in >> vtheta >> sep >> vphi;
+ in >> vtheta >> sep >> vphi;
- if (in.eof())
- break;
+ if (in.eof())
+ break;
- (*velocity_values)[i%n_theta][i/n_theta][0] = vtheta/cmyr_si;
- (*velocity_values)[i%n_theta][i/n_theta][1] = vphi/cmyr_si;
+ (*velocity_values)[i%n_theta][i/n_theta][0] = vtheta/cmyr_si;
+ (*velocity_values)[i%n_theta][i/n_theta][1] = vphi/cmyr_si;
- i++;
- }
+ i++;
}
+ }
- template <int dim>
- Tensor<1,dim>
- GPlatesLookup::surface_velocity(const Point<dim> &position_, const double time_weight) const
- {
+ template <int dim>
+ Tensor<1,dim>
+ GPlatesLookup::surface_velocity(const Point<dim> &position_, const double time_weight) const
+ {
- Tensor<1,dim> tensor_position;
- for (unsigned int i = 0 ; i < dim; i++) tensor_position[i] = position_[i];
+ Tensor<1,dim> tensor_position;
+ for (unsigned int i = 0 ; i < dim; i++) tensor_position[i] = position_[i];
- Tensor<1,3> position;
- if (dim == 2)
- {
- position = rotate(convert_tensor<dim,3>(tensor_position),rotation_axis,rotation_angle);
- }
- else
- position = convert_tensor<dim,3>(tensor_position);
+ Tensor<1,3> position;
+ if (dim == 2)
+ {
+ position = rotate(convert_tensor<dim,3>(tensor_position),rotation_axis,rotation_angle);
+ }
+ else
+ position = convert_tensor<dim,3>(tensor_position);
- const Tensor<1,3> scoord = spherical_surface_coordinates(position);
+ const Tensor<1,3> scoord = spherical_surface_coordinates(position);
- const double idtheta = get_idtheta(scoord[0]);
- const unsigned int idxtheta = static_cast<unsigned int>(idtheta);
- const unsigned int nexttheta = idxtheta + 1;
+ const double idtheta = get_idtheta(scoord[0]);
+ const unsigned int idxtheta = static_cast<unsigned int>(idtheta);
+ const unsigned int nexttheta = idxtheta + 1;
- const double idphi = get_idphi(scoord[1]);
- const unsigned int idxphi = static_cast<unsigned int>(idphi);
- const unsigned int nextphi = (idxphi+1) % velocity_values->n_cols();
+ const double idphi = get_idphi(scoord[1]);
+ const unsigned int idxphi = static_cast<unsigned int>(idphi);
+ const unsigned int nextphi = (idxphi+1) % velocity_values->n_cols();
- Assert(idxtheta<velocity_values->n_rows(), ExcMessage("not in range"));
- Assert(idxphi <velocity_values->n_cols(), ExcMessage("not in range"));
+ Assert(idxtheta<velocity_values->n_rows(), ExcMessage("not in range"));
+ Assert(idxphi <velocity_values->n_cols(), ExcMessage("not in range"));
- Assert(nexttheta<velocity_values->n_rows(), ExcMessage("not in range"));
- Assert(nextphi <velocity_values->n_cols(), ExcMessage("not in range"));
+ Assert(nexttheta<velocity_values->n_rows(), ExcMessage("not in range"));
+ Assert(nextphi <velocity_values->n_cols(), ExcMessage("not in range"));
- // compute the coordinates of this point in the
- // reference cell between the data points
- const double xi = idtheta-static_cast<double>(idxtheta);
- const double eta = idphi-static_cast<double>(idxphi);
+ // compute the coordinates of this point in the
+ // reference cell between the data points
+ const double xi = idtheta-static_cast<double>(idxtheta);
+ const double eta = idphi-static_cast<double>(idxphi);
- Assert ((0 <= xi) && (xi <= 1), ExcInternalError());
- Assert ((0 <= eta) && (eta <= 1), ExcInternalError());
- Assert ((0 <= time_weight) && (time_weight <= 1), ExcInternalError());
+ Assert ((0 <= xi) && (xi <= 1), ExcInternalError());
+ Assert ((0 <= eta) && (eta <= 1), ExcInternalError());
+ Assert ((0 <= time_weight) && (time_weight <= 1), ExcInternalError());
- // use xi, eta and time_weight for a trilinear interpolation
- // TODO: interpolation in cartesian probably more accurate
+ // use xi, eta and time_weight for a trilinear interpolation
+ // TODO: interpolation in cartesian probably more accurate
- Tensor<1,2> surf_vel;
+ Tensor<1,2> surf_vel;
- surf_vel[0] = time_weight *
- ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][0] +
- xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][0] +
- (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][0] +
- xi *eta *(*velocity_values)[nexttheta][nextphi][0]);
+ surf_vel[0] = time_weight *
+ ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][0] +
+ xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][0] +
+ (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][0] +
+ xi *eta *(*velocity_values)[nexttheta][nextphi][0]);
- surf_vel[1] = time_weight *
- ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][1] +
- xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][1] +
- (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][1] +
- xi *eta *(*velocity_values)[nexttheta][nextphi][1]);
+ surf_vel[1] = time_weight *
+ ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][1] +
+ xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][1] +
+ (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][1] +
+ xi *eta *(*velocity_values)[nexttheta][nextphi][1]);
- surf_vel[0] += (1-time_weight) *
- ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][0] +
- xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][0] +
- (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][0] +
- xi *eta *(*old_velocity_values)[nexttheta][nextphi][0]);
+ surf_vel[0] += (1-time_weight) *
+ ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][0] +
+ xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][0] +
+ (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][0] +
+ xi *eta *(*old_velocity_values)[nexttheta][nextphi][0]);
- surf_vel[1] += (1-time_weight) *
- ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][1] +
- xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][1] +
- (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][1] +
- xi *eta *(*old_velocity_values)[nexttheta][nextphi][1]);
+ surf_vel[1] += (1-time_weight) *
+ ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][1] +
+ xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][1] +
+ (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][1] +
+ xi *eta *(*old_velocity_values)[nexttheta][nextphi][1]);
- const Tensor<1,3> cart_velo = sphere_to_cart_velocity(surf_vel,scoord);
- Tensor<1,dim> velos;
+ const Tensor<1,3> cart_velo = sphere_to_cart_velocity(surf_vel,scoord);
+ Tensor<1,dim> velos;
- if (dim == 2)
- velos = convert_tensor<3,dim>(rotate(cart_velo,rotation_axis,-1.0*rotation_angle));
- else
- velos = convert_tensor<3,dim>(cart_velo);
+ if (dim == 2)
+ velos = convert_tensor<3,dim>(rotate(cart_velo,rotation_axis,-1.0*rotation_angle));
+ else
+ velos = convert_tensor<3,dim>(cart_velo);
- return velos;
- }
+ return velos;
+ }
- Tensor<1,3>
- GPlatesLookup::rotate (const Tensor<1,3> &position, const Tensor<1,3> &rotation_axis, const double angle) const
- {
- Tensor<1,3> cross;
- cross_product(cross,position,rotation_axis);
- Tensor<1,3> newpos = (1-std::cos(angle)) * rotation_axis*(rotation_axis*position) +
- std::cos(angle) * position + std::sin(angle) * cross;
- return newpos;
- }
+ Tensor<1,3>
+ GPlatesLookup::rotate (const Tensor<1,3> &position, const Tensor<1,3> &rotation_axis, const double angle) const
+ {
+ Tensor<1,3> cross;
+ cross_product(cross,position,rotation_axis);
+ Tensor<1,3> newpos = (1-std::cos(angle)) * rotation_axis*(rotation_axis*position) +
+ std::cos(angle) * position + std::sin(angle) * cross;
+ return newpos;
+ }
- Tensor<1,3>
- GPlatesLookup::spherical_surface_coordinates(const Tensor<1,3> &position) const
- {
- Tensor<1,3> scoord;
+ Tensor<1,3>
+ GPlatesLookup::spherical_surface_coordinates(const Tensor<1,3> &position) const
+ {
+ Tensor<1,3> scoord;
- scoord[0] = std::acos(position[2]/std::sqrt(position.norm_square())); // Theta
- scoord[1] = std::atan2(position[1],position[0]); // Phi
- if (scoord[1] < 0.0) scoord[1] = 2*numbers::PI + scoord[1]; // correct phi to [0,2*pi]
- scoord[2] = std::sqrt(position.norm_square()); // R
- return scoord;
- }
+ scoord[0] = std::acos(position[2]/std::sqrt(position.norm_square())); // Theta
+ scoord[1] = std::atan2(position[1],position[0]); // Phi
+ if (scoord[1] < 0.0) scoord[1] = 2*numbers::PI + scoord[1]; // correct phi to [0,2*pi]
+ scoord[2] = std::sqrt(position.norm_square()); // R
+ return scoord;
+ }
- Tensor<1,3>
- GPlatesLookup::cartesian_surface_coordinates(const Tensor<1,2> &sposition) const
- {
- Tensor<1,3> ccoord;
+ Tensor<1,3>
+ GPlatesLookup::cartesian_surface_coordinates(const Tensor<1,2> &sposition) const
+ {
+ Tensor<1,3> ccoord;
- ccoord[0] = std::sin(sposition[0] * std::cos(sposition[1])); // X
- ccoord[1] = std::sin(sposition[0] * std::sin(sposition[1])); // Y
- ccoord[2] = std::cos(sposition[0]); // Z
- return ccoord;
- }
+ ccoord[0] = std::sin(sposition[0] * std::cos(sposition[1])); // X
+ ccoord[1] = std::sin(sposition[0] * std::sin(sposition[1])); // Y
+ ccoord[2] = std::cos(sposition[0]); // Z
+ return ccoord;
+ }
- Tensor<1,3>
- GPlatesLookup::sphere_to_cart_velocity(const Tensor<1,2> &s_velocities, const Tensor<1,3> &s_position) const
- {
- Tensor<1,3> velocity;
+ Tensor<1,3>
+ GPlatesLookup::sphere_to_cart_velocity(const Tensor<1,2> &s_velocities, const Tensor<1,3> &s_position) const
+ {
+ Tensor<1,3> velocity;
- velocity[0] = -1.0 * s_velocities[1] * std::sin(s_position[1]) + s_velocities[0]*std::cos(s_position[0])*std::cos(s_position[1]);
- velocity[1] = s_velocities[1]*std::cos(s_position[1])+s_velocities[0]*std::cos(s_position[0])*std::sin(s_position[1]);
- velocity[2] = -1.0*s_velocities[0]*std::sin(s_position[0]);
- return velocity;
- }
+ velocity[0] = -1.0 * s_velocities[1] * std::sin(s_position[1]) + s_velocities[0]*std::cos(s_position[0])*std::cos(s_position[1]);
+ velocity[1] = s_velocities[1]*std::cos(s_position[1])+s_velocities[0]*std::cos(s_position[0])*std::sin(s_position[1]);
+ velocity[2] = -1.0*s_velocities[0]*std::sin(s_position[0]);
+ return velocity;
+ }
- template <int in, int out>
- Tensor<1,out>
- GPlatesLookup::convert_tensor (Tensor<1,in> old_tensor) const
- {
- Tensor<1,out> new_tensor;
- for (unsigned int i = 0; i < out; i++)
- if (i < in) new_tensor[i] = old_tensor[i];
- else new_tensor[i] = 0.0;
+ template <int in, int out>
+ Tensor<1,out>
+ GPlatesLookup::convert_tensor (Tensor<1,in> old_tensor) const
+ {
+ Tensor<1,out> new_tensor;
+ for (unsigned int i = 0; i < out; i++)
+ if (i < in) new_tensor[i] = old_tensor[i];
+ else new_tensor[i] = 0.0;
- return new_tensor;
- }
+ return new_tensor;
+ }
- double
- GPlatesLookup::get_idphi(const double phi_) const
- {
- const double phi = std::max(std::min(phi_,2*numbers::PI-1e-7),0.0);
+ double
+ GPlatesLookup::get_idphi(const double phi_) const
+ {
+ const double phi = std::max(std::min(phi_,2*numbers::PI-1e-7),0.0);
- Assert(phi>=0, ExcMessage("not in range"));
- Assert(phi<=2*numbers::PI, ExcMessage("not in range"));
- return phi/delta_phi;
- }
+ Assert(phi>=0, ExcMessage("not in range"));
+ Assert(phi<=2*numbers::PI, ExcMessage("not in range"));
+ return phi/delta_phi;
+ }
- double
- GPlatesLookup::get_idtheta(const double theta_) const
- {
- double theta = std::max(std::min(theta_,numbers::PI-1e-7),0.0);
+ double
+ GPlatesLookup::get_idtheta(const double theta_) const
+ {
+ double theta = std::max(std::min(theta_,numbers::PI-1e-7),0.0);
- Assert(theta>=0, ExcMessage("not in range"));
- Assert(theta<=numbers::PI, ExcMessage("not in range"));
- return theta/delta_theta;
- }
+ Assert(theta>=0, ExcMessage("not in range"));
+ Assert(theta<=numbers::PI, ExcMessage("not in range"));
+ return theta/delta_theta;
+ }
}
template <int dim>
GPlates<dim>::GPlates ()
- :
- current_time(0.0),
- current_time_step(0),
- time_step(0.0),
- time_weight(0.0),
- time_dependent(true),
- point1("0.0,0.0"),
- point2("0.0,0.0"),
- lookup(NULL)
+ :
+ current_time(0.0),
+ current_time_step(0),
+ time_step(0.0),
+ time_weight(0.0),
+ time_dependent(true),
+ point1("0.0,0.0"),
+ point2("0.0,0.0"),
+ lookup(NULL)
{}
@@ -397,11 +397,11 @@
// Load next velocity file for interpolation
try
{
- // If the time step was large enough to move forward more
- // then one velocity file, we need to load both current files
- // to stay accurate in interpolation
- if (current_time_step > old_time_step + 1)
- lookup->load_file(create_filename(current_time_step));
+ // If the time step was large enough to move forward more
+ // then one velocity file, we need to load both current files
+ // to stay accurate in interpolation
+ if (current_time_step > old_time_step + 1)
+ lookup->load_file(create_filename(current_time_step));
lookup->load_file(create_filename(current_time_step+1));
}
More information about the CIG-COMMITS
mailing list