Does the following work for you?

f1[x_, y_] := Sin[x y]
vecF = Normalize[{D[f1[x, y], x], D[f1[x, y], y], -1}];
scalarR = f1[x, y] - z;
vecData = Table[vecF, {x, 0, 3, 0.5}, {y, 0, 3, 0.5}, {z, -1, 1, 0.5}];

For this, I created a regularly spaced array of vectors in 3D. The iteration on z is for that purpose. The scalar region function, scalarR, is used as the rule for RegionFunction. BTW, it seems like you have a redundant Normalize in your original code, so I left it out.

vG = ListVectorPlot3D[vecData, DataRange -> {{0, 3}, {0, 3}, {-1, 1}}, 
  VectorPoints -> 20, 
  RegionFunction -> Function[{x, y, z}, -0.2 <= scalarR <= 0.2]]

cG = ContourPlot3D[scalarR == 0, {x, 0, 3}, {y, 0, 3}, {z, -1, 1}, 
  Mesh -> None, ContourStyle -> Opacity[0.5, Green]]

Show[vG, cG]

Dear David,

I think I understand better now what you mean. That behaviour is odd and has been discussed before:

http://forums.wolfram.com/mathgroup/archive/2009/Feb/msg00770.html

http://mathematica.stackexchange.com/questions/26519/making-a-listvectorplot3d-from-data-in-an-external-file

http://mathematica.stackexchange.com/questions/27921/problem-with-listvectorplot3d

Is it true that you want something like this:

Show[Graphics3D[{Thick, Arrow[{#[[1]], #[[1]] + #[[2]]}]} & /@ nvf], 
 ListPlot3D[nvf[[All, 1]]]]

I couldn't make it work with ListVectorPlot3D. This, however, is similar and does work

    scalarField = x^2 - y^2 - z;
    vectorField = D[scalarField, {{x, y, z}}];
    data = Table[vectorField, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}];

v = ListVectorPlot3D[data, DataRange -> {{-2, 2}, {-2, 2}, {-2, 2}}, 
  VectorPoints -> 25, VectorScale -> {0.1, Scaled[0.5]}, 
  RegionFunction -> Function[{x, y, z}, -0.1 <= scalarField <= 0.1]]

c = ContourPlot3D[
  scalarField == 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None, 
  ContourStyle -> Opacity[0.5, Green]]

Show[v, c]

I know that this does not really answer the question of why ListVectorPlot3D does not do it.

Cheers, Marco

Here is the image. The plot plots the surface on which the 3D points for the vector field are generated, and the points as extracted from the field list.

Thanks to you both. But I still find this peculiar. I think my nvf has the format specified in the second entry of the docs for ListVectorPlot3D, which is { {{x1,y1,y2},{vx1,vy1,vz1}}, . . . }, and has rank 3. -- which is just a list of field points and vector values.

The unflattened nvecs does plot as you show, Marco, but I think all the field points should lie on the surface f[x,y] = Sin[x y]. Here I plot the field points of nvf on the surface:

Show[Plot3D[
  f[x, y], {x, 0, 3}, {y, 0, 3}], {PointSize[0.02], 
   Point[nvf[[All, 1]]]} // Graphics3D]

But the field points for the unflattened nvecs don't lie on the surface where they were generated, but seem to be on a grid, as though Mathemetica were interpreting it as some array, which is like the first entry in the documentation for ListVectorPlot3D.