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.

Grading

This assignment is worth 6.5% of your grade.

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

    1. 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!

  1. 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. 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.

  2. 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. 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. Thanks for the reply.

        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. 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. 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. 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 :)).

          3. 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.

          4. 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.

          5. 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.

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

          7. 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?

          8. 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.

          9. 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?

          10. “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.

          11. 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.

  3. 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. 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. 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.

  4. 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. 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.

      1. 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.

    1. 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.

  5. 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. 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. 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.

  6. 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. 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. 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 🙁

  7. 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. 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.

  8. 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. 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. 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];

  9. Hi, 2 questions about halfedge.
    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. 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. 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. 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?

  10. 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. 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. 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. 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!

Leave a Reply