/* phi_tri_cm.c  compute coefficients of phi_tri.c functions
 *               for a specific triangle
 *               call three times:
 *               phi_tri_cm22(x1, y1, x2, y2, x3, y3, debug, cm1);
 *               phi_tri_cm22(x2, y2, x3, y3, x1, y1, debug, cm2);
 *               phi_tri_cm22(x3, y3, x1, y1, x2, y2, debug, cm3);
 *
 *   see phi_tri.h for numbering convention cm12 cm22 cm32
 */

#include <stdio.h>
#include <stdlib.h>
#include "phi_tri.h"
void simeq(int n, double A[], double y[], double x[]); /* solve  Ax=y */

void phi_tri_cm12(double x1, double y1, double x2, double y2,
                  double x3, double y3, int kdebug, double cm[])
{
  /* specific to first order in two dimensions for triangles */
  /* phi12(x,y)=cm[0]+cm[1]*x+cm[2]*y */ 
  double am[9], ym[3];
  double sum;
  int i, j;

  am[0]=1.0;
  am[1]=x1;
  am[2]=y1;
  am[3]=1.0;
  am[4]=x2;
  am[5]=y2;
  am[6]=1.0;
  am[7]=x3;
  am[8]=y3;

  ym[0]=1.0; /* 1.0 for first, all others zero */
  ym[1]=0.0;
  ym[2]=0.0;
  simeq(3,am,ym,cm);
  if(kdebug)
  {
    printf("cm matrix  of  am cm = ym   solve for cm \n");
    for(i=0; i<3; i++) printf(" %8.2f",cm[i]);
    printf("\n");
  }
  /* check solution */
  for(i=0; i<3; i++)
  {
    sum=0.0;
    for(j=0; j<3; j++)
      sum=sum+am[i*3+j]*cm[j];
    if(i==0 && abs(sum-1.0)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
    if(i!=0 && abs(sum)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
   }
   if(abs(phi12(cm,x1,y1)-1.0)>1.0e-6)
      printf("BAD cm x1=%f, y1=%f \n", x1, y1);
   if(abs(phi12(cm,x2,y2))>1.0e-6)
      printf("BAD cm x2=%f, y2=%f \n", x2, y2);
   if(abs(phi12(cm,x3,y3))>1.0e-6)
      printf("BAD cm x3=%f, y3=%f \n", x3, y3);
} /* end phi_tri_cm12 */

void phi_tri_cm22(double x1, double y1, double x2, double y2,
              double x3, double y3, int kdebug, double cm[])
{
  /* specific to second order in two dimensions for triangles */
  /* phi22(x,y)=cm[0]+cm[1]*x+cm[2]*y+cm[3]*x*x+cm[4]*x*y+cm[5]*y*y */ 
  double x12, y12, x13, y13, x23, y23;
  double am[36], ym[6];
  double sum;
  int i, j;

  am[0]=1.0;
  am[1]=x1;
  am[2]=y1;
  am[3]=am[1]*am[1];
  am[4]=am[1]*am[2];
  am[5]=am[2]*am[2];
  am[6]=1.0;
  am[7]=x2;
  am[8]=y2;
  am[9]=am[7]*am[7];
  am[10]=am[7]*am[8];
  am[11]=am[8]*am[8];
  am[12]=1.0;
  am[13]=x3;
  am[14]=y3;
  am[15]=am[13]*am[13];
  am[16]=am[13]*am[14];
  am[17]=am[14]*am[14];
  x12=(x1+x2)/2.0;
  y12=(y1+y2)/2.0;
  am[18]=1.0;
  am[19]=x12;
  am[20]=y12;
  am[21]=am[19]*am[19];
  am[22]=am[19]*am[20];
  am[23]=am[20]*am[20];
  x13=(x1+x3)/2.0;
  y13=(y1+y3)/2.0;
  am[24]=1.0;
  am[25]=x13;
  am[26]=y13;
  am[27]=am[25]*am[25];
  am[28]=am[25]*am[26];
  am[29]=am[26]*am[26];
  x23=(x2+x3)/2.0;
  y23=(y2+y3)/2.0;
  am[30]=1.0;
  am[31]=x23;
  am[32]=y23;
  am[33]=am[31]*am[31];
  am[34]=am[31]*am[32];
  am[35]=am[32]*am[32];

  ym[0]=1.0; /* 1.0 for first, all others zero */
  ym[1]=0.0;
  ym[2]=0.0;
  ym[3]=0.0;
  ym[4]=0.0;
  ym[5]=0.0;
  simeq(6,am,ym,cm);
  if(kdebug)
  {
    printf("cm matrix  of  am cm = ym   solve for cm \n");
    for(i=0; i<6; i++) printf(" %8.2f",cm[i]);
    printf("\n");
  }
  /* check solution */
  for(i=0; i<6; i++)
  {
    sum=0.0;
    for(j=0; j<6; j++)
      sum=sum+am[i*6+j]*cm[j];
    if(i==0 && abs(sum-1.0)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
    if(i!=0 && abs(sum)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
   }
   if(abs(phi22(cm,x1,y1)-1.0)>1.0e-6)
      printf("BAD cm x1=%f, y1=%f \n", x1, y1);
   if(abs(phi22(cm,x2,y2))>1.0e-6)
      printf("BAD cm x2=%f, y2=%f \n", x2, y2);
   if(abs(phi22(cm,x3,y3))>1.0e-6)
      printf("BAD cm x3=%f, y3=%f \n", x3, y3);
   if(abs(phi22(cm,x12,y12))>1.0e-6)
      printf("BAD cm x12=%f, y12=%f \n", x12, y12);
   if(abs(phi22(cm,x13,y13))>1.0e-6)
      printf("BAD cm x13=%f, y13=%f \n", x13, y13);
   if(abs(phi22(cm,x23,y23))>1.0e-6)
      printf("BAD cm x23=%f, y23=%f \n", x23, y23);
} /* end phi_tri_cm22 */

void phi_tri_cm32(double x1, double y1, double x2, double y2,
              double x3, double y3, int kdebug, double cm[])
{
  /* specific to third order in two dimensions for triangles */
  /* phi32(x,y)=cm[0]+cm[1]*x+cm[2]*y+cm[3]*x*x+cm[4]*x*y+cm[5]*y*y+ */
  /*            cm[6]*x*x*x+cm[7]*x*x*y+cm[8]*x*y*y+cm[9]*y*y*y */ 
  double x112, y112, x122, y122, x113, y113,
               x133, y133, x223, y223, x233, y233, x123, y123;
  double am[100], ym[10];
  double sum;
  int i, j;

  am[0]=1.0;
  am[1]=x1;
  am[2]=y1;
  am[3]=am[1]*am[1];
  am[4]=am[1]*am[2];
  am[5]=am[2]*am[2];
  am[6]=am[3]*am[1];
  am[7]=am[3]*am[2];
  am[8]=am[1]*am[5];
  am[9]=am[5]*am[2];

  am[10]=1.0;
  am[11]=x2;
  am[12]=y2;
  am[13]=am[11]*am[11];
  am[14]=am[11]*am[12];
  am[15]=am[12]*am[12];
  am[16]=am[13]*am[11];
  am[17]=am[13]*am[12];
  am[18]=am[11]*am[15];
  am[19]=am[15]*am[12];

  am[20]=1.0;
  am[21]=x3;
  am[22]=y3;
  am[23]=am[21]*am[21];
  am[24]=am[21]*am[22];
  am[25]=am[22]*am[22];
  am[26]=am[23]*am[21];
  am[27]=am[23]*am[22];
  am[28]=am[21]*am[25];
  am[29]=am[25]*am[22];

  x112=(2.0*x1+x2)/3.0;
  y112=(2.0*y1+y2)/3.0;
  am[30]=1.0;
  am[31]=x112;
  am[32]=y112;
  am[33]=am[31]*am[31];
  am[34]=am[31]*am[32];
  am[35]=am[32]*am[32];
  am[36]=am[33]*am[31];
  am[37]=am[33]*am[32];
  am[38]=am[31]*am[35];
  am[39]=am[35]*am[32];

  x122=(x1+2.0*x2)/3.0;
  y122=(y1+2.0*y2)/3.0;
  am[40]=1.0;
  am[41]=x122;
  am[42]=y122;
  am[43]=am[41]*am[41];
  am[44]=am[41]*am[42];
  am[45]=am[42]*am[42];
  am[46]=am[43]*am[41];
  am[47]=am[43]*am[42];
  am[48]=am[41]*am[45];
  am[49]=am[45]*am[42];

  x113=(2.0*x1+x3)/3.0;
  y113=(2.0*y1+y3)/3.0;
  am[50]=1.0;
  am[51]=x113;
  am[52]=y113;
  am[53]=am[51]*am[51];
  am[54]=am[51]*am[52];
  am[55]=am[52]*am[52];
  am[56]=am[53]*am[51];
  am[57]=am[53]*am[52];
  am[58]=am[51]*am[55];
  am[59]=am[55]*am[52];

  x133=(x1+2.0*x3)/3.0;
  y133=(y1+2.0*y3)/3.0;
  am[60]=1.0;
  am[61]=x133;
  am[62]=y133;
  am[63]=am[61]*am[61];
  am[64]=am[61]*am[62];
  am[65]=am[62]*am[62];
  am[66]=am[63]*am[61];
  am[67]=am[63]*am[62];
  am[68]=am[61]*am[65];
  am[69]=am[65]*am[62];

  x223=(2.0*x2+x3)/3.0;
  y223=(2.0*y2+y3)/3.0;
  am[70]=1.0;
  am[71]=x223;
  am[72]=y223;
  am[73]=am[71]*am[71];
  am[74]=am[71]*am[72];
  am[75]=am[72]*am[72];
  am[76]=am[73]*am[71];
  am[77]=am[73]*am[72];
  am[78]=am[71]*am[75];
  am[79]=am[75]*am[72];

  x233=(x2+2.0*x3)/3.0;
  y233=(y2+2.0*y3)/3.0;
  am[80]=1.0;
  am[81]=x233;
  am[82]=y233;
  am[83]=am[81]*am[81];
  am[84]=am[81]*am[82];
  am[85]=am[82]*am[82];
  am[86]=am[83]*am[81];
  am[87]=am[83]*am[82];
  am[88]=am[81]*am[85];
  am[89]=am[85]*am[82];

  x123=(x1+x2+x3)/3.0;
  y123=(y1+y2+y3)/3.0;
  am[90]=1.0;
  am[91]=x123;
  am[92]=y123;
  am[93]=am[91]*am[91];
  am[94]=am[91]*am[92];
  am[95]=am[92]*am[92];
  am[96]=am[93]*am[91];
  am[97]=am[93]*am[92];
  am[98]=am[91]*am[95];
  am[99]=am[95]*am[92];

  ym[0]=1.0; /* 1.0 for first, all others zero */
  ym[1]=0.0;
  ym[2]=0.0;
  ym[3]=0.0;
  ym[4]=0.0;
  ym[5]=0.0;
  ym[6]=0.0;
  ym[7]=0.0;
  ym[8]=0.0;
  ym[9]=0.0;
  simeq(10,am,ym,cm);
  if(kdebug)
  {
    printf("cm matrix  of  am cm = ym   solve for cm \n");
    for(i=0; i<5; i++) printf(" %8.2f",cm[i]);
    printf("\n");
    for(i=5; i<10; i++) printf(" %8.2f",cm[i]);
    printf("\n");
  }
  /* check solution */
  for(i=0; i<10; i++)
  {
    sum=0.0;
    for(j=0; j<10; j++)
      sum=sum+am[i*10+j]*cm[j];
    if(i==0 && abs(sum-1.0)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
    if(i!=0 && abs(sum)>1.0e-6)
        printf("BAD cm i=%d, sum=%f\n", i, sum);
   }
   if(abs(phi32(cm,x1,y1)-1.0)>1.0e-6)
      printf("BAD cm x1=%f, y1=%f \n", x1, y1);
   if(abs(phi32(cm,x2,y2))>1.0e-6)
      printf("BAD cm x2=%f, y2=%f \n", x2, y2);
   if(abs(phi32(cm,x3,y3))>1.0e-6)
      printf("BAD cm x3=%f, y3=%f \n", x3, y3);

   if(abs(phi32(cm,x112,y112))>1.0e-6)
      printf("BAD cm x112=%f, y112=%f \n", x112, y112);
   if(abs(phi32(cm,x122,y122))>1.0e-6)
      printf("BAD cm x122=%f, y122=%f \n", x122, y122);
   if(abs(phi32(cm,x113,y113))>1.0e-6)
      printf("BAD cm x113=%f, y113=%f \n", x113, y113);
   if(abs(phi32(cm,x133,y133))>1.0e-6)
      printf("BAD cm x133=%f, y133=%f \n", x133, y133);
   if(abs(phi32(cm,x223,y223))>1.0e-6)
      printf("BAD cm x223=%f, y223=%f \n", x223, y223);
   if(abs(phi32(cm,x233,y233))>1.0e-6)
      printf("BAD cm x233=%f, y233=%f \n", x233, y233);
   if(abs(phi32(cm,x123,y123))>1.0e-6)
      printf("BAD cm x123=%f, y123=%f \n", x123, y123);
} /* end phi_tri_cm32 */