Enhancing 3D Models:
Animation and Product Designs with Subdivision Surfaces

Alec Shedelbower, Wolfram Research

#WolframTechConf

Talk Overview

Challenges of Modeling Smooth Surfaces

How Subdivision Surfaces Work

SubdivisionRegion Representation and Usage

Workflows

Applications

Challenges of Smooth
Surface Modeling

The Polygon Paradigm

Polygons balance computational simplicity with flexibility.

images=Image[#,ImageSize->200]&/@{,};
Grid[{images,Style[#,"Text",12]&/@{"GPU Hardware","Modeling Software"}},Spacings->{3, 0}]

Reality is Smooth, Polygons are Not

Most objects contain some continuous curvature, while polygons are fundamentally discrete.

mesh=ResourceData["Sample 3D Model: Fan Disk"];

smoothModels={
SubdivisionRegion[ResourceData["Sample 3D Model: Seashell"],2],
Module[{a=,b=},CSGRegion["Difference",{Style[a,White],Style[b,White]}]],
SubdivisionRegion[RepairMesh[ResourceData["Sample 3D Model: Hammerhead Shark"]],1]
};
smoothModels[[{1,3}]]=Table[MeshRegion[m,VertexNormals->Region`Mesh`MeshCellNormals[m,0]],{m,smoothModels[[{1,3}]]}];
smoothModels=Graphics3D[{EdgeForm[],#},ImageSize->{Automatic,100},Boxed->False,Method->{"ShrinkWrap" -> True}]&/@smoothModels;
Grid[{smoothModels}]

Solution — Use More Polygons

Smooth surfaces are approximated by using more polygons along with some shading tricks.

Clear[renderSmoothNormals]
renderSmoothNormals[mesh_,nt_]:=
Block[{faces,vn,fn,coords,flat,smooth,interp,icoords,ifaces,fcount,temp,gc,hash},
faces=Join@@(MeshCells[mesh,2,"Multicells"->True][[All,1]]);
vn=Region`Mesh`MeshCellNormals[mesh,0];
fn=Region`Mesh`MeshCellNormals[mesh,2];
coords=MeshCoordinates[mesh];
hash=Hash[mesh];
If[ListQ[$meshNormalsCache[hash]],
{flat,smooth}=$meshNormalsCache[hash],
flat=Join@@Table[ConstantArray[fn[[fi]],Length[faces[[fi]]]],{fi,Length[faces]}];
smooth=Join@@Table[vn[[faces[[fi]]]],{fi,Length[faces]}];
$meshNormalsCache[hash]={flat,smooth}
];

interp=nt*smooth+flat*(1-nt);

icoords=Join@@Table[coords[[faces[[fi]]]],{fi,Length[faces]}];
fcount=1;
ifaces=Table[temp=Range[fcount,fcount+Length[faces[[fi]]]-1];fcount+=Length[faces[[fi]]];temp,{fi,Length[faces]}];

gc=GraphicsComplex[icoords,Polygon[ifaces],VertexNormals->interp];

Graphics3D[{EdgeForm[],gc},Method->{"ShrinkWrap" -> True},Boxed->False,ImageSize->{Automatic,170}]
]
$normalLODs:=
Block[{base,lods},
base=ResourceData["Stanford Bunny"];
lods=SimplifyMesh[base,{"MaxVertexCount",#}]&/@{100,200,400,800,4000};
lods[[3;;-1]]=Remesh/@lods[[3;;-1]];
$normalLODs=lods;
$normalLODs
]

Manipulate[renderSmoothNormals[$normalLODs[[li]],nt],
Style["Polygons",Bold],
{{li,1,""},1,Length[$normalLODs],1},
Style["Normal Smoothing",Bold],
{{nt,0,""},0,1},
ControlPlacement->Left
]

A Job for Machines

Manually adding polygons, even with editor tooling, is tedious and time consuming.

meshes=SubdivisionRegion[,{0,1,2,3},ImageSize->120];
labels=Table[ToString[MeshCellCount[m,2]],{m,meshes}];
Grid[{meshes,Style[#<>" Triangles","Text",12]&/@labels}]

Is there a way to automate it?

Splines — Smooth by Design

Splines are a fundamentally smooth representation for curves/surfaces.

cpts=BlockRandom[Table[{i,j,RandomReal[{-1,1}]},{i,5},{j,5}],RandomSeeding->10];
surf=BSplineSurface[cpts];
font=Graphics[{Haloing[Black,1],Text[Style["Curve",80,,FontFamily->"Charm"]]},ImageSize->{200,100}];
surf=Graphics3D[{EdgeForm[AbsoluteThickness[1]],surf},Boxed->False,ImageSize->200];
Grid[{{font,surf}},Spacings->{3, 1}]

Splines — High Level Control

Define surfaces through a structured grid of control points. Each point on the surface is a weighted average of those control points.

locatorSpline=DynamicModule[{pt={{1,1},{3,5},{5,0},{7,5},{9,0}}},LocatorPane[Dynamic[pt],Dynamic@Graphics[{{GrayLevel[0.9],Line[pt]},BSplineCurve[pt,SplineDegree->2],Red,Point[pt]},PlotRange->{{0,10},{0,5}},PlotRangePadding->0.2],Appearance->None,ImageSize->200]];
locatorSegment=DynamicModule[{pt=N@Table[{t,4(Sin[t]+1)/2},{t,0,9,0.5}]},LocatorPane[Dynamic[pt],Dynamic@Graphics[{{GrayLevel[0.9],Line[pt]},Line[pt],Red,Point[pt]},PlotRange->{{0,10},{0,5}},PlotRangePadding->0.4],Appearance->None,ImageSize->200]];
Grid[{{locatorSegment,locatorSpline},Style[#,"Text",12]&/@{"Traditional","Spline"}},Spacings->2]

Splines — Conversion to Polygons

Spline regions are well-parameterized, making it easy to get polygons.

Manipulate[
DynamicModule[{cpts,func,reg,coords,polygons},
cpts=BlockRandom[Table[{i,j,RandomReal[{-1,1}]},{i,5},{j,5}],RandomSeeding->34];
func=BSplineFunction[cpts,SplineDegree->2];
reg=BSplineSurface[cpts,SplineDegree->2];
coords=Table[func[u,v],{u,0,1,1/n},{v,0,1,1/n}];
polygons=Table[Polygon[{coords[[i,j]],coords[[i+1,j]],coords[[i+1,j+1]],coords[[i,j+1]]}],{j,Length[coords]-1},{i,Length[coords]-1}];
Graphics3D[{polygons},Boxed->False,ImageSize->{Automatic,170}]
],
Style["Polygons per Side",Bold],
{{n,1,""},1,20,1},
ControlPlacement->Left
]

Splines — The Price of Perfection

Complex Shape Difficulty: Surfaces constructed from rectangular spline patches, which is challenging for some topologies.

Patch Alignment Issues: Adjacent patches must connect seamlessly to avoid discontinuities.

Is there a way to combine the robustness of polygons with the automated smoothness and high-level control of splines?

Subdivision Regions:
The Best of Both Worlds

What are Subdivision Regions?

Regions constructed by an iterative process that creates smooth meshes from coarse ones.

base=;
sub=SubdivisionRegion[base,2];
sub=MeshRegion[MeshCoordinates[sub],MeshCells[sub,2],VertexNormals->EstimatedPointNormals[sub],Method->{"CoplanarityTolerance"->Infinity}];
imgSize=150;
sub=Rasterize[Show[sub,Method->{"ShrinkWrap"->True}],RasterSize->4*imgSize,ImageSize->imgSize];
base=Show[base,ImageSize->imgSize,Method->{"ShrinkWrap"->True}];
row={base,sub};
Grid[{Riffle[row,$arrow]},Spacings->{2, 1}]

Subdivision Steps

For each polygon in our mesh:

Subdivide – Split the polygon into smaller polygons (often four)

Smooth – Adjust vertex positions based on local neighborhood

z=0;
cell=Polygon[{{-1,-1,z},{1,-1,z},{1,1,z},{-1,1,z}}];
scells=Table[Polygon[{{pr[[1,1]],pr[[2,1]],z},{pr[[1,2]],pr[[2,1]],z},{pr[[1,2]],pr[[2,2]],z},{pr[[1,1]],pr[[2,2]],z}}],{pr,{{{-1,0},{-1,0}},{{0,1},{-1,0}},{{0,1},{0,1}},{{-1,0},{0,1}}}}];
acells=MeshPrimitives[SubdivisionRegion[Region`LinearToMeshRegion[cell],1],2];

cell={cell};
Framed[#,FrameMargins->{{0,0},{0,0}},FrameStyle->None]&@Grid[{Riffle[Graphics3D[{EdgeForm[{,AbsoluteThickness[1.5]}],,#},ViewPoint->{1,-2.4`,3.`},Lighting->,Boxed->False,ImageSize->120]&/@{cell,scells,acells},$arrow],Riffle[Style[#,"Text",12]&/@{"Original","Subdivide","Smooth"},""]},Spacings->{0.75,1}]

Subdivision Scheme

Determines exactly how polygons are subdivided and smoothed.

style={BaseStyle->EdgeForm[{Black,AbsoluteThickness[1.25]}],ImageSize->80};

pentagon=quad=MeshRegion[,style];
quad=MeshRegion[,style];
tri=MeshRegion[,style];
columns={
{tri,quad,pentagon},
ConstantArray[$arrow,3],
{SubdivisionRegion[tri,1,"CatmullClark",style],SubdivisionRegion[quad,1,"CatmullClark",style],SubdivisionRegion[pentagon,1,"CatmullClark",style]},
{SubdivisionRegion[tri,1,"QuadTri",style],SubdivisionRegion[quad,1,"QuadTri",style],SubdivisionRegion[pentagon,1,"QuadTri",style]},
{SubdivisionRegion[tri,1,"Loop",style],Style["–",Bold,LightGray,32],Style["–",Bold,LightGray,32]}
};

headers=Style[#,"Text",12]&/@{"Original","","Catmull-Clark","QuadTri","Loop"};

Grid[Prepend[Transpose[columns],headers],Spacings->{1,0.5}]

Subdivision Level

Describes how many iterations have been applied to the base mesh.

levels=SubdivisionRegion[Cube[],{0,1,2,0},ImageSize->Tiny];
labels=Style[Row[{"Level ",#}],"Text",12]&/@{0,1,2,1};
labels[[4]]="";
levels[[4]]="...";
Grid[{Riffle[levels,$arrow],Riffle[labels,""]},Spacings->{0.75, 1}]

Limit Surface

Surface converges as the level approaches infinity.

levels=SubdivisionRegion[Cube[],{3,4,0,Infinity},ImageSize->Tiny];
labels=Style[Row[{"Level ",#}],"Text",12]&/@{3,4,0,Style["∞",20]};
levels[[3]]="...";
labels[[3]]="";
Grid[{Riffle[levels,$arrow],Riffle[labels,""]},Spacings->{0.75, 1}]

Limit surface is mathematically smooth ((C^(1+))-continuous*).

*Continuity varies between subdivision schemes, and while they have at least (C^1) across most of their surfaces, there may be a finite number extraordinary points of lower continuity.

Subdivision Regions as Generalized Splines

columns={
{Style["Splines","Subsection"],
Style[Row[{Style["  ✓ ",],"Mathematically Smooth"}],"Text"],
Style[Row[{Style["  ✓ ",],"High-Level Shape Control"}],"Text"],
Style[Row[{Style["  ✓ ",],"Multiple Polygon Resolution"}],"Text"],
Style[Row[{Style["  ✗ ",],"Topologically Difficult"}],"Text"]
},
{Style["Subdivision Regions","Subsection"],
Style[Row[{Style["  ✓ ",],"Mathematically Smooth"}],"Text"],
Style[Row[{Style["  ✓ ",],"High-Level Shape Control"}],"Text"],
Style[Row[{Style["  ✓ ",],"Multiple Polygon Resolution"}],"Text"],
Style[Row[{Style["  ✓ ",],"Topologically Robust"}],"Text"]
}
};
Grid[{Grid[List/@#,Alignment->Left,Spacings->{1, 1}]&/@columns},Spacings->{2, 1}]

SubdivisionRegion
Representation & Usage

Usage — Subdivision Level

({{, SubdivisionRegion[mesh,n]
gives the n(^th) subdivision level of mesh.}})

SubdivisionRegion[,1]

Usage — Subdivision Level

({{, SubdivisionRegion[mesh,n]
gives the n(^th) subdivision level of mesh.}})

Manipulate[SubdivisionRegion[,level],{level,0,3,1}]

Usage — Limit Surface

({{, SubdivisionRegion[mesh] 
represents the limit region of applying a subdivision scheme to the mesh region mesh.}})

SubdivisionRegion[]

Usage — Limit Surface

({{, SubdivisionRegion[mesh] 
represents the limit region of applying a subdivision scheme to the mesh region mesh.}})

RegionEqual[
SubdivisionRegion[],
SubdivisionRegion[,Infinity]
]

Usage — Schemes

({{, SubdivisionRegion[…,proc]
uses the subdivision refinement scheme proc.}})

SubdivisionRegion[,2,"Loop"]

Usage — Schemes

({{, SubdivisionRegion[…,proc]
uses the subdivision refinement scheme proc.}})

SubdivisionRegion[,2,#,ImageSize->100]&/@{"CatmullClark","Loop","QuadTri"}

Usage — Schemes

({{, SubdivisionRegion[…,proc]
uses the subdivision refinement scheme proc.}})

base=;
size=80;
regions={
Region[base,ImageSize->size],
SubdivisionRegion[base,"CatmullClark",ImageSize->size],
SubdivisionRegion[base,"Loop",ImageSize->size],
SubdivisionRegion[base,"QuadTri",ImageSize->size]
};
regions=Show[#,Method->{"ShrinkWrap"->True}]&/@regions;
regions[[1]]=Framed[regions[[1]],FrameStyle->None,FrameMargins->{{0,40},{0,0}}];
labels=Style[#,"Text",12]&/@{"Base","Catmull-Clark","Loop","QuadTri"};
Grid[{regions},Spacings->{3, 1}]

base=;
size=80;
regions={
Region[base,ImageSize->size],
SubdivisionRegion[base,"CatmullClark",ImageSize->size],
SubdivisionRegion[base,"Loop",ImageSize->size],
SubdivisionRegion[base,"QuadTri",ImageSize->size]
};
regions=Show[#,Method->{"ShrinkWrap"->True}]&/@regions;
regions[[1]]=Framed[regions[[1]],FrameStyle->None,FrameMargins->{{0,40},{0,0}}];
labels=Style[#,"Text",12]&/@{"Base","Catmull-Clark","Loop","QuadTri"};
Grid[{regions,labels},Spacings->{3, 1}]

Subdivision Workflows

Subdivision Workflows

There are two main workflows for using subdivision surfaces:

Mesh Processing Operation

Control Mesh Modeling

Workflow — Mesh Processing Operation

Treat it as another mesh processing step:

SimplifyMesh ⇒ Reduce polygon count

Remesh ⇒ Conform polygons

SmoothMesh ⇒ Reduce high frequency noise

SubdivisionRegion ⇒ Round off sharp corners

Workflow — Mesh Processing Operation

SubdivisionRegion[,3]

Workflow — Control Mesh Modeling

Treat the base mesh like the control cage in spline modeling.

Save and edit control cage

Apply surface subdivision as a post processing effect

Workflow — Control Mesh Modeling

regs={,,};
subs=SubdivisionRegion[#,3]&/@regs;

regLabels=Table[
cage=MeshCellCount[regs[[i]],2];
Style[Row[{"Polygons: ",ToString[cage],Style[If[i>1," (+"<>ToString[cage-MeshCellCount[regs[[i-1]],2]]<>")",""],]}],"Text",12]
,{i,Length[regs]}];

subLabels=Table[
sub=MeshCellCount[subs[[i]],2];
Style[Row[{"Polygons: ",ToString[sub],Style[If[i>1," (+"<>ToString[sub-MeshCellCount[subs[[i-1]],2]]<>")",""],]}],"Text",12]
,{i,Length[subs]}];

imgSize={Automatic,120};
regs=Show[#,ImageSize->imgSize]&/@regs;
subs=Show[Graphics3D[{EdgeForm[],#},Boxed->False],ImageSize->imgSize]&/@subs;

rows={
regs,regLabels,
subs,subLabels
};
rows=Table[Riffle[rows[[ri]],If[EvenQ[ri],"",$arrow]],{ri,Length@rows}];
Grid[rows,Spacings->{1, 0}]

Subdivision Applications

Applications — Product Rendering

Produce high-fidelity meshes that hold up under close scrutiny.

base=;
subdivided=SubdivisionRegion[base,2];
smooth=MeshRegion[subdivided,VertexNormals->EstimatedPointNormals[subdivided,Method->"FaceWeighted"],];

connectedMeshComponents[mesh_]:=
Block[{a2,comps,faces,coords,res},
a2=MeshConnectivityGraph[mesh,2];
faces=MeshCells[mesh,2];
comps=ConnectedComponents[a2][[All,All,2]];
coords=MeshCoordinates[mesh];
res=Table[MeshRegion[coords,faces[[comp]],Method->{"CoplanarityTolerance"->Infinity}],{comp,comps}];
res
]
components=connectedMeshComponents[smooth];

smoothComponents=Table[MeshRegion[comp,VertexNormals->EstimatedPointNormals[comp,Method->"FaceWeighted"],],{comp,components}];
separated={
MaterialShading[{"Rubber",}],smoothComponents[[1]],
MaterialShading[{"Rubber",}],smoothComponents[[3]],
MaterialShading[{"Plastic",}],smoothComponents[[{2,4,5,6}]],
MaterialShading[{"Plastic",}],smoothComponents[[{7}]],
MaterialShading[{"Plastic",}],smoothComponents[[{8}]]
};

graphics=Table[Graphics3D[{EdgeForm[],object},Boxed->False,Lighting->"ThreePoint",ImageSize->140,Method->{"ShrinkWrap" -> True}],{object,{{EdgeForm[{Black,AbsoluteThickness[0.25]}],base},smooth,separated}}];

Framed[Grid[{Riffle[graphics,$arrow]},Spacings->1],FrameMargins->{{0,0},{30,30}},FrameStyle->None]

Applications — Product Rendering

Grid[{Image[#,ImageSize->220]&/@{,}},Spacings->{3, 1}]

Applications — Level of Detail (LOD) Generation

Optimize render times by reducing triangle counts for distant objects.

levels=SubdivisionRegion[,{0,1,2},]

Grid[]

Applications — Animation and Simulation

Render deformed meshes without issue.

simulateAndRenderCloth[]

Big Picture — When to use Subdivision Surfaces?

columns={
{Style["Ideal For...","Subsection"],
Style[Row[{Style["  ✓ ",],"Modeling Smooth Surfaces"}],"Text"],
Style[Row[{Style["  ✓ ",],"Adaptive Fidelity (LOD)"}],"Text"],
Style[Row[{Style["  ✓ ",],"Animated/Deformed Surfaces"}],"Text"],
Style[Row[{Style["  ✓ ",],"Rounding Off Simple Meshes"}],"Text"]
},
{Style["Not Ideal For...","Subsection"],
Style[Row[{Style["  ✗ ",],"Precise Measurements (CAD)"}],"Text"],
Style[Row[{Style["  ✗ ",],"Non-Manifold Meshes"}],"Text"],
Style[Row[{Style["  ✗ ",],"Smoothing Large Meshes"}],"Text"]
}
};
Grid[{Grid[List/@#,Alignment->Left,Spacings->{1, 1}]&/@columns},Spacings->{2, 1},Alignment->Top]

Optimizing 3D Models: Techniques and Best Practices for Efficient Mesh Processing

Jaebum Jung, Thursday 3:00 PM CDT

Wolfram R&D YouTube Channel
@WolframRD

#WolframTechConf

Initialization Code

Graphics

Arrow

$arrow=;

Cloth Simulation

constructVerletDataForRectangularCloth

Clear[constructVerletDataForRectangularCloth]
constructVerletDataForRectangularCloth[{nx_,ny_},size:{sx_,sy_}]:=
Block[{neighbors,positions,adjacencyList,adj,dists,indices,nindex,index,corners,neighborOffsets,neighborCounts,accumAdjacentCount,accum},
positions=Join@@Table[{Rescale[xi,{1,nx},{-sx/2.0,sx/2.0}],Rescale[yi,{1,ny},{-sy/2.0,sy/2.0}],0.0},{xi,nx},{yi,ny}];
neighborOffsets=Join[
{{-1,-1},{-1,0},{-1,1},{0,-1},{0,1},{1,-1},{1,0},{1,1}},
{{-2,0},{2,0},{0,-2},{0,2}},
{{-2,-2},{2,-2},{2,2},{-2,2}}
];
neighborCounts=CreateDataStructure["DynamicArray"];
adjacencyList=Join@@Table[
neighbors=Table[{xi,yi}+offset,{offset,neighborOffsets}];
neighbors=Table[If[neighbor[[1]]<=0||neighbor[[2]]<=0||neighbor[[1]]>nx||neighbor[[2]]>ny,Nothing,(neighbor[[1]]-1)*ny+neighbor[[2]]],{neighbor,neighbors}];
neighborCounts["Append",Length[neighbors]];
neighbors,
{xi,nx},{yi,ny}
];
dists=Table[Table[
EuclideanDistance[positions[[i]],positions[[ni]]],
{ni,adjacencyList[[i]]}
],
{i,Length[adjacencyList]}
];
index=0;
indices=Join@@Table[
corners={xi,yi}+#&/@{{0,0},{1,0},{1,1},{0,1}};
Table[(c[[1]]-1)*ny+c[[2]],{c,corners}]
,{xi,nx-1},{yi,ny-1}
];
accum=0;
accumAdjacentCount=Table[accum+=neighborCounts[[i]];accum,{i,Length[positions]}];
{positions,Join@@adjacencyList,Join@@dists,accumAdjacentCount,indices}
]
(*{positions,adjacency,dists,accumAdjacentCount,indices}=constructVerletDataForRectangularCloth[{3,3},{1,1}];
{positions,adjacency,dists,accumAdjacentCount,indices}//Column*)

verletStep

func=Function[{
Typed[prev,"PackedArray"::["MachineReal",2]],
Typed[icurr,"PackedArray"::["MachineReal",2]],
Typed[movability,"PackedArray"::["MachineReal",1]],
Typed[accelerations,"PackedArray"::["MachineReal",2]],
Typed[frictions,"PackedArray"::["MachineReal",1]],
Typed[distances,"PackedArray"::["MachineReal",1]],
Typed[adjacency,"PackedArray"::["MachineInteger",1]],
Typed[adjacencyCount,"PackedArray"::["MachineInteger",1]],
Typed[accumAdjacentCount,"PackedArray"::["MachineInteger",1]],
Typed[ballCenters,"PackedArray"::["MachineReal",2]],
Typed[ballRadii,"PackedArray"::["MachineReal",1]],
Typed[floorHeight,"MachineReal"],
Typed[Δt,"MachineReal"],
Typed[maxIterations,"MachineInteger"],
Typed[distanceThreshold,"MachineReal"],
Typed[defaultFriction,"MachineReal"],
Typed[collisionFriction,"MachineReal"],
Typed[staticFrictionMinSpeed,"MachineReal"]
},
Module[{n,next,curr,velocities,index,collidedQ,ai,pos1,pos2,ratio,delta,dist,targetDist,diff,corr,center,radius,iters,length,collisionBuffer,maxAdjustedDist,nextFrictions,vel,res1},
n=Length[prev];
curr=icurr;

collisionBuffer=0.02;

velocities=(curr-prev)/Δt;
velocities+=accelerations;
Do[
vel=velocities[[i]];
velocities[[i]]=If[
frictions[[i]]<=collisionFriction&&Dot[vel,vel]<=staticFrictionMinSpeed*staticFrictionMinSpeed,{0.,0.,0.},vel],
{i,n}
];
velocities*=movability;
curr+=Δt*(velocities*frictions);

nextFrictions=ConstantArray[defaultFriction,n];

(* Solve Iterations *)
iters=0;
maxAdjustedDist=100000.0;
Do[
If[maxAdjustedDist<=distanceThreshold,Break[]];

maxAdjustedDist=0.0;

iters++;
next=ConstantArray[0.0,{n,3}];

(* Distance Constraint *)
Do[
index=If[i===1,1,accumAdjacentCount[[i-1]]]+j;
ai=adjacency[[index]];


If[i<ai,Continue[]];
If[movability[[ai]]+movability[[i]]<=0.0,Continue[]];

targetDist=distances[[index]];

pos1=curr[[i]];
pos2=curr[[ai]];
delta=pos1-pos2;
dist=Norm[delta];
If[dist<0.00001,
Continue[];
];
diff=(dist-targetDist);
ratio=movability[[ai]]/(movability[[ai]]+movability[[i]]);

corr=delta*diff/dist;
pos1=pos1-corr*ratio;
pos2=pos2+corr*(1.0-ratio);

next[[i]]+=pos1;
next[[ai]]+=pos2;

,{i,n},{j,adjacencyCount[[i]]}
];
next/=adjacencyCount;

Do[
pos1=next[[i]];
Do[
center=ballCenters[[bi]];
radius=ballRadii[[bi]];
radius+=collisionBuffer;
length=Norm[pos1-center];
If[length>radius,Continue[]];
(*continueIterating=True;*)
pos1=center+((pos1-center)/length)*(radius);
nextFrictions[[i]]=collisionFriction,
{bi,Length[ballCenters]}
];
next[[i]]=pos1,
{i,n}
];

Do[
pos1=next[[i]];
If[pos1[[3]]>=floorHeight+collisionBuffer,Continue[]];
pos1[[3]]=floorHeight+collisionBuffer;
next[[i]]=pos1;
nextFrictions[[i]]=collisionFriction;
,{i,n}
];

Do[dist=Norm[curr[[i]]-next[[i]]];maxAdjustedDist=Max[maxAdjustedDist,dist],{i,n}];

curr=next;
,maxIterations
];

res1={curr,nextFrictions};
{res1,iters}
]
];

(*verletStep:=(verletStep=FunctionCompile[func])*)

verletStep=CompiledCodeFunction[];

runVerletSimulation

Clear[runVerletSimulation]
runVerletSimulation[data:{positions_,adjacency_,distances_,accumAdjacentCount_,indices_},pinned_,{ballFunc_Function,floorHeight_},Δt_,steps_,maxIterations_,distanceThreshold_,defaultFriction_,collisionFriction_,staticFrictionMinSpeed_]/;(Head[verletStep]===CompiledCodeFunction||Head[verletStep]===CompiledFunction):=
Block[{n,prev,curr,next,accelerations,frictions,adjacencyCount,balls,ballCenters,ballRadii,results,iters,percentage,elapsedTime,movability},
n=Length[positions];

prev=positions;
curr=positions;

accelerations=ConstantArray[{0,0,-9.8}*Δt^2,n];
frictions=ConstantArray[defaultFriction,n];
adjacencyCount=Table[accumAdjacentCount[[i]]-If[i===1,0,accumAdjacentCount[[i-1]]],{i,n}];

movability=If[pinned===None,ConstantArray[1.0,Length[positions]],Map[If[#,0.0,1.0]&,pinned]];

$progress=0.0;
$iterationCounts=ConstantArray[0,steps];
PrintTemporary[Dynamic[Grid[{{solverIterationPlot[$iterationCounts]},{ProgressIndicator[$progress]}}]]];

results=Table[
percentage=(si-1)/(steps-1);
elapsedTime=si*Δt;

balls=ballFunc[percentage,elapsedTime];
If[balls==={},balls={Ball[10000{1,1,1},0.0]}];
ballCenters=N[balls[[All,1]]];
ballRadii=N[balls[[All,-1]]];

{{next,frictions},iters}=verletStep[prev,curr,movability,accelerations,frictions,distances,adjacency,adjacencyCount,accumAdjacentCount,ballCenters,ballRadii,floorHeight,Δt,maxIterations,distanceThreshold,defaultFriction,collisionFriction,staticFrictionMinSpeed];
{prev,curr}={curr,next};
$progress+=1.0/steps;
$iterationCounts[[si]]=iters;
curr,
{si,steps}
];

results/;ListQ[results]
]

Clear[solverIterationPlot]
solverIterationPlot[iterations_]:=Framed[ListPlot[iterations,Joined->True,Filling->Axis,Ticks->None,Axes->None,ImageMargins->0,PlotRangePadding->None,ImageSize->{200},PlotStyle->{AbsoluteThickness[1]},PlotRange->{Automatic,{0,Max[iterations]+1}}],FrameMargins->0]

renderSimulationAtStep

Clear[renderSimulationAtStep]
renderSimulationAtStep[sim_,data_,obstacles_,step_,opts___]:=
Block[{indices},
indices=data[[5]];
DynamicModule[{t=(step-1)/(Length[sim]-1)},
Graphics3D[{{EdgeForm[Black],FaceForm[],GraphicsComplex[sim[[step]],Polygon[indices]],AbsolutePointSize[2],Point[sim[[step]]]},obstacles[t,step*Δt]},Lighting->"Accent",Boxed->False,PlotRange->{{-1,1},{-1,1},{-2,0.1}},opts]
]
]

subdividedRenderSimulationAtStep

Clear[subdividedRenderSimulationAtStep]
subdividedRenderSimulationAtStep[sim_,data_,obstacles_,step_,opts___]:=
Block[{indices,mesh,sub},
indices=data[[5]];
mesh=MeshRegion[sim[[step]],Polygon[indices],Method->{"CoplanarityTolerance"->Infinity}];
sub=SubdivisionRegion[mesh];
DynamicModule[{t=(step-1)/(Length[sim]-1)},
Graphics3D[{{MaterialShading[{"Satin",Red}],sub},obstacles[t,step*Δt]},Lighting->"ThreePoint",Boxed->False,PlotRange->{{-1,1},{-1,1},{-2,0.1}},opts]
]
]

simulateAndRenderCloth

Clear[simulateAndRenderCloth]
simulateAndRenderCloth[]:=
Block[{data,obstacles,sim},
data=constructVerletDataForRectangularCloth[{24,24},{2,2}];
obstacles=Function[{t,time},{Ball[{0,0,-1},0.5]}];
sim=runVerletSimulation[data,None,{obstacles,-2},0.05,250,50,0.005,0.95,0.6,0.1];
With[{isim=sim,idata=data,iobstacles=obstacles},
Row[{
Manipulate[renderSimulationAtStep[isim,idata,iobstacles,step,ImageSize->{180,180}],{step,1,Length[isim]-1,1}],
Manipulate[subdividedRenderSimulationAtStep[isim,idata,iobstacles,step,ImageSize->{180,180}],{step,1,Length[isim]-1,1}]
}]
]
]