The first aspect we will have to take a closer look at, is how we can generate our starting point, i.e. the initial state before the simulation is run.
We’ve already discussed the first half of that in a previous post, where we’ve generated a simple triangle mesh of a sphere. So, I’ll first recap the most important part of that and then dive into how we can populate that mesh with the information we need for the actual simulation, like plate-IDs, plate-velocities and elevations.
Creating the mesh
To recap, what we need is a set of random points, that are uniformly distributed on the surface of a sphere, where each point represents part of a tectonic plate. To generate uniformly distributed points[1], we utilize the golden angle to create a Fibonacci Sphere, by placing points walking downwards from one of the poles and rotating radians[2] each step (image A).
Once we have a set of points, the next step is to generate a triangle mesh from them. There are many possible ways to triangulate a set of points, but the one we’ll use is the Delaunay triangulation, because it avoids long thin triangles, that could be problematic during the simulation, and also because its dual-mesh is the Voronoi diagram of the mesh, which we’ll utilize for other parts of the simulation. While there are more efficient algorithms to generate a Delaunay triangulation of a set of points, since our points all lay on the surface of a sphere, the easiest approach is to generate the convex hull of the points (using the QuickHull algorithm), which in this case coincides with the Delaunay triangulation of the points (image B).
The algorithm described above distributes the points completely uniformly. But for a more natural look, it’s actually preferable to sacrifice some uniformity to achieve a less predictable structure. Luckily, all that’s required to achieve this, is to offset each point by a random amount (image C). And as long as that offset is smaller than our step-size, the points are still guaranteed to be far enough apart.
This gives use the following code to generate our mesh:
// Define the data layer that contains the vertex positions in 3D Cartesian coordinates
constexpr auto position_info = Layer_info<Vec3, Layer_type::vertex>("position");
// Acquire all necessary resources from the world structure
auto mesh = world.lock_mesh();
auto [positions] = world.lock_layer(position_info);
auto rand = world.lock_random();
// Create N vertices
mesh->add_vertices(vertex_count);
// Pre-calculate the golden angle and the step-size for calculating y,
// so we just need to multiply them with i in the loop
constexpr auto golden_angle = 2.f * std::numbers::pi_v<float> * (2.f - std::numbers::phi_v<float>);
const auto step_size = 2.f / (vertices - 1);
for(std::int32_t i = 0; i < vertices; i++) {
const auto offset = perturbation > 0 ? rand->uniform(0.f, perturbation) : 0.f;
auto ii = std::min(i + offset, vertices - 1.f);
// Calculate the x/y/z position of the current vertex on the unit sphere
const auto y = 1.f - step_size * ii;
const auto r = std::sqrt(1 - y * y);
const auto theta = golden_angle * ii;
const auto x = std::cos(theta) * r;
const auto z = std::sin(theta) * r;
// Set the vertex position (multiplying with the radius of our sphere)
positions[Vertex(i)] = Vec3{x, y, z} * radius;
}
// Triangulate the generated points
auto qh = quickhull::QuickHull<float>{};
auto hull = qh.getConvexHull(&positions.begin()->x, positions.size(), false, true);
auto& indices = hull.getIndexBuffer();
for(std::size_t i = 0; i < indices.size(); i += 3) {
mesh->add_face(indices[i], indices[i+1], indices[i+2]);
}
Generating the first tectonic plates
Besides the mesh itself, there are also a couple of values the simulation uses, that we need to initialize to create the initial set of tectonic plates. As already mentioned, these values are attached to the meshes vertices to describe the properties of the area around the vertex and are defined as follows in the code:
enum class Crust_type : int8_t { none, oceanic, continental };
// Position of the vertex in space (initialized by the process above)
constexpr auto layer_position = Layer_info<Vec3, Layer_type::vertex>("position");
// Integer ID of the plate, the vertex belongs to
constexpr auto layer_plate_id = Layer_info<int32_t, Layer_type::vertex>("plate_id");
// Dominant type of the crust in this area
constexpr auto layer_plate_type = Layer_info<Plate_type, Layer_type::vertex>("crust_type");
// Velocity of this piece of crust
constexpr auto layer_velocity = Layer_info<Vec3, Layer_type::vertex>("velocity");
// Height of this piece of crust above/below the (arbitrary) sea level
constexpr auto layer_elevation = Layer_info<float, Layer_type::vertex>("elevation");
// The timestamp at which the crust has formed
constexpr auto layer_created = Layer_info<float, Layer_type::vertex>("crust_created");
The process of generating the values itself is actually incredibly simple. First we pick random vertices as seed-points to create tectonic plates, assign them each a random elevation/velocity/… and finally use a flood-fill type algorithm to grow them in all directions, until all vertices are assigned to a plate.
Initialize seed vertices
First we’ll need to pick a couple of random points that will be the seeds from which we then grow our tectonic plates. The number of plates we want to generate, as well as what percentage of these will be oceanic/continental crust, is defined by a set of parameters. Then we just loop times, choosing a random vertex[3] in each iteration and assigning it a random initial state (ID of the plate, type of crust, age, elevation and its initial velocity).
const auto count = rand.uniform(plates_min, plates_max);
const auto ocean_count = static_cast<int>(std::floor(count * ocean_prop));
auto next_id = std::int32_t(1);
auto open_vertices = std::vector<Vertex>();
open_vertices.reserve(static_cast<std::size_t>(count) * 10u);
for(auto i = 0; i < count; i++) {
// Try to find a free vertex to use as a seed point
if(auto seed_vertex = random_vertex(mesh, rng, [&ids = ids](auto v) { return ids[v] == 0; })) {
const auto ocean = i < ocean_count; // First X plates are oceanic, rest is continental
ids[*seed_vertex] = next_id++;
types[*seed_vertex] = ocean ? Crust_type::oceanic : Crust_type::continental;
// Assign other properties based on crust type
// ...
// Remember the seed vertices, so we can use them for the flood-fill
open_vertices.emplace_back(*seed_vertex);
}
}
Most of the properties are pretty straightforward. The ID and crust type are directly defined in the code above. And the crusts age and elevation are then randomly chosen from a range of acceptable values for the given type of crust. The only property that is a bit more complicated to compute is the crusts initial velocity. Because I’ve opted to use a simple Cartesian coordinate system, the velocity is a 3D vector that could point in any direction. But, since the vertices are restrained to the sphere’s surface, the velocity should ideally be locally flat, i.e. perpendicular to the surface normal at any given position. Since our mesh is a sphere, we can derive the surface normal by normalizing the position vector of the current vertex. Based on this vector, we can then construct an orthonormal basis at this point, as shown in the code below[4]. We can then generate a random 2D vector from an angle and a target velocity and transform it into the coordinate basis we just defined.
// Generate random 2D velocity vector
const auto angle = rng.uniform(0.2f, 2.f * 3.141f);
const auto speed = rng.uniform(0.01f, 1.f) * (ocean ? max_ocean_speed : max_continent_speed);
const auto vel_2d = Vec2{std::sin(angle), std::cos(angle)} * speed;
// Build 3D coordinate system at current vertex
const auto normal = normalized(positions[*seed_vertex]);
const auto tangent = cross(normal, Vec3(-normal.z, normal.x, normal.y));
const auto bitangent = cross(normal, tangent);
// Transform the 2D vector to lay in the 3D plane defined by the tangent and bitangent
velocities[*seed_vertex] = vel_2d.x*tangent + vel_2d.y*bitangent;
Flood Fill
The final step of the process is then to iteratively grow each plate from its initial seed, until all vertices are assigned to a tectonic plate.
We already have a list of the seed vertices, that we memorized in the previous step. So, all we have to do now is iterate over each of these vertices neighbors, copy their properties to these neighbors, if they are still uninitialized, and repeat that recursively until there are no uninitialized vertices left.
auto new_opened_vertices = std::vector<Vertex>();
new_opened_vertices.reserve(open_vertices.size());
// Iterate until the last iteration didn't contain any new, formally unassigned, vertices
while(!open_vertices.empty()) {
for(auto& seed : open_vertices) {
for(Vertex n : v.neighbors(mesh)) {
if(types[n] == Crust_type::none) {
ids[n] = ids[seed];
types[n] = types[seed];
elevations[n] = elevations[seed] + /* ... */;
created_times[n] = created_times[seed] + /* ... */;
velocities[n] = normalized_velocity(positions[n], velocities[seed]);
// Remember as starting point for next iteration of outer loop
new_opened_vertices.emplace_back(n);
}
}
}
// Iterate again with the new points,
// but in a random order to avoid growing the first plates larger than the rest
rng.shuffle(new_opened_vertices);
std::swap(new_opened_vertices, open_vertices);
new_opened_vertices.clear();
}
The implementation is again relatively straightforward. We are using a secondary vertex-list, that we populate with the new vertices we initialized in this loop. Iterate over the previous list of seeds, to add each of their neighbors to their plate and remember them for the next iteration. And after each iteration of the outer loop, we swap the two vertex-lists and continue until we had one iteration where no new vertices have been assigned, so we know we are done. The only aspect we have to consider, is that this algorithm favors plates that are visited earlier because they tend to have more unassigned neighbors to “convert” to their plate. Since the first seed also adds its converted neighbors to the front of the next-iterations list, this effect could disproportionately grow the first couple of plates. To counteract this, we need to randomize the order of the vertices in each iteration, so it cancels out over multiple iterations. And as an added bonus, this also tends to generate less circular and more natural-looking shapes.
The process of assigning a new vertex to a seed’s plate is, as can be seen above, again pretty simple. Most of the properties are copied directly from the original plate, applying small perturbations to both the crust formation time[5] and its elevation[6]. The only property that needs additional consideration is (again) the velocity. Because the velocity vector should still be perpendicular to the new vertices surface normal, we can’t just copy it from the original vertex. Instead, we need to re-normalize it using the following function:
Vec3 normalized_velocity(const Vec3& position, const Vec3& velocity) {
constexpr auto dt = 100.f;
auto next_position = position + velocity * dt; // Move position one step with given velocity
next_position = normalized(next_position) * length(position); // Normalize new position back to sphere surface
return (next_position - position) / dt; // Calculate actual velocity
}
With that out of the way, in the following posts we’ll look into the simulation algorithm, that iterates from this initial state to generate our terrain.