[cig-commits] commit 2289 by gassmoeller to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 2289

Updated the gplates plugin to use the official gplates 1.3 format. Also fixed a small bug with prescribing the velocities in the first timestep.

U   trunk/aspect/source/velocity_boundary_conditions/gplates.cc


http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2289&peg=2289

Diff:
Modified: trunk/aspect/source/velocity_boundary_conditions/gplates.cc
===================================================================
--- trunk/aspect/source/velocity_boundary_conditions/gplates.cc	2014-02-03 11:36:29 UTC (rev 2288)
+++ trunk/aspect/source/velocity_boundary_conditions/gplates.cc	2014-02-03 16:02:16 UTC (rev 2289)
@@ -104,9 +104,9 @@
         // 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 unsigned int n_points = pt.get_child("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:domainSet.gml:MultiPoint").size();
+        const unsigned int n_theta = static_cast<unsigned int>(0.5 + std::sqrt(0.25 + n_points/2));
+        const unsigned int n_phi = 2 * (n_theta - 1);
 
         if ((delta_theta != 0.0) || (delta_phi != 0.0))
           {
@@ -140,8 +140,8 @@
             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_phi][i%n_phi][0] = vtheta/cmyr_si;
+            (*velocity_values)[i/n_phi][i%n_phi][1] = vphi/cmyr_si;
 
             i++;
           }
@@ -434,7 +434,7 @@
     GPlates<dim>::
     boundary_velocity (const Point<dim> &position) const
     {
-      if (current_time > 0.0)
+      if (current_time >= 0.0)
         return lookup->surface_velocity(position,time_weight);
       else
         return Tensor<1,dim> ();