From 7715b83e98b0b4df267d44333b9e55f7c0b19949 Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Thu, 5 Dec 2013 11:41:25 +0000
Subject: [PATCH] fix tangentil jump

git-svn-id: https://svn.dealii.org/trunk@31899 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/integrators/laplace.h | 34 +++++++++++++------
 1 file changed, 23 insertions(+), 11 deletions(-)

diff --git a/deal.II/include/deal.II/integrators/laplace.h b/deal.II/include/deal.II/integrators/laplace.h
index 61a58afbb8..0148bbcaa5 100644
--- a/deal.II/include/deal.II/integrators/laplace.h
+++ b/deal.II/include/deal.II/integrators/laplace.h
@@ -405,28 +405,40 @@ namespace LocalIntegrators
 		  double v1dotn = 0.;
 		  double u2dotn = 0.;
 		  double v2dotn = 0.;
-		  
+
+		  double ngradu1n = 0.;
+		  double ngradv1n = 0.;
+		  double ngradu2n = 0.;
+		  double ngradv2n = 0.;
+
 		  for (unsigned int d=0; d<dim; ++d)
 		    {
 		      u1dotn += n(d)*fe1.shape_value_component(j,k,d);
 		      v1dotn += n(d)*fe1.shape_value_component(i,k,d);
 		      u2dotn += n(d)*fe2.shape_value_component(j,k,d);
 		      v2dotn += n(d)*fe2.shape_value_component(i,k,d);
+
+		      ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
+		      ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
+		      ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
+		      ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
 		    }
 			  
 		  for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
 		    {
+                      const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n(d);
+		      const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
 		      
-						       // multiply by 
-                      const double vi = fe1.shape_value_component(i,k,d)*(1-v1dotn);
-                      const double dnvi = n * fe1.shape_grad_component(i,k,d)*(1-v1dotn);
-                      const double ve = fe2.shape_value_component(i,k,d)*(1-v2dotn);
-                      const double dnve = n * fe2.shape_grad_component(i,k,d)*(1-v2dotn);
-                      const double ui = fe1.shape_value_component(j,k,d)*(1-u1dotn);
-                      const double dnui = n * fe1.shape_grad_component(j,k,d)*(1-u1dotn);
-                      const double ue = fe2.shape_value_component(j,k,d)*(1-u2dotn);
-                      const double dnue = n * fe2.shape_grad_component(j,k,d)*(1-u2dotn);
-                      M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
+                      const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n(d);
+                      const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n(d);
+		      
+                      const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n(d);
+                      const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n(d);
+		      
+                      const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n(d);
+                      const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n(d);
+
+		      M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
                       M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
                       M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
                       M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
-- 
2.39.5