How to discretize the boundary of two regions?

Question

Problem

How to discretize a surface constructed with multiple components?

Unsuccessful attempts include combining the objects top and hull using the Mathematica commands Graphics3D, Union, and RegionUnion.

Other posts investigated

Combine regions?

Boundary discretize region of ellipsoid returns a three dimensional region

How to combine regions of two 3D plots

RegionUnion issues with many Regions

Example

Combine and mesh the surface defined by these two blocks:

{length, beam, draft} = {50, 3, 4}
pmin={0, 0, 0};
pmax={length, beam, draft};
hull = Cuboid[pmin, pmax];
{topLength, height} = {30, 3};
pmin = {10, 0, draft};
pmax = pmin + {topLength, beam, height};
top = Cuboid[pmin, pmax];
ohp = RegionUnion[top, hull]

Fails to discretize:

BoundaryDiscretizeRegion[ohp, MaxCellMeasure -> {"Length" -> 5}]
BoundaryDiscretizeRegion: A non-degenerate region is expected at position 1

Answer 1

Here is an option using OpenCascadeLink. OpenCascade is an open source 3D CAD package that often does a better job retaining sharp features with boolean operations and seems to be fairly robust.

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
{length, beam, draft} = {50, 3, 4};
pmin = {0, 0, 0};
pmax = {length, beam, draft};
hull = Cuboid[pmin, pmax];
{topLength, height} = {30, 3};
pmin = {10, 0, draft};
pmax = pmin + {topLength, beam, height};
top = Cuboid[pmin, pmax];
shape1 = OpenCascadeShape[hull];
shape2 = OpenCascadeShape[top];
union = OpenCascadeShapeUnion[shape1, shape2];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
BoundaryDiscretizeRegion[
 MeshRegion[
  MeshOrderAlteration[
   ToElementMesh[bmesh, MaxCellMeasure -> {"Length" -> 5}], 1]]]

Answer 2

Here is a slightly different approach also using OpenCascadeLink

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[ohp, "BoundaryMeshGenerator" -> "OpenCasdade"]
MeshRegion[bmesh]

Note, however, there is a slight difference in the result compared to Tim's answer. In this case the union is created. I.e. no subdivision between the two cuboids. Tim's answer is a more general approach.