# Assignment 1 (Coding): Discrete Exterior Calculus

For the coding portion of your first assignment, you will implement the discrete exterior calculus (DEC) operators $\star_0$, $\star_1$, $\star_2$, $d_0$ and $d_1$. Once implemented, you will be able to apply these operators to a scalar function (as depicted above) by pressing the “*” and “d” button in the viewer. The diagram shown above will be updated to indicate what kind of differential k-form is currently displayed. These basic operations will be the starting point for many of the algorithms we will implement throughout the rest of the class; the visualization (and implementation!) should help you build further intuition about what these operators mean and how they work.

Getting Started

• Please clone this repository. It contains a fast and flexible framework for 3D geometry processing implemented in Javascript. Over the course of the semester, you will implement all of your coding assignments here. Please note: If you already cloned the repository during recitation, clone again!
• For this assignment, you need to implement the following routines:
1. in core/geometry.js
1. cotan
2. barycentricDualArea
2. in core/discrete-exterior-calculus.js
1. buildHodgeStar0Form
2. buildHodgeStar1Form
3. buildHodgeStar2Form
4. buildExteriorDerivative0Form
5. buildExteriorDerivative1Form

In practice, a simple and efficient way to compute the cotangent of the angle $\theta$ between two vectors $u$ and $v$ is to use the cross product and the dot product rather than calling any trigonometric functions directly; we ask that you implement your solution this way. (Hint: how are the dot and cross product of two vectors related to the cosine and sine of the angle between them?)

In case we have not yet covered it in class, the barycentric dual area associated with a vertex $i$ is equal to one-third the area of all triangles $ijk$ touching $i$.

The discrete Hodge star and discrete exterior derivatives are introduced in Section 3.8 of the course notes; the matrix representation of these operators (which you need to implement!) will be discussed in class. They were also basically covered already in our discussion of signed incidence matrices, in the lecture on the simplicial complex.

Notes

• This assignment comes with a viewer (projects/discrete-exterior-calculus/index.html) which lets you apply your operators on random k-forms and visualize the results.
• This assignment also comes with a grading script (tests/discrete-exterior-calculus/test.html) which you can use to verify the correctness of your operators.
• The code framework is implemented in Javascript, which means no compilation or installation is necessary on any platform. You can simply get started by opening the index.html file in projects/discrete-exterior-calculus/ in a web browser. We recommend using Chrome or Firefox. Safari has poor WebGL performance.
• If you do not have prior experience with Javascript, do not worry! You should be able to get a handle on Javascript syntax by reading through some of the code in the framework (a good place to start might be core/geometry.js). The framework also contains extensive documentation (see docs/index.html) with examples on how to use the halfedge data structure and the linear algebra classes.
• All browsers come with tools for debugging (for instance the JavaScript Console in Chrome).

Submission Instructions

Please rename your geometry.js and discrete-exterior-calculus.js files to geometry.txt and discrete-exterior-calculus.txt (respectively) and submit them in a single zip file called solution.zip by email to Geometry.Collective@gmail.com.

## 68 thoughts on “Assignment 1 (Coding): Discrete Exterior Calculus”

1. kye says:

Are both parts of Assignment 1 due at the same time (5:59 pm on 9/27)?

2. yongOne says:

We can always assume we are working with a triangular mesh? And that it is always properly connected (or on a boundary)?

1. Keenan says:

Yes, you can assume you always have a triangle mesh (because otherwise things like the Hodge star do not have an immediate definition) and you can also assume the mesh is manifold, possibly with boundary (since that’s all a halfedge data structure can represent). In real code, it’s also a great practice to check whether the input satisfies the requirements of the algorithm, e.g., whether all faces have three sides. We don’t require that you make this check, but you might want to try it just for your own edification!

3. peiya_wang says:

How do I test to see if my cotan function is correct? The test.html doesn’t include information on whether this function is correct or not. Do I need to manually create a random triangle somehow and give cotan a half edge?

1. rohansawhney says:

Yes, the test.html files in tests/geometry and tests/discrete-exterior-calculus do not check the correctness the cotan function directly. However, you will have to call the cotan function when you are constructing one of the DEC operators. You can check this operator’s correctness (and hence the cotan function’s) with the test.html file in tests/discrete-exterior-calculus.

1. apoorva says:

We’re supposed to build all operators for general triangle meshes, right? Not necessarily equilateral triangles?

1. rohansawhney says:

Yes, don’t make any assumptions about the triangles in the mesh.

4. plav says:

I am trying to finish the part of the code, but I have few questions:

1) Shouldn’t we implement circumscentricDualArea in order to implement buildHodgeStar0Form?
2) Is the following interpretation correct? The indices we use to access the elements in the array of geometry.mesh are different from the id indices of the elements (represented in the dictionaries, vertexIndex, edgeIndex and faceIndex). The positions in the matrices for the hodge star are related to the id indices of the vertices, edges and faces.

1. rohansawhney says:

1) We use barycentric dual areas in our implementation because they are always positive, irrespective of the mesh they are being calculated on. This isn’t generally true for circumcentric dual areas.
2) Yes, assume that the indices of the elements in the mesh lists (such as vertices, edges, faces etc) are different from the indices stored as values in the vertexIndex, edgeIndex, faceIndex etc. dictionaries. You should construct your matrices using the latter.

1. plav says:

1) I think my hodge star functions are running, but somehow returning the wrong result as I get the following output: “AssertionError: expected false to equal true”. Is there a test case, that I can use to check my result?

2) In the buildHodgeStar2Form I also got an error “TypeError: f.isBoundaryLoop is not a function”. Maybe the object f I am passing to geometry.area(f) might not be a face, so it does not have a method f.isBoundaryLoop(). However, f is an object from the array geometry.mesh.faces, so it should be a face.

(I assume that I can comment about these topics here, but if I am wrong, tell me and delete my comment. In this case I will send the questions by email or I will present in the OH).

1. rohansawhney says:

1) You should see a green checkmark against the functions that have been implemented correctly, otherwise you’ll see an error in red like the one you mentioned. Other than the test file, you could visualize the result of applying your operators on differential k forms with the index.html file in projects/discrete-exterior-calculus. Ask yourself what you should see if say you apply the discrete hodge 0 star to a differential 0 form, or a discrete hodge 1 star to a differential 1 form?

2) If you are iterating over objects in geometry.mesh.faces, then you should not get that error. Might be best to discuss this one in office hours.

1. plav says:

I am implementing the functions and I was able to pass the tests for the buildHodgeStar0Form. For the buildHodgeStar1Form, I used the following methods to calculate the length of the dual:
-distances between centroids and midpoint
-distances between circumcenters and midpoint

Both failed. I checked in the slides that it is possible to use the opposite angles to calculate it. Our data structures support this except that the function Geometry.angle() is not implemented.

Could you give me a hint about the right method to calculate the length of the dual or the whole expression for buildHodgeStar1Form?

2. rohansawhney says:

Take at look at the slide titled “Computing Volumes” in Discrete Exterior Calculus slides!

3. plav says:

But in this case, shouldn’t I have a function angle implemented?

4. rohansawhney says:

Try implementing it in the way its described in the assignment writeup:

“In practice, a simple and efficient way to compute the cotangent of the angle θ between two vectors u and v is to use the cross product and the dot product rather than calling any trigonometric functions directly; we ask that you implement your solution this way. (Hint: how are the dot and cross product of two vectors related to the cosine and sine of the angle between them?)”

You could also implement it using the angle function (which you’ll have to implement on the next assignment in any case :)).

5. plav says:

You are right, I have done all that and used the equations you cited. After studying the last topics, I completely forgot them. Thanks argain, Rohan.

6. plav says:

My function for discrete hodge 1 is not working at all. I am using precisely the formula indicated. It might be something with my cotan function. However, I tried everything and I talked with Josh yesterday… and I found no solution here. Do you have any hint or test case?
Thanks again.

7. rohansawhney says:

Does your cotan function return 0 when the input halfedge lies on the boundary? If you have the exterior derivative 0 form (i.e., d0) working, check visually whether your discrete hodge 1 star implementation rotates the vector field by 90 degrees.

8. plav says:

Yes it is returning 0 for half-edges on the boundary.
I will implement the derivative now to try to visualize.

9. plav says:

1) My functions for hodge 1 and derivatives do not work at all… and I checked the logic manually

2) Can I assume that face.adjacentEdges() and face.adjacentHalfedges() return objects in the same position (i.e, the first halfedge will be the halfedge or the twin of the first edge and so on)?

3) Is there any specification for numbers in the discrete calculus functions? Ex: should 0 be a float or an integer?

10. rohansawhney says:

2) No, you shouldn’t make that assumption.

3) Using 1 vs 1.0 or 0 vs 0.0 shouldn’t make a difference in the test script.

11. plav says:

One last thing. I am assuming that the orientation of the edge is equivalent to the half edge it references in its attributes. Is this correct?

12. rohansawhney says:

You don’t have to, but unfortunately the test script makes that assumption.

13. plav says:

“2) No, you shouldn’t make that assumption.”

So the Class face does not return edges, vertices (and halfedges aligned with the vertices) starting from the same relative position and order? I feel that this should be the necessary in order to compare orientation.

14. plav says:

… figured out another method.

15. rohansawhney says:

The orientation of a face on a halfedge mesh is determined only by the orientation of the halfedges in each face. It doesn’t matter which halfedge face.adjacentHalfedges() returns first, or more generally in what order the halfedges are returned (this is an implementation detail about which no assumptions should be made), the choice is arbitrary. Given a halfedge h in the interior of a triangle mesh, the only conditions that need to hold are that h.next.next.next == h and h.twin.twin == h. Similarly, the order in which vertices and edges are returned by face.adjacentVertices() and face.adjacentEdges() depends on the implementation, there is no hard condition on what order they should be returned in.

5. ahagerty says:

I guess I’m a little confused by our implementation of the Discrete Exterior Calculus. We want, for instance, the Hodge star on 0 forms to be encoded by a matrix- What would each entry of this matrix represent? Is this the matrix which acts on a 0-form (represented as a vector of size Number-Of-Vertices) and spits out a 2-form (also represented as a vector of size Number-Of-Vertices)? Or am I misinterpreting something?

1. Keenan says:

Sounds like you’re on the right track—we’ll cover all this stuff in class on Tue/Thu. (Just in case you didn’t get the announcement, the final due date for the homework is now FRIDAY :-))

Pretty much everything you need for the assignment is also spelled out in the slides on discrete exterior calculus. Look especially for slides about (i) assigning indices to mesh elements, and (ii) representing the discrete exterior derivative and discrete Hodge star as matrices.

2. rohansawhney says:

Thats correct, the discrete hodge star matrix applied to a primal 0-form (of size Number-Of-Vertices) on the primal mesh spits out a dual 2 form (also of size Number-Of-Vertices) on the dual mesh. Take a look at the formal definition of the discrete hodge star in section 3.3.4 of the written homework, that should give you an idea of what each entry of this matrix should be.

6. aishwaryapawar says:

In the computation of the Discrete Exterior Derivative in the function buildExteriorDerivative1Form(geometry, faceIndex, edgeIndex), is there a way to check if the relative orientation is taken into account accurately? I am getting a bit confused in assigning orientations for the face and the edges. Thank you!

1. rohansawhney says:

If your orientations are correct, then d1 d0 should equal a |F| x |V| zero matrix. I suggest taking a look at the slide titled “Discrete Exterior Derivative d1—Example” in the Discrete Exterior Calculus slides. It builds d1 on a mesh with only two faces that share a common edge. That might help clear things up.

7. peiya_wang says:

I don’t believe gmail will let you send a .zip file? How should I get around this?

1. yangxio says:

Now that js files are blocked, would you mind if we include our own github repository link as a way of submission? Thanks.

1. rohansawhney says:

We’d like you to rename your geometry.js and discrete-exterior-calculus.js files to geometry.txt and discrete-exterior-calculus.txt (respectively) and submit them in a single zip file. This makes it easier for us to batch download all the assignments.

2. rohansawhney says:

We’d like you to rename your geometry.js and discrete-exterior-calculus.js files to geometry.txt and discrete-exterior-calculus.txt (respectively) and submit them in a single zip file. This makes it easier for us to batch download all the assignments.

8. apoorva says:

I’m not sure if I’m reading something wrong, but it seems like the last line of the formal definitions of the discrete hodge star(s) doesn’t typecheck. It says Area(f*), and f* is a dual vertex. What’s the area of a vertex?

1. Keenan says:

It’s 1, always. 🙂

(Unless you’re asking about the area of the dual cell associated with a vertex, which is different.)

1. apoorva says:

The definition as given reads $(\star_2\alpha_2)(f^*) = \frac{1}{\text{Area}(f^*)}\alpha_2(f)$. If the area of a vertex is always 1, this would be $(\star_2\alpha_2)(f^*) = \alpha_2(f)$, but that doesn’t feel right… Should it be $\frac{1}{\text{Area}(f)}\alpha_2(f)$?

1. Keenan says:

In 2D, to go from a primal 2-form (on triangles) to a dual 0-form (on triangle centers), you should divide by triangle area and multiply by 1, i.e., divide by the primal volume and multiply by the dual volume. Likewise, to go from a primal 0-form to a dual 2-form, you should divide by 1 and then multiply by the dual cell area.

2. Josh says:

It looks like you are looking at the definition from an old version of the homework writeup. The latest version “V9” has the typo fixed.

1. apoorva says:

I see. Thanks!

9. Yifanh says:

Maybe I am too ignorant about Javascript, but is there any way to set/get an element of a SparseMatrix? For a SparseMatrix m , both m and m.data do not have set() or get(); m[0,0], m(0,0), m.data(0,0), m.data[0,0] are all undefined.

1. apoorva says:

I think there isn’t a way to directly modify entries of a SparseMatrix. We have to prepare the data separately (either as a DenseMatrix representing the diagonal of our matrix or a list of Triplets specifying the entries of the matrix), and pass it into one of the SparseMatrix constructors.

1. rohansawhney says:

Yes, sparse matrices are generally constructed “all at once” through the Triplet object. This is a consequence of the data structure (Compressed Sparse Row/Column format) used internally to represent sparse matrices, which allows for efficient multiplication and addition between sparse matrices, but makes access expensive. Setting and getting are O(n) operations, where n is the number of entries in the matrix.

As for the brackets, operator overloading is not supported in Javascript, which means that its not possible to use [] or () with matrix classes. This also means that we can’t use operators such as +, -, *, /, +=, -=, *=, /*, % etc on matrices (or on any Javascript object other than primitives such as number and string). As a result, we are stuck with get, set, plus, minus, multiply, increment, decrement etc 🙁

1. Yifanh says:

Great explanation, thanks guys. I guess I need to do more real coding instead of being spoiled by Matlab…

10. Yifanh says:

I think there is a specific rule used in the test for determining the orientation of edges/faces. Just being consistent with myself doesn’t work. If I use “edges always points to vertex with higher index”, it doesn’t work. If I use “edge direction agrees with the halfedge property in it” then it works. It might be a good idea for the TAs to make a clearance on this.

1. rohansawhney says:

Yes, you are right. Good catch. I forgot to take this into consideration when designing the test script. I apologize. Unfortunately, modifying the test script to track arbitrary user defined orientations is going to be difficult at this point. There will be a manual check of anyone’s discrete exterior derivative 1 form implementation if it fails in the test script.

11. ahagerty says:

I am having trouble with some practical aspects of the coding… I am not sure how to manipulate dictionaries in Javascript, and google isn’t being my friend.

Is the dictionary something that takes in a vertex and returns and index or is it something that takes in an index and returns a vertex? How do I use it to map back and forth?

1. yangxio says:

the “vertexIndex” eats a vertex in geometry.mesh.vertices and spits out an indexing number. should use “dictionary[key]” to access the “value”.

there is actually documentations in the entire DDG js folder. just look for those auto-generated html files. most have sample codes.

2. rohansawhney says:

Dictionaries in the DEC class take a mesh element such as a vertex, edge or face as input and return an index that specifies its position in a matrix. So for example, if you are trying to access the index of a vertex v in the matrix, you would do so in code as follows:

let index = vertexIndex[v];

12. intrepidowl says:

Is there a good way to print output for debugging? Simply using console.log doesn’t seem to work in the test cases when I’m failing.

1. rohansawhney says:

Depends, what are you trying to print?

2. kye says:

Make sure your Chrome inspector has debug output turned on. “Default levels > Info”

13. Leslie says:

1. The halfedge in our code is defined with its own orientation or we need to set for each one? (which means we need to get all the edges, use command e.halfedge?)
2. I am not sure it is correct or not. Let an edge $e$ = $|ab|$, the orientation of halfedge is defined from point $a$ to point $b$, when we use command h = e.halfedge, $h$ is $|ab|$ with orientation from a to b. If we use command h.vertex, it will return vertex $a$ or $b$?
Thanks!

1. rohansawhney says:

Each halfedge comes with a predefined orientation, from the vertex at its base (h.vertex) to the vertex it points to (h.next.vertex or h.twin.vertex).

1. rohansawhney says:

14. kye says:

Can we pick all faces to be oriented one way (e.g. CCW)?

1. kye says:

Or, in general, pick whatever orientations we want for the faces?

1. rohansawhney says:

The orientation of all the faces is CCW and is determined by the orientation of its halfedges. The latter is set automatically when the mesh is loaded.

1. Leslie says:

Does this mean f.halfedge and e.halfedge (for same edge) can have different orientations? And what is f.halfedge? Thanks!

1. Josh says:

Yes, f.halfedge and e.halfedge do not have the same orientation always.

f.halfedge is one of the halfedges with the same orientation as the face (this could also been answered by checking the “docs” folder in your repository)

2. rohansawhney says:

Assuming f.halfedge.edge and e are the same edge, its possible that f.halfedge and e.halfedge are the same halfedge (i.e., they represent the same object in memory), but this is not guaranteed. e.halfedge can be any one of the two halfedges incident on the edge e. So, either f.halfedge and e.halfedge will be the same halfedge or f.halfedge and e.halfedge.twin will be the same halfedge. Does that clear things up?

15. Leslie says:

I think I understand this point. Correct me if I am wrong.
Basically, there are two kinds of halfedge, f.halfedge and e.halfedge. Both of them are predesigned and is set automatically when the mesh is loaded.
However, f.halfedge are decided to make the orientation of the face CCW. e.halfedge may be set randomly or based on some other theory. Therefore, when e.halfedge and f.halfedge are same halfedge, the corresponding element $d_{ij}$ of derivative matrix $d_1$ is 1, otherwise will be -1.
When we use command e.halfedge, we will get an edge with specific orientation. When we use command f.halfedge, it will return maybe more than one halfedges of more than one edges, and these halfedges are associated with this face.
Are all of these correct?
Thanks!

1. rohansawhney says:

Parts of what you mentioned are correct, but I think two separate ideas are being conflated here. The first is how connectivity information is encoded in the halfedge mesh data structure, and the second is how the the exterior derivative matrix d1 is constructed. I recommend looking at the description of the halfedge data structure in the slides, and in module-Core.html in the docs folder. Let me try to walk through each of the points you brought up.

1) “Basically, there are two kinds of halfedge, f.halfedge and e.halfedge.”

There is only one kind of halfedge. f.halfedge and e.halfedge could be different halfedges (i.e., different objects in memory), but they aren’t different kinds of halfedges. Both are objects instantiated from the same Halfedge class. The mesh class stores a list of halfedges, and each face and edge in the mesh has a reference to one of the halfedges stored in this list.

2) “However, f.halfedge are decided to make the orientation of the face CCW. e.halfedge may be set randomly or based on some other theory.”

Yes, by convention, the orientation of the faces in the mesh is CCW. This orientation (determined by the 3 halfedges of the face) doesn’t necessarily have anything to do with how the exterior derivative matrix d1 is constructed. Furthermore, the choice of which halfedge an edge stores as a reference is arbitrary. It can be any one of the two halfedges that are incident on the edge. This too doesn’t necessarily have anything to do with how d1 is constructed.

3) “Therefore, when e.halfedge and f.halfedge are same halfedge, the corresponding element dij
of derivative matrix d1 is 1, otherwise will be -1. When we use command e.halfedge, we will get an edge with specific orientation.”

The edge by itself has no orientation. But to construct the matrix d1, you need to assign each edge in the mesh an orientation. And you have absolute freedom is deciding how this is done, independent of how the halfedge mesh is constructed. The only condition is that for an edge shared by two faces, it should have an opposite orientation for each face, i.e. +1 for one face and -1 for the other. One way to assign an orientation to an edge could be to give it the same orientation as the halfedge is stores a reference to. Another way could be to assign it the orientation of its halfedge’s twin. Its up to you.

Unfortunately, as I mentioned earlier today in the email I sent out, the grading script incorrectly assumes that all edges are assigned an orientation that agrees with either their halfedge or their halfedge’s twin. Put another way, if you decide to give one edge the same orientation as its halfedge, then all the other edges should also have the same orientation as their own halfedge.

1. rohansawhney says:

One other thing. Your construction of d1 should be consistent with your construction d0 – edges should have the same orientation in both cases. d1 times d0 should equal a |F| x |V| zero matrix. Why?

1. Leslie says:

Thank you so much!
d1 times d0 should equal a |F|times |V|zero matrix if we keep the consistency of halfedge, since we apply discrete exterior derivative twice.
Thanks again!