1
$\begingroup$

A friend of mine brought this problem to my attention. I've been playing with it and I think I'm on the right path- I'm not sure what's wrong with my code- Mathematica says that there is a missing brace, but I can't find it. Cross product method:

Clear[AA, BB, CC, P, AB, AC, AP];
AA = {0, 0, 0}; BB = {0, 1, 0}; CC = {1, 0, 0};
PP = {1/4, 1/4, 0};
AB = BB - AA;
AP = PP - AA;
AC = CC - AA;
ZP = Cross[AB, AP]
ZC = Cross[AB, AC]
If[ZP[[3]]*ZC[[3]] > 0, Print["P and C are on the same side of AB"], 
Print["P and C are not on the same side of AB"]]

Dot Product method: AP.AB=uAB.AB+vAC.AB AP.AC=uAB.AC+vAC.AC

Finite Element Soln using Cross product:

gridsize = 0.05;
grid = Flatten[
Join[Transpose[
Outer[List, Range[0, 3, gridsize], Range[0, 1, gridsize], {0}]], 
Transpose[
Outer[List, Range[2, 3, gridsize], Range[1, 3, gridsize], {0}]]],
2];

nodes = {{0, 0, 0}, {2, 0, 0}, {3, 0, 0}, {0, 1, 0}, {2, 1, 0}, {3, 1,
0}, {2, 3, 0}, {3, 3, 0}};
elements = {{1, 2, 4}, {2, 5, 4}, {2, 6, 5}, {2, 3, 6}, {5, 6, 7}, {6,
8, 7}};
Do[(* This loop goes through each point in the grid with index i *)
point = grid[[i]];
Do[(*This loop goes through each triangle with index j *)

vertices = nodes[[elements[[j]]]];
Flag = True;
Do[(* This loop goes through each side of the triangle with index \k *)
AB = vertices[[Mod[k + 1, 3, 1]]] - vertices[[Mod[k, 3, 1]]];(* 
vector from k to k+1 *)

AP = point - vertices[[Mod[k, 3, 1]]]; (* vector from k to point *)

AC = vertices[[Mod[k + 2, 3, 1]]] - vertices[[Mod[k, 3, 1]]];(* 
vector from k to k+2 *)
Z1 = Cross[AB, AP];
Z2 = Cross[AB, AC];
If[Z1[[3]]*Z2[[3]] >= 0, , Flag = False],
{k, 1, 3}]; 

(* If the point passes the test, then we should set the value of the       solution with the appropriate value for that triangle. *)

If[Flag, 
grid[[i, 3]] = a[[j]]*point[[1]] + b[[j]]*point[[2]] + c[[j]]],
{j, 1, Length[elements]}],
{i, 1, Length[grid]}];

My Soln:

Do[(* This loop goes through each point in the grid with index i *)

point = grid[[i]];
Do[(*This loop goes through each triangle with index j *)

vertices = nodes[[elements[[j]]]];
Flag = True;
Do[(* This loop goes through each side of the triangle with index k \*)
AB = vertices[[Mod[k + 1, 3, 1]]] - vertices[[Mod[k, 3, 1]]];(* 
vector from k to k+1 *)
AP = point - vertices[[Mod[k, 3, 1]]]; (* 
vector from k to point *)

AC = vertices[[Mod[k + 2, 3, 1]]] - vertices[[Mod[k, 3, 1]]];(* 
vector from k to k+2 *)
(* this is where I try to implement Dot Product in my soln *)

soln = Solve[
AP.AB == u*AB.AB + v*AC.AB && AP.AC == u*AB.AC + v*AC.AC, {u, v}]
{u, v} = {u, v} /. soln[[1]];
If[u >= 0 && v >= 0 && u + v <= 1, 
Print["The point ", PP, " is inside the triangle"], 
Print["The point ", PP, " is outside the triangle"]];
(* If the point passes the test, 
then we should set the value of the solution with the appropriate \
value for that triangle. *)

If[Flag, 
grid[[i, 3]] = a[[j]]*point[[1]] + b[[j]]*point[[2]] + c[[j]]],
{j, 1, Length[elements]}],
{i, 1, Length[grid]}];

Thank you for your time, any advice would be greatly appreciated!

$\endgroup$
  • 1
    $\begingroup$ the line soln = Solve... needs a semicolon. Is there a specific question here? "Please debug my code" is overly broad. $\endgroup$ – george2079 Oct 5 '16 at 11:18