How to calculate interior control points for continuous bicubic Bezier surfaces for C2

I am trying to calculate control points for a bicubic bezier surface of smoothly connecting tubes. This example gives 16 breakpoints for the corresponding Bézier patch (in OpenSCAD syntax):

[[[ 2 ,  2 , 0], [2.5, 1.5, 0], [3.5,  1 , 0], [4,  1 , 0]],
 [[ 2 ,  2 , 1], [2.5, 1.5, 1], [3.5,  1 , 1], [4,  1 , 1]],
 [[1.3, 1.3, 2], [1.8,  0 , 2], [ 3 , 0.5, 2], [4, 0.5, 2]],
 [[ 0 ,  0 , 2], [1.8,  0 , 2], [ 3 ,  0 , 2], [4,  0 , 2]]]

      

        
        
        
      

    

The corner points are given; control points (red) at the edges were computed tangentially with explicit directions: The patch can be mirrored on the z = 0, y = 0, x = 4 and x = y planes to generate a larger folded surface that is C1 continuous along all four seams. This surface is also C2 continuous along the seams in the z = 0, y = 0 and x = 4 planes, but not so along the seam in the x = y plane (see V-shaped yellow normal vectors):

For the desired C2 continuity, the surface normals along the breaking weld must be coplanar to the x = y plane.

Is it possible to set the four inner controls such that normals along all seams come out at will?

If so, what formula or algorithm should the internal control points be calculated to ensure a completely C2 continuous folded surface?

If this is not the case, what is the approach to create the required patch?