[cig-commits] [commit] master: Bugfix and additional rotation angle output. (6b56730)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Aug 18 03:47:30 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/3adfd87dec1546e9cd43c7d0a85230cc6d56c458...43f425f57ef5959f8b3b406b08aea3cf203cd979
>---------------------------------------------------------------
commit 6b5673034a9230ca4842e2487db478c2c2ec30c4
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date: Tue Aug 12 16:05:19 2014 +0200
Bugfix and additional rotation angle output.
>---------------------------------------------------------------
6b5673034a9230ca4842e2487db478c2c2ec30c4
doc/modules/changes.h | 10 +
.../aspect/velocity_boundary_conditions/gplates.h | 39 +++-
source/velocity_boundary_conditions/gplates.cc | 207 +++++++++++++++++++--
tests/gplates-1-3/screen-output | 28 +--
tests/gplates-1-4/screen-output | 28 +--
tests/{gplates-1-3.prm => gplates-rotation.prm} | 7 +-
tests/gplates-rotation/screen-output | 44 +++++
tests/{gplates-1-3 => gplates-rotation}/statistics | 2 +-
8 files changed, 293 insertions(+), 72 deletions(-)
diff --git a/doc/modules/changes.h b/doc/modules/changes.h
index c783964..fcb193c 100644
--- a/doc/modules/changes.h
+++ b/doc/modules/changes.h
@@ -6,6 +6,16 @@
*
*
* <ol>
+ * <li> Changed: The GPlates velocity boundary plugin prescribed only the normal
+ * vector of the model plane of 2D models. Therefore it was very hard to rotate
+ * the model into its actual orientation, when comparing it with other datasets.
+ * Now the first of the user input points is always rotated to a known position
+ * (0,1,0), and the plugin outputs the according rotation angles to rotate the
+ * model into its proper orientation and the inverse angles to rotate other
+ * datasets into the model plane.
+ * <br>
+ * (Rene Gassmoeller, 2014/08/15)
+ *
* <li> New: There are numerous places in the input file where one can or
* has to input boundary indicators. These indicators identify individual
* parts of the boundary of the domain, such as the inner or outer surfaces
diff --git a/include/aspect/velocity_boundary_conditions/gplates.h b/include/aspect/velocity_boundary_conditions/gplates.h
index a6be311..0b2cf37 100644
--- a/include/aspect/velocity_boundary_conditions/gplates.h
+++ b/include/aspect/velocity_boundary_conditions/gplates.h
@@ -23,6 +23,7 @@
#define __aspect__velocity_boundary_conditions_gplates_h
#include <aspect/velocity_boundary_conditions/interface.h>
+#include <deal.II/base/std_cxx1x/array.h>
#include <aspect/simulator_access.h>
@@ -117,12 +118,11 @@ namespace aspect
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 matrix, which describes the rotation by which a 2D model needs
+ * 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;
+ Tensor<2,3> rotation_matrix;
/**
* Determines the width of the velocity interpolation zone around
@@ -140,9 +140,34 @@ namespace aspect
* defined angle
*/
Tensor<1,3>
- rotate (const Tensor<1,3> &position,
- const Tensor<1,3> &rotation_axis,
- const double angle) const;
+ rotate_around_axis (const Tensor<1,3> &position,
+ const Tensor<1,3> &rotation_axis,
+ const double angle) const;
+
+ /**
+ * A function that returns the corresponding paraview angles for a
+ * rotation described by a rotation matrix. These differ from the
+ * usually used euler angles by assuming a rotation around the coordinate
+ * axes in the order y-x-z (instead of the often used z-x-z)
+ */
+ std_cxx1x::array<double,3>
+ angles_from_matrix (const Tensor<2,3> &rotation_matrix) const;
+
+ /**
+ * A function that returns the corresponding rotation axis/angle for a
+ * rotation described by a rotation matrix.
+ */
+ double
+ rotation_axis_from_matrix (Tensor<1,3> &rotation_axis,
+ const Tensor<2,3> &rotation_matrix) const;
+
+ /**
+ * A function that returns the corresponding euler angles for a
+ * rotation described by rotation axis and angle.
+ */
+ Tensor<2,3>
+ rotation_matrix_from_axis (const Tensor<1,3> &rotation_axis,
+ const double rotation_angle) const;
/**
* Convert a tensor of rank 1 and dimension in to rank 1 and
diff --git a/source/velocity_boundary_conditions/gplates.cc b/source/velocity_boundary_conditions/gplates.cc
index 76d4aab..a783c91 100644
--- a/source/velocity_boundary_conditions/gplates.cc
+++ b/source/velocity_boundary_conditions/gplates.cc
@@ -53,13 +53,13 @@ namespace aspect
delta_theta(0.0),
interpolation_width(interpolation_width_)
{
- // get the Cartesian coordinates of the points the 2D model shall lie in
+ // get the Cartesian coordinates of the points the 2D model will 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(convert_tensor<2,3>(surface_point_one));
const Tensor<1,3> point_two = cartesian_surface_coordinates(convert_tensor<2,3>(surface_point_two));
- //Set up the normal vector of an unrotated 2D spherical shell
+ // 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);
@@ -68,23 +68,49 @@ namespace aspect
// 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();
if ((rotated_normal_vector - unrotated_normal_vector).norm() > 1e-3)
{
// Calculate the crossing line of the two normals,
- // which will be the rotation axis to transform them
- // into each other
+ // which will be the rotation axis to transform the one
+ // normal into the other
+ Tensor<1,3> rotation_axis;
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);
+ const double rotation_angle = std::acos(rotated_normal_vector*unrotated_normal_vector);
+
+ rotation_matrix = rotation_matrix_from_axis(rotation_axis,rotation_angle);
+
+ // Now apply the rotation that will project point_one onto the known point
+ // (0,1,0).
+ const Tensor<1,3> rotated_point_one = transpose(rotation_matrix) * point_one;
+
+ const double point_one_coords[3] = {0.0,1.0,0.0};
+ const Tensor<1,3> final_point_one (point_one_coords);
+
+ const double second_rotation_angle = std::acos(rotated_point_one*final_point_one);
+ Tensor<1,3> second_rotation_axis;
+ cross_product(second_rotation_axis,final_point_one,rotated_point_one);
+ second_rotation_axis /= second_rotation_axis.norm();
+
+ const Tensor<2,3> second_rotation_matrix = rotation_matrix_from_axis(second_rotation_axis,second_rotation_angle);
+
+ // The final rotation used for the model will be the combined
+ // rotation of the two operation above. This is achieved by a
+ // matrix multiplication of the rotation matrices.
+ // This concatenation of rotations is the reason for using a
+ // rotation matrix instead of a combined rotation_axis + angle
+ rotation_matrix = rotation_matrix * second_rotation_matrix;
}
else
{
- rotation_axis = unrotated_normal_vector;
- rotation_angle = 0.0;
+ rotation_matrix[0][0] = 1.0;
+ rotation_matrix[1][1] = 1.0;
+ rotation_matrix[2][2] = 1.0;
}
}
@@ -103,13 +129,24 @@ namespace aspect
<< std::endl;
if (dim == 2)
{
- output << " Input point 1 spherical coordinates: " << surface_point_one << std::endl
- << " Input point 1 Cartesian output coordinates: " << rotate(point_one,rotation_axis,-rotation_angle) << std::endl
- << " Input point 2 spherical coordinates: " << surface_point_two << std::endl
- << " Input point 2 Cartesian output coordinates: " << rotate(point_two,rotation_axis,-rotation_angle) << std::endl
+ Tensor<1,3> rotation_axis;
+ const double rotation_angle = rotation_axis_from_matrix(rotation_axis,rotation_matrix);
+
+ std_cxx1x::array<double,3> angles = angles_from_matrix(rotation_matrix);
+ std_cxx1x::array<double,3> back_angles = angles_from_matrix(transpose(rotation_matrix));
+
+ output << " Input point 1 spherical coordinates: " << surface_point_one << std::endl
+ << " Input point 1 normalized cartesian coordinates: " << point_one << std::endl
+ << " Input point 1 rotated model coordinates: " << transpose(rotation_matrix) * point_one << std::endl
+ << " Input point 2 spherical coordinates: " << surface_point_two << std::endl
+ << " Input point 2 normalized cartesian coordinates: " << point_two << std::endl
+ << " Input point 2 rotated model coordinates: " << transpose(rotation_matrix) * point_two << std::endl
<< std::endl << std::setprecision(2)
<< " Model will be rotated by " << -rotation_angle*180/numbers::PI
<< " degrees around axis " << rotation_axis << std::endl
+ << " The ParaView rotation angles are: " << angles[0] << " " << angles [1] << " " << angles[2] << std::endl
+ << " The inverse ParaView rotation angles are: " << back_angles[0] << " " << back_angles [1] << " " << back_angles[2]
+
<< std::endl;
}
@@ -221,7 +258,7 @@ namespace aspect
{
Tensor<1,3> internal_position;
if (dim == 2)
- internal_position = rotate(convert_tensor<dim,3>(position),rotation_axis,rotation_angle);
+ internal_position = rotation_matrix * convert_tensor<dim,3>(position);
else
internal_position = convert_tensor<dim,3>(position);
@@ -233,7 +270,7 @@ namespace aspect
if (dim == 2)
// convert_tensor conveniently also handles the projection to the 2D plane by
// omitting the z-component of velocity (since the 2D model lies in the x-y plane).
- output_boundary_velocity = convert_tensor<3,dim>(rotate(interpolated_velocity,rotation_axis,-1.0*rotation_angle));
+ output_boundary_velocity = convert_tensor<3,dim>(transpose(rotation_matrix) * interpolated_velocity);
else
output_boundary_velocity = convert_tensor<3,dim>(interpolated_velocity);
@@ -397,7 +434,7 @@ namespace aspect
// Calculate the rotation angle from the inner product rule
const double local_rotation_angle = std::acos(data_position*point_position);
- const Tensor<1,3> point_velocity = rotate(data_velocity,local_rotation_axis,-local_rotation_angle);
+ const Tensor<1,3> point_velocity = rotate_around_axis(data_velocity,local_rotation_axis,local_rotation_angle);
return point_velocity;
}
@@ -459,10 +496,10 @@ namespace aspect
}
Tensor<1,3>
- GPlatesLookup::rotate (const Tensor<1,3> &position, const Tensor<1,3> &rotation_axis, const double angle) const
+ GPlatesLookup::rotate_around_axis (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);
+ cross_product(cross,rotation_axis,position);
const Tensor<1,3> newpos = (1-std::cos(angle)) * rotation_axis*(rotation_axis*position) +
std::cos(angle) * position + std::sin(angle) * cross;
return newpos;
@@ -490,6 +527,144 @@ namespace aspect
return velocity;
}
+ Tensor<2,3>
+ GPlatesLookup::rotation_matrix_from_axis (const Tensor<1,3> &rotation_axis,
+ const double rotation_angle) const
+ {
+ Tensor<2,3> rotation_matrix;
+ rotation_matrix[0][0] = (1-std::cos(rotation_angle)) * rotation_axis[0]*rotation_axis[0] + std::cos(rotation_angle);
+ rotation_matrix[0][1] = (1-std::cos(rotation_angle)) * rotation_axis[0]*rotation_axis[1] - rotation_axis[2] * std::sin(rotation_angle);
+ rotation_matrix[0][2] = (1-std::cos(rotation_angle)) * rotation_axis[0]*rotation_axis[2] + rotation_axis[1] * std::sin(rotation_angle);
+ rotation_matrix[1][0] = (1-std::cos(rotation_angle)) * rotation_axis[1]*rotation_axis[0] + rotation_axis[2] * std::sin(rotation_angle);
+ rotation_matrix[1][1] = (1-std::cos(rotation_angle)) * rotation_axis[1]*rotation_axis[1] + std::cos(rotation_angle);
+ rotation_matrix[1][2] = (1-std::cos(rotation_angle)) * rotation_axis[1]*rotation_axis[2] - rotation_axis[0] * std::sin(rotation_angle);
+ rotation_matrix[2][0] = (1-std::cos(rotation_angle)) * rotation_axis[2]*rotation_axis[0] - rotation_axis[1] * std::sin(rotation_angle);
+ rotation_matrix[2][1] = (1-std::cos(rotation_angle)) * rotation_axis[2]*rotation_axis[1] + rotation_axis[0] * std::sin(rotation_angle);
+ rotation_matrix[2][2] = (1-std::cos(rotation_angle)) * rotation_axis[2]*rotation_axis[2] + std::cos(rotation_angle);
+ return rotation_matrix;
+ }
+
+ double
+ GPlatesLookup::rotation_axis_from_matrix (Tensor<1,3> &rotation_axis,
+ const Tensor<2,3> &rotation_matrix) const
+ {
+ double rotation_angle = std::acos(0.5 * (rotation_matrix[0][0] + rotation_matrix[1][1] + rotation_matrix[2][2] - 1));
+
+ if (rotation_angle > std::numeric_limits<double>::min())
+ {
+ rotation_axis[0] = (rotation_matrix[2][1] - rotation_matrix[1][2]) / (2*std::sin(rotation_angle));
+ rotation_axis[1] = (rotation_matrix[0][2] - rotation_matrix[2][0]) / (2*std::sin(rotation_angle));
+ rotation_axis[2] = (rotation_matrix[1][0] - rotation_matrix[0][1]) / (2*std::sin(rotation_angle));
+ }
+ else
+ {
+ rotation_axis[0] = 0.0;
+ rotation_axis[1] = 0.0;
+ rotation_axis[2] = 1.0;
+ }
+
+ return rotation_angle;
+ }
+
+ std_cxx1x::array<double,3>
+ GPlatesLookup::angles_from_matrix(const Tensor<2,3> &rotation_matrix) const
+ {
+ std_cxx1x::array<double,3> orientation;
+
+ /*
+ * The following code is part of the VTK project and copied here for
+ * compatibility to the paraview rotation formalism. It is protected by
+ * the following license:
+ *
+ *
+ * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted under certain conditions. See
+ * http://www.kitware.com/Copyright.htm for details.
+
+ * This software is distributed WITHOUT ANY WARRANTY; without even
+ * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ * PURPOSE. See the above copyright notice for more information.
+ *
+ *
+ * The original code in the VTK-6.0.0 source folder is found in:
+ * /Common/Transforms/vtkTransform.cxx in the function GetOrientation()
+ *
+ */
+
+ // first rotate about y axis
+ const double x2 = rotation_matrix[2][0];
+ const double y2 = rotation_matrix[2][1];
+ const double z2 = rotation_matrix[2][2];
+
+ const double x3 = rotation_matrix[1][0];
+ const double y3 = rotation_matrix[1][1];
+ const double z3 = rotation_matrix[1][2];
+
+ double d1 = sqrt(x2*x2 + z2*z2);
+
+ double cosTheta, sinTheta;
+ if (d1 < std::numeric_limits<double>::min())
+ {
+ cosTheta = 1.0;
+ sinTheta = 0.0;
+ }
+ else
+ {
+ cosTheta = z2/d1;
+ sinTheta = x2/d1;
+ }
+
+ double theta = atan2(sinTheta, cosTheta);
+ orientation[1] = - theta * 180 / numbers::PI;
+
+ // now rotate about x axis
+ double d = sqrt(x2*x2 + y2*y2 + z2*z2);
+
+ double sinPhi, cosPhi;
+ if (d < std::numeric_limits<double>::min())
+ {
+ sinPhi = 0.0;
+ cosPhi = 1.0;
+ }
+ else if (d1 < std::numeric_limits<double>::min())
+ {
+ sinPhi = y2/d;
+ cosPhi = z2/d;
+ }
+ else
+ {
+ sinPhi = y2/d;
+ cosPhi = (x2*x2 + z2*z2)/(d1*d);
+ }
+
+ double phi = atan2(sinPhi, cosPhi);
+ orientation[0] = phi * 180 / numbers::PI;
+
+ // finally, rotate about z
+ double x3p = x3*cosTheta - z3*sinTheta;
+ double y3p = - sinPhi*sinTheta*x3 + cosPhi*y3 - sinPhi*cosTheta*z3;
+ double d2 = sqrt(x3p*x3p + y3p*y3p);
+
+ double cosAlpha, sinAlpha;
+ if (d2 < std::numeric_limits<double>::min())
+ {
+ cosAlpha = 1.0;
+ sinAlpha = 0.0;
+ }
+ else
+ {
+ cosAlpha = y3p/d2;
+ sinAlpha = x3p/d2;
+ }
+
+ double alpha = atan2(sinAlpha, cosAlpha);
+ orientation[2] = alpha * 180 / numbers::PI;
+ return orientation;
+ }
+
template <int in, int out>
Tensor<1,out>
GPlatesLookup::convert_tensor (const Tensor<1,in> &old_tensor) const
diff --git a/tests/gplates-1-3/screen-output b/tests/gplates-1-3/screen-output
index 5e9ecd4..b4de6f6 100644
--- a/tests/gplates-1-3/screen-output
+++ b/tests/gplates-1-3/screen-output
@@ -1,10 +1,3 @@
------------------------------------------------------------------------------
--- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
--- . version 1.1.pre
--- . running in DEBUG mode
--- . running with 1 MPI process
--- . using Trilinos
------------------------------------------------------------------------------
Number of active cells: 768 (on 4 levels)
Number of degrees of freedom: 10,656 (6,528+864+3,264)
@@ -13,12 +6,15 @@ Number of degrees of freedom: 10,656 (6,528+864+3,264)
Setting up GPlates boundary velocity plugin.
Input point 1 spherical coordinates: 1.571 4.870
- Input point 1 Cartesian output coordinates: 0.157 -0.988 0.000
+ Input point 1 normalized cartesian coordinates: 0.157 -0.988 0.000
+ Input point 1 rotated model coordinates: 0.157 -0.988 0.000
Input point 2 spherical coordinates: 1.571 5.240
- Input point 2 Cartesian output coordinates: 0.503 -0.864 0.000
+ Input point 2 normalized cartesian coordinates: 0.503 -0.864 0.000
+ Input point 2 rotated model coordinates: 0.503 -0.864 0.000
Model will be rotated by 0.00 degrees around axis 0.00 0.00 1.00
-
+ The ParaView rotation angles are: 0.00 0.00 0.00
+ The inverse ParaView rotation angles are: 0.00 0.00 0.00
Loading GPlates boundary velocity file ASPECT_DIR/data/velocity-boundary-conditions/gplates/current_day.gpml.
@@ -43,18 +39,6 @@ Termination requested by criterion: end time
+---------------------------------------------+------------+------------+
-| Total wallclock time elapsed since start | 1.58s | |
-| | | |
-| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
-| Assemble Stokes system | 1 | 0.225s | 14% |
-| Assemble temperature system | 1 | 0.24s | 15% |
-| Build Stokes preconditioner | 1 | 0.243s | 15% |
-| Build temperature preconditioner| 1 | 0.0112s | 0.71% |
-| Solve Stokes system | 1 | 0.486s | 31% |
-| Solve temperature system | 1 | 0.0014s | 0.089% |
-| Initialization | 2 | 0.122s | 7.8% |
-| Postprocessing | 1 | 0.031s | 2% |
-| Setup dof systems | 1 | 0.101s | 6.4% |
+---------------------------------+-----------+------------+------------+
diff --git a/tests/gplates-1-4/screen-output b/tests/gplates-1-4/screen-output
index 5ae0392..88a61ac 100644
--- a/tests/gplates-1-4/screen-output
+++ b/tests/gplates-1-4/screen-output
@@ -1,10 +1,3 @@
------------------------------------------------------------------------------
--- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
--- . version 1.1.pre
--- . running in DEBUG mode
--- . running with 1 MPI process
--- . using Trilinos
------------------------------------------------------------------------------
Number of active cells: 768 (on 4 levels)
Number of degrees of freedom: 10,656 (6,528+864+3,264)
@@ -13,12 +6,15 @@ Number of degrees of freedom: 10,656 (6,528+864+3,264)
Setting up GPlates boundary velocity plugin.
Input point 1 spherical coordinates: 1.571 4.870
- Input point 1 Cartesian output coordinates: 0.157 -0.988 0.000
+ Input point 1 normalized cartesian coordinates: 0.157 -0.988 0.000
+ Input point 1 rotated model coordinates: 0.157 -0.988 0.000
Input point 2 spherical coordinates: 1.571 5.240
- Input point 2 Cartesian output coordinates: 0.503 -0.864 0.000
+ Input point 2 normalized cartesian coordinates: 0.503 -0.864 0.000
+ Input point 2 rotated model coordinates: 0.503 -0.864 0.000
Model will be rotated by 0.00 degrees around axis 0.00 0.00 1.00
-
+ The ParaView rotation angles are: 0.00 0.00 0.00
+ The inverse ParaView rotation angles are: 0.00 0.00 0.00
Loading GPlates boundary velocity file ASPECT_DIR/data/velocity-boundary-conditions/gplates/current_day_1.4.gpml.
@@ -43,18 +39,6 @@ Termination requested by criterion: end time
+---------------------------------------------+------------+------------+
-| Total wallclock time elapsed since start | 1.93s | |
-| | | |
-| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
-| Assemble Stokes system | 1 | 0.404s | 21% |
-| Assemble temperature system | 1 | 0.3s | 16% |
-| Build Stokes preconditioner | 1 | 0.308s | 16% |
-| Build temperature preconditioner| 1 | 0.0107s | 0.55% |
-| Solve Stokes system | 1 | 0.494s | 26% |
-| Solve temperature system | 1 | 0.0134s | 0.7% |
-| Initialization | 2 | 0.158s | 8.2% |
-| Postprocessing | 1 | 0.031s | 1.6% |
-| Setup dof systems | 1 | 0.0995s | 5.2% |
+---------------------------------+-----------+------------+------------+
diff --git a/tests/gplates-1-3.prm b/tests/gplates-rotation.prm
similarity index 90%
copy from tests/gplates-1-3.prm
copy to tests/gplates-rotation.prm
index cb86c9c..817a983 100644
--- a/tests/gplates-1-3.prm
+++ b/tests/gplates-rotation.prm
@@ -1,5 +1,5 @@
# A simple setup for for using the GPlates interface in a 2d shell.
-# This test uses the grid produced by GPlates 1.3
+# This tests the rotation functionality of a 2D model.
set Dimension = 2
set Use years in output instead of seconds = true
@@ -41,11 +41,10 @@ end
subsection Boundary velocity model
subsection GPlates model
- set Data directory = $ASPECT_SOURCE_DIR/data/velocity-boundary-conditions/gplates/
set Velocity file name = current_day.gpml
set Time step = 1e6
- set Point one = 1.5708,4.87
- set Point two = 1.5708,5.24
+ set Point one = 1.0,1.87
+ set Point two = 2.0,5.24
set Interpolation width = 2000000
end
end
diff --git a/tests/gplates-rotation/screen-output b/tests/gplates-rotation/screen-output
new file mode 100644
index 0000000..9688d23
--- /dev/null
+++ b/tests/gplates-rotation/screen-output
@@ -0,0 +1,44 @@
+
+Number of active cells: 768 (on 4 levels)
+Number of degrees of freedom: 10,656 (6,528+864+3,264)
+
+
+ Setting up GPlates boundary velocity plugin.
+
+ Input point 1 spherical coordinates: 1.000 1.870
+ Input point 1 normalized cartesian coordinates: -0.248 0.804 0.540
+ Input point 1 rotated model coordinates: 0.000 1.000 0.000
+ Input point 2 spherical coordinates: 2.000 5.240
+ Input point 2 normalized cartesian coordinates: 0.458 -0.786 -0.416
+ Input point 2 rotated model coordinates: -0.243 -0.970 0.000
+
+ Model will be rotated by -154.50 degrees around axis -0.06 0.95 0.31
+ The ParaView rotation angles are: 32.70 148.06 17.14
+ The inverse ParaView rotation angles are: 36.45 -152.58 -1.64
+
+ Loading GPlates boundary velocity file ASPECT_DIR/data/velocity-boundary-conditions/gplates/current_day.gpml.
+
+
+ Loading GPlates boundary velocity file ASPECT_DIR/data/velocity-boundary-conditions/gplates/current_day.gpml.
+
+
+ Loading new velocity file did not succeed.
+ Assuming constant boundary conditions for rest of model run.
+
+*** Timestep 0: t=0 years
+ Solving temperature system... 0 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 30+15 iterations.
+
+ Postprocessing:
+ RMS, max velocity: 0.0209 m/year, 0.1 m/year
+ Temperature min/avg/max: 273 K, 1575 K, 2600 K
+ Heat fluxes through boundary parts: -9.361e+05 W, 1.919e+06 W
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
++---------------------------------+-----------+------------+------------+
++---------------------------------+-----------+------------+------------+
+
diff --git a/tests/gplates-1-3/statistics b/tests/gplates-rotation/statistics
similarity index 85%
copy from tests/gplates-1-3/statistics
copy to tests/gplates-rotation/statistics
index 3d59c62..feff876 100644
--- a/tests/gplates-1-3/statistics
+++ b/tests/gplates-rotation/statistics
@@ -16,4 +16,4 @@
# 16: Average nondimensional temperature (K)
# 17: Outward heat flux through boundary with indicator 0 ("inner") (W)
# 18: Outward heat flux through boundary with indicator 1 ("outer") (W)
-0 0.0000e+00 768 7392 3264 0 45 108 64 1.6661e+06 2.56804983e-02 1.01393827e-01 2.73000000e+02 1.57484593e+03 2.60000000e+03 5.59452482e-01 -9.36071065e+05 1.91866167e+06
+0 0.0000e+00 768 7392 3264 0 45 108 64 1.7217e+06 2.09138819e-02 1.00395547e-01 2.73000000e+02 1.57484593e+03 2.60000000e+03 5.59452482e-01 -9.36071065e+05 1.91866167e+06
More information about the CIG-COMMITS
mailing list