Finding the largest jacobian in a finite element mesh

Question

So I'm trying to implement the following problem in a couple of different languages to get a rough idea of how different languages work & feel. I'm completely new to C, and this is the first code I've cobbled together using the C reference and Google/SO.

The problem is this:

• Generate a list of n random doubles in (0, 1) and sort it
• Create n "nodes" or points, so a Node is just a type with three doubles x,y,z
• set x of each node to the corresponding entry in the list
• generate a list of "elements", called a "mesh", where each element has two nodes
• compute the Jacobian for each element (see jacobian below)
• print the largest Jacobian of all elements

The following code works, I think. Since I have no experience in C, I want to know:

• How can the code be made more idiomatic?
• How can the code be made more performant?

#include <stdlib.h>
#include <math.h>
#include <stdio.h>

struct Node
{
    double x, y, z;
};

struct Element
{
    struct Node node1, node2;
};

struct Mesh
{
    struct Element* elements;
    int num_elements;
};

int compare_doubles(const void* a, const void* b)
{
    double arg1 = *(const double*)a;
    double arg2 = *(const double*)b;

    if (arg1 < arg2) return -1;
    if (arg1 > arg2) return 1;
    return 0;
}

double* generate_coordinates(int num_nodes)
{
    double* coordinates = malloc(num_nodes * sizeof(double));
    for (int i=0; i < num_nodes; i++)
    {
        coordinates[i] = (double) rand() / RAND_MAX;
    }

    qsort(coordinates, num_nodes, sizeof(double), compare_doubles);

    return coordinates;
}

struct Mesh generate_mesh(double* coordinates, int num_nodes, int num_elements)
{
    struct Node* nodes = malloc(num_nodes * sizeof *nodes);
    for (int i=0; i < num_nodes; i++)
    {
        nodes[i].x = coordinates[i];
        nodes[i].y = 0.0;
        nodes[i].z = 0.0;
    }

    struct Element* elements = malloc(num_elements * sizeof *elements);
    for (int i=0; i < num_nodes; i++)
    {
        elements[i].node1 = nodes[i];
        elements[i].node2 = nodes[i+1];
    }

    struct Mesh mesh;
    mesh.elements = elements;
    mesh.num_elements = num_elements;

    return mesh;
}

double jacobian(struct Element elem)
{
    double dx = elem.node2.x - elem.node1.x;
    double dy = elem.node2.y - elem.node1.y;
    double dz = elem.node2.z - elem.node1.z;
    return 0.5 * sqrt(dx*dx + dy*dy + dz*dz);
}

double find_max_jacobian(struct Mesh mesh)
{
    double max = 0.0;
    for (int i=0; i < mesh.num_elements; i++)
    {
        double jac = jacobian(mesh.elements[i]);
        if (jac > max)
            max = jac;
    }

    return max;
}

int main()
{
    int num_elements = 5000000;
    int num_nodes = num_elements + 1;

    double* coordinates = generate_coordinates(num_nodes);
    struct Mesh mesh = generate_mesh(coordinates, num_nodes, num_elements);

    printf("Max jacobian: %e\n", find_max_jacobian(mesh));
    return 0;
}

---

Edit: I just realised that I'm storing the nodes twice (right?). So I've made the nodes in the element just pointers, and modified the access in jacobian, so

struct Element
{
    struct Node *node1, *node2;
};
// ...
    elements[i].node1 = &nodes[i];
    elements[i].node2 = &nodes[i+1];
// ...
    double dx = elem.node2->x - elem.node1->x;

Answer 1

Welcome to the wonderful world of C programming where you have to do everything yourself.

The C programming language is basically a high level assembly language, there is no garbage collection and the programmer has to manage the memory themselves.

The Good

Before I start discussing what can be improved, I'll mention the good things:
- The code is quite readable
- No casts from malloc() to the receiving variable
- The code follows the Single Responsibility and KISS principles, lots of nice small functions, although the function struct Mesh generate_mesh(double* coordinates, int num_nodes, int num_elements) could have 2 sub-functions to implement it, each of the for loops could be a function.
- There is only one Magic Number, and that one is forgivable, the 5000000 used to initalize the variable num_elements

Error Checking

The allocation functions malloc(), calloc() and realloc() have been known to fail. This happens less often now than historically because computers generally have more memory than they used to, but there are still cases of limited memory such as in embedded programming.

When an allocation function fails it returns NULL. Reference through a NULL pointer causes unknown behavior, and generally the program will terminate without an explanation. It is a good habit to test the results of the allocation before using the pointer in code:

    double* coordinates = malloc(num_nodes * sizeof(double));
    if (!coordinates)
    {
        fprintf(stderr, "Allocation of coordinates failed, exiting program\n");
        exit(EXIT_FAILURE);
    }

    for (int i=0; i < num_nodes; i++)
    {
        coordinates[i] = (double) rand() / RAND_MAX;
    }

There are other ways to handle this error, you might want to learn about setjmp and longjmp. There is a discussion about setjmp and longjmp in this stackoverflow question. There are problems using exit() in some cases.

Prevent Memory Leaks

As mentioned above, C doesn't have garbage collection, the programmer is responsible for de-allocating memory through the use of free(). If this program was part of a larger more complex program there would be memory leaks because none of the allocated memory is ever freed. This could lead to other memory allocation failing.

In reference to the edit, you could free the nodes array after the content of the nodes were copied into the elements.

Prefer calloc Over malloc for Arrays

The function calloc(size_t num, size_t size) is generally preferred over malloc(size_t size) for arrays. The first parameter is the number of items to allocate, the second number is the size of the items. The parameters obviously make sense for arrays. One of the benefits of calloc() is that it zeros out each element of the array when the memory is allocated.

When allocating memory the code become more maintainable when the sizeof() function contains the variable the memory is being assigned to rather than the actual type. This allows the programmer or maintainer t change the type of the variable without having to change the code of the allocation itself. An example is :

    double* coordinates = calloc(num_nodes, sizeof(*coordinates));
    if (!coordinates)
    {
        fprintf(stderr, "Allocation of coordinates failed, exiting program\n");
        exit(EXIT_FAILURE);
    }

In the code above changing the type of the variable coordinates automatically changes the size of the items being allocated.