I am using Dunavant quadrature in my software, specifically this file by John Burkardt. Recently, I wanted to convert the code into a constexpr code in C++. But doing so requires me understanding what the code does. I have read the code and the paper and I am confused about the suborders. I have specifically three questions:

  1. What are the suborders, how are they related to the list of points and weights?

  2. Why are the suborders apprarently calculated in three dimensions and then only two of them are used?

  3. Why does the author use a mod function to calculate some suborders like in this code?

    for (k = 0; k < 3; k++) {
       xy[0 + o * 2] = suborder_xyz[i4_wrap(k, 0, 2) + s * 3];
       xy[1 + o * 2] = suborder_xyz[i4_wrap(k + 1, 0, 2) + s * 3];
       w[o] = suborder_w[s];
       o = o + 1;
  • $\begingroup$ This is a question about a specific software package. You will be better off asking the author(s) themselves since only they can help you with your questions. $\endgroup$ – Wolfgang Bangerth Nov 5 '17 at 18:24
  • $\begingroup$ I have actually mailed the author and have got a response from them. I am wondering now if I should delete the question or post an answer to my own question. $\endgroup$ – user1800 Nov 5 '17 at 19:41
  • 5
    $\begingroup$ @JayeshBadwaik If you can focus the answer (and the question by editing it) on the mathematical/algorithmical aspects rather than on this specific implementation (which you can of course still refer to, but then just for the sake of illustration), it would be on-topic in my opinion. (Self-answering is explicitly encouraged on the Stack Exchange network; the next person struggling with this issue will appreciate it.) $\endgroup$ – Christian Clason Nov 5 '17 at 20:49

The integration formulas for triangles are defined in terms of the three area coordinates of the triangle and a weight. Since the sum of the three area coordinates must equal one, only two of them are independent (possibly this is what you mean when you say "calculated in three dimensions and then only two of them are used").

It is important that an integration formula return the same value irrespective of how the vertices of the triangle are ordered, i.e. the coordinates of the integration points should not be biased toward a particular vertex. The centroid of the triangle is one such point where all three area coordinates are equal. This would be a "suborder" of size one. If two of the area coordinates at a particular integration point are the same, then there should be two additional integration points that are permutations of these coordinates; this would be a suborder of size three. Finally, if all three area coordinates at an integration point are different, there would be six permutations of these; i.e. a suborder of size six.

The mod function is a clever way of calculating these different permutations. But, in this case, since there are only three or six of them, simply listing them explicitly is somewhat more obvious and not an unreasonable way to implement the routine.


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.