@@ -591,7 +591,6 @@ class SmokeSolver {
591591 bool isDirichletPressure = (flag == 4 );
592592 auto vNgbr = Vec3s::zero (); // dirichletVelocity->tree().getValue(neighbor);
593593
594- // TODO: Double check this:
595594 if (isNeumannPressure) {
596595 double delta = 0.0 ;
597596 // Neumann pressure from bbox
@@ -620,7 +619,7 @@ class SmokeSolver {
620619 } else if (isDirichletPressure) {
621620 diagonal -= 1.0 ;
622621 source -= dirichletBC;
623- #if 0 // supposedly the same as the two lines above--checked on Friday.
622+ #if 0 // supposedly the same as the two lines above
624623 // Dirichlet pressure
625624 if (neighbor.x() + 1 == ijk.x() /* left x-face */) {
626625 diagonal -= 1.0;
@@ -704,24 +703,6 @@ class SmokeSolver {
704703 mPressure ->setName (" pressure" );
705704 }
706705
707- void render ()
708- {
709- if (mVERBOSE ) printRelevantVelocity (" velocity init" );
710-
711- float divBefore = computeDivergence (mDivBefore , mVCurr , " before" );
712- if (mVERBOSE ) printGrid (*mDivBefore );
713-
714- // Make the velocity divergence free by solving Poisson Equation and subtracting the pressure gradient
715- pressureProjection ();
716- subtractPressureGradFromVel ();
717-
718- float divAfter = computeDivergence (mDivAfter , mVCurr , " after" );
719- if (mVERBOSE ) printGrid (*mDivAfter );
720-
721- writeVDBsDebug (1 );
722- }
723-
724-
725706 template <class GridType >
726707 typename GridType::Ptr
727708 initGridBgAndName (typename GridType::ValueType background, std::string name)
@@ -732,44 +713,6 @@ class SmokeSolver {
732713 return grid;
733714 }
734715
735- template <class GridType >
736- void printGrid (const GridType& grid, std::string nameFromUser = " " ) {
737- using ValueType = typename GridType::ValueType;
738- auto name = nameFromUser != " " ? nameFromUser : grid.getName ();
739- std::cout << " printGrid::Printing grid " << name << std::endl;
740- auto acc = grid.getAccessor ();
741- for (auto iter = grid.beginValueOn (); iter; ++iter) {
742- math::Coord ijk = iter.getCoord ();
743- std::cout << " val" << ijk << " = " << acc.getValue (ijk) << std::endl;
744- }
745- std::cout << std::endl;
746- }
747-
748- void printRelevantVelocity (std::string nameFromUser = " " ) {
749- std::cout << " printRelevantVelocity::printing " << nameFromUser << std::endl;
750- auto flagsAcc = mFlags ->getAccessor ();
751- auto vAcc = mVCurr ->getAccessor ();
752- for (auto iter = mFlags ->beginValueOn (); iter; ++iter) {
753- math::Coord ijk = iter.getCoord ();
754- math::Coord im1jk = ijk.offsetBy (-1 , 0 , 0 );
755- math::Coord ijm1k = ijk.offsetBy (0 , -1 , 0 );
756- math::Coord ijkm1 = ijk.offsetBy (0 , 0 , -1 );
757-
758- int flag = flagsAcc.getValue (ijk);
759- int flagim1jk = flagsAcc.getValue (im1jk);
760- int flagijm1k = flagsAcc.getValue (ijm1k);
761- int flagijkm1 = flagsAcc.getValue (ijkm1);
762-
763- if (flag == 1 ) {
764- std::cout << " vel" << ijk << " = " << vAcc.getValue (ijk) << std::endl;
765- } else {
766- if (flagim1jk == 1 || flagijm1k == 1 || flagijkm1 == 1 ) {
767- std::cout << " vel" << ijk << " = " << vAcc.getValue (ijk) << std::endl;
768- }
769- }
770- }
771- }
772-
773716 void initFlags ()
774717 {
775718 mFlags = initGridBgAndName<Int32Grid>(0 , " flags" );
@@ -814,7 +757,6 @@ class SmokeSolver {
814757 divGrid = tools::divergence (*vecGrid);
815758 divGrid->tree ().topologyIntersection (mInteriorPressure ->tree ());
816759 float div = computeLInfinity (*divGrid);
817- std::cout << " Divergence " << suffix.c_str () << " = " << div << std::endl;
818760 return div;
819761 }
820762
@@ -855,24 +797,6 @@ class SmokeSolver {
855797 }
856798 }
857799
858- void writeVDBsDebug (int const frame) {
859- std::ostringstream ss;
860- ss << " INIT_DEBUG" << std::setw (3 ) << std::setfill (' 0' ) << frame << " .vdb" ;
861- std::string fileName (ss.str ());
862- io::File file (fileName.c_str ());
863-
864- openvdb::GridPtrVec grids;
865- grids.push_back (mFlags );
866- grids.push_back (mInteriorPressure );
867- grids.push_back (mVCurr );
868- grids.push_back (mDivBefore );
869- grids.push_back (mDivAfter );
870- grids.push_back (mPressure );
871-
872- file.write (grids);
873- file.close ();
874- }
875-
876800 bool mVERBOSE = false ;
877801
878802 float mVoxelSize = 0 .1f ;
@@ -900,6 +824,7 @@ TEST_F(TestPoissonSolver, testRemoveDivergence)
900824 SmokeSolver smoke (0 .1f );
901825
902826 float divBefore = smoke.computeDivergence (smoke.mDivBefore , smoke.mVCurr , " before" );
827+ EXPECT_NEAR (divBefore, -39.425 , 1 .e -4f );
903828
904829 // Make the velocity divergence free by solving Poisson Equation and subtracting the pressure gradient
905830 smoke.pressureProjection ();
@@ -913,5 +838,4 @@ TEST_F(TestPoissonSolver, testRemoveDivergence)
913838
914839 float divAfter = smoke.computeDivergence (smoke.mDivAfter , smoke.mVCurr , " after" );
915840 EXPECT_TRUE (divAfter < 1 .e -3f );
916- smoke.writeVDBsDebug (1 /* frame */ );
917841}
0 commit comments