[WIP] Pbl?

.travis.yml

@@ -55,7 +55,7 @@ install:
 
     # hdf5
     - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then sudo $apt_get_install libhdf5-7; fi
-    - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then sudo $apt_get_install libpango-1.0-0 libpangocairo-1.0-0 libhdf5-dev; fi
+    - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then sudo $apt_get_install  -o Dpkg::Options::="--force-confdef" -o Dpkg::Options::="--force-confold" libpango-1.0-0 libpangocairo-1.0-0 libhdf5-dev; fi
     - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then sudo $apt_get_install libhdf5-serial-dev hdf5-tools; fi
     - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then brew install hdf5 --with-cxx; fi
 

libmpdata++/output/gnuplot.hpp

@@ -138,7 +138,7 @@ namespace libmpdataxx
 
       // helper constructs to make it compilable for both 1D and 2D versions
       std::string binfmt(blitz::Array<typename parent_t::real_t, 1>) { throw std::logic_error("binfmt() only for 2D!"); }
-      std::string binfmt(blitz::Array<typename parent_t::real_t, 2> a) { return gp->binfmt(a) + " scan=yx "; }
+      std::string binfmt(blitz::Array<typename parent_t::real_t, 2> a) { return gp->binfmt(a.transpose(blitz::secondDim, blitz::firstDim)) + " scan=yx "; }
 
       void record(const int var)
       {

libmpdata++/solvers/detail/mpdata_rhs_vip_prs_2d_common.hpp

@@ -60,11 +60,17 @@ namespace libmpdataxx
             lap_tmp1(this->ijk) /= (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
             lap_tmp2(this->ijk) /= (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
           }
+          if (this->mem->G)
+          {
+            lap_tmp1(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp2(this->ijk) *= (*this->mem->G)(this->ijk);
+          }
           this->set_edges(lap_tmp1, lap_tmp2, i, j, 0);
           this->xchng_pres(lap_tmp1, i, j);
           this->xchng_pres(lap_tmp2, i, j);
           ,
           formulae::nabla::div(lap_tmp1, lap_tmp2, i, j, dx, dy)
+          / formulae::G<ct_params_t::opts BOOST_PP_COMMA() 0>(*this->mem->G, i, j) 
         );
 
         auto err_init(
@@ -84,11 +90,17 @@ namespace libmpdataxx
             lap_tmp1(this->ijk) /= (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
             lap_tmp2(this->ijk) /= (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
           }
+          if (this->mem->G)
+          {
+            lap_tmp1(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp2(this->ijk) *= (*this->mem->G)(this->ijk);
+          }
           this->set_edges(lap_tmp1, lap_tmp2, i, j, -1);
           this->xchng_pres(lap_tmp1, i, j);
           this->xchng_pres(lap_tmp2, i, j);
           ,
           formulae::nabla::div(lap_tmp1, lap_tmp2, i, j, dx, dy)
+          / formulae::G<ct_params_t::opts BOOST_PP_COMMA() 0>(*this->mem->G, i, j) 
         );
 
 	void ini_pressure()

libmpdata++/solvers/detail/mpdata_rhs_vip_prs_3d_common.hpp

@@ -58,12 +58,19 @@ namespace libmpdataxx
           lap_tmp1(this->ijk) = formulae::nabla::grad<0>(arr, i, j, k, dx);
           lap_tmp2(this->ijk) = formulae::nabla::grad<1>(arr, j, k, i, dy);
           lap_tmp3(this->ijk) = formulae::nabla::grad<2>(arr, k, i, j, dz);
+          if (this->mem->G)
+          {
+            lap_tmp1(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp2(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp3(this->ijk) *= (*this->mem->G)(this->ijk);
+          }
           this->set_edges(lap_tmp1, lap_tmp2, lap_tmp3, i, j, k, 0);
           this->xchng_pres(lap_tmp1, i, j, k);
           this->xchng_pres(lap_tmp2, i, j, k);
           this->xchng_pres(lap_tmp3, i, j, k);
           ,
           formulae::nabla::div(lap_tmp1, lap_tmp2, lap_tmp3, i, j, k, dx, dy, dz)
+          / formulae::G<ct_params_t::opts BOOST_PP_COMMA() 0>(*this->mem->G, i, j, k)
         );
 
         auto err_init(
@@ -82,12 +89,19 @@ namespace libmpdataxx
           lap_tmp1(this->ijk) = formulae::nabla::grad<0>(arr, i, j, k, dx) - v1(this->ijk);
           lap_tmp2(this->ijk) = formulae::nabla::grad<1>(arr, j, k, i, dy) - v2(this->ijk);
           lap_tmp3(this->ijk) = formulae::nabla::grad<2>(arr, k, i, j, dz) - v3(this->ijk);
+          if (this->mem->G)
+          {
+            lap_tmp1(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp2(this->ijk) *= (*this->mem->G)(this->ijk);
+            lap_tmp3(this->ijk) *= (*this->mem->G)(this->ijk);
+          }
           this->set_edges(lap_tmp1, lap_tmp2, lap_tmp3, i, j, k, -1);
           this->xchng_pres(lap_tmp1, i, j, k);
           this->xchng_pres(lap_tmp2, i, j, k);
           this->xchng_pres(lap_tmp3, i, j, k);
           ,
           formulae::nabla::div(lap_tmp1, lap_tmp2, lap_tmp3, i, j, k, dx, dy, dz)
+          / formulae::G<ct_params_t::opts BOOST_PP_COMMA() 0>(*this->mem->G, i, j, k)
         );
 
 	void ini_pressure()

libmpdata++/solvers/mpdata_rhs_vip.hpp

@@ -120,7 +120,10 @@ namespace libmpdataxx
 	} 
 	else
 	{ 
-	  assert(false); // TODO (perhaps better change definition of stash?)
+	  this->mem->GC[d](pi<d>(i+h,j)) = this->dt / di * .5 * (
+	    (*this->mem->G)(pi<d>(i,    j)) * psi(pi<d>(i,    j)) + 
+	    (*this->mem->G)(pi<d>(i + 1,j)) * psi(pi<d>(i + 1,j))
+	  );
 	}
       }  
 
@@ -227,7 +230,10 @@ namespace libmpdataxx
 	} 
 	else
 	{ 
-	  assert(false); // TODO (perhaps better change definition of stash?)
+	  this->mem->GC[d](pi<d>(i+h, j, k)) = this->dt / di * .5 * (
+	    (*this->mem->G)(pi<d>(i,     j, k)) * psi(pi<d>(i,     j, k)) +
+	    (*this->mem->G)(pi<d>(i + 1, j, k)) * psi(pi<d>(i + 1, j, k))
+	  );
 	}
       }