Shader: Rounded edge between the two intersecting faces?
There's an OSL bevel Shader for Cycles that I believe was originally written by the person who built Cycles.  Several people have been making adjustments to it and have been posting the results to the following thread on BlenderArtists.org  http://www.blenderartists.org/forum/showthread.php?329295-Bevel-shader
OSL Code
Paste the following code into a new Text Datablock in Blender, and use it in the script node of an OSL node tree:
void rng_seed(output int rng, int seed)
{
  int chash = seed;
  if (chash == 0) chash = 1;
  rng = chash * 30391861;
}

float rng_uniform(output int rng)
{
  float res = rng / float(2137483647) * 0.5 + 0.5;
  rng *= 30391861;
  return res;
}

void to_unit_disk(float x, float y, output float x_out, output float y_out)
{
  float r, phi;
  float a = 2.0 * x - 1.0;
  float b = 2.0 * y - 1.0;
    
  if(a &gt; -b) 
  { if(a &gt; b) 
    { r = a;
      phi = M_PI_4 *(b/a);
    }
    else 
    { r = b;
      phi = M_PI_4 *(2.0 - a/b);
  } }
  else 
  { if(a &lt; b) 
    { r = -a;
      phi = M_PI_4 *(4.0 + b/a);
    }
    else 
    { r = -b;
      if(b != 0.0) phi = M_PI_4 *(6.0 - a/b);
      else phi = 0.0;
  } }
  x_out = r * cos(phi);
  y_out = r * sin(phi);
}

void make_orthonormals(vector N, output vector a, output vector b)
{
  if(N[0] != N[1] || N[0] != N[2]) a = cross(vector(1, 1, 1), N);
  else a = cross(vector(-1, 1, 1), N);
  
  a = normalize(a);
  b = cross(N, a);
}

vector sample_cos_hemisphere(vector N, float randu, float randv)
{
  vector T, B;
    
  make_orthonormals(N, T, B);
  to_unit_disk(randu, randv, randu, randv);
  float costheta = sqrt(max(1.0 - randu * randu - randv * randv, 0.0));

  return randu * T + randv * B + costheta * N;
}

shader node_occlusion2(
  color Effect = color(0),
  color No_Effect = color(1),
  int Mode = 0, /* 0: Concave (AO) 1:Convex (Wear) 2:Both */
  int InvertEffect = 0,
  float Distance = 0.2,
  int Samples = 1,
  output color Color = 0,
  output float Fac = 0,
  output normal outNormal = N
)
{
  int i, rng;
  float f, randu, randv, ray_t, hits = 0;
  vector ray_P, ray_R;
  normal hit_normal = N;
  float hit_dist;

  f = fmod(cellnoise(P*123456.0), 1.0);
  rng_seed(rng, int(f * 21374647));
  
  for(i = 0; i &lt; Samples; i++) 
  { randu = rng_uniform(rng);
    randv = rng_uniform(rng);
       
    ray_P = P;
    ray_R = sample_cos_hemisphere(-N, randu, randv);
    ray_t = Distance;
    
        
    if (!Mode)
    { if(trace(ray_P, -ray_R, &quot;maxdist&quot;, ray_t)) {
            hits += 1.0;
            int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
            outNormal = outNormal + (hit_normal);
        }
    }
    else if (Mode == 1)
    { if(trace(ray_P, ray_R, &quot;maxdist&quot;, ray_t)) {
           hits += 1.0;
           int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
           outNormal = outNormal - (hit_normal);
        }
    }
    else { 
        if(trace(ray_P, -ray_R, &quot;maxdist&quot;, ray_t)) {
            hits += 1.0;
            int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
            outNormal = outNormal + (hit_normal);
        }
        if(trace(ray_P, ray_R, &quot;maxdist&quot;, ray_t)) {
           hits += 1.0;
           int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
           outNormal = outNormal - (hit_normal);
        }
    } 
  }
  Fac = 1.0 - (hits/Samples);
  if(InvertEffect) Color = mix(No_Effect, Effect, Fac);
  else Color = mix(Effect, No_Effect, Fac);
  outNormal = normalize(outNormal);
}

Here's a render of many spheres that are all separate objects.  It works well in some cases but not always.