[cig-commits] [commit] master: Now always remove the radial component of velocity from the GPlates plugin. (01d4773)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/aspect

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/e2d2428a203331b4832dbc4e5f4d716934c22694...f68c1b584dc0211691228ef830ea1d2538da258d

>---------------------------------------------------------------

commit 01d4773941dcf119e62e30b9379be02a51ed4793
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date:   Fri May 23 18:47:00 2014 -0500

    Now always remove the radial component of velocity from the GPlates plugin.
    
    Before the removal we check that this component is small, to ensure the rotation and projection is working properly. If it is small and we do not want to have it in the solution we might as well remove it completely.


>---------------------------------------------------------------

01d4773941dcf119e62e30b9379be02a51ed4793
 source/velocity_boundary_conditions/gplates.cc | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/source/velocity_boundary_conditions/gplates.cc b/source/velocity_boundary_conditions/gplates.cc
index ca190f9..02b2f35 100644
--- a/source/velocity_boundary_conditions/gplates.cc
+++ b/source/velocity_boundary_conditions/gplates.cc
@@ -264,6 +264,12 @@ namespace aspect
                                         ExcMessage("Error in velocity boundary module interpolation. "
                                                    "Radial component of velocity should be zero, but is not."));
 
+                            // After checking that the residual normal velocity is small
+                            // we might as well remove it to increase compatibility of the
+                            // pressure right hand side.
+                            const Tensor<1,3> normalized_position = position / position.norm();
+                            surf_vel -=  (surf_vel * normalized_position) * normalized_position;
+
                             return surf_vel / n_interpolation_weight;
                           }
                         else
@@ -292,6 +298,12 @@ namespace aspect
                      ExcMessage("Error in velocity boundary module interpolation. "
                                 "Radial component of velocity should be zero, but is not."));
 
+        // After checking that the residual normal velocity is small
+        // we might as well remove it to increase compatibility of the
+        // pressure right hand side.
+        const Tensor<1,3> normalized_position = position / position.norm();
+        surf_vel -=  (surf_vel * normalized_position) * normalized_position;
+
         // We have interpolated over the whole dataset of the provided velocities, which is ok.
         // Velocity will be constant for all evaluation points in this case.
         return surf_vel / n_interpolation_weight;