I wish to find the surface normal, and intersection from any arbitrary shape. I can do this using the Gaphics3D shapes as follows:

I have a constraint function which finds logic to give the allowed values which are along the line and is on one of the surfaces, but not within any two shapes. So, the function to do create the logic is called constraints[] and intersection[] which then just uses NSolve[] to find the solution with the line, i.e.

NSolve[{x, y, z} [Element] line && (surfaceequations1 || surfaceequations2 || surfaceequations3) && ! (volumeequation1 && volumeequation2) && ! (volumeequation2 && volumeequation3) && ! (volumeequation1 && volumeequation3), {x, y, z}]

constraints[shapes__] := 
 And[## & @@ (Not /@ 
      Through[(RegionMember[RegionIntersection@##] & @@@ 
          Subsets[{shapes}, {2}])@#]), 
   RegionMember[RegionUnion @@ (RegionBoundary /@ {shapes})]@#] &

intersections[l_, s__] := 
 NSolve[{x, y, z} \[Element] l && constraints[s][{x, y, z}], {x, y, z}]

shape1 = Cuboid[{5, 5, 1}, {6, 25, 3}];
shape2 = Cuboid[{6, 5, 1}, {21, 6, 3}];

shapes = {shape1, shape2};

line = InfiniteLine[{{5.5, 10.5, 0.5}, {5.5, 10.5, 1}}];

intersection = intersections[line, ##] & @@ shapes;
intersection = DeleteDuplicates[intersection];
points = {x, y, z} /. intersection;
Graphics3D[{line, shapes, {Red, PointSize[0.05], Point[points]}}]

This works

I can also find the surface normals of the shapes as follows

shape1 =  Cuboid[{5, 5, 1}, {15, 15, 4/3}];
shapes = {shape1};

standardize[a_ <= b_] := a - b;
standardize[a_ >= b_] := b - a;

regnormal[reg_, {x_, y_, z_}] := 
  Module[{impl}, 
   impl = LogicalExpand@
     Simplify[
      RegionMember[reg, {x, y, z}], {x, y, z} \[Element] Reals];
   If[Head@impl === Or, impl = List @@ impl, impl = List@impl];
   impl = Replace[impl, {Verbatim[And][a___] :> {a}, e_ :> {e}}, 1];
   Piecewise[
    Flatten[Function[{component}, 
       Table[{D[standardize[component[[i]]], {{x, y, z}}], 
         Simplify[(And @@ Drop[component, {i}] /. {LessEqual -> Less, 
              GreaterEqual -> 
               Greater}) && (component[[i]] /. {LessEqual -> Equal, 
              GreaterEqual -> Equal}), 
          TransformationFunctions -> {Automatic, 
            Reduce[#, {}, Reals] &}]}, {i, Length@component}]] /@ 
      impl, 1], Indeterminate]];

normals = regnormal[RegionUnion[##] & @@ shapes, {x, y, z}]

But I wish to do it with more complex shapes, i.e. ones I may have designed in a modelling package, for instance, the Example Data helicopter.

image = Import["ExampleData/helicopter.dxf.gz"]

How can I find the surface normals, and the intersection with a line with shapes like this? Mathematica must be doing it for the Graphics plots. Any ideas?