debugging uniform lattice

8 files changed, 243 insertions, 173 deletions

diff --git a/extern/softbody/src/admmpd_energy.cpp b/extern/softbody/src/admmpd_energy.cpp
index 36f0bfa1478..41aebbe0739 100644
--- a/extern/softbody/src/admmpd_energy.cpp
+++ b/extern/softbody/src/admmpd_energy.cpp
@@ -94,7 +94,7 @@ int EnergyTerm::init_tet(
 	volume = edges.determinant() / 6.0f;
 	if( volume < 0 )
 	{
-		printf("**EnergyTerm::init_tet: Inverted initial tet");
+		printf("**admmpd::EnergyTerm Error: Inverted initial tet: %f\n",volume);
 		return 0;
 	}
 	double k = lame.m_bulk_mod;
diff --git a/extern/softbody/src/admmpd_lattice.cpp b/extern/softbody/src/admmpd_lattice.cpp
index be7a5bc63d2..b0c8f6ca556 100644
--- a/extern/softbody/src/admmpd_lattice.cpp
+++ b/extern/softbody/src/admmpd_lattice.cpp
@@ -4,6 +4,7 @@
 #include "admmpd_lattice.h"
 #include "admmpd_math.h"
 #include <iostream>
+#include <unordered_map>
 
 //#include "vdb.h"
 
@@ -52,94 +53,114 @@ inline void map_to_x4i(const std::vector<Vector4i> &x_vec, Eigen::MatrixXi *x)
 
 bool Lattice::generate(
 	const Eigen::MatrixXd &V,
-    Eigen::MatrixXd *x, // lattice vertices, n x 3
-    Eigen::MatrixXi *tets) // lattice elements, m x 4
+    Eigen::MatrixXd *X, // lattice vertices, n x 3
+    Eigen::MatrixXi *T) // lattice elements, m x 4
 {
-	// All vertices enclosed in 5 tets!
-	// Will generate actual lattice in the future.
-	AlignedBox<double,3> box;
-	vtx = V;
-	int nv = vtx.rows();
-	for (int i=0; i<nv; ++i)
+	AlignedBox<double,3> aabb;
+	int nv = V.rows();
+	for(int i=0; i<nv; ++i)
+		aabb.extend(V.row(i).transpose());
+	
+	aabb.extend(aabb.min()-Vector3d::Ones()*1e-6);
+	aabb.extend(aabb.max()+Vector3d::Ones()*1e-6);
+	Vector3d mesh_boxmin = aabb.min();
+	Vector3d sizes = aabb.sizes();
+	double grid_dx = sizes.maxCoeff() * 0.2;
+
+	// Generate vertices and tets
+	std::vector<Vector3d> verts;
+	std::vector<Vector4i> tets;
 	{
-		Vector3d v = vtx.row(i);
-		box.extend(v);
-	}
-	box.extend(box.min() - Vector3d::Ones() * 1e-12);
-	box.extend(box.max() + Vector3d::Ones() * 1e-12);
+		std::unordered_map<std::string, int> vertex_map; // [i,j,k]->index
 
-	std::vector<Vector3d> x_vec;
-	std::vector<Vector4i> t_vec;
-	create_packed_tets(box.min(),box.max(),x_vec,t_vec);
+		auto grid_hash = [&](const Vector3d &x)
+		{
+			Vector3i ll = Vector3i( // nearest gridcell
+				std::round((x[0]-mesh_boxmin[0])/grid_dx),
+				std::round((x[1]-mesh_boxmin[1])/grid_dx),
+				std::round((x[2]-mesh_boxmin[2])/grid_dx));
+			return std::to_string(ll[0])+' '+std::to_string(ll[1])+' '+std::to_string(ll[2]);
+		};
+
+		auto add_box = [&](int i_, int j_, int k_)
+		{
+			Vector3d min = mesh_boxmin + Vector3d(i_*grid_dx, j_*grid_dx, k_*grid_dx);
+			Vector3d max = mesh_boxmin + Vector3d((i_+1)*grid_dx, (j_+1)*grid_dx, (k_+1)*grid_dx);
+			std::vector<Vector3d> v = {
+				// Top plane, clockwise looking down
+				max,
+				Vector3d(min[0], max[1], max[2]),
+				Vector3d(min[0], max[1], min[2]),
+				Vector3d(max[0], max[1], min[2]),
+				// Bottom plan, clockwise looking down
+				Vector3d(max[0], min[1], max[2]),
+				Vector3d(min[0], min[1], max[2]),
+				min,
+				Vector3d(max[0], min[1], min[2])
+			};
+			// Add vertices and get indices of the box
+			std::vector<int> b;
+			for(int i=0; i<8; ++i)
+			{
+				std::string hash = grid_hash(v[i]);
+				if( vertex_map.count(hash)==0 )
+				{
+					vertex_map[hash] = verts.size();
+					verts.emplace_back(v[i]);
+				}
+				b.emplace_back(vertex_map[hash]);
+			}
+			// From the box, create five new tets
+			std::vector<Vector4i> new_tets = {
+				Vector4i( b[0], b[5], b[7], b[4] ),
+				Vector4i( b[5], b[7], b[2], b[0] ),
+				Vector4i( b[5], b[0], b[2], b[1] ),
+				Vector4i( b[7], b[2], b[0], b[3] ),
+				Vector4i( b[5], b[2], b[7], b[6] )
+			};
+			for(int i=0; i<5; ++i)
+				tets.emplace_back(new_tets[i]);
+		};
+
+		int ni = std::ceil(sizes[0]/grid_dx);
+		int nj = std::ceil(sizes[1]/grid_dx);
+		int nk = std::ceil(sizes[2]/grid_dx);
+		for(int i=0; i<ni; ++i)
+		{
+			for(int j=0; j<nj; ++j)
+			{
+				for(int k=0; k<nk; ++k)
+				{
+					add_box(i,j,k);
+				}
+			}
+		}
 
-	// Copy vector data to output
-	map_to_x3d(x_vec, x);
-	map_to_x4i(t_vec, tets);
-	return compute_vtx_tet_mapping(&vtx, &vtx_to_tet, &barys, x, tets);
+	} // end create boxes and vertices
 
-} // end gen lattice
+	// Copy data into matrices
+	{
+		nv = verts.size();
+		X->resize(nv,3);
+		for(int i=0; i<nv; ++i){
+			for(int j=0; j<3; ++j){
+				X->operator()(i,j) = verts[i][j];
+			}
+		}
+		int nt = tets.size();
+		T->resize(nt,4);
+		for(int i=0; i<nt; ++i){
+			for(int j=0; j<4; ++j){
+				T->operator()(i,j) = tets[i][j];
+			}
+		}
+	}
 
-//
-// Original source (BSD-2):
-// github.com/mattoverby/mclscene/blob/master/include/MCL/EmbeddedMesh.hpp
-//
-void Lattice::create_packed_tets(
-	const Eigen::Vector3d &min,
-	const Eigen::Vector3d &max,
-	std::vector<Vector3d> &verts,
-	std::vector<Vector4i> &tets )
-{
+	return compute_vtx_tet_mapping(
+		&V, &vtx_to_tet, &barys,
+		X, T);
 
-	// Top plane, clockwise looking down
-	Vector3d a = max;
-	Vector3d b( min[0], max[1], max[2] );
-	Vector3d c( min[0], max[1], min[2] );
-	Vector3d d( max[0], max[1], min[2] );
-
-	// Bottom plan, clockwise looking down
-	Vector3d e( max[0], min[1], max[2] );
-	Vector3d f( min[0], min[1], max[2] );
-	Vector3d g( min[0], min[1], min[2] );
-	Vector3d h( max[0], min[1], min[2] );
-
-	// Add the verts
-	int nv = verts.size();
-	verts.emplace_back( a ); // 0
-	verts.emplace_back( b ); // 1
-	verts.emplace_back( c ); // 2
-	verts.emplace_back( d ); // 3
-	verts.emplace_back( e ); // 4
-	verts.emplace_back( f ); // 5
-	verts.emplace_back( g ); // 6
-	verts.emplace_back( h ); // 7
-
-	// Pack 5 tets into the cube
-	Vector4i t0( 0, 5, 7, 4 );
-	Vector4i t1( 5, 7, 2, 0 );
-	Vector4i t2( 5, 0, 2, 1 );
-	Vector4i t3( 7, 2, 0, 3 );
-	Vector4i t4( 5, 2, 7, 6 );
-	Vector4i offset(nv,nv,nv,nv);
-
-	// Add the tets
-	tets.emplace_back( t0+offset );
-	tets.emplace_back( t1+offset );
-	tets.emplace_back( t2+offset );
-	tets.emplace_back( t3+offset );
-	tets.emplace_back( t4+offset );
-
-//vdb_point(a[0],a[1],a[2]);
-//vdb_point(b[0],b[1],b[2]);
-//vdb_point(c[0],c[1],c[2]);
-//vdb_point(d[0],d[1],d[2]);
-//vdb_point(e[0],e[1],e[2]);
-//vdb_point(f[0],f[1],f[2]);
-//vdb_point(g[0],g[1],g[2]);
-//vdb_point(h[0],h[1],h[2]);
-//vdb_line(float x0, float y0, float z0, 
-//             float x1, float y1, float z1);
-
-} // end create packed tets
+} // end gen lattice
 
 bool Lattice::compute_vtx_tet_mapping(
 	const Eigen::MatrixXd *vtx_, // embedded vertices, p x 3
@@ -204,27 +225,5 @@ Eigen::Vector3d Lattice::get_mapped_vertex(
 		x_or_v->row(tet[2]) * b[2] +
 		x_or_v->row(tet[3]) * b[3]);
 }
-/*
-void Lattice::map_to_object(
-	const Eigen::MatrixXd *x,
-	const Eigen::MatrixXi *tets,
-	float (*vertexCos)[3])
-{
-  int nv = vtx.rows();
-  for (int i=0; i<nv; ++i)
-  {
-    int t_idx = vtx_to_tet[i];
-    RowVector4i tet = tets->row(t_idx);
-    RowVector4d b = barys.row(i);
-    vtx.row(i) =
-		x->row(tet[0]) * b[0] +
-		x->row(tet[1]) * b[1] +
-		x->row(tet[2]) * b[2] +
-		x->row(tet[3]) * b[3];
-	vertexCos[i][0] = vtx(i,0);
-	vertexCos[i][1] = vtx(i,1);
-	vertexCos[i][2] = vtx(i,2);
-  }
-} // end map to object
-*/
+
 } // namespace admmpd
\ No newline at end of file
diff --git a/extern/softbody/src/admmpd_lattice.h b/extern/softbody/src/admmpd_lattice.h
index b341de636ec..15e2d391d11 100644
--- a/extern/softbody/src/admmpd_lattice.h
+++ b/extern/softbody/src/admmpd_lattice.h
@@ -11,8 +11,6 @@ namespace admmpd {
 
 class Lattice {
 public:
-
-    Eigen::MatrixXd vtx; // rest embedded vertices, p x 3
     Eigen::VectorXi vtx_to_tet; // what tet vtx is embedded in, p x 1
     Eigen::MatrixXd barys; // barycoords of the embedding, p x 4
 
@@ -22,17 +20,15 @@ public:
         Eigen::MatrixXd *x, // lattice vertices, n x 3
         Eigen::MatrixXi *tets); // lattice elements, m x 4
 
-    // Given boxmin and boxmax, adds
-    // 5 tets (and verts) to the vectors
-    void create_packed_tets(
-        const Eigen::Vector3d &min,
-        const Eigen::Vector3d &max,
-        std::vector<Eigen::Vector3d> &verts,
-        std::vector<Eigen::Vector4i> &tets );
+    // Returns the vtx mapped from x/v and tets
+    Eigen::Vector3d get_mapped_vertex(
+        int idx,
+        const Eigen::MatrixXd *x_or_v,
+        const Eigen::MatrixXi *tets );
+
+protected:
 
     // Returns true on success
-    // Works on ptrs intead of local vars in case
-    // I want to use it for something else.
     bool compute_vtx_tet_mapping(
         const Eigen::MatrixXd *vtx_, // embedded vertices, p x 3
         Eigen::VectorXi *vtx_to_tet_, // what tet vtx is embedded in, p x 1
@@ -40,16 +36,6 @@ public:
         const Eigen::MatrixXd *x_, // lattice vertices, n x 3
         const Eigen::MatrixXi *tets_); // lattice elements, m x 4
 
-    // Returns the vtx mapped from x/v and tets
-    Eigen::Vector3d get_mapped_vertex(
-        int idx,
-        const Eigen::MatrixXd *x_or_v,
-        const Eigen::MatrixXi *tets );
-//    void map_to_object(
-//        const Eigen::MatrixXd *x,
-//        const Eigen::MatrixXi *tets,
-//        float (*vertexCos)[3]);
-
 }; // class Lattice
 
 } // namespace admmpd
diff --git a/extern/softbody/src/admmpd_solver.cpp b/extern/softbody/src/admmpd_solver.cpp
index 4756cb76fac..298ba7ff173 100644
--- a/extern/softbody/src/admmpd_solver.cpp
+++ b/extern/softbody/src/admmpd_solver.cpp
@@ -32,9 +32,16 @@ bool Solver::init(
 	if (!data || !options)
 		throw std::runtime_error("init: data/options null");
 
+	if (V.rows()==0)
+		throw std::runtime_error("init: no input vertices");
+
 	data->x = V;
+	data->v.resize(V.rows(),3);
+	data->v.setZero();
 	data->tets = T;
-	compute_matrices(options,data);
+	if (!compute_matrices(options,data))
+		return false;
+
 	return true;
 } // end init
 
@@ -42,6 +49,8 @@ int Solver::solve(
 	const Options *options,
 	Data *data)
 {
+	if ((int)data->A.nonZeros()==0)
+		return 0;
 
 	// Init the solve which computes
 	// quantaties like M_xbar and makes sure
@@ -276,12 +285,15 @@ void Solver::solve_conjugate_gradients(
 
 } // end solve conjugate gradients
 
-void Solver::compute_matrices(
+bool Solver::compute_matrices(
 	const Options *options,
 	Data *data)
 {
 	// Allocate per-vertex data
 	int nx = data->x.rows();
+	if (nx==0)
+		return false;
+
 	data->x_start = data->x;
 	data->M_xbar.resize(nx,3);
 	data->M_xbar.setZero();
@@ -298,6 +310,11 @@ void Solver::compute_matrices(
 	// Add per-element energies to data
 	std::vector<Triplet<double> > trips;
 	append_energies(options,data,trips);
+	if (trips.size()==0)
+	{
+		printf("**admmpd::Solver Error: No reduction coeffs\n");
+		return false;
+	}
 	int n_row_D = trips.back().row()+1;
 	double dt2 = options->timestep_s * options->timestep_s;
 	if (options->timestep_s <= 0)
@@ -338,6 +355,8 @@ void Solver::compute_matrices(
 	data->u.resize(n_row_D,3);
 	data->u.setZero();
 
+	return true;
+
 } // end compute matrices
 
 void Solver::compute_masses(
diff --git a/extern/softbody/src/admmpd_solver.h b/extern/softbody/src/admmpd_solver.h
index 0305d535f39..c690bae189c 100644
--- a/extern/softbody/src/admmpd_solver.h
+++ b/extern/softbody/src/admmpd_solver.h
@@ -107,11 +107,7 @@ protected:
         const Options *options,
         Data *data);
 
-    void compute_lattice(
-        const Options *options,
-        Data *data);
-
-    void compute_matrices(
+    bool compute_matrices(
         const Options *options,
         Data *data);
 
diff --git a/intern/softbody/admmpd_api.cpp b/intern/softbody/admmpd_api.cpp
index ade895a1971..f60097666c5 100644
--- a/intern/softbody/admmpd_api.cpp
+++ b/intern/softbody/admmpd_api.cpp
@@ -38,7 +38,7 @@
 struct ADMMPDInternalData {
   admmpd::Options *options;
   admmpd::Data *data;
-//  admmpd::Lattice *lattice;
+  admmpd::Lattice *lattice;
   int in_totverts; // number of input verts
 };
 
@@ -57,23 +57,18 @@ void admmpd_dealloc(ADMMPDInterfaceData *iface)
         delete iface->data->options;
     if(iface->data->data)
         delete iface->data->data;
+    if(iface->data->lattice)
+        delete iface->data->lattice;
     delete iface->data;
   }
 
   iface->data = NULL;
 }
 
-int admmpd_init(ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_faces)
+static int admmpd_init_with_tetgen(
+  ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_faces,
+  Eigen::MatrixXd *V, Eigen::MatrixXi *T)
 {
-
-  if (iface==NULL)
-    return 0;
-  if (in_verts==NULL || in_faces==NULL)
-    return 0;
-  if (iface->mesh_totverts<=0 || iface->mesh_totfaces<=0)
-    return 0;
-
-  // Generate tets
   TetGenRemeshData tg;
   init_tetgenremeshdata(&tg);
   tg.in_verts = in_verts;
@@ -83,29 +78,70 @@ int admmpd_init(ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_fa
   bool success = tetgen_resmesh(&tg);
   if (!success || tg.out_tottets==0)
     return 0;
-
+  
   // Create initializer for ADMMPD
   iface->totverts = tg.out_totverts;
   int nv = tg.out_totverts;
   int nt = tg.out_tottets;
-  Eigen::MatrixXd V(nv,3);
-  Eigen::MatrixXi T(nt,4);
-  V.setZero();
+  V->resize(nv,3);
+  T->resize(nt,4);
+  V->setZero();
   for (int i=0; i<nv; ++i)
   {
     for (int j=0; j<3; ++j)
     {
-      V(i,j) = tg.out_verts[i*3+j];
+      V->operator()(i,j) = tg.out_verts[i*3+j];
     }
   }
-  T.setZero();
+  T->setZero();
   for (int i=0; i<nt; ++i)
   {
-    T(i,0) = tg.out_tets[i*4+0];
-    T(i,1) = tg.out_tets[i*4+1];
-    T(i,2) = tg.out_tets[i*4+2];
-    T(i,3) = tg.out_tets[i*4+3];
+    T->operator()(i,0) = tg.out_tets[i*4+0];
+    T->operator()(i,1) = tg.out_tets[i*4+1];
+    T->operator()(i,2) = tg.out_tets[i*4+2];
+    T->operator()(i,3) = tg.out_tets[i*4+3];
   }
+  // Clean up tetgen data
+  MEM_freeN(tg.out_tets);
+  MEM_freeN(tg.out_facets);
+  MEM_freeN(tg.out_verts);
+  return 1;
+}
+
+static int admmpd_init_with_lattice(
+  ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_faces,
+  Eigen::MatrixXd *V, Eigen::MatrixXi *T)
+{
+  (void)(in_faces);
+  int nv = iface->mesh_totverts;
+  Eigen::MatrixXd in_V(nv,3);
+  for (int i=0; i<nv; ++i)
+  {
+    for (int j=0; j<3; ++j)
+    {
+      in_V(i,j) = in_verts[i*3+j];
+    }
+  }
+
+  iface->totverts = 0;
+  bool success = iface->data->lattice->generate(in_V,V,T);
+  if (success)
+  {
+    iface->totverts = V->rows();
+    return 1;
+  }
+  return 0;
+}
+
+int admmpd_init(ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_faces)
+{
+
+  if (iface==NULL)
+    return 0;
+  if (in_verts==NULL || in_faces==NULL)
+    return 0;
+  if (iface->mesh_totverts<=0 || iface->mesh_totfaces<=0)
+    return 0;
 
   // Generate solver data
   iface->data = new ADMMPDInternalData();
@@ -113,6 +149,27 @@ int admmpd_init(ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_fa
   admmpd::Options *options = iface->data->options;
   iface->data->data = new admmpd::Data();
   admmpd::Data *data = iface->data->data;
+  iface->data->lattice = new admmpd::Lattice();
+
+  // Generate tets and vertices
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi T;
+  int gen_success = 0;
+  switch (iface->init_mode)
+  {
+    default:
+    case 0:
+      gen_success = admmpd_init_with_tetgen(iface,in_verts,in_faces,&V,&T);
+      break;
+    case 1:
+      gen_success = admmpd_init_with_lattice(iface,in_verts,in_faces,&V,&T);
+      break;
+  }
+  if (!gen_success)
+  {
+    printf("**ADMMPD Failed to generate tets\n");
+    return 0;
+  } 
 
   // Initialize
   bool init_success = false;
@@ -125,11 +182,13 @@ int admmpd_init(ADMMPDInterfaceData *iface, float *in_verts, unsigned int *in_fa
     printf("**ADMMPD Error on init: %s\n", e.what());
   }
 
-  // Clean up tetgen data
-  MEM_freeN(tg.out_tets);
-  MEM_freeN(tg.out_facets);
-  MEM_freeN(tg.out_verts);
-  return int(init_success);
+  if (!init_success)
+  {
+    printf("**ADMMPD Failed to initialize\n");
+    return 0;
+  }
+
+  return 1;
 }
 
 void admmpd_copy_from_bodypoint(ADMMPDInterfaceData *iface, const BodyPoint *pts)
@@ -165,17 +224,24 @@ void admmpd_copy_to_bodypoint_and_object(ADMMPDInterfaceData *iface, BodyPoint *
       }
     }
 
-    // TODO eventually replace with mapping. For now
-    // we assume TetGen, which the first N vertices are
-    // mapped one-to-one.
-    if (vertexCos != NULL && i<iface->mesh_totverts)
+    if (vertexCos != NULL && iface->init_mode==0 && i<iface->mesh_totverts)
     {
       vertexCos[i][0] = iface->data->data->x(i,0);
       vertexCos[i][1] = iface->data->data->x(i,1);
       vertexCos[i][2] = iface->data->data->x(i,2);
     }
+  } // end loop all verts
+
+  if (vertexCos != NULL && iface->init_mode==1)
+  {
+//    Eigen::Vector3d xi = iface->data->lattice->get_mapped_vertex(
+//      i, &iface->data->data->x, &iface->data->data->tets);
+//    vertexCos[i][0] = xi[0];
+//    vertexCos[i][1] = xi[1];
+//    vertexCos[i][2] = xi[2];
   }
-}
+
+} // end map ADMMPD to bodypoint and object
 
 void admmpd_solve(ADMMPDInterfaceData *iface)
 {
@@ -183,6 +249,9 @@ void admmpd_solve(ADMMPDInterfaceData *iface)
   if (iface == NULL)
     return;
 
+  if (iface->data == NULL || iface->data->options == NULL || iface->data->data == NULL)
+    return;
+
   try
   {
     admmpd::Solver().solve(iface->data->options,iface->data->data);
diff --git a/intern/softbody/admmpd_api.h b/intern/softbody/admmpd_api.h
index 36cbe684e9e..b90a89360bd 100644
--- a/intern/softbody/admmpd_api.h
+++ b/intern/softbody/admmpd_api.h
@@ -36,6 +36,7 @@ typedef struct ADMMPDInterfaceData {
     int totverts; // number of deformable verts (output)
     int mesh_totverts; // number of surface mesh vertices (input)
     int mesh_totfaces; // number of surface mesh faces (input)
+    int init_mode; // 0=tetgen, 1=lattice
     // Solver data used internally
     struct ADMMPDInternalData *data;
 } ADMMPDInterfaceData;
diff --git a/source/blender/blenkernel/intern/softbody.c b/source/blender/blenkernel/intern/softbody.c
index 29e37c67a9c..38cb9f0db08 100644
--- a/source/blender/blenkernel/intern/softbody.c
+++ b/source/blender/blenkernel/intern/softbody.c
@@ -949,7 +949,7 @@ static void free_softbody_intern(SoftBody *sb)
     sb->bpoint = NULL;
     sb->bspring = NULL;
 
-    if (sb->admmpd) {
+    if (sb->admmpd != NULL) {
       admmpd_dealloc(sb->admmpd);
       MEM_freeN(sb->admmpd);
       sb->admmpd = NULL;
@@ -3579,6 +3579,12 @@ static void init_admmpd_interface(Object *ob, float (*vertexCos)[3])
   Mesh *me = ob->data;
   SoftBody *sb = ob->soft;
 
+  if (sb->admmpd != NULL)
+  {
+    admmpd_dealloc(sb->admmpd);
+    MEM_freeN(sb->admmpd);
+    sb->admmpd = NULL;
+  }
   sb->admmpd = (ADMMPDInterfaceData*)MEM_callocN(sizeof(ADMMPDInterfaceData), "SoftBody_ADMMPD");
 
   // Resize data
@@ -3615,6 +3621,7 @@ static void init_admmpd_interface(Object *ob, float (*vertexCos)[3])
   looptri = NULL;
 
   // Initalize solver and free tmp data
+  sb->admmpd->init_mode = 1;
   admmpd_init(sb->admmpd, in_verts, in_faces);
   MEM_freeN(in_verts);
   MEM_freeN(in_faces);
@@ -3686,7 +3693,7 @@ void sbObjectStep_admmpd(
   }
 
   // If it's the first frame we reset the cache and data
-  if (framenr == startframe) {
+  if (framenr == startframe || sb->admmpd == NULL) {
     BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
     init_admmpd_interface(ob,vertexCos);
     BKE_ptcache_validate(cache, framenr);
@@ -3696,13 +3703,6 @@ void sbObjectStep_admmpd(
     return;
   }
 
-  // Initialize simulation data if needed.
-  if (sb->admmpd == NULL)
-  {
-    printf("\n\n\n\n\n**ADMMPD data null!\n\n\n\n\n");
-    return;    
-  }
-
   // Read from cache
   {
     bool can_write_cache = DEG_is_active(depsgraph);