9
$\begingroup$

I want to work on subdivision surfaces. Unfortunately, I don’t have any source code to start with. I need some Mathematica codes for applying Catmull-Clark and Doo-Sabin methods. I would like to request anyone on this platform who have such codes to kindly share them with me.

The codes I have from: http://hakenberg.de/subdivision/enclosed_volume.htm are giving the error code:

Part::partd: Part specification DooSabinListAlt[[1]] is longer than depth of object.

I will appreciate your help.

$\endgroup$
  • $\begingroup$ There is an implementation of Catmull-Clark on Rosetta Code $\endgroup$ – C. E. Apr 18 at 20:25
19
$\begingroup$

Catmull-Clark Subdivision

Indeed, I have some code for Catmull-Clark subdivision and I planned to post it here for quite some time. This seems to be a good opportunity.

The code is optimized for performance, so it involves a lot of CompiledFunction and SparseArray hacks. I am sorry if you find it somewhat unidiomatic.

CatmullClarkSubdivisionMatrix creates the, well, subdivision matrix while the actual subdivision is performed by CatmullClarkSubdivide. The code below assumes that the surface is a manifold-like polyhedral mesh in $\mathbb{R}^3$, possibly with boundary and not necessarily orientable.

Before we start, you might be also interested in Loop subdivision; that one is implemented here.

Application

First, we have the load the code from the section "Implementation" below. Then we can employ the function CatmullClarkSubdivide:

pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1}, {1, 2, 6, 5}, {2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};

M = polymesh[pts, polys];

MList = NestList[CatmullClarkSubdivide, M, 4];

GraphicsRow[
 GridMeshPlot /@ MList,
 ImageSize -> Full
 ]

enter image description here

A somewhat more complex example:

R = ExampleData[{"Geometry3D", "Triceratops"}, "MeshRegion"];
M = polymesh[MeshCoordinates[R], MeshCells[R, 2, "Multicells" -> True][[1, 1]]];
MList = NestList[CatmullClarkSubdivide, M, 3];
GraphicsRow[
 MeshPlot /@ MList, 
 ImageSize -> Full
 ]

enter image description here

Two ways to subdivide the boundary:

pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1},(*{1,2,6,5},*){2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};
M = polymesh[pts, polys];

With averaging, applying the standard $(1/8, 1/2, 1/8)$-subdivision along the boundary curves:

MList = NestList[
   CatmullClarkSubdivide[#, "AverageBoundaryPoints" -> True] &, 
   M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]

enter image description here

And without averaging:

enter image description here

Implementation

getEdgesFromPolygons = Compile[{{f, _Integer, 1}},
   Table[
    {
     Min[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]], 
     Max[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]]
     },
    {i, 1, Length[f]}
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

getSubdividedPolygons = 
  Compile[{{qq, _Integer, 1}, {ee, _Integer, 1}, {n, _Integer}},
   Table[
    {
     Compile`GetElement[qq, i],
     Compile`GetElement[ee, i],
     n,
     Compile`GetElement[ee, Mod[i - 1, Length[qq], 1]]
     },
    {i, 1, Length[qq]}],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

AccumulateIntegerList = Compile[{{list, _Integer, 1}},
   Block[{c = 0, r = 0},
    Table[
     If[i <= Length[list],
      r = c; c += Compile`GetElement[list, i]; r,
      c
      ]
     , {i, 1, Length[list] + 1}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

cExtractIntegerFromSparseMatrix = Compile[
   {{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer, 
     1}, {background, _Integer},
    {i, _Integer}, {j, _Integer}},
   Block[{k, c},
    k = Compile`GetElement[rp, i] + 1;
    c = Compile`GetElement[rp, i + 1];
    While[k < c + 1 && Compile`GetElement[ci, k] != j, ++k];
    If[k == c + 1, background, Compile`GetElement[vals, k]]
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

ToPack = Developer`ToPackedArray;

polymesh::usage = "";
polymesh /: polymesh[points0_, polygons0_] :=
 Module[{polygons},
  polygons = ToPack[polygons0];
  polymesh[
   Association[
    "MeshCoordinates" -> ToPack[N[points0]],
    "MeshCells" -> Association[
      0 -> Partition[Range[Length[points0]], 1],
      1 -> DeleteDuplicates[ToPack[Flatten[getEdgesFromPolygons[polygons], 1]]],
      2 -> polygons
      ]
    ]
   ]
  ];

polymesh /: MeshCoordinates[M_polymesh] := M[[1]][["MeshCoordinates"]];

polymesh /: MeshCells[M_polymesh, d_Integer] := M[[1]][["MeshCells", Key[d]]];

polymesh /: MeshCellCount[M_polymesh, d_Integer] := Length[MeshCells[M, d]];

GridMeshPlot::usage = "";
polymesh /: GridMeshPlot[M_polymesh] := Graphics3D[{
   ColorData[97][1], Specularity[White, 30], EdgeForm[{Thin, Black}];
   GraphicsComplex[MeshCoordinates[M], Polygon[MeshCells[M, 2]]],
   },
  Lighting -> "Neutral",
  Boxed -> False
  ];

MeshPlot::usage = "";
polymesh /: MeshPlot[M_polymesh] := Graphics3D[{
   ColorData[97][1], Specularity[White, 30], EdgeForm[],
   GraphicsComplex[MeshCoordinates[M], Polygon[MeshCells[M, 2]]],
   },
  Lighting -> "Neutral",
  Boxed -> False
  ];

SignedPolygonsNeighEdges::usage = "";
polymesh /: SignedPolygonsNeighEdges[M_polymesh] := 
 Module[{edges, n, A00, i, j},
  edges = MeshCells[M, 1];
  n = MeshCellCount[M, 0];
  A00 = SparseArray`SparseArraySort@SparseArray[
     Rule[
      Join[edges, Transpose[Transpose[edges][[{2, 1}]]]],
      Join[Range[1, Length[edges]], Range[-1, -Length[edges], -1]]
      ],
     {n, n}
     ];
  {i, j} = Transpose[Join @@ With[{cf = Compile[{{p, _Integer, 1}},
         Transpose[{p, RotateLeft[p]}],
         RuntimeAttributes -> {Listable},
         Parallelization -> True
         ]},
      cf[MeshCells[M, 2]]
      ]];
  Internal`PartitionRagged[
   cExtractIntegerFromSparseMatrix[
    A00["NonzeroValues"], A00["RowPointers"], 
    Flatten[A00["ColumnIndices"]], 0, i, j
    ],
   Length /@ MeshCells[M, 2]
   ]
  ];

SubdividedPolygons::usage = "";
polymesh /: SubdividedPolygons[M_polymesh] :=
  With[{
    n0 = MeshCellCount[M, 0],
    n1 = MeshCellCount[M, 1],
    n2 = MeshCellCount[M, 2]
    },
   Flatten[getSubdividedPolygons[
     MeshCells[M, 2],
     Abs[SignedPolygonsNeighEdges[M]] + n0,
     Range[1 + n0 + n1, n0 + n1 + n2]
     ], 1]
   ];

getConnectivityMatrix::usage = "";
getConnectivityMatrix[n_Integer, cells_List] :=
  With[{m = Length[cells]},
   If[m > 0,
    Module[{A, lens, nn, rp},
     lens = Compile[{{cell, _Integer, 1}},
        Length[cell],
        CompilationTarget -> "WVM",
        RuntimeAttributes -> {Listable},
        Parallelization -> True
        ][cells];
     rp = AccumulateIntegerList[lens];
     nn = rp[[-1]];
     A = SparseArray @@ {Automatic, {m, n}, 0, {1, {rp, Partition[Flatten[cells], 1]}, ConstantArray[1, nn]}}]
    ,
    {}
    ]
   ];

getMeshCellAdjacencyMatrix::usage = "";
getMeshCellAdjacencyMatrix[A_?MatrixQ, d_Integer] :=  
  If[Length[A] > 0,
   With[{B = A.A\[Transpose]},
    SparseArray[UnitStep[B - DiagonalMatrix[Diagonal[B]] - d]]
    ],
   {}
   ];

getMeshCellAdjacencyMatrix[Ad10_?MatrixQ, A0d2_?MatrixQ, d1_Integer, 
   d2_Integer] := If[(Length[Ad10] > 0) && (Length[A0d2] > 0),
   With[{B = Ad10.A0d2}, SparseArray[
     If[d1 == d2,
      UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1],
      UnitStep[B - (Min[d1, d2] + 1)]]
     ]
    ],
   {}
   ];

MeshCellAdjacencyMatrix::usage = "";

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, 0, 0 _] := SparseArray[
  Join[MeshCells[M, 1], 
    Transpose[Reverse[Transpose[MeshCells[M, 1]]]]] -> 1,
  {1, 1} MeshCellCount[M, 0]
  ];

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, 0 _, d_Integer] := 
  With[{cells = MeshCells[M, d]},
   If[Length[cells] > 0,
    Transpose[getConnectivityMatrix[MeshCellCount[M, 0], MeshCells[M, d]]],
    {}
    ]
   ];

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer, 0 _] :=
  With[{A = MeshCellAdjacencyMatrix[M, 0, d]},
   If[Length[A] > 0,
    Transpose[MeshCellAdjacencyMatrix[M, 0, d]],
    {}
    ]
   ];

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer, d_Integer] := 
 getMeshCellAdjacencyMatrix[MeshCellAdjacencyMatrix[M, d, 0], d]

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer] := 
  MeshCellAdjacencyMatrix[M, d, d];

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d1_Integer, d2_Integer] := 
  Module[{r, m1, m2},
   {m1, m2} = MinMax[{d1, d2}];
   r = getMeshCellAdjacencyMatrix[
     MeshCellAdjacencyMatrix[M, m1, 0],
     MeshCellAdjacencyMatrix[M, 0, m2],
     m1,
     m2
     ];
   If[d1 < d2, r, If[Length[r] > 0, Transpose[r], {}]]
   ];

CatmullClarkSubdivisionMatrix::usage = "";
polymesh /: CatmullClarkSubdivisionMatrix[M_polymesh, OptionsPattern[{"AverageBoundaryPoints" -> True}]] := 
  Module[{avgbndpQ, A02, A01, A10, valences, bplist, edgevalencelist, χbndp, χbndpcomp, χbnde, χbndecomp, belist, A20, A12, vB, eB, pB, n0, n1},      
   avgbndpQ = OptionValue["AverageBoundaryPoints"];       
   n0 = MeshCellCount[M, 0];
   n1 = MeshCellCount[M, 1];
   A02 = MeshCellAdjacencyMatrix[M, 0, 2];
   A20 = Transpose[A02];
   A01 = MeshCellAdjacencyMatrix[M, 0, 1];
   A10 = Transpose[A01];
   A12 = getMeshCellAdjacencyMatrix[A10, A02, 1, 2];
   valences = N[Total[A01, {2}]];
   belist = Random`Private`PositionsOf[Total[A12, {2}], 1];

   χbnde = SparseArray[Transpose[{belist}] -> 1, {n1}, 0];
   χbndecomp = (1. - Normal[χbnde]);
   bplist = Union @@ MeshCells[M, 1][[belist]];
   χbndp = SparseArray[Partition[bplist, 1] -> 1, {n0}, 0];
   χbndpcomp = (1. - Normal[χbndp]);

   pB = A20 SparseArray[1./(Length /@ MeshCells[M, 2])];
   eB = (0.5 χbnde) Transpose[χbndp A01] + SparseArray[0.25 χbndecomp] (A10 + A12.pB);
   vB = Plus[
      SparseArray[χbndpcomp/valences^2] (A02.pB + A01.A10), 
      DiagonalMatrix[SparseArray[χbndpcomp (1. - 3./valences) + If[avgbndpQ, 0.75, 1.] Normal[χbndp]]]
      ];
   If[avgbndpQ,
     vB += (0.125 χbndp) Transpose[χbndp MeshCellAdjacencyMatrix[M, 0, 0]]
     ];
   Join[vB, eB, pB]
   ];

CatmullClarkSubdivide::usage = ""; 
polymesh /: CatmullClarkSubdivide[M0_polymesh, OptionsPattern[{
    "Subdivisions" -> 1,
    "AverageBoundaryPoints" -> True
    }]
  ] :=
 Module[{t, M, A},
  M = M0;
  If[OptionValue["Subdivisions"] > 0,
   PrintTemporary["Subdividing..."];
   t = AbsoluteTiming[
      A = CatmullClarkSubdivisionMatrix[M, "AverageBoundaryPoints" -> OptionValue["AverageBoundaryPoints"]];
      M = polymesh[A.MeshCoordinates[M], SubdividedPolygons[M]];
      ][[1]];
   PrintTemporary["Subdivision done. Time elapsed: ", ToString[t]];
   ];
  If[OptionValue["Subdivisions"] > 1,
   M = CatmullClarkSubdivide[M,
     "Subdivisions" -> OptionValue["Subdivisions"] - 1,
     "AverageBoundaryPoints" -> OptionValue["AverageBoundaryPoints"]
     ]
   ];
  M
  ];
$\endgroup$
  • $\begingroup$ (+1) thanks for posting. One thing I do not get is the purpose of such sub divisions. What are they used for? $\endgroup$ – user21 Apr 18 at 11:04
  • $\begingroup$ Dear Henrik Schumacher, I appreciate your help since this is exactly what i was looking for Catmull-Clark . Anyone else with Doo-Sabin subdivision surface codes can also share them with us. $\endgroup$ – Rong Jin Apr 18 at 11:08
  • 1
    $\begingroup$ @user21 Finite elements! ;) The limiting surfaces are proven to be C^2 apart from "exceptional" points where the initial mesh does not have the generic valence (for triangle meshes, the generic valence is 6 and for quad meshes, the generic valuence is 4). So, you get $C^1$-finite elements for surfaces. But the curvature aorund the exceptional vertices blows up so much that subdivision surface are actually not overly useful for minimizing curvature energies... (you have to add many quadrature points around the curvature singularieties. $\endgroup$ – Henrik Schumacher Apr 18 at 11:08
  • 1
    $\begingroup$ @user21 The fine thing about the limiting surfaces: You can evaluate also exact positions and normals with finite cost from the initial mesh (the so-called "control mesh"). $\endgroup$ – Henrik Schumacher Apr 18 at 11:10
  • 2
    $\begingroup$ @user21 So, subdivision surfaces are not that much used for finite element analysis, although there has been a small hype in the graphics processing community a few years ago. Mostly, subdivision surfaces are used as a flexible tool for surface design quite similar how splines and NURBS surfaces are used: The vertices of the initial mesh provide you with controls to deform the limiting surface in a surprisingly intuitive and predictable way. And you have also a hierarchy of mesh resolutions, so you can generate quick previews and high quality meshes with the same pipeline. $\endgroup$ – Henrik Schumacher Apr 18 at 11:13
13
$\begingroup$

Doo-Sabin Subdivision

To my own surprise, Doo-Sabin subdivision is in many ways much easier to implement than Catmull-Clark subdivision. The only real problem I met was to compute the faces created at vertices correctly. The method I use for this is reasonable fast for meshes with not-to-high vertex degrees, but it is not guaranteed that the subdivision of an oriented mesh is again oriented.

Application

a = ConstantArray[1, {3, 3, 3}];
Do[a[[i, j, k]] = 0;, {i, {1, 3}}, {j, {1, 3}}, {k, {1, 3}}];
R = RegionBoundary[ArrayMesh[a]];
M = polymesh[MeshCoordinates[R], MeshCells[R, 2, "Multicells" -> True][[1, 1]]];

MList = NestList[DooSabinSubdivide, M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]

enter image description here

An example with nontrivial boundary:

pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1},(*{1,2,6,5},*){2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};
M = polymesh[pts, polys];

MList = NestList[DooSabinSubdivide, M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]

enter image description here

Implementation

In addition to the code from my other post, one requires the following. This version also features correct(?) handling of boundaries.

getDooSabinSubdivisionMasks = Compile[{{n, _Real}},
   Block[{ω, a},
    ω = 2. Pi/n;
    a = (n + 5.)/(4. n);
    Flatten@Table[If[i == j, a, (3. + 2. Cos[ω (i - j)])/(4. n)], {i, 1, n}, {j, 1, n}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

getDooSabinCombinatorics = Compile[{{face, _Integer, 1}, {idx, _Integer}},
   Flatten[
    Table[
     {i, Compile`GetElement[face, j]},
     {i, idx - Length[face] + 1, idx}, {j, 1, Length[face]}],
    1
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

DooSabinSubdivisionMatrix::usage = "";
polymesh /: DooSabinSubdivisionMatrix[M_polymesh] := 
  Module[{lens, acc, n0, n2, polys, L, A, A12, bndelist, edgesneighpolys, e1, e2, i, j, newbndplist, bndn1}, 
   n0 = MeshCellCount[M, 0];
   n2 = MeshCellCount[M, 2];
   polys = MeshCells[M, 2];
   lens = ToPack[Length /@ polys];
   acc = AccumulateIntegerList[lens];
   L = SparseArray[
     Rule[
      Join @@ getDooSabinCombinatorics[MeshCells[M, 2], Rest[acc]], 
      Join @@ getDooSabinSubdivisionMasks[Range[Max[lens]]][[lens]]], {acc[[-1]], n0}];

   A12 = MeshCellAdjacencyMatrix[M, 1, 2];
   bndelist = Random`Private`PositionsOf[Total[A12, {2}], 1];
   If[Length[bndelist] > 0,
    A = SparseArray @@ {Automatic, {n2, n0}, 0, {1, {acc, Partition[Join @@ polys, 1]}, Join @@ (Most[acc] + Range[lens])}};
    edgesneighpolys = A12[[bndelist]]["AdjacencyLists"];
    {e1, e2} = Transpose[MeshCells[M, 1][[bndelist]]];
    {i, j} = Transpose[Riffle[Transpose[{Flatten[edgesneighpolys], e1}], Transpose[{Flatten[edgesneighpolys], e2}]]];
    newbndplist = cExtractIntegerFromSparseMatrix[A["NonzeroValues"], A["RowPointers"], Flatten[A["ColumnIndices"]], 0, i, j];
    bndn1 = Length[bndelist];
    L[[newbndplist]] = 
     SparseArray[
      Rule[
        Transpose[{Join[Range[2 bndn1], Range[2 bndn1]], Join[Riffle[e1, e2], Riffle[e2, e1]]}], Join[ConstantArray[3./4., 2 bndn1], 
        ConstantArray[1./4., 2 bndn1]]], {2 bndn1, n0}];
    ];
   L
   ];

getDooSabinEdgeQuads = Compile[{{edge, _Integer, 1}, {polyidx, _Integer, 1}},
   {
    {Compile`GetElement[polyidx, 1], Compile`GetElement[edge, 1]},
    {Compile`GetElement[polyidx, 1], Compile`GetElement[edge, 2]},
    {Compile`GetElement[polyidx, 2], Compile`GetElement[edge, 2]},
    {Compile`GetElement[polyidx, 2], Compile`GetElement[edge, 1]}
    },
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

getDooSabinVertexFaces = Compile[{{dualedges, _Integer, 2}, {idx, _Integer, 1}},
   Block[{n, q, p, i, k, c, j},
    n = Length[idx];
    q = dualedges[[idx]];
    p = Table[0, n];
    p[[1]] = Compile`GetElement[q, 1, 1];
    p[[2]] = c = Compile`GetElement[q, 1, 2];
    i = 1;
    k = 2;
    While[k < n,
     j = 1;
     While[
      And[j < n, Or[i == j, c != Compile`GetElement[q, j, 1] && c != Compile`GetElement[q, j, 2]]],
      j++
      ];
     k++;
     i = j;
     p[[k]] = c = If[c == Compile`GetElement[q, j, 1],
        Compile`GetElement[q, j, 2],
        Compile`GetElement[q, j, 1]
        ];
     ];
    p
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

DooSabinSubdividedPolygons::usage = "";
polymesh /: DooSabinSubdividedPolygons[M_polymesh] := 
  Module[{lens, acc, polys, edges, edgesneighpolys, n2, n0, A, B, A12, edgevalences, vertexfaces, edgefaces, facefaces, i, j, t, intelist, bndelist, intplist, intedges},
   n0 = MeshCellCount[M, 0];
   polys = MeshCells[M, 2];
   n2 = Length[polys];
   lens = ToPack[Length /@ polys];
   acc = AccumulateIntegerList[lens];

   facefaces = Internal`PartitionRagged[Range[acc[[-1]]], lens];

   A12 = MeshCellAdjacencyMatrix[M, 1, 2];
   edgevalences = Total[A12, {2}];
   intelist = Random`Private`PositionsOf[edgevalences, 2];
   bndelist = Random`Private`PositionsOf[edgevalences, 1];
   intedges = MeshCells[M, 1][[intelist]];
   intplist = Complement[Range[n0], Flatten[MeshCells[M, 1][[bndelist]]]];
   edgesneighpolys = A12[[intelist]]["AdjacencyLists"];
   {i, j} = 
    Transpose[Join @@ getDooSabinEdgeQuads[intedges, edgesneighpolys]];
   A = SparseArray @@ {Automatic, {n2, n0}, 0, {1, {acc, Partition[Join @@ polys, 1]}, Join @@ (Most[acc] + Range[lens])}};
   edgefaces = Partition[
     cExtractIntegerFromSparseMatrix[
      A["NonzeroValues"],
      A["RowPointers"],
      Flatten[A["ColumnIndices"]], 0, i, j
      ],
     4
     ];

   B = SparseArray[
     Transpose[{Flatten[intedges], Range[2 Length[intelist]]}] -> 1,
     {n0, 2 Length[intelist]}
     ];
   vertexfaces = getDooSabinVertexFaces[
     ArrayReshape[edgefaces[[All, {1, 4, 2, 3}]], {2 Length[edgefaces], 2}],
     B[[intplist]]["AdjacencyLists"]
     ];
   Join[facefaces, edgefaces, vertexfaces]
   ];

DooSabinSubdivide::usage = "";
polymesh /: DooSabinSubdivide[M0_polymesh, OptionsPattern[{"Subdivisions" -> 1}]
   ] := Module[{lens, acc, pat, polys, edges, edgesneighpolys, m, n, nn, A, B, A02, valences, vertexfaces, edgefaces, facefaces, i, j, t, M},
   PrintTemporary["Subdividing..."];
   t = AbsoluteTiming[
      If[OptionValue["Subdivisions"] > 0,
       M = M0;
       M = polymesh[
         DooSabinSubdivisionMatrix[M].MeshCoordinates[M],
         DooSabinSubdividedPolygons[M]
         ];
       ,
       M = M0;
       ]
      ][[1]];
   PrintTemporary["Doo-Sabin subdivision done. Time elapsed: ", ToString[t]];

   If[OptionValue["Subdivisions"] > 1,
    M = DooSabinSubdivide[M, "Subdivisions" -> OptionValue["Subdivisions"] - 1]
    ];
   M
   ];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.