(* 
      Expand for functions for dealing with x/
          y coeffs for triangles (Eqn1 in Igarashi05) *)

Rot90 = {{0, 1}, {-1, 
            0}};

    FindXY[v0_,  v1_,  v2_]  :=  
        Module[
{v01,  v01N,  vLocal, x, 
            y}, 
v01  =  
            v1 - v0; 
v01N  =  {{0, 1}, {-1, 0}} . 
              v01; 
vLocal  =  
            v2  -  v0; 
  
x  =  (vLocal . 
                  v01) / (v01 . 
                  v01); 
y  =  (vLocal . 
                  v01N)/(v01N . v01N); 
{x, 
            y}
];

    ExpandXY[v0_, v1_, v2_, xy_] :=  
        Module [
{},  
          v0  +  xy[[1]]*(v1 - v0)  +  
            xy[[2]]*
              Rot90 . (v1 - v0)
];

    CheckXY[v0_, v1_, v2_, xy_] :=  
        Module [
{},  
          v2  ==  
            v0  +  xy[[1]]*(v1 - v0)  +  
              xy[[2]]*
                Rot90 . (v1 - v0)
];

\(\(\( (*\ 
      declare\ relevant\ variables . \ 
          here\ v\_i\ is\ what\ we\ are\ solving\ for\ *) \)\(
\)\(\(v\_0 = {v\_\(0, x\), 
          v\_\(0, y\)};\)

    \(v\_1 = {v\_\(1, x\), v\_\(1, y\)};\)

    \(v\_2 = {v\_\(2, x\), v\_\(2, y\)};\)

    \(\(v\^d\)\_0 = {\(v\^d\)\_\(0, x\), \(v\^d\)\_\(0, y\)};\)

    \(\(v\^d\)\_1 = {\(v\^d\)\_\(1, x\), \(v\^d\)\_\(1, y\)};\)

    \(\(v\^d\)\_2 = {\(v\^d\)\_\(2, x\), \(v\^d\)\_\(2, y\)};\)

    \(v\&^\)\_0\  = \ {\(v\&^\)\_\(0, x\), \(v\&^\)\_\(0, y\)}; \ \ \ \ \ \ \ \ \  (*\ these\ are\ the\ points\ from\ step\ 1\ of\ deformation, \ 
      which\ we\ know\ *) 
\(v\&^\)\_1\  = \ {\(v\&^\)\_\(1, x\), \(v\&^\)\_\(1, y\)};

    \(\(v\&^\)\_2\  = \ {\(v\&^\)\_\(2, x\), \(v\&^\)\_\(2, y\)};\)

 (*\ 
      Equations\ for\ the\ 3\ vertices\ of\ the\ triangle\ *) 

    \(\(v\^d\)\_2\  = \ 
        v\_0\  + \ \(x\_\(0, 1\)\) \((v\_1 - v\_0)\) + \(y\_\(0, 1\)\) \((\ 
              Rot90\  . \ \((v\_1 - v\_0)\))\);\)

    \(\(v\^d\)\_1\  = \ 
        v\_2\  + \ \(x\_\(2, 0\)\) \((v\_0 - v\_2)\) + \(y\_\(2, 0\)\) \((\ 
              Rot90\  . \ \((v\_0 - v\_2)\))\);\)

    \(\(v\^d\)\_0\  = \ 
        v\_1\  + \ \(x\_\(1, 2\)\) \((v\_2 - v\_1)\) + \(y\_\(1, 2\)\) \((\ 
              Rot90\  . \ \((v\_2 - 
                    v\_1)\))\);\)

 (*\ 
      error\ equations\ for\ each\ point\ in\ triangle . \ 
          rewrite\ v2\ in\ terms\ of\ v0\ and\ v1\ \ *) 

    \(v2Rules\  = \ {\ 
          v\_\(2, x\) -> \ 
            v\_\(0, x\) + \((\(-v\_\(0, x\)\) + 
                    v\_\(1, x\))\)\ x\_\(0, 1\) + \((\(-v\_\(0, y\)\) + 
                    v\_\(1, y\))\)\ y\_\(0, 1\), \ \ v\_\(2, y\) -> \ 
            v\_\(0, y\) + \((\(-v\_\(0, y\)\) + 
                    v\_\(1, y\))\)\ x\_\(0, 1\) + \((v\_\(0, x\) - 
                    v\_\(1, x\))\)\ y\_\(0, 1\)};\)

    \(Ev2\  = \ \((\(v\^d\)\_2 - \(v\&^\)\_2)\) . \((\(v\^d\)\_2 - \(v\&^\)\_2)\);\)

    \(Ev1\  = \ \((\(v\^d\)\_1 - \(v\&^\)\_1)\) . \((\(v\^d\)\_1 - \(v\&^\)\_1)\)\  /. \ v2Rules;\)

    \(Ev0\  = \ \((\(v\^d\)\_0 - \(v\&^\)\_0)\) . \((\(v\^d\)\_0 - \(v\&^\)\_0)\)\  /. \ v2Rules;\)\)\)\)

\(\(\(
\)\( (*\ 
      Now\ set\ up\ a\ test\ case\ *) \)\(
\)\(
\)\(P0 = {\(-1\), \(-1\)}; \ \ P1\  = \ {0, 
        0}; \ \ P2\  = \ {0, \(-1\)};

    \(XY\_\(0, 1\) = \ FindXY[P0, P1, P2];\)

    \(XY\_\(2, 0\) = \ FindXY[P2, P0, \ P1];\)

    \(XY\_\(1, 2\) = \ 
        FindXY[P1, P2, P0];\)
 (*Print["XYChecks: ", 
        CheckXY[\ P0, \ P1, \ P2, \ XY\_\(0, 1\)], 
        CheckXY[\ P2, \ P0, \ P1, \ XY\_\(2, 0\)], \ \ \ CheckXY[P1, \ P2, \ 
          P0, \ XY\_\(1, 2\)]]*) 


    \(FullRules\  = \ {v\_\(0, x\) -> \ P0[\([1]\)], \ 
          v\_\(0, y\) -> \ P0[\([2]\)], \ 
          v\_\(1, x\) -> \ P1[\([1]\)], 
          v\_\(1, y\) -> \ P1[\([2]\)], \ 
          v\_\(2, x\) -> \ P2[\([1]\)], 
          v\_\(2, y\) -> \ 
            P2[\([2]\)], 
x\_\(0, 1\) -> \ 
            XY\_\(0, 1\)[\([1]\)], \ 
          y\_\(0, 1\) -> XY\_\(0, 1\)[\([2]\)], \ 
          x\_\(2, 0\) -> XY\_\(2, 0\)[\([1]\)], \ 
          y\_\(2, 0\) -> \ XY\_\(2, 0\)[\([2]\)], \ 
          x\_\(1, 2\) -> \ XY\_\(1, 2\)[\([1]\)], \ 
          y\_\(1, 2\) -> \ XY\_\(1, 2\)[\([2]\)]};\)

    \(XYRules\  = \ {x\_\(0, 1\) -> \ XY\_\(0, 1\)[\([1]\)], \ 
          y\_\(0, 1\) -> XY\_\(0, 1\)[\([2]\)], \ 
          x\_\(2, 0\) -> XY\_\(2, 0\)[\([1]\)], \ 
          y\_\(2, 0\) -> \ XY\_\(2, 0\)[\([2]\)], \ 
          x\_\(1, 2\) -> \ XY\_\(1, 2\)[\([1]\)], \ 
          y\_\(1, 2\) -> \ XY\_\(1, 2\)[\([2]\)]};\)

    \(HatRules\  = \ {\(v\&^\)\_\(0, x\) -> \ 
            P0[\([1]\)], \(v\&^\)\_\(0, y\) -> \ 
            P0[\([2]\)], \ \(v\&^\)\_\(1, x\) -> \ 
            P1[\([1]\)], \(v\&^\)\_\(1, y\) -> \ 
            P1[\([2]\)], \ \(v\&^\)\_\(2, x\) -> \ 
            P2[\([1]\)], \(v\&^\)\_\(2, y\) -> \ 
            P2[\([2]\)]};\)

    Print["EVChecks ShouldBeZero: ", \ 
\(Ev2 /. 
          FullRules\) /. HatRules, " ", \ \(Ev1 /. FullRules\) /. 
        HatRules, " ", \ \(Ev2\  /. \ FullRules\)\  /. \ 
        HatRules]\)\)\)

\(\(\( (*\ Ok, \ 
      now\ solve\ the\ linear\ system\ of\ equations\ *) \)\(
\)\(
\)\(\(EF\  = \ \ Ev0\  + \ 
          Ev1\  + \ Ev2;\)

    Print["ShouldBeZeroish: ", \(EF\  /. \ FullRules\)\  /. \ 
        HatRules]


    \(EqnSet\  = \ {
D[\(EF\  /. \ HatRules\)\  /. \ 
                XYRules, \ v\_\(0, x\)]\  == \ 
            0, 
D[\(EF\  /. \ HatRules\)\  /. \ XYRules, \ 
              v\_\(0, y\)]\  == \ 
            0, 
D[\(EF\  /. \ HatRules\)\  /. \ XYRules, \ 
              v\_\(1, x\)]\  == \ 
            0, 
D[\(EF\  /. \ HatRules\)\  /. \ XYRules, \ 
              v\_\(1, y\)]\  == \ 0\ \ };\)

    Solve[
      EqnSet, {v\_\(0, x\), \ v\_\(0, y\), \ \ v\_\(1, x\), \ 
        v\_\(1, y\)}]

    \)\)\)

(* 
      Sanity check  -   should get same result when differentiating component - 
        wise *) 
EqnSet  =  {
(D[Ev0, (Subscript[v, 0, x])]  +   D[Ev1, (Subscript[v, 0, x])]  +  
                    D[Ev2, (Subscript[v, 0, x])] /.  HatRules  /.  
                XYRules)  ==  
            0, 
(D[Ev0, (Subscript[v, 0, y])]  +   D[Ev1, 
                      (Subscript[v, 0, y])]  +  D[Ev2, (Subscript[v, 0, y])] /.  
                  HatRules  /.  XYRules)  ==  
            0, 
(D[Ev0, (Subscript[v, 1, x])]  +   D[Ev1, 
                      (Subscript[v, 1, x])]  +  D[Ev2, (Subscript[v, 1, x])] /.  
                  HatRules  /.  XYRules)  ==  
            0, 
(D[Ev0, (Subscript[v, 1, y])]  +   D[Ev1, 
                      (Subscript[v, 1, y])]  +  D[Ev2, (Subscript[v, 1, y])] /.  
                  HatRules  /.  XYRules)  ==  
            0  };

    Solve[
      EqnSet, {(Subscript[v, 0, x]),  (Subscript[v, 0, y]),   (Subscript[v, 1, x]),  
        (Subscript[v, 1, y])}]

\(\(\(
\)\( (*\ Ok, \ 
      here\ is\ a\ useful\ test\ for\ coding . \ 
          This\ will\ produce\ the\ F\ and\ C\ matrices\ that
one\ would\ use\ in\ the\ code . \ So, \ by\ setting\ P0, 
      P1, P2\ above, \ 
      you\ can\ check\ that\ your\ matrices
and\ vectors\ are\ correct\ \((\ 
          they\ \ printed\ below\ )\)\ *) \)\(
\)\(
\)\(Print["This is the F matrix (before inverting): "]

    FMatrix\  = \ {
Coefficient[\ 
          D[EF\  /. \ XYRules, \ v\_\(0, x\)]\ , \ {v\_\(0, x\), \ 
            v\_\(0, y\), \ \ v\_\(1, x\), \ 
            v\_\(1, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ v\_\(0, y\)]\ , \ {v\_\(0, x\), \ 
            v\_\(0, y\), \ \ v\_\(1, x\), \ 
            v\_\(1, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ v\_\(1, x\)]\ , \ {v\_\(0, x\), \ 
            v\_\(0, y\), \ \ v\_\(1, x\), \ 
            v\_\(1, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ v\_\(1, y\)]\ , \ {v\_\(0, x\), \ 
            v\_\(0, y\), \ \ v\_\(1, x\), \ 
            v\_\(1, y\)}]\ }

    Print["This is the C matrix: "]

    CVec\  = \ {
Coefficient[\ 
          D[EF\  /. \ XYRules, \ 
            v\_\(0, x\)]\ , \ {\(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ 
            v\_\(0, y\)]\ , \ {\(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ 
            v\_\(1, x\)]\ , \ {\(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)}], 
Coefficient[\ 
          D[EF\  /. \ XYRules, \ 
            v\_\(1, y\)]\ , \ {\(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)}]\ }

    Print["This is C multiplied by the deformed points:"]

    CVec2\  = \ 
      CVec . {\(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)}\  /. \ 
        HatRules

    Print["And here is the solution (Should be {P0,P1}):"]

    LinearSolve[FMatrix, \(-CVec2\)]\)\)\)

\(\(\(
\)\(
\)\( (*\ 
      this\ is\ the\ nasty\ bit\ where\ I\ simplified\ the\ equations\ so\ they\ are\ reasonable\ to\ type\ into\ code . \ 
          Perhaps\ there\ is\ a\ better\ way\ to\ do\ this, \ 
      but\ I\ don' 
          t\ know\ what\ it\ is ... *) \)\(
\)\(
\)\(\(FormulaRules\  = \ {

\((x\_\(1, 2\)\ y\_\(0, 1\) + \((\(-1\) + 
                      x\_\(0, 1\))\)\ y\_\(1, 2\))\) -> \ 
            k\_1, \ \ \ \  (*\ 
            a\ *) 
\((x\_\(1, 2\)\ y\_\(0, 1\) - 
                y\_\(1, 2\) + x\_\(0, 1\)\ y\_\(1, 2\))\) -> \ 
            k\_1, 
\((\(-x\_\(1, 2\)\)\ y\_\(0, 1\) + 
                y\_\(1, 2\) - 
                x\_\(0, 1\)\ y\_\(1, 2\))\) -> \ \(-k\_1\), 
\((\(-x\_\(1, 2\)\)\ y\_\(0, 1\) - \((\(-1\) + 
                      x\_\(0, 1\))\)\ y\_\(1, 2\))\) -> \ \(-k\_1\), \ \  (*\ a\ *) 

\((\(-x\_\(1, 2\)\) + 
                x\_\(0, 1\)\ x\_\(1, 2\) - 
                y\_\(0, 1\)\ y\_\(1, 2\))\) -> \ 
            k\_2, 
\((\((\(-1\) + 
                      x\_\(0, 1\))\)\ x\_\(1, 2\) - 
                y\_\(0, 1\)\ y\_\(1, 2\))\) -> \ 
            k\_2, 
\((\(-1\) - \((\(-1\) + 
                      x\_\(0, 1\))\)\ x\_\(1, 2\) + 
                y\_\(0, 1\)\ y\_\(1, 2\))\) -> \ \((\(-1\)\  - \ 
                k\_2)\), 

\((\(-y\_\(0, 1\)\) + x\_\(2, 0\)\ y\_\(0, 1\) + x\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            k\_3, 
\((y\_\(0, 1\) - 
                x\_\(2, 0\)\ y\_\(0, 1\) - 
                x\_\(0, 1\)\ y\_\(2, 0\))\) -> \ \(-k\_3\), 
\((\((\(-1\) + x\_\(2, 0\))\)\ y\_\(0, 1\) + 
                x\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            k\_3, 

\((\((\(-1\) + 
                      x\_\(0, 1\))\)\ y\_\(0, 1\) + 
                x\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            k\_4, 
\((\((\(-1\) + 
                      x\_\(2, 0\))\)\ x\_\(0, 1\) - 
                y\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            k\_5, 
\((\(-1\) - \((\(-1\) + 
                      x\_\(2, 0\))\)\ x\_\(0, 1\) + 
                y\_\(0, 1\)\ y\_\(2, 0\))\) -> \ \(-1\) - 
              k\_5, 
\((\(-x\_\(0, 1\)\) + 
                x\_\(0, 1\)\ x\_\(2, 0\) - 
                y\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            k\_5, 

\((\(-1\) + 
                x\_\(0, 1\))\) -> \ 
            a, 
\((\((\(-1\) + x\_\(0, 1\))\)\^2 + 
                y\_\(0, 1\)\%2)\) -> \ a\_1, \ \ \  (*\ 
            this\ is\ a\^2\  + \ 
              y\ *) 
\((x\_\(0, 1\)\%2 + 
                y\_\(0, 1\)\%2)\) -> \ 
            a\_2, 
\((\(-x\_\(0, 1\)\%2\) - 
                y\_\(0, 1\)\%2)\) -> \ \(-a\_2\), 

\((\(-1\) + x\_\(2, 0\))\) -> \ 
            b, 
\((\((\(-1\) + x\_\(2, 0\))\)\^2 + 
                y\_\(2, 0\)\%2)\) -> \ 
            b\_1, 

\((x\_\(1, 2\)\%2 + 
                y\_\(1, 2\)\%2)\) -> \ 
            c\_2

};\)


    \(RulesRound2\  = \ {
\((1 + 2\ a\ x\_\(1, 2\) + 
                a\_1\ x\_\(1, 2\)\%2 - 2\ y\_\(0, 1\)\ y\_\(1, 2\) + 
                a\_1\ y\_\(1, 2\)\%2)\) -> \ 
            r\_1, 
\((\(-b\)\ x\_\(0, 1\) - 
                b\_1\ x\_\(0, 1\)\%2 + 
                y\_\(0, 1\)\ \((\(-b\_1\)\ y\_\(0, 1\) + 
                      y\_\(2, 0\))\))\) -> \ 
            r\_2, 
\((\(-a\)\ x\_\(1, 2\) - 
                a\_1\ x\_\(1, 2\)\%2 + 
                y\_\(1, 2\)\ \((y\_\(0, 1\) - 
                      a\_1\ y\_\(1, 2\))\))\) -> \ 
            r\_3, 
\((a\ x\_\(0, 1\) + 
                y\_\(0, 1\)\%2)\) -> \ 
            r\_5, 
\((\(-a\)\ x\_\(0, 1\) - 
                y\_\(0, 1\)\%2)\) -> \ \(-r\_5\), 
\((\(-b\)\ y\_\(0, 1\) - x\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            r\_6, 
\((1 + 2\ b\ x\_\(0, 1\) + 
                b\_1\ x\_\(0, 1\)\%2 + b\_1\ y\_\(0, 1\)\%2 - 
                2\ y\_\(0, 1\)\ y\_\(2, 0\))\) -> \ 
            r\_7, 
\((1 + 2\ b\ x\_\(0, 1\) + 
                b\_1\ x\_\(0, 1\)\%2 + 
                y\_\(0, 1\)\ \((b\_1\ y\_\(0, 1\) - 
                      2\ y\_\(2, 0\))\))\) -> \ 
            r\_7
};\)

    \(CoeffVec\  = \ {v\_\(0, x\), v\_\(0, y\), v\_\(1, x\), 
          v\_\(1, y\), \(v\&^\)\_\(0, x\), \ \(v\&^\)\_\(0, y\), \(v\&^\)\_\(1, x\), \ \(v\&^\)\_\(1, y\), \(v\&^\)\_\(2, x\), \ \(v\&^\)\_\(2, y\)};\)


    Print\ ["First  Row:"]

    Row1A\  = \ \(Coefficient[FullSimplify[D[Ev0, v\_\(0, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row1B\  = \ \(Coefficient[FullSimplify[D[Ev1, v\_\(0, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row1C\  = \ \(Coefficient[FullSimplify[D[Ev2, v\_\(0, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Print\ ["Second  Row:"]

    Row2A\  = \ \(Coefficient[FullSimplify[D[Ev0, v\_\(0, y\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row2B\  = \ \(Coefficient[FullSimplify[D[Ev1, v\_\(0, y\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row2C\  = \ \(Coefficient[FullSimplify[D[Ev2, v\_\(0, y\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Print\ ["Third  Row:"]

    Row3A\  = \ \(Coefficient[FullSimplify[D[Ev0, v\_\(1, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row3B\  = \ \(Coefficient[FullSimplify[D[Ev1, v\_\(1, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row3C\  = \ \(Coefficient[FullSimplify[D[Ev2, v\_\(1, x\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Print\ ["Fourth  Row:"]

    Row4A\  = \ \(Coefficient[FullSimplify[D[Ev0, v\_\(1, y\)]], 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row4B\  = \ \(Coefficient[FullSimplify[D[Ev1, v\_\(1, y\)]], \ 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2

    Row4C\  = \ \(Coefficient[FullSimplify[D[Ev2, v\_\(1, y\)]], 
            CoeffVec]\  /. \ FormulaRules\) /. RulesRound2\)\)\)