Function findQuadraticPolynomial

plugin/seq/plotPDF.cpp:1753–1798

//----------------------------------------------------------------------
// P2 Finite Element
//----------------------------------------------------------------------

#define P2_BARYCENTER
#define P2_EDGE

void findQuadraticPolynomial(double *const phi, const double *const vx, const double *const vy, const double *const f_P2)
{
    const int NQ = 6; // number of unknowns in quadratic polynomial

    double *a[NQ];
    a[0] = new double [NQ*(NQ+1)];  // NQ*(NQ+1) matrix for Gauss Elimination
    for(int i = 1; i < NQ; i++)
	a[i] = a[i-1] + (NQ+1);

#ifdef P2_BARYCENTER
    const double cx = (vx[0]+vx[1]+vx[2]) / 3;
    const double cy = (vy[0]+vy[1]+vy[2]) / 3;
#endif

#if defined(P2_BARYCENTER) && defined(P2_EDGE)
    const double mx[] = { (vx[1]+vx[2])/2, (vx[2]+vx[0])/2, (vx[0]+vx[1])/2 };
    const double my[] = { (vy[1]+vy[2])/2, (vy[2]+vy[0])/2, (vy[0]+vy[1])/2 };
    const double ex[] = { 0.99*mx[0]+0.01*cx, 0.99*mx[1]+0.01*cx, 0.99*mx[2]+0.01*cx };
    const double ey[] = { 0.99*my[0]+0.01*cy, 0.99*my[1]+0.01*cy, 0.99*my[2]+0.01*cy };
#elif defined(P2_BARYCENTER)
    const double ex[] = { (vx[0]+cx)/2, (vx[1]+cx)/2, (vx[2]+cx)/2 };
    const double ey[] = { (vy[0]+cy)/2, (vy[1]+cy)/2, (vy[2]+cy)/2 };
#elif defined(P2_EDGE)
    const double ex[] = { (vx[1]+vx[2])/2, (vx[2]+vx[0])/2, (vx[0]+vx[1])/2 };
    const double ey[] = { (vy[1]+vy[2])/2, (vy[2]+vy[0])/2, (vy[0]+vy[1])/2 };
#endif

    const double x[] = { vx[0], vx[1], vx[2], ex[0], ex[1], ex[2] };
    const double y[] = { vy[0], vy[1], vy[2], ey[0], ey[1], ey[2] };

    for(int i = 0; i < NQ; i++){
	a[i][0] = x[i] * x[i];
	a[i][1] = x[i] * y[i];
	a[i][2] = y[i] * y[i];
	a[i][3] = x[i];
	a[i][4] = y[i];
	a[i][5] = 1;
	a[i][6] = f_P2[i];
    }

    GaussElimination(phi, a, NQ);

    delete [] a[0];