From bd79240e200d8c6819262d56f1877768643f2494 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Wed, 9 May 2018 11:31:48 +0200
Subject: [PATCH] correction des approximations de u et k

---
 src/mesh/Mesh.hpp                     |  8 +++---
 src/scheme/FiniteVolumesDiffusion.hpp | 35 +++++++++++++++++----------
 2 files changed, 26 insertions(+), 17 deletions(-)

diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 2a9bb022d..58b7a0520 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -50,7 +50,7 @@ public:
 
   // pas constant
 
-
+  /*
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes())
@@ -60,11 +60,11 @@ public:
   	m_xr[r][0] = r*delta_x;
       });
   }
-
+*/
 
   // pas non constant
 
-  /*    
+      
   Mesh(const Connectivity& connectivity)
   : m_connectivity(connectivity),
     m_xr("xr", connectivity.numberOfNodes())
@@ -85,7 +85,7 @@ public:
   	}
       });
   }
-  */
+  
 
   ~Mesh()
   {
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 6e2215256..760a28350 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -92,32 +92,38 @@ private:
       = m_mesh.connectivity().faceCellLocalFace();
 
     const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
+    const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
 
     Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
         Rdd sum = zero;
 	double sum2 = 0.;
+	double sum3 = 0.;
 	for (int j=0; j<face_nb_cells(l); ++j) {
 	    int cell_here = face_cells(l,j);
 	    int local_face_number_in_cell = face_cell_local_face(l,j);
 	    sum -= tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell));
-	    sum2 += kj(cell_here); 
+	    sum2 += kj(cell_here)*Vj(cell_here);
+	    sum3 += Vj(cell_here);
 	}
 
-	m_Fl(l) = ((sum2/face_nb_cells(l))/Vl(l))*sum;
-
+	//m_Fl(l) = ((sum2/face_nb_cells(l))/Vl(l))*sum;
+	m_Fl(l) = ((sum2/sum3)/Vl(l))*sum;
 	
 	// Conditions aux bords 
 	
+	/*
+	// Essai 1
 	int cell_here = face_cells(0,0);
 	int local_face_number_in_cell = face_cell_local_face(0,0);
-	//m_Fl(0) = -(kL(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
-	m_Fl(0) = -(kL(0) + kj(0))*0.5*(tensorProduct(uj(0), Cjr(0, 0)));
-
+	m_Fl(0) = -(kL(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
 	cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
 	local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
-	//m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+	m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+
+	// Essai 2
+	m_Fl(0) = -(kL(0) + kj(0))*0.5*(tensorProduct(uj(0), Cjr(0, 0)));
 	m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(m_mesh.numberOfCells()-1))*0.5*(tensorProduct(uj(m_mesh.numberOfCells()-1), Cjr(m_mesh.numberOfCells()-1,1))); 
-	
+	*/
 	
       });
 
@@ -137,21 +143,24 @@ private:
     const Kokkos::View<const unsigned short**> face_cell_local_face
       = m_mesh.connectivity().faceCellLocalFace();
 
+    const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+
     Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
         Rd sum = zero;
+	double sum2 = 0.;
 	for (int j=0; j<face_nb_cells(l); ++j) {
 	  int cell_here = face_cells(l,j);
-	  sum += uj(cell_here);
+	  sum += Vj(cell_here)*uj(cell_here);
+	  sum2 += Vj(cell_here);
 	}                
 	 
-	m_Gl(l) = (1./face_nb_cells(l))*Fl(l)*sum;
+	//m_Gl(l) = (1./face_nb_cells(l))*Fl(l)*sum;
+	m_Gl(l) = (1./sum2)*Fl(l)*sum;
       });
 
-    // C'est juste un test, mais ca devrait marcher
+    // Provisoire
     m_Gl(0) = 0;
     m_Gl(m_mesh.numberOfFaces()-1) = 0;
-    // m_Gl(0) = 0.5*m_Gl(0);
-    // m_Gl(m_mesh.numberOfFaces()-1) = 0.5*m_Gl(m_mesh.numberOfFaces()-1);
 
     return m_Gl ;
   }
-- 
GitLab