Here is a short example. I assume that by overlay you are looking for intersects; namely to dissolve all polygons that intersect from both layers.
library(sp)
library(rgeos)
# Create a dataset
poly <- SpatialPolygons(list(
Polygons(list(Polygon(coords = matrix(c(1, 1, 4, 3, 4, 2, 1, 1), ncol = 2, byrow = TRUE))), ID = "1"),
Polygons(list(Polygon(coords = matrix(c(3, 1.5, 3.5, 2, 4, 2, 3, 1.5), ncol = 2, byrow = TRUE))), ID = "2"),
Polygons(list(Polygon(coords = matrix(c(4, 1, 5, 2, 5, 1, 4, 1), ncol = 2, byrow = TRUE))), ID = "3")
))
# Split dataset to two SpatialPolygon objects
polyA <- poly[1, ]
polyB <- poly[2:3, ]
# Show dataset
plot(poly, axes = TRUE)
plot(polyA, axes = TRUE, add = TRUE)
plot(polyB, axes = TRUE, add =TRUE, col = "red")
You can see the sample data set below. PolyB has two polygons from which one intersects the polygon in PolyA.
Using gUnion will result in one polygon of all polygons in both objects, as you have suggested: plot(gUnion(polyA, polyB)).
Yet, if you select only those polygons that intersects polyB[polyA, ], dissolve will give you the expected result:
plot(gUnion(polyA, polyB[polyA, ]))
You can and should subset the first layer as well, gUnion(polyA[polyB, ], polyB[polyA, ]), if it has more than one feature.
edit
If you want to disaggregate the multi-polygon feature into single-polygon features afterwards you can simply use the disaggregate function from the raster package.
