TODO

Author: ESeNonFossiIo

include/interfaces/ALE_navier_stokes.h

@@ -265,9 +265,12 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
 // variables:
 // velocity:
       const ResidualType &div_u = div_us[quad];
+      const Tensor<1, dim, ResidualType> &u = us[quad];
       const Tensor<1, dim, ResidualType> &u_dot = us_dot[quad];
       const Tensor<2, dim, ResidualType> &grad_u = grad_us[quad];
       const Tensor <2, dim, ResidualType> &sym_grad_u = sym_grad_us[quad];
+      
+
 // displacement
       const Tensor<1, dim, ResidualType> &d_dot = ds_dot[quad];
       const Tensor<2, dim, ResidualType> &grad_d = grad_ds[quad];
@@ -308,7 +311,7 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
 // pressure:
           auto p_test = fev[pressure].value(i,quad);
           auto grad_p_test = fev[pressure].gradient(i,quad);
-
+          
           // Preconditioner:
           residual[1][i] +=
             (
@@ -321,22 +324,22 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
             )*JxW[quad];
 
           // Navier Stokes:
-          res_u =
+          res_u = J_ale * (
 // time derivative term
-            rho*scalar_product( u_dot  , u_test )
-//
-            + scalar_product( grad_u * ( F_inv * ( u_old - d_dot ) )  , u_test )
+            + rho * scalar_product( 
+              u_dot + grad_u * ( F_inv * ( u_old - d_dot ) ), u_test )
 //
-            + scalar_product(sigma * Ft_inv, grad_u_test )
+            + scalar_product(sigma * Ft_inv, grad_u_test ))
 // divergence free constriant
-            - div_u  * p_test;
+            - div_u * p_test;
+            // - scalar_product(J_ale * F_inv * u, grad_p_test); // TODO: muovere la divergenza a sinistra! se no c'é anche l'integrale sul bordo.
 
           // ALE:
           res_d =
 // Impose armonicity of d
             scalar_product( grad_d , grad_d_test );
 
-          residual[0][i] += (res_u * J_ale + res_d)*JxW[quad];
+          residual[0][i] += (res_u + res_d)*JxW[quad];
         }
     }