auto void add_edge(toposlic_vert_id_t kv0, toposlic_vert_id_t kv1);
      /* Add arcs {Arc[2*ke]} and {Arc[2*ke+1]} from vertices {Vpos[kv0]} and {Vpos[kv1]}.
        Then increments {*ke}. */
    
    /* Add primary vertices of wall, and longitudinal edges: */
    for (int32_t r = 0; r <= NR; r++)
      { toposlic_vert_id_t kv0 = add_main_vert(&kv, 0, r); /* Main vertex {s,r}. */
        for (int32_t s = 0; s < NS; s++)
          { toposlic_vert_id_t kv1; /* Main vertex {s+1,r}. */
            if (s == NS-1)
              { kv1 = kv - NS; }
            else
              { kv1 = add_main_vert(&kv, s+1, r); }
            /* Horizontal edge from {s,r} to {s+1,r}: */
            add_edge(kv0, kv1);
            kv0 = kv1;
          }
      }

    /* Add secondary vertices and meridional edges: */
    for (int32_t r = 0; r < NR; r++)
      { for (int32_t s = 0; s < NS; s++)
          { toposlic_vert_id_t kv0 = r*NS + s;
            toposlic_vert_id_t kv1 = kv0 + NS;
            toposlic_vert_id_t kvp = kv0;
            for (int32_t b = 0; b < NB; b++)
              { /* Add extra vertex {b} along {L}. */
                toposlic_vert_id_t kvb = add_interp_vert(s, r, b);
                /* Add meridional arc to {vb}: */ 
                add_edge(kvp, kvb);
                kvp = kvb;
              }
            /* Add final edge from {kvp} to {kv1}: */
            add_edge(kvp, kv1);
          }
       }

    assert(mesh.NE == NE);
    assert(mesh.NV == NV);
    
    /* Set the {.skip} links to be {.lnext}: */
    for (int32_t r = 0; r < NR; r++)
      { for (int32_t s = 0; s < NS; s++)
          { toposlic_arc_id_t kapo = 2*(r*NS + s); /* Arc from {s,r} to {s+1,r}. */
            toposlic_arc_id_t kamo = kapo + (s == 0 ? 2*NS : 0) - 2; /* Arc from {s-1,r} to {s,r}. */
            toposlic_arc_id_t kapp = kapo + 2*NS; /* Arc from {s,r+1} to {s+1,r+1}. */
            toposlic_arc_id_t kamp = kamo + 2*NS; /* Arc from {s-1,r+1} to {s,r+1}. */
            
            if (r == 0) { Arc[kamo].skip = kapo; }
            if (r == NR-1) { Arc[kapp+1].skip = kamp+1; }
            
            toposlic_arc_id_t kpu = kamo;
            toposlic_arc_id_t kpv = kapo;
            for (int32_t b = 0; b <= NB; b++)
              { toposlic_arc_id_t kab = 2*(NS*(NR+1) + (r*NS + s)*(NB+1)) + 2*b; /* First arc up. */
                Arc[kpu].skip = kab;
                Arc[kab+1].skip = kpv;
                kpu = kab;
                kpv = kab+1;
              }
            Arc[kpu].skip = kamp + 1;
            Arc[kpu+1].skip = kapp;
          }
      }

    (*Arc_P) = Arc;
    (*Vpos_P) = Vpos;
    (*Vlab_P) = Vlab;
    
    void add_edge(toposlic_vert_id_t kv0, toposlic_vert_id_t kv1)
      {
        int32_t ke = mesh->NE;
        assert(ke < NE);
        toposlic_arc_id_t ka = 2*ke;
        Arc[ka] = (toposlic_arc_t){ .skip = -1, .ivorg = kv0, .pred = ka, .succ = ka };
        ka++;
        Arc[ka] = (toposlic_arc_t){ .skip = -1, .ivorg = kv1, .pred = ka, .succ = ka };
        mesh->NE = ke + 1;
      }

toposlic_coord_t *pick_planes(int32_t NP, toposlic_coord_t Zmax, int8_t debug)
  { 
    toposlic_coord_t *Zplane = malloc(NP*sizeof(toposlic_coord_t));
    assert(Zplane != NULL);
    double w = 6*M_PI; /* Freq of smooth variation (must be even multiple of {M_PI}). */
    double a = 0.5/w;  /* Amplitude of smooth variation. */
    double d = ((double)random())/((double)INT32_MAX)*2*M_PI; /* Phase of smooth variation. */
    for (toposlic_plane_id_t ip = 0; ip < NP; ip++)
      { /* Generate a step with gradual and random variation: */
        double r = ((double)random())/((double)INT32_MAX);
        double f = ((double)ip + 0.25*r + 0.5)/((double)NP);
        f = f + a*(sin(w*f + d) - sin(d));  /* Assumes that {w} is even multiple of {M_PI}. */
        toposlic_coord_t Zp = (toposlic_coord_t)floor((1-f)*Zmin + f*Zmax + 0.5);
        /* Ensure {Zp} is odd: */
        Zp = 2*(Zp/2) + 1;
        /* Ensure strictly increasing and between {Zmin..Zmax}: */
        assert((Zmin < Zp) && (Zp < Zmax));
        if (ip > 0) { assert(Zp > Zplane[ip-1]); }
        Zplane[ip] = Zp;
      }
    if (debug) { fprintf(stderr, "  plane Z range = {%+d ..%+d}\n", Zplane[0], Zplane[NP-1]); } 
    assert(Zplane[0] >= Zbot-2);
    assert(Zplane[NP-1] <= Ztop+2);
    return Zplane;
  }