This is my best attempt using python and not using geo-nodes: if you create two separate meshes and the combine them by either using a boolean union modifier or by simply joining them together, the resulting mesh contains the right angles. Just un-triple-quote which ever option you prefer.
import bpy
vertices=[[0,0,1],[0,1,1]]
edges=[[0,1]]
faces = []
name = "my_bmesh"
mesh = bpy.data.meshes.new(name)
obj1 = bpy.data.objects.new(name, mesh)
coll = bpy.data.collections.get("Collection")
coll.objects.link(obj1)
bpy.context.view_layer.objects.active = obj1
mesh.from_pydata(vertices, edges, faces)
dim = (0.3, 0.3)
vertices=[[0,0,0],[0,0,1 + dim[0]]]
edges=[[0, 1]]
faces = []
name = "my_bmesh"
mesh = bpy.data.meshes.new(name)
obj2 = bpy.data.objects.new(name, mesh)
coll.objects.link(obj2)
bpy.context.view_layer.objects.active = obj2
mesh.from_pydata(vertices, edges, faces)
Add skin modifier
mod_skin = obj1.modifiers.new('Skin', 'SKIN')
Add radius
for v in obj1.data.skin_vertices[0].data:
v.radius = dim
mod_skin = obj2.modifiers.new('Skin', 'SKIN')
for v in obj2.data.skin_vertices[0].data:
v.radius = dim
OPTION 1 --- BOOLEAN ROUTE (NO INNER GEOMETRY)
"""
bpy.ops.object.modifier_add(type='BOOLEAN')
bool = bpy.context.object.modifiers["Boolean"]
bool.object = obj1
bool.operation = 'UNION'
bpy.context.view_layer.objects.active = obj1
bpy.ops.object.modifier_apply(modifier="Skin")
bpy.context.view_layer.objects.active = obj2
bpy.ops.object.modifier_apply(modifier="Skin")
bpy.ops.object.modifier_apply(modifier="Boolean")
bpy.data.objects.remove(obj1)
"""
OPTION 2 --- JOIN ROUTE (CONTAINS INNER GEOMETRY)
"""
obj1.select_set(True)
obj2.select_set(True)
bpy.ops.object.join()
can uncomment the following line if you only want the surrounding mesh
#bpy.ops.object.modifier_apply(modifier="Skin")
"""
If you want to go the geometry node route, this is the best node set-up I could come up with:

On top of the geometry nodes, be sure to add the solidify modifier as follows

Here's the geometry nodes script:
import bpy
vertices=[[0,0,0],[0,0,1],[0,1,1]]
edges=[[0,1],[1,2]]
faces = []
name = "my_bmesh"
mesh = bpy.data.meshes.new(name)
obj = bpy.data.objects.new(name, mesh)
coll = bpy.data.collections.get("Collection")
coll.objects.link(obj)
bpy.context.view_layer.objects.active = obj
mesh.from_pydata(vertices, edges, faces)
add GeometryNodes modifier
bpy.ops.object.modifier_add(type="NODES")
bpy.ops.node.new_geometry_node_group_assign()
access active object node_group
node_group = bpy.context.object.modifiers["GeometryNodes"].node_group
add socket
inputs = node_group.inputs
add nodes
nodes = node_group.nodes
mesh_to_curve = nodes.new(type = "GeometryNodeMeshToCurve")
mesh_to_curve .location.x -= 150
mesh_to_curve .location.y -= 0
curve_to_mesh = nodes.new(type="GeometryNodeCurveToMesh")
curve_to_mesh.location.x += 20
curve_to_mesh.location.y += 0
bez = nodes.new(type="GeometryNodeCurvePrimitiveBezierSegment")
bez.location.x += -150
bez.location.y -= 150
bez.inputs[1].default_value[0] = .3
bez.inputs[2].default_value[0] = 0
bez.inputs[2].default_value[1] = 0
bez.inputs[4].default_value[0] = -.3
connect
links = node_group.links
links.new(nodes["Group Input"].outputs["Geometry"], mesh_to_curve.inputs["Mesh"])
links.new(mesh_to_curve.outputs["Curve"], curve_to_mesh.inputs["Curve"])
links.new(curve_to_mesh.outputs["Mesh"], nodes["Group Output"].inputs["Geometry"])
links.new(bez.outputs["Curve"], curve_to_mesh.inputs["Profile Curve"])
bpy.ops.object.modifier_add(type='SOLIDIFY')
solidify = bpy.context.object.modifiers["Solidify"]
solidify.thickness = 0.6
solidify.use_even_offset = True
solidify.offset = 0
Also, for this geometry nodes route, if you don't like the shading\smoothness, you can apply the geo-nodes and then auto-smooth the mesh or you can apply a Remesh modifier and set it to sharp or you can bump up the resolution of the bezier segment