Fixes to assemblers and test.

Author: MFraters

source/simulator/assemblers/newton_stokes.cc

@@ -464,31 +464,13 @@ namespace aspect
               if (enable_prescribed_dilation)
                 {
                   const unsigned int index_direction=fe.system_to_component_index(i).first;
-                  // Only want the velocity components and not the pressure one (which is the last one), so add 1
-                  if (introspection.is_stokes_component(index_direction+1))
-                    data.local_rhs(i) += (
-                                           // RHS of - (div u,q) = - (R,q)
-                                           - pressure_scaling
-                                           * (prescribed_dilation->dilation_rhs_term[index_direction][q] -
-                                              prescribed_dilation->dilation_lhs_term[q] *
-                                              scratch.material_model_inputs.pressure[q])
-                                           * scratch.phi_p[i]
-                                         ) * JxW;
-                }
-
-              // Only assemble this term if we are running incompressible, otherwise this term
-              // is already included on the LHS of the equation.
-              if (enable_prescribed_dilation && !material_model_is_compressible)
-                {
-                  const unsigned int index_direction=fe.system_to_component_index(i).first;
-                  // Only want the velocity components and not the pressure one (which is the last one), so add 1
-                  if (introspection.is_stokes_component(index_direction+1))
-                    data.local_rhs(i) += (
-                                           // RHS of momentum eqn: - \int 2/3 eta R, div v
-                                           - 2.0 / 3.0 * eta
-                                           * prescribed_dilation->dilation_rhs_term[index_direction][q]
-                                           * scratch.div_phi_u[i]
-                                         ) * JxW;
+                  data.local_rhs(i) += (
+                                         - pressure_scaling
+                                         * (prescribed_dilation->dilation_rhs_term[index_direction][q] -
+                                            prescribed_dilation->dilation_lhs_term[q] *
+                                            scratch.material_model_inputs.pressure[q])
+                                         * scratch.phi_p[i]
+                                       ) * JxW;
                 }
             }
 

source/simulator/assemblers/stokes.cc

@@ -456,12 +456,16 @@ namespace aspect
                   const unsigned int index_direction=fe.system_to_component_index(i).first;
                   // Only want the velocity components and not the pressure one (which is the last one), so add 1
                   if (introspection.is_stokes_component(index_direction+1))
-                    data.local_rhs(i) += (
-                                           // RHS of - (div u,q) = - (R,q)
-                                           - pressure_scaling
-                                           * prescribed_dilation->dilation_rhs_term[index_direction][q]
-                                           * scratch.phi_p[i]
-                                         ) * JxW;
+                    {
+                      const unsigned int index_direction=fe.system_to_component_index(i).first;
+                      data.local_rhs(i) += (
+                                             - pressure_scaling
+                                             * (prescribed_dilation->dilation_rhs_term[index_direction][q] -
+                                                prescribed_dilation->dilation_lhs_term[q] *
+                                                scratch.material_model_inputs.pressure[q])
+                                             * scratch.phi_p[i]
+                                           ) * JxW;
+                    }
                 }
 
               // Only assemble this term if we are running incompressible, otherwise this term

tests/prescribed_dilation.cc

@@ -207,11 +207,11 @@ namespace aspect
                 }
               if (prescribed_dilation)
                 {
-                  prescribed_dilation->dilation_rhs_term[i] = x;
-                  prescribed_dilation->dilation_lhs_term[0][i] = 0.;
-                  prescribed_dilation->dilation_lhs_term[1][i] = 0.;
+                  prescribed_dilation->dilation_lhs_term[i] = 0;
+                  prescribed_dilation->dilation_rhs_term[0][i] = x;
+                  prescribed_dilation->dilation_rhs_term[1][i] = x;
                   if (dim == 3)
-                    prescribed_dilation->dilation_lhs_term[2][i] = 0.;
+                    prescribed_dilation->dilation_rhs_term[2][i] = x;
                 }
 
             }