Function plot_P2_isoline_body

plugin/seq/plotPDF.cpp:2527–2698

	// X' = - sqrt( a(Y')^2 + b )
	trackHyperbolaCore( Cy, Cx, -1, a, b, Zy_minus, Vy, Vx );

    }
    return;
}

void plot_P2_isoline_body( std::stringstream &Content, const Fem2D::Mesh &Th, const KN<double> &f_P2_,
			   const std::vector<double> &isoline_val, const double cmin, const double cmax,
			   const KNM<double> &palette, const double scale, const double ar,
			   const double x0, const double y0,
			   const int marginl, const int marginb,
			   const bool monochrome, const bool logscale,
			   const int NISOLINES, const double linewidth )
{
    const double EPS = 1e-10;

    std::stringstream &st = Content;

    const int NQ = 6; // number of coefficients in quadratic polynomial
    const int &nTriangles = Th.nt;

    st << "q\n";
    st << linewidth << " w\n"; // setlinewidth
    st << "1 0 0 1 " << PLOTPDFVAR::PADDING+marginl << " " << PLOTPDFVAR::PADDING+marginb << " cm\n";
	
    for(int k = 0; k < nTriangles; k++){

	const int &v0 = Th(k,0); const int &v1 = Th(k,1); const int &v2 = Th(k,2);

	const double vx[] = { Th(v0).x,  Th(v1).x,  Th(v2).x };
	const double vy[] = { Th(v0).y,  Th(v1).y,  Th(v2).y };

	// f_P2[i] = phi(v[i]), f_P2[i+3] = phi(e[i]), i=0,1,2.
	const double *const f_P2 = f_P2_ + k*NQ;
	// in vector2flow
	//const double f_P2[NQ] = { f2[3*k], f2[3*k+1], f2[3*k+2], f2m[3*k], f2m[3*k+1], f2m[3*k+2] };

	// quadratic polynomials : ax^2 + bxy + cy^2 + dx + ey + f
	double phi[ NQ ];
	const double &a = phi[0]; const double &b = phi[1]; const double &c = phi[2];
	const double &d = phi[3]; const double &e = phi[4]; const double &f = phi[5];

	findQuadraticPolynomial(phi, vx, vy, f_P2);
	
	const bool isLinear = (fabs(a) + fabs(b) + fabs(c) < EPS * (fabs(d)+fabs(e) + fabs(f)));

	// normalize Quadratic Polynomial
	double PHI[9];
	findCanonicalForm( PHI, phi );

	// PHI(X,Y) = lambda1*X*X + lambda2*Y*Y + D*X + E*Y + f
        //          = lambda1 * (X + D/(2*lambda1))^2 + lambda2 * (Y + E/(2*lambda2))^2
        //            + ( -D*D/(4*lambda1) - E*E/(4*lambda2) + f)
        //          = lambda1*(X + D/(2*lambda1))^2 + lambda2*(Y + E/(2*lambda2))^2 + F
	const double &lambda1 = PHI[0]; const double &lambda2 = PHI[1];
	const double &D = PHI[6];   const double &E = PHI[7];   double &F = PHI[8];

	// lambda1*X^2 + lambda2*Y^2 + D*X + E*Y + f = value,
	// if lambda2 == E == 0, then lambda1*X^2 + D*X + f = value, i.e. X=const