3

In How to obtain the cell-adjacency graph of a mesh? there are some great answers that show how to obtain the cell adjacency matrix of a mesh. How would one generalize these to incorporate the orientations of the cells? That is, I want to compute the boundary operators, $\partial_i:C_i(M)\to C_{i-1}(M)$ on the mesh that show up, for example, if you want to compute the simplicial homology...

I've been pretty stuck even getting the induced orientation of edges and vertices

user51761
  • 255
  • 1
  • 3
  • 2
    On the ref page of ElementMesh in the scope section there are a BoundaryConnectivity and a VertexBoundaryConnectivity section. Perhaps useful. – user21 Jul 16 '19 at 06:34
  • 1
    Do you look for an answer for dimension 0 -- 3 or for general dimension? – Henrik Schumacher Jul 16 '19 at 12:51
  • I have an implementation for quite general simplicial manifolds (or rather for finite-dimensional simplicial complexes that are uniquely defined by their top-dimensional simplices). If you are interested, I could make it stand-alone and post it here. – Henrik Schumacher Jul 16 '19 at 15:04
  • 1
    I'm mostly interested in these low dimensions if there's an easier answer there – Yousuf Soliman Jul 16 '19 at 15:09

1 Answers1

3

Okay here is my code. Unfortunately, it is machine-generated code and also not well documented. But it does the job.

Code dump:

SetAttributes[simplicialmanifold, HoldAll];

$simplicialmanifoldCounter = 0;

ToPack = Developer`ToPackedArray;

simplicialmanifold /: 
  ClearAllCache[x_simplicialmanifold, 
   arg_: All] := (ClearCache[x, arg]; ClearPersistentCache[x, arg]; );

simplicialmanifold /: 
  Cache[x : simplicialmanifold[$object_], 
   args___] := $object[["Cache", args]];

simplicialmanifold /: CacheKeys[x : simplicialmanifold[$object_]] := 
  Keys[$object[["Cache"]]];

simplicialmanifold /: 
  ClearCache[
   x : simplicialmanifold[$object_]] := ($object[["Cache"]] = 
     Association[]; );

simplicialmanifold /: 
  ClearCache[x : simplicialmanifold[$object_], 
   All] := ($object[["Cache"]] = Association[]; );

simplicialmanifold /: 
  ClearCache[x : simplicialmanifold[$object_], s_String] := 
  KeyDropFrom[$object[["Cache"]], s];

simplicialmanifold /: 
  ClearCache[x : simplicialmanifold[$object_], p_] := 
  KeyDropFrom[$object[["Cache"]], ToString[p, InputForm]];

simplicialmanifold /: 
  ClearCache[x : simplicialmanifold[$object_], {s___String}] := 
  KeyDropFrom[$object[["Cache"]], {s}];

simplicialmanifold /: 
  ClearCache[x : simplicialmanifold[$object_], {p__}] := 
  KeyDropFrom[$object[["Cache"]], 
   Function[z, If[StringQ[z], z, ToString[z, InputForm]]] /@ 
    Flatten[{p}]];

simplicialmanifold /: 
  SetCache[x : simplicialmanifold[$object_], pos_, 
   val_] := ($object[["Cache", pos]] = val; );

simplicialmanifold /: 
  PersistentCache[x : simplicialmanifold[$object_], 
   args___] := $object[["PersistentCache", args]];

simplicialmanifold /: 
  PersistentCacheKeys[x : simplicialmanifold[$object_]] := 
  Keys[$object[["PersistentCache"]]];

simplicialmanifold /: 
  ClearPersistentCache[
   x : simplicialmanifold[$object_]] := ($object[["PersistentCache"]] \
= Association[]; );

simplicialmanifold /: 
  ClearPersistentCache[x : simplicialmanifold[$object_], 
   All] := ($object[["PersistentCache"]] = Association[]; );

simplicialmanifold /: 
  ClearPersistentCache[x : simplicialmanifold[$object_], s_String] := 
  KeyDropFrom[$object[["PersistentCache"]], s];

simplicialmanifold /: 
  ClearPersistentCache[x : simplicialmanifold[$object_], p_] := 
  KeyDropFrom[$object[["PersistentCache"]], ToString[p, InputForm]];

simplicialmanifold /: 
  ClearPersistentCache[
   x : simplicialmanifold[$object_], {s___String}] := 
  KeyDropFrom[$object[["PersistentCache"]], {s}];

simplicialmanifold /: 
  ClearPersistentCache[x : simplicialmanifold[$object_], {p__}] := 
  KeyDropFrom[$object[["PersistentCache"]], 
   Function[z, If[StringQ[z], z, ToString[z, InputForm]]] /@ 
    Flatten[{p}]];

simplicialmanifold /: 
  ClearPersistentCacheDependingOn[x_simplicialmanifold, s_] := 
  ClearPersistentCache[x, 
   VertexInComponent[CallGraph[$PM], ToString /@ Flatten[{s}]]];

simplicialmanifold /: 
  SetPersistentCache[x : simplicialmanifold[$object_], pos_, 
   val_] := ($object[["PersistentCache", pos]] = val; );

simplicialmanifold /: DeepCopy[x$ : simplicialmanifold[$nobject$_]] := 
  Module[{$object$ = DeepCopy /@ $nobject$}, 
   simplicialmanifold[$object$]];

simplicialmanifold /: (y_) \[LeftArrow] (x : 
     simplicialmanifold[$nobject_]) := y = DeepCopy[x];

simplicialmanifold /: (x : simplicialmanifold[$object_])[[1, 
    args__]] := $object[[args]];

simplicialmanifold /: (x : simplicialmanifold[$object_])[[s_String, 
    args___]] := $object[[s, args]];

simplicialmanifold /: (x_simplicialmanifold)[s_String, args___] := 
  x[[1]][[s, args]];

simplicialmanifold /: 
  Initialize[simplicialmanifold, data0$_Association] := 
  With[{c$ = ++$simplicialmanifoldCounter},
    ToExpression[
    StringJoin["$", "simplicialmanifold", "$", ToString[c$]],
     InputForm, 
    Function[data$, SetAttributes[data$, Temporary]; data$ = data0$; 
     If[ ! KeyExistsQ[data$, "Cache"], 
      AppendTo[data$, "Cache" -> Association[]]]; 
     If[ ! KeyExistsQ[data$, "PersistentCache"], 
      AppendTo[data$, "PersistentCache" -> Association[]]]; 
     If[ ! KeyExistsQ[data$, "Settings"], 
      AppendTo[data$, "Settings" -> Association[]]]; 
     simplicialmanifold[data$], HoldAll]
    ]
   ];

simplicialmanifold /: Equal[a__simplicialmanifold] := 
  Equal @@ KeyDrop[
    Flatten[{a}][[All, 1]], {"Cache", "PersistentCache"}];

simplicialmanifold /: Identifier[x : simplicialmanifold[$object_]] := 
  ToString[Unevaluated[$object]];

simplicialmanifold /: IDNumber[x : simplicialmanifold[$object_]] := 
  With[{s = ToString[Unevaluated[$object]]}, 
   ToExpression[
    StringTake[s, {StringPosition[s, "$"][[-1, -1]] + 1, -1}]]];

simplicialmanifold /: ObjectQ[simplicialmanifold] := True;

simplicialmanifold /: 
  Settings[x : simplicialmanifold[$object_], 
   args___] := $object[["Settings", args]];

simplicialmanifold /: SettingsKeys[x : simplicialmanifold[$object_]] := 
  Keys[$object[["Settings"]]];

simplicialmanifold /: 
  ClearSettings[
   x : simplicialmanifold[$object_]] := ($object[["Settings"]] = 
     Association[]; );

simplicialmanifold /: 
  ClearSettings[x : simplicialmanifold[$object_], 
   All] := ($object[["Settings"]] = Association[]; );

simplicialmanifold /: 
  ClearSettings[x : simplicialmanifold[$object_], s_String] := 
  KeyDropFrom[$object[["Settings"]], s];

simplicialmanifold /: 
  ClearSettings[x : simplicialmanifold[$object_], p_] := 
  KeyDropFrom[$object[["Settings"]], ToString[p, InputForm]];

simplicialmanifold /: 
  ClearSettings[x : simplicialmanifold[$object_], {s___String}] := 
  KeyDropFrom[$object[["Settings"]], {s}];

simplicialmanifold /: 
  ClearSettings[x : simplicialmanifold[$object_], {p__}] := 
  KeyDropFrom[$object[["Settings"]], 
   Function[z, If[StringQ[z], z, ToString[z, InputForm]]] /@ 
    Flatten[{p}]];

simplicialmanifold /: 
  SetSettings[x : simplicialmanifold[$object_], pos_, 
   val_] := ($object[["Settings", pos]] = val; );

simplicialmanifold /: 
  Init[simplicialmanifold, points0_, simplices0_] := 
  Initialize[simplicialmanifold, 
   Association["VertexCoordinates" -> ToPack[N[points0]], 
    "TopSimplices" -> ToPack[simplices0], 
    "Dimension" -> Length[points0[[1]]]]];

getSimplexFaces = Compile[{{simplex, _Integer, 1}},
   Table[Delete[simplex, i], {i, 1, Length[simplex]}],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

SortThread = Compile[{{list, _Integer, 1}},
   Sort[list],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

simplicialmanifold /: 
  VertexCoordinates[M : simplicialmanifold[$simplicialmanifold_], 
   args___] := M[[1]][["VertexCoordinates", args]];

simplicialmanifold /: 
  DofCount[M : simplicialmanifold[$simplicialmanifold_]] := If[
   KeyExistsQ[$simplicialmanifold[["Cache"]], "DofCount"], 
   $simplicialmanifold[["Cache", "DofCount"]], 
   $simplicialmanifold[["Cache", "DofCount"]] = 
    Times @@ Dimensions[M[[1]][["VertexCoordinates"]]]
   ];

simplicialmanifold /: 
  VertexCount[M : simplicialmanifold[$simplicialmanifold_]] := 
  If[KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
    "VertexCount"], $simplicialmanifold[["PersistentCache", 
     "VertexCount"]], $simplicialmanifold[["PersistentCache", 
      "VertexCount"]] = Length[M[[1]][["VertexCoordinates"]]]];

simplicialmanifold /: 
  TopSimplices[M : simplicialmanifold[$simplicialmanifold_], 
   args___] := M[[1]][["TopSimplices", args]];

simplicialmanifold /: 
  TopSimplexCount[M : simplicialmanifold[$simplicialmanifold_]] := 
  If[KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
    "TopSimplexCount"], 
   $simplicialmanifold[["PersistentCache", "TopSimplexCount"]], 
   $simplicialmanifold[["PersistentCache", "TopSimplexCount"]] = 
    Length[M[[1]][["TopSimplices"]]]
   ];

simplicialmanifold /: 
  TopSimplexData[M : simplicialmanifold[$simplicialmanifold_]] := 
  If[KeyExistsQ[$simplicialmanifold[["Cache"]], 
    "TopSimplexData"], $simplicialmanifold[["Cache", 
     "TopSimplexData"]], $simplicialmanifold[["Cache", 
      "TopSimplexData"]] = 
    Partition[VertexCoordinates[M][[Flatten[TopSimplices[M]]]], 
     IntrinsicDimension[M] + 1]
   ];

simplicialmanifold /: 
  IntrinsicDimension[M : simplicialmanifold[$simplicialmanifold_]] := 
  If[
   KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
    "IntrinsicDimension"], $simplicialmanifold[["PersistentCache", 
     "IntrinsicDimension"]],
   $simplicialmanifold[["PersistentCache", "IntrinsicDimension"]] = 
    Dimensions[TopSimplices[M]][[2]] - 1
   ];

simplicialmanifold /: 
  SimplexCombinatorics[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (
   If[ ! KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
      "SimplexCombinatorics"], $simplicialmanifold[["PersistentCache",
        "SimplexCombinatorics"]] = Association[]
    ]; 
   If[
    AssociationQ[$simplicialmanifold[["PersistentCache", 
       "SimplexCombinatorics"]]],
    If[KeyExistsQ[
      $simplicialmanifold[["PersistentCache", 
        "SimplexCombinatorics"]], {(First @* List)[dim, _Integer]}
      ],
     $simplicialmanifold[["PersistentCache", "SimplexCombinatorics", 
       Key[{(First @* List)[dim, _Integer]}]]]
     , $simplicialmanifold[["PersistentCache", "SimplexCombinatorics",
         Key[{(First @* List)[dim, _Integer]}]]] = 
      Module[{simplices, idx, faces, sortedfaces}, 
       If[TrueQ[dim == IntrinsicDimension[M]], 
        With[{topSimplices = SortThread[TopSimplices[M]]}, 
         AssociationThread[topSimplices, 
          Range[Length[topSimplices]]]], 
        If[Inequality[0, LessEqual, dim, Less, 
          IntrinsicDimension[M]], 
         idx = Subsets[Range[dim + 2], {dim + 1}]; 
         faces = Partition[
           Flatten[
            Keys[SimplexCombinatorics[M, dim + 1]][[All, 
              Flatten[idx]]]], dim + 1]; 
         If[dim == IntrinsicDimension[M] - 1, 
          faces = SortThread[faces]]; 
         sortedfaces = DeleteDuplicates[Sort[faces]]; 
         AssociationThread[sortedfaces, Range[Length[sortedfaces]]], 
         Print["Error"]; $Failed]]
       ]
     ], 
    Print[
     StringJoin["Call to ", "SimplexCombinatorics", 
      " with arguments ", {"dim_Integer"}, " : ", "PersistentCache", 
      " with key ", "SimplexCombinatorics", 
      " is not an Association. Caching ignored."]]; 
    Module[{simplices, idx, faces, sortedfaces}, 
     If[TrueQ[dim == IntrinsicDimension[M]], 
      With[{topSimplices = SortThread[TopSimplices[M]]}, 
       AssociationThread[topSimplices, Range[Length[topSimplices]]]], 
      If[Inequality[0, LessEqual, dim, Less, IntrinsicDimension[M]], 
       idx = Subsets[Range[dim + 2], {dim + 1}]; 
       faces = Partition[
         Flatten[Keys[SimplexCombinatorics[M, dim + 1]][[All, 
            Flatten[idx]]]], dim + 1]; 
       If[dim == IntrinsicDimension[M] - 1, 
        faces = SortThread[faces]]; 
       sortedfaces = DeleteDuplicates[Sort[faces]]; 
       AssociationThread[sortedfaces, Range[Length[sortedfaces]]], 
       Print["Error"]; $Failed]]]]
   );

simplicialmanifold /: 
  Simplices[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (If[ ! 
     KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
      "Simplices"], $simplicialmanifold[["PersistentCache", 
       "Simplices"]] = Association[]]; 
   If[AssociationQ[$simplicialmanifold[["PersistentCache", 
       "Simplices"]]], 
    If[KeyExistsQ[$simplicialmanifold[["PersistentCache", 
        "Simplices"]], {(First @* List)[
        dim, _Integer]}], $simplicialmanifold[["PersistentCache", 
       "Simplices", 
       Key[{(First @* List)[
          dim, _Integer]}]]], $simplicialmanifold[["PersistentCache", 
        "Simplices", Key[{(First @* List)[dim, _Integer]}]]] = 
      If[dim == IntrinsicDimension[M], TopSimplices[M], 
       Keys[SimplexCombinatorics[M, dim]]]], 
    Print[StringJoin["Call to ", "Simplices", 
      " with arguments ", {"dim_Integer"}, " : ", "PersistentCache", 
      " with key ", "Simplices", 
      " is not an Association. Caching ignored."]]; 
    If[dim == IntrinsicDimension[M], TopSimplices[M], 
     Keys[SimplexCombinatorics[M, dim]]]]);

simplicialmanifold /: 
  SimplexData[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (Partition[
    VertexCoordinates[M][[Flatten[Simplices[M, dim]]]], dim + 1]);

simplicialmanifold /: 
  SimplexCount[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (If[ ! 
     KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
      "SimplexCount"], $simplicialmanifold[["PersistentCache", 
       "SimplexCount"]] = Association[]]; 
   If[AssociationQ[$simplicialmanifold[["PersistentCache", 
       "SimplexCount"]]], 
    If[KeyExistsQ[$simplicialmanifold[["PersistentCache", 
        "SimplexCount"]], {(First @* List)[
        dim, _Integer]}], $simplicialmanifold[["PersistentCache", 
       "SimplexCount", 
       Key[{(First @* List)[
          dim, _Integer]}]]], $simplicialmanifold[["PersistentCache", 
        "SimplexCount", Key[{(First @* List)[dim, _Integer]}]]] = 
      Length[Keys[SimplexCombinatorics[M, dim]]]], 
    Print[StringJoin["Call to ", "SimplexCount", 
      " with arguments ", {"dim_Integer"}, " : ", "PersistentCache", 
      " with key ", "SimplexCount", 
      " is not an Association. Caching ignored."]]; 
    Length[Keys[SimplexCombinatorics[M, dim]]]]);

simplicialmanifold /: 
  ToTopSimplexMatrix[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (If[ ! 
     KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
      "ToTopSimplexMatrix"], $simplicialmanifold[["PersistentCache", 
       "ToTopSimplexMatrix"]] = Association[]]; 
   If[AssociationQ[$simplicialmanifold[["PersistentCache", 
       "ToTopSimplexMatrix"]]], 
    If[KeyExistsQ[$simplicialmanifold[["PersistentCache", 
        "ToTopSimplexMatrix"]], {(First @* List)[
        dim, _Integer]}], $simplicialmanifold[["PersistentCache", 
       "ToTopSimplexMatrix", 
       Key[{(First @* List)[
          dim, _Integer]}]]], $simplicialmanifold[["PersistentCache", 
        "ToTopSimplexMatrix", 
        Key[{(First @* List)[dim, _Integer]}]]] = 
      Module[{simplices, faces, m, idx}, simplices = TopSimplices[M]; 
       idx = Subsets[Range[IntrinsicDimension[M] + 1], {dim + 1}]; 
       faces = Partition[Flatten[simplices[[All, Flatten[idx]]]], 
         dim + 1]; m = Binomial[IntrinsicDimension[M] + 1, dim + 1]; 
       SparseArray[
        Transpose[{Range[m*Length[simplices]], 
           SimplexLookup[M, faces]}] -> 
         Signature /@ faces, {m*Length[simplices], 
         SimplexCount[M, dim]}]]], 
    Print[StringJoin["Call to ", "ToTopSimplexMatrix", 
      " with arguments ", {"dim_Integer"}, " : ", "PersistentCache", 
      " with key ", "ToTopSimplexMatrix", 
      " is not an Association. Caching ignored."]]; 
    Module[{simplices, faces, m, idx}, simplices = TopSimplices[M]; 
     idx = Subsets[Range[IntrinsicDimension[M] + 1], {dim + 1}]; 
     faces = Partition[Flatten[simplices[[All, Flatten[idx]]]], 
       dim + 1]; m = Binomial[IntrinsicDimension[M] + 1, dim + 1]; 
     SparseArray[
      Transpose[{Range[m*Length[simplices]], 
         SimplexLookup[M, faces]}] -> 
       Signature /@ faces, {m*Length[simplices], 
       SimplexCount[M, dim]}]]]);

simplicialmanifold /: 
  SimplexLookup[
   M : simplicialmanifold[$simplicialmanifold_], (simplices_)?
    MatrixQ] := 
  Lookup[SimplexCombinatorics[M, Dimensions[simplices][[2]] - 1], 
   SortThread[simplices]];

simplicialmanifold /: 
  BoundaryOperator[M : simplicialmanifold[$simplicialmanifold_], 
   dim_Integer] := (If[ ! 
     KeyExistsQ[$simplicialmanifold[["PersistentCache"]], 
      "BoundaryOperator"], $simplicialmanifold[["PersistentCache", 
       "BoundaryOperator"]] = Association[]]; 
   If[AssociationQ[$simplicialmanifold[["PersistentCache", 
       "BoundaryOperator"]]], 
    If[KeyExistsQ[$simplicialmanifold[["PersistentCache", 
        "BoundaryOperator"]], {(First @* List)[
        dim, _Integer]}], $simplicialmanifold[["PersistentCache", 
       "BoundaryOperator", 
       Key[{(First @* List)[
          dim, _Integer]}]]], $simplicialmanifold[["PersistentCache", 
        "BoundaryOperator", Key[{(First @* List)[dim, _Integer]}]]] = 
      If[1 <= dim <= IntrinsicDimension[M], 
       Module[{simplices, idx, faces, m}, 
        simplices = Simplices[M, dim]; m = Length[simplices]; 
        idx = (Delete[Range[dim + 1], #1] & ) /@ Range[dim + 1]; 
        faces = Partition[Flatten[simplices[[All, Flatten[idx]]]], 
          dim]; SparseArray[
         Transpose[{SimplexLookup[M, faces], 
            Flatten[
             Transpose[ConstantArray[Range[m], Length[idx]]]]}] -> 
          If[dim == IntrinsicDimension[M], 
           Flatten[ConstantArray[(-1)^Range[0, dim], m]]*
            Signature /@ faces, 
           Flatten[
            ConstantArray[(-1)^Range[0, dim], m]]], {SimplexCount[M, 
           dim - 1], m}]], {}]], 
    Print[StringJoin["Call to ", "BoundaryOperator", 
      " with arguments ", {"dim_Integer"}, " : ", "PersistentCache", 
      " with key ", "BoundaryOperator", 
      " is not an Association. Caching ignored."]]; 
    If[1 <= dim <= IntrinsicDimension[M], 
     Module[{simplices, idx, faces, m}, 
      simplices = Simplices[M, dim]; m = Length[simplices]; 
      idx = (Delete[Range[dim + 1], #1] & ) /@ Range[dim + 1]; 
      faces = Partition[Flatten[simplices[[All, Flatten[idx]]]], 
        dim]; SparseArray[
       Transpose[{SimplexLookup[M, faces], 
          Flatten[Transpose[ConstantArray[Range[m], Length[idx]]]]}] -> 
        If[dim == IntrinsicDimension[M], 
         Flatten[ConstantArray[(-1)^Range[0, dim], m]]*
          Signature /@ faces, 
         Flatten[ConstantArray[(-1)^Range[0, dim], 
           m]]], {SimplexCount[M, dim - 1], m}]], {}]]);

simplicialmanifold /: 
  BoundaryOperators[M : simplicialmanifold[$simplicialmanifold_]] := 
  Association[
   Table[i -> BoundaryOperator[M, i], {i, 1, IntrinsicDimension[M]}]];

ToSimplicialManifold[M_MeshRegion] := 
  Init[simplicialmanifold, MeshCoordinates[M], 
   MeshCells[M, RegionDimension[M], "Multicells" -> True][[1, 1]]];

Here is a short usage example:

R = DiscretizeRegion[Ball[]];
M = ToSimplicialManifold[R];
A = BoundaryOperators[M]

What is returned is an association A containing the boundary operators, where A[i] is a matrix representation of $\partial_i \colon C_i(M) \to C_{i-1}(M)$.

This list of $i$-dimensional oriented simplices that were used as basis can be obtained with

Simplices[M, i]
Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • Hm, I am curious, what does it mean "machine-generated" code? – Pinti Jul 17 '19 at 04:41
  • 1
    @Pinti There is a lot of boiler-plate code, e.g., for the management of caches (the simplicialmanifold data type is able to remember computed results; this speeds up execution for data that is frequently needed). So I wrote tools that generate the final code from a few lines of code that contain the actual behavior. But unfortunately, this code is hradly readable, in particular because it is not formatted... – Henrik Schumacher Jul 17 '19 at 12:03