Ok, everything I know about this matter comes from this Gamasutra article. To interpolate positions over the triangle, you use this formula:

To calculate two tangent vectors in order to calculate the surface normal. You need to use two partial derivatives of the above formula given as so:

How this looks in code for a 3rd order bezier triangle, or how I think this looks in code, can be seen in my domain shader.

#define p300 0
#define p030 1
#define p003 2
#define p210 3
#define p201 4
#define p120 5
#define p021 6
#define p012 7
#define p102 8
#define p111 9

#define U bary.x
#define V bary.y
#define W bary.z

cbuffer PerFrameBuffer : register(b0)
{
	//float4x4 World; This is just the identity matrix, so not needed
	float4x4 ViewProjection;
	float4 vecEye; 
	float4 LightDirection;
	float4 LightColor; 
};

struct PatchTess
{
	float EdgeTess[3]	: SV_TessFactor;
	float InsideTess	: SV_InsideTessFactor;
};

struct HullOut
{
	float3 PositionWorld	: POSITION;
	float2 MainTexCoord		: TEXCOORD;
};

struct DomainOut
{
	float4 PositionH		: SV_Position;
	float3 Normal			: NORMAL0;
	float3 Tangent			: TANGENT0;
	float3 Bitangent		: BITANGENT0;
	float3 View				: NORMAL1;
	float2 MainTexCoord		: TEXCOORD0;
};

[domain("tri")]
DomainOut DS(PatchTess patchTess, float3 bary : SV_DomainLocation, const OutputPatch<HullOut,10> cp)
{
	DomainOut dout;

	float3 position =
		cp[p300].PositionWorld * pow(U,3) +
		cp[p030].PositionWorld * pow(V,3) +
		cp[p003].PositionWorld * pow(W,3) +
		cp[p210].PositionWorld * 3 * pow(U,2) * V +
		cp[p201].PositionWorld * 3 * pow(U,2) * W +
		cp[p120].PositionWorld * 3 * U * pow(V,2) +
		cp[p021].PositionWorld * 3 * pow(V,2) * W +
		cp[p012].PositionWorld * 3 * V * pow(W,2) +
		cp[p102].PositionWorld * 3 * U * pow(W,2) +
		cp[p111].PositionWorld * 6 * U * V * W;

	float3 tangent =
		cp[p300].PositionWorld * pow(U,2) +
		cp[p120].PositionWorld * pow(V,2) +
		cp[p102].PositionWorld * pow(W,2) +
		cp[p201].PositionWorld * 2 * U * W +
		cp[p210].PositionWorld * 2 * U * V +
		cp[p111].PositionWorld * 2 * V * W;

	float3 bitangent = 
		cp[p030].PositionWorld * pow(V,2) +
		cp[p012].PositionWorld * pow(W,2) +
		cp[p210].PositionWorld * pow(U,2) +
		cp[p120].PositionWorld * 2 * U * V +
		cp[p021].PositionWorld * 2 * V * W +
		cp[p111].PositionWorld * 2 * U * W;

	tangent = normalize(tangent);
	bitangent = normalize(bitangent);
	float3 normal = normalize(cross(tangent,bitangent));

	dout.View = vecEye.xyz - position.xyz;
	dout.PositionH = mul(float4(position,1.0f),ViewProjection);
	dout.Normal = normal;
	dout.Tangent = tangent;
	dout.Bitangent = bitangent;

	dout.MainTexCoord = cp[p300].MainTexCoord * U + cp[p030].MainTexCoord * V + cp[p003].MainTexCoord * W;

	return dout;
}

However, the results I am getting are not quite right, and I wonder if I have some misunderstanding of how this is supposed to work. The positions calculated across each bezier triangle are exactly correct, but I am having trouble calculating a usable and continuous normal from the tangents. The following screenshot uses the above domain shader and is basically a series of instanced bezier triangles defining a height map drawn in two draw calls. One with triangles pointing down and another with triangles pointing up.

[attachment=14966:Bad Normals.jpg]

You can see here some strange effects, the directional light source comes from the top, and there is some sense that some triangles have proper lighting while others do not. I originally thought that the angle between the two tangents was unpredictable, so some of the normals were pointing the opposite direction, but if I brute force the normals to point up with the following code the situation is improved but still not right.

	float3 normal = normalize(cross(tangent,bitangent));
	if(normal.y < 0.0)
		normal = normalize(cross(bitangent,tangent));

[attachment=14967:Bad Normals2.jpg]

Notice that each tri instance has one continuous edge with its neighbor, but the other two edges are wrong. Also I have no idea what causes the diagonal artifact in both images other than my calculations being just wrong. Ultimately I would like a heighmap made from bezier triangles with continuous normals across the whole field. Does anyone here have experience working with bezier triangles who can tell me where my assumptions are wrong?

Ok, I did some googling, and I found something called 'de Casteljau's algorithm', which seems to be a lot more convenient for determining the surface normals, and is just as fast or faster for calculating the point on the triangle. The penultimate step in the algorithm are three points that define the tangent plane of the target point, which is a big bonus. The only problem now is that my normals are not continuous from tri to tri, which surprises me because I thought that because the shared edges are identical, the normals would be also, but I guess not. Here are the new domain shader and screenshot:

#define p300 0
#define p030 1
#define p003 2
#define p210 3
#define p201 4
#define p120 5
#define p021 6
#define p012 7
#define p102 8
#define p111 9

#define U bary.x
#define V bary.y
#define W bary.z

cbuffer PerFrameBuffer : register(b0)
{
	//float4x4 World; This is just the identity matrix, so not needed
	float4x4 ViewProjection;
	float4 vecEye; 
	float4 LightDirection;
	float4 LightColor; 
};

struct PatchTess
{
	float EdgeTess[3]	: SV_TessFactor;
	float InsideTess	: SV_InsideTessFactor;
};

struct HullOut
{
	float3 PositionWorld	: POSITION;
	float2 MainTexCoord		: TEXCOORD;
};

struct DomainOut
{
	float4 PositionH		: SV_Position;
	float3 Normal			: NORMAL0;
	float3 Tangent			: TANGENT0;
	float3 Bitangent		: BITANGENT0;
	float3 View				: NORMAL1;
	float2 MainTexCoord		: TEXCOORD0;
};

[domain("tri")]
DomainOut DS(PatchTess patchTess, float3 bary : SV_DomainLocation, const OutputPatch<HullOut,10> cp)
{
	DomainOut dout;

	float3 p1_200 = U * cp[p300].PositionWorld + V * cp[p210].PositionWorld + W * cp[p201].PositionWorld;
	float3 p1_110 = U * cp[p210].PositionWorld + V * cp[p120].PositionWorld + W * cp[p111].PositionWorld;
	float3 p1_101 = U * cp[p201].PositionWorld + V * cp[p111].PositionWorld + W * cp[p102].PositionWorld;
	float3 p1_020 = U * cp[p120].PositionWorld + V * cp[p030].PositionWorld + W * cp[p021].PositionWorld;
	float3 p1_011 = U * cp[p111].PositionWorld + V * cp[p021].PositionWorld + W * cp[p012].PositionWorld;
	float3 p1_002 = U * cp[p102].PositionWorld + V * cp[p012].PositionWorld + W * cp[p003].PositionWorld;

	float3 p2_100 = U * p1_200 + V * p1_110 + W * p1_101;
	float3 p2_010 = U * p1_110 + V * p1_020 + W * p1_011;
	float3 p2_001 = U * p1_101 + V * p1_011 + W * p1_002;

	float3 position = U * p2_100 + V * p2_010 + W * p2_001;

	float3 tangent = p2_010 - p2_100;
	float3 bitangent = p2_001 - p2_100;

	tangent = normalize(tangent);
	bitangent = normalize(bitangent);
	float3 normal = normalize(cross(bitangent,tangent));

	dout.View = vecEye.xyz - position.xyz;
	dout.PositionH = mul(float4(position,1.0f),ViewProjection);
	dout.Normal = normal;
	dout.Tangent = tangent;
	dout.Bitangent = bitangent;

	dout.MainTexCoord = cp[p300].MainTexCoord * U + cp[p030].MainTexCoord * V + cp[p003].MainTexCoord * W;

	return dout;
}

[attachment=14979:Bad Normals3.jpg]

I've done some more research on the matter, and my second example is actually right. What was wrong was my assumption that because bezier triangles are easy to line up positionally, that they would also have continuous normals. That is not the case, and if you think about it, it makes sense. To get continuous normals, you really need neighbor information. I'll have to come up with something else.

I may be able to generate normals on the edge only, which are identical with their neighbor edge, and then interpolate across the triangle in a way that looks plausible. Hopefully there won't be too much computation involved.

I had researched PN-triangles a bit, but they seem to be very specific for exploding flat surfaces into curved ones on a convex object. It may still be useful for me, and I'll probably try that first, but I have some other ideas I'd like to try too. Now that I can see something, I'm a lot more comfortable tinkering with it.

EDIT: One of the things that gives me doubts about PN-triangles in this case is that my corner normals are always the up vector (0,1,0), so there's nothing to mirror. Although, that might just make things easier, I'm not sure yet.

Wow, I did a little experiment this morning I wanted to show. I think it really explains what I'm trying to do. I adjusted my control points to make the hex shapes more apparent, because previously they were too smooth. You can see the value of getting the normals right across the triangle. I would sure hate to lose that in the quest for continuous normals!

EDIT: I know exactly what to do. I'll report back when I have a chance to implement it.

[attachment=15125:Bad Normals4.jpg]

Well, I tried something similar to PN-Triangles, in that I am interpolating normals from the edge, and here is the result. It's sort of continuous across triangles, but not very well defined like the above example. I may get better results if I increase the order of control points to interpolate, but then we are getting into slow territory. It's a tough problem!

[attachment=15126:Bad Normals5.jpg]

Here is my new hull shader and domain shader. In the hull shader I am using de Casteljau's on the curve made by the edge in order to get the n110, n011 and n101 normals for interpolation. This is a bit different than with standard PN-Triangles:

#define p300 0
#define p030 1
#define p003 2
#define p210 3
#define p201 4
#define p120 5
#define p021 6
#define p012 7
#define p102 8
#define p111 9

#define n200 0
#define n020 1
#define n002 2
#define n110 3
#define n011 4
#define n101 5

struct VertexOutput
{
	float3 PositionWorld	: POSITION;
	float2 MainTexCoord		: TEXCOORD;
};

struct PatchConstants
{
	float EdgeTess[3]	: SV_TessFactor;
	float InsideTess	: SV_InsideTessFactor;

	float3 NormalCP[6]	: NORMAL;
};

PatchConstants PatchHS(InputPatch<VertexOutput,10> cp, uint patchID : SV_PrimitiveID)
{
	PatchConstants pt;

	pt.EdgeTess[0] = 12;
	pt.EdgeTess[1] = 12;
	pt.EdgeTess[2] = 12;
	pt.InsideTess = 12;

	pt.NormalCP[n200] = float3(0,1,0);
	pt.NormalCP[n020] = float3(0,1,0);
	pt.NormalCP[n002] = float3(0,1,0);

	float t = 0.5f;
	float3 bitangent;

	float3 p1_0 = (1.0f - t) * cp[p300].PositionWorld + t * cp[p210].PositionWorld;
	float3 p1_1 = (1.0f - t) * cp[p210].PositionWorld + t * cp[p120].PositionWorld;
	float3 p1_2 = (1.0f - t) * cp[p120].PositionWorld + t * cp[p030].PositionWorld;

	float3 p2_0 = (1.0f - t) * p1_0 + t * p1_1;
	float3 p2_1 = (1.0f - t) * p1_1 + t * p1_2;
	
	float3 tangent = normalize(p2_1 - p2_0);
	float3 offline = normalize(cp[p030].PositionWorld - cp[p300].PositionWorld);
	if(cp[p030].PositionWorld.y > cp[p300].PositionWorld.y)
	{
		bitangent = normalize(cross(offline,tangent));
	}
	else
	{
		bitangent = normalize(cross(tangent,offline));
	}
	pt.NormalCP[n110] = normalize(cross(bitangent, tangent));

	p1_0 = (1.0f - t) * cp[p030].PositionWorld + t * cp[p021].PositionWorld;
	p1_1 = (1.0f - t) * cp[p021].PositionWorld + t * cp[p012].PositionWorld;
	p1_2 = (1.0f - t) * cp[p012].PositionWorld + t * cp[p003].PositionWorld;

	p2_0 = (1.0f - t) * p1_0 + t * p1_1;
	p2_1 = (1.0f - t) * p1_1 + t * p1_2;
	
	tangent = normalize(p2_1 - p2_0);
	offline = normalize(cp[p003].PositionWorld - cp[p030].PositionWorld);
	if(cp[p003].PositionWorld.y > cp[p030].PositionWorld.y)
	{
		bitangent = normalize(cross(offline,tangent));
	}
	else
	{
		bitangent = normalize(cross(tangent,offline));
	}
	pt.NormalCP[n011] = normalize(cross(bitangent, tangent));

	p1_0 = (1.0f - t) * cp[p003].PositionWorld + t * cp[p102].PositionWorld;
	p1_1 = (1.0f - t) * cp[p102].PositionWorld + t * cp[p201].PositionWorld;
	p1_2 = (1.0f - t) * cp[p201].PositionWorld + t * cp[p300].PositionWorld;

	p2_0 = (1.0f - t) * p1_0 + t * p1_1;
	p2_1 = (1.0f - t) * p1_1 + t * p1_2;
	
	tangent = normalize(p2_1 - p2_0);
	offline = normalize(cp[p300].PositionWorld - cp[p003].PositionWorld);
	if(cp[p300].PositionWorld.y > cp[p003].PositionWorld.y)
	{
		bitangent = normalize(cross(offline,tangent));
	}
	else
	{
		bitangent = normalize(cross(tangent,offline));
	}
	pt.NormalCP[n101] = normalize(cross(bitangent, tangent));

	return pt;
}

struct HullOut
{
	float3 PositionWorld	: POSITION;
	float2 MainTexCoord		: TEXCOORD;
};

[domain("tri")]
[partitioning("fractional_even")]
[outputtopology("triangle_cw")]
[outputcontrolpoints(10)]
[patchconstantfunc("PatchHS")]
HullOut HS(InputPatch<VertexOutput,10> p, uint i : SV_OutputControlPointID, uint patchId : SV_PrimitiveID)
{
	HullOut hout;

	hout = p;

	return hout;
}

And here is the interpolation in the domain shader:

#define p300 0
#define p030 1
#define p003 2
#define p210 3
#define p201 4
#define p120 5
#define p021 6
#define p012 7
#define p102 8
#define p111 9

#define n200 0
#define n020 1
#define n002 2
#define n110 3
#define n011 4
#define n101 5

#define U bary.x
#define V bary.y
#define W bary.z

cbuffer PerFrameBuffer : register(b0)
{
	//float4x4 World; This is just the identity matrix, so not needed
	float4x4 ViewProjection;
	float4 vecEye; 
	float4 LightDirection;
	float4 LightColor; 
};

struct PatchConstants
{
	float EdgeTess[3]	: SV_TessFactor;
	float InsideTess	: SV_InsideTessFactor;

	float3 NormalCP[6]	: NORMAL;
};

struct HullOut
{
	float3 PositionWorld	: POSITION;
	float2 MainTexCoord		: TEXCOORD;
};

struct DomainOut
{
	float4 PositionH		: SV_Position;
	float3 Normal			: NORMAL0;
	float3 Tangent			: TANGENT0;
	float3 Bitangent		: BITANGENT0;
	float3 View				: NORMAL1;
	float2 MainTexCoord		: TEXCOORD0;
};

[domain("tri")]
DomainOut DS(PatchConstants pc, float3 bary : SV_DomainLocation, const OutputPatch<HullOut,10> cp)
{
	DomainOut dout;

	float3 p1_200 = U * cp[p300].PositionWorld + V * cp[p210].PositionWorld + W * cp[p201].PositionWorld;
	float3 p1_110 = U * cp[p210].PositionWorld + V * cp[p120].PositionWorld + W * cp[p111].PositionWorld;
	float3 p1_101 = U * cp[p201].PositionWorld + V * cp[p111].PositionWorld + W * cp[p102].PositionWorld;
	float3 p1_020 = U * cp[p120].PositionWorld + V * cp[p030].PositionWorld + W * cp[p021].PositionWorld;
	float3 p1_011 = U * cp[p111].PositionWorld + V * cp[p021].PositionWorld + W * cp[p012].PositionWorld;
	float3 p1_002 = U * cp[p102].PositionWorld + V * cp[p012].PositionWorld + W * cp[p003].PositionWorld;

	float3 p2_100 = U * p1_200 + V * p1_110 + W * p1_101;
	float3 p2_010 = U * p1_110 + V * p1_020 + W * p1_011;
	float3 p2_001 = U * p1_101 + V * p1_011 + W * p1_002;

	float3 position = U * p2_100 + V * p2_010 + W * p2_001;

	//float3 tangent = p2_010 - p2_100;
	//float3 bitangent = p2_001 - p2_100;

	//tangent = normalize(tangent);
	//bitangent = normalize(bitangent);
	float3 normal = pow(U,2) * pc.NormalCP[n200] +
		pow(V,2) * pc.NormalCP[n020] +
		pow(W,2) * pc.NormalCP[n002] + 
		U * V * pc.NormalCP[n110] + 
		U * W * pc.NormalCP[n101] +
		V * W * pc.NormalCP[n011];

	dout.View = vecEye.xyz - position.xyz;
	dout.PositionH = mul(float4(position,1.0f),ViewProjection);
	dout.Normal = normalize(normal);
	dout.Tangent = float3(0,1,0);
	dout.Bitangent = float3(0,1,0);

	dout.MainTexCoord = cp[p300].MainTexCoord * U + cp[p030].MainTexCoord * V + cp[p003].MainTexCoord * W;

	return dout;
}

Off topic: I often actually marvel at the output of erronuous renderings, e.g. the third screenshot. Maybe I should start a thread about "The beauty of rendering glitches" .

Heh, I guess cold hard mathematics brings all the 'honesty' a piece of art could possibly require. :)

What I'm going to do next is go ahead and calculate the 3 full edge normals in the domain shader, and interpolate those somehow with the 'real' normals weighted toward the center and the edge normals taking over as you approach an edge, just to see if I can solve the problem.

The only problem is that the domain shader is eventually going to have to do a lot of heavy lifting yet with normal map, displacement map and color texturing, in addition to carving river channels and road bumps. I hope the bottleneck isn't to tight.