I have two sets of inputs that are to be fed to a simple function. The inputs are essentially a set of triangles and the normal for each triangle. The function will compute the gradient about the points (shown in blue). The notebook can be downloaded from the following Dropbox
Link
https://www.dropbox.com/s/yas32nfccd2dzj4/debug%20code.nb?dl=0
The two sets of triangles are supposed to be similar (only separated by a translation in space). I compute the area gradient about point pt1
grad1 = Block[{ptTri, normal, cross, target, facept,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, source = pt1},
Do[
ptTri = opentr1[[i]];
normal = normOpentr1[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross), {i, 1, Length@normOpentr1}];
Do[
ptTri = closedtri1[[j]];
normal = normClosedtr1[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross), {j, 1, Length@normClosedtr1}];
0.7*closedS + 1*openS
]
(*{0.110728, 0.0466838, 0.752509}*)
Likewise if I compute the area gradient about the second point pt2
I get the same answer.
grad2 = Block[{ptTri, normal, cross, target, facept,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, source = pt2},
Do[
ptTri = opentr2[[i]];
normal = normOpentr2[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross);, {i, 1, Length@normOpentr2}];
Do[
ptTri = closedtr2[[j]];
normal = normClosedtr2[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross),{j, 1, Length@normClosedtr2}];
0.7*closedS + 1*openS
]
(*{0.110728, 0.0466838, 0.752509}*)
The two answers are quite similar as the difference between grad1
and grad2
is {1.83187*10^-15, -5.55112*10^-16,4.44089*10^-16}
HOWEVER
the same code when under Compile
does not give the same answers for the two datasets.
With[{epcc = 0.7, epco = 1.},
surfaceGradFn2 =
Compile[{{point, _Real, 1}, {opentr, _Real, 3}, {normalO, _Real,2}, {closedtr, _Real, 3}, {normalC,
_Real, 2}},
Block[{ptTri, source = point, normal, target, facept, cross,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, nO = normalO,
nC = normalC, OT = opentr, CT = closedtr},
Do[
ptTri = OT[[i]];
normal = nO[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross);
, {i, 1, Length@nO}];
Do[
ptTri = CT[[j]];
normal = nC[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross);
, {j, 1, Length@nC}];
epcc*closedS + epco*openS
], CompilationTarget -> "C"]
];
Though the structure of the code is the same I guess, the outputs are completely different.
surfaceGradFn2[pt1, opentr1, normOpentr1, closedtri1, normClosedtr1]
(*{-0.0157582, 0.386426, 0.582118} -> this answer is now different but I gave the same inputs as before *)
surfaceGradFn2[pt2, opentr2, normOpentr2, closedtr2, normClosedtr2]
(* {0.110728, 0.0466838, 0.752509} -> answer is same as before with the same inputs *)
What am I doing wrong here? Why does not the compiled code work while the non-compiled version works? I wish to have the compiled version of the code working.
also, please ignore what the code intends to do. If you may, kindly look at why the compiled output is different.Thanks and grateful for your help.