# Signed cell adjacency matrices (boundary operators)

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

• On the ref page of ElementMesh in the scope section there are a BoundaryConnectivity and a VertexBoundaryConnectivity section. Perhaps useful. Jul 16 '19 at 6:34
• Do you look for an answer for dimension 0 -- 3 or for general dimension? 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. Jul 16 '19 at 15:04
• I'm mostly interested in these low dimensions if there's an easier answer there Jul 16 '19 at 15:09

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 = DeveloperToPackedArray; 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]

• Hm, I am curious, what does it mean "machine-generated" code? Jul 17 '19 at 4:41
• @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... Jul 17 '19 at 12:03