/*
 $Id: feplot2.cc,v 1.3 1996/11/20 10:00:06 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "feplot2.h"
#include "triang2.h"

#include "utils.h"
#include "numerics.h"
#include "triangtempl.h"

#include "cmdpars.h"
extern CmdPars Cmd;

static const Real EPS = 1e-6;

//-------------------------------------------------------------------------


FEPlotMESH2:: FEPlotMESH2(MESH* t, int plotType, char* caption, float size)
{
    backCol 	 = WHITE;
    dirichletCol = BLUE; 
    neumannCol 	 = GREEN; 
    cauchyCol 	 = CYAN;
    triFineCol	 = BLACK;
    triCoarseCol = MAGENTA;
    levelCol     = RED;
    percentage   = False; 

    mesh = t->castToMESH2();

    noOfLevels = 10;   Cmd.get("PlotLevels", &noOfLevels);
    --noOfLevels;

    const char *dummyCaption = "FePlot",  *cp;

    if (caption==0)  cp = dummyCaption;
    else		cp = caption;

    plotSize = size;
    plot = new Plot(plotType, size, cp);
    setMinMax();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: updateMesh(MESH* t) { mesh = t->castToMESH2(); }

//-------------------------------------------------------------------------

void FEPlotMESH2:: compLineLevels(Vector<Real>& x)
{
    Real xmax = RMIN, xmin = RMAX, diff, xVal;
    int k;
    PT2 *p;

    DListIter<PT2>  iterPT(mesh->points());

    while ((p =iterPT.all()))
    {
	if (x[p->node] > xmax) xmax = x[p->node];
	if (x[p->node] < xmin) xmin = x[p->node];
    }

    // if (percentage) xmin = 0.0;
    diff = (xmax-xmin) / (noOfLevels + 1);

    xVal = xmin + diff; // * ((percentage)? 1.0 : 0.5);

    for (k=1; k<=noOfLevels; k++)
    {
	levelsAt[k] = xVal;
	xVal += diff;
    }
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotContourLine(Real xh, EDG2 *ed1, EDG2 *ed2, Vector<Real>& x)
{
    double xa[2], ya[2], fac, quot;
    PT2 *p11 = ed1->p1, *p12 = ed1->p2, *p21 = ed2->p1, *p22 = ed2->p2;
    
    quot = x[p12->node] - x[p11->node];
    if (fabs(quot) < SMALL) quot = 0.5;
    fac = (xh-x[p11->node]) / quot;

    xa[0] = (p11->x)+fac*((p12->x)-(p11->x));
    ya[0] = (p11->y)+fac*((p12->y)-(p11->y));

    quot = x[p22->node] - x[p21->node];
    if (fabs(quot) < SMALL) quot = 0.5;
    fac = (xh-x[p21->node]) / quot;

    xa[1] = (p21->x)+fac*((p22->x)-(p21->x));
    ya[1] = (p21->y)+fac*((p22->y)-(p21->y));

    if      (equal(xh,0.0)) plot->set(PENCOL, GREEN);
    else if (xh > 0.0)      plot->set(PENCOL, RED);
    else	   	    plot->set(PENCOL, BLUE);

    plot->line(xa, ya, 2);
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotSolTriangle(TR2 *t, Vector<Real>& x)
{
    int  dec, i;
    Real zl, z1, z2, z3;
    
    z1 = x[t->p1->node]; z2 = x[t->p2->node]; z3 = x[t->p3->node];

    for (i=1; i<=noOfLevels; i++) 	
    {
        dec = 0;
        zl = levelsAt[i];
        if (z1 > zl) dec += 1;
        if (z2 > zl) dec += 2;
        if (z3 > zl) dec += 4;
	
        if (dec == 0)  break;
        else if ((dec == 1) || (dec == 6))  plotContourLine(zl,t->e3,t->e2,x);
        else if ((dec == 2) || (dec == 5))  plotContourLine(zl,t->e3,t->e1,x);
        else if ((dec == 3) || (dec == 4))  plotContourLine(zl,t->e1,t->e2,x);
    }
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotSolution(Vector<Real>& u)
{
    int  node;
    TR2* t;
    PT2* pt;

    levelsAt.resize(noOfLevels);
    compLineLevels(u);

    plot->set(PENCOL, levelCol);
    plot->set(PENSIZE, MEDIUM);

    DListIter<TR2>  iter(mesh->triangles());
    while ((t = iter.all())) plotSolTriangle(t,u);


    plot->set(FONTSIZE,SMALL);

    char line[100];
    Real uMin = machMax(Real(0)), uMax = -machMax(Real(0));

    DListIter<PT2> ptIter(mesh->points());
    while ((pt = ptIter.all()))
    {
	node = pt->node;
	if (u[node] < uMin) uMin = u[node];
	if (u[node] > uMax) uMax = u[node];
    }

    Real xdiff = (xMax-xMin)/10;
    Real ydiff = (yMax-yMin)/10.;

    sprintf(line,"Min: %1.4g   Max: %1.4g", uMin, uMax);
    plot->text(xMin+xdiff, yMax-0.5*ydiff, line);

    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotPointNodes()
{
    char line[10];

    plotElements(); 

    double xdiff = (xMax-xMin)/50.;
    double ydiff = (yMax-yMin)/50.;

    if (plotSize > 0.5) plot->set(FONTSIZE, MEDIUM);
    else     		plot->set(FONTSIZE, SMALL);

    plot->set(FONTCOL,  BLACK);

    PT2* pt;
    DListIter<PT2>  iter(mesh->points());

    while ((pt = iter.all())) 
    {
	sprintf (line, "%1d(%1d)", pt->node, pt->depth);
	plot->text ((pt->x)+xdiff, (pt->y)+ydiff, line);
    }

    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotEdgeNodes()
{
    char line[10];

    plotElements(); 

    double xdiff = (xMax-xMin)/50.;

    if (plotSize > 0.5) plot->set(FONTSIZE, MEDIUM);
    else     		plot->set(FONTSIZE, SMALL);

    plot->set(FONTCOL,  BLACK);

    EDG2 *ed;
    DListIter<EDG2> iter(mesh->edges());

    while ((ed = iter.all()))
    {
	// sprintf (line, "%1d(%1d)", edg->node, pt->depth);
	sprintf (line, "%1d", ed->node);
	plot->text (0.5*(ed->p1->x + ed->p2->x) + xdiff, 
		    0.5*(ed->p1->y + ed->p2->y), line);
    }
    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotTriangleNodes()
{
    TR2* t;
    char line[10];

    plotElements(); 

    if (plotSize > 0.5) plot->set(FONTSIZE, MEDIUM);
    else     		plot->set(FONTSIZE, SMALL);

    plot->set(FONTCOL, BLACK);

    DListIter<TR2>  iter(mesh->triangles());
    Vector<Real> xG(2);

    while ((t = iter.all()))
    {
	t->centerOfGravity(xG);
	// sprintf (line, "%1d(%1d)", ++no, t->depth);
	sprintf (line, "%1d", t->depth);
	plot->text (xG[1], xG[2], line);
    }
    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotElements()
{
    int depth;
    EDG2* ed;
    plot->set(PENSIZE, SMALL);
    plot->set(PENCOL,  triFineCol);

    DListIter<EDG2> iter(mesh->edges());
    while ((ed = iter.allOfHistory())) { if (!ed->refined()) plotEdge(ed); }

    plot->set(PENSIZE, MEDIUM);
    plot->set(PENCOL,  triCoarseCol);

    depth = 0;
    DListIter<EDG2> iter2(mesh->edges(), depth);
    while ((ed = iter2.allOfHistory())) { if (!ed->refined()) plotEdge(ed); }

    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotBoundary()
{
    // int depth=0;
    EDG2* ed;
    // DListIter<EDG2> iter(mesh->edges(), depth);
    DListIter<EDG2> iter(mesh->edges());

    plot->set(PENSIZE, BIG);

    plot->set(PENCOL, dirichletCol);
    //  while (ed = iter.allOfHistory()) plotEdge(ed, DIRICHLET);
    while ((ed = iter.all())) plotEdge(ed, DIRICHLET);

    plot->set(PENCOL, neumannCol);
    iter.reset();
    // while (ed = iter.allOfHistory()) plotEdge(ed, NEUMANN);

    while ((ed = iter.all())) plotEdge(ed, NEUMANN);

    plot->set(PENCOL, cauchyCol);
    iter.reset();
    // while (ed = iter.allOfHistory()) plotEdge(ed, CAUCHY);

    while ((ed = iter.all())) plotEdge(ed, CAUCHY);

    plot->flush();
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: plotEdge(EDG2* ed, int type)
{
    Real x[2], y[2];

    if (type > -999) { if (ed->boundP != type) return; }

    x[0] = ed->p1->x; x[1] = ed->p2->x;
    y[0] = ed->p1->y; y[1] = ed->p2->y;
    plot->line(x, y, 2);
}
/*-------------------------------------------------------------------------*/

/*
void FEPlotMESH2:: plotInitialEdges ()
{
    Real x[2], y[2];

    EDG2* ed;
    for (ed=mesh->FirstEdgOfInitial(); ed; ed=mesh->NextEdgOfInitial())
    {
	if (ed->son != nil) continue;
	x[0] = ed->p1->x; x[1] = ed->p2->x;
	y[0] = ed->p1->y; y[1] = ed->p2->y;
	plot->line (x, y, 2);
	
    }
}

void FEPlotMESH2:: plotActiveEdges (int type)
{
    Real x[2], y[2];

    EDG2* ed;
    for (ed=mesh->FirstEdg(); ed; d=mesh->NextEdg())
        if (ed->boundP == type)
	{
	    if (ed->son != nil) continue;
	    x[0] = ed->p1->x; x[1] = ed->p2->x;
	    y[0] = ed->p1->y; y[1] = ed->p2->y;
	    plot->line (x, y, 2);
	    
	}
}
*/
//-------------------------------------------------------------------------


void FEPlotMESH2:: setMinMax()
{
    int depth = 0;
    double dx, dy;
    PT2* p;
    xMin= 1e20; xMax= -1e20; yMin= 1e20, yMax= -1e20;

    DListIter<PT2>  iter(mesh->points(), depth);

    while ((p = iter.all())) 
    {
	if (p->x < xMin) xMin = p->x;
	if (p->x > xMax) xMax = p->x;
	if (p->y < yMin) yMin = p->y;
	if (p->y > yMax) yMax = p->y;
    }
 
    dx = (xMax-xMin)/10.0;
    dy = (yMax-yMin)/10.0;
    xMin -= dx; xMax += dx;
    yMin -= dy; yMax += dy;
    plot->set (MINX, xMin);
    plot->set (MINY, yMin);
    plot->set (MAXX, xMax);
    plot->set (MAXY, yMax);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

struct L_Tr_E
{
    TR2 *t;
    Real FarZ;
};

struct iIE
{
    EDG *e;
    Real x, y;
};

static L_Tr_E *TrL  = 0;

static Real bCubeX[8] = { -1.0, -1.0,  1.0,  1.0, -1.0, -1.0,  1.0,  1.0 };
static Real bCubeY[8] = { -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0 };
static Real bCubeZ[8] = { -1.0,  1.0,  1.0, -1.0, -1.0,  1.0,  1.0, -1.0 };
static Real tr_bCubeX[8], tr_bCubeY[8], tr_bCubeZ[8];
static int  bCubeSurf[6][5] = {
    {0, 1, 2, 3, 0},
    {0, 1, 5, 4, 0},
    {1, 2, 6, 5, 1},
    {2, 3, 7, 6, 2},
    {3, 0, 4, 7, 3},
    {4, 5, 6, 7, 4}
};
static Real bCubeSZ[8];

static int  len_TrL, No_L_Tr_E;

static Real Xoff, Yoff, Zoff, scalX, scalY, scalZ;
static Real V11, V12, V13, V21, V22, V23, V31, V32, V33, V34;
static Real SinAng, CosAng;
static Real Minjac, Maxjac, Totjac, Meanjac;
static int InnerCol;
static int DrawatF;

//-------------------------------------------------------------------------

static Real GrJacobi(Real p1x, Real p1y, Real p2x, Real p2y,
		      Real p3x, Real p3y)
{
    return (p2x - p1x) * (p3y - p1y) - (p3x - p1x) * (p2y - p1y);
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: VTr_bCube()
{
    int i, j, k;
    Real factor = REALPI / 180.;
    Real aux;
    double angle, dx, dy;
    
    for (i = 0; i < 8; i++)
    {
        tr_bCubeX[i] = V11*bCubeX[i] + V12*bCubeY[i] + V13*bCubeZ[i];
        tr_bCubeY[i] = V21*bCubeX[i] + V22*bCubeY[i] + V23*bCubeZ[i];
        tr_bCubeZ[i] = V31*bCubeX[i] + V32*bCubeY[i] + V33*bCubeZ[i] + V34;
    }
    dx = (double) (tr_bCubeX[4] - tr_bCubeX[0]);
    dy = (double) (tr_bCubeY[4] - tr_bCubeY[0]);
    angle = Abs(dy) > EPS ? (-atan2(dx,dy)) : (0.0);
    SinAng = (Real) sin(angle); CosAng = (Real) cos(angle);
    angle = angle / factor;
    for (i = 0; i < 8; i++)
    {
        aux = tr_bCubeX[i];
        tr_bCubeX[i] = aux *  CosAng + tr_bCubeY[i] *  SinAng;
        tr_bCubeY[i] = aux * -SinAng + tr_bCubeY[i] *  CosAng;
    }
    for (i = 0; i < 6; i++)
    {
        j = bCubeSurf[i][0]; k = bCubeSurf[i][2];
        bCubeSZ[i] = (tr_bCubeZ[j] + tr_bCubeZ[k]) / 2.0;
    }
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: ViewTransPoint(PT2 *p, Real z, Real& rTrx, Real& rTry, 
			      Real& rTrz)
{
    Real xn, yn, zn, auxx, auxy;
    
    xn =  (p->x - Xoff)/scalX - 2.0;
    zn = -(p->y - Yoff)/scalY - 2.0;
    yn =     (z - Zoff)/scalZ - 1.7321;

    auxx = V11 * xn + V12 * yn + V13 * zn;
    auxy = V21 * xn + V22 * yn + V23 * zn;
    rTrx = auxx *  CosAng + auxy *  SinAng; 
    rTry = auxx * -SinAng + auxy *  CosAng;
    rTrz = V31 * xn + V32 * yn + V33 * zn + V34;
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: Dr_bCube(int background)
{
    int i, j, k, l;
    Real x[2], y[2];
    
    if (background)
    {
        for (i = 0; i < 6; i++)
	{
            if (bCubeSZ[i] < 0.0)
	    {
                for (j = 0; j < 4; j++)
		{
                    k = bCubeSurf[i][j]; l = bCubeSurf[i][j+1];
                    if ((k == 1) || (l == 1)) plot->set(PENSIZE, BIG);
                    else		      plot->set(PENSIZE, SMALL);

                    x[0] = tr_bCubeX[k];  x[1] = tr_bCubeX[l];
                    y[0] = tr_bCubeY[k];  y[1] = tr_bCubeY[l];
                    plot->line(x, y, 2);
		}
	    }
	}
    }
    else
    {
        for (i = 0; i < 6; i++)
	{
            if (bCubeSZ[i] > 0.0)
	    {
                for (j = 0; j < 4; j++)
		{
                    k = bCubeSurf[i][j]; l = bCubeSurf[i][j+1];
                    if ((k == 1) || (l == 1)) plot->set(PENSIZE, BIG);
                    else		      plot->set(PENSIZE, SMALL);

                    x[0] = tr_bCubeX[k]; x[1] = tr_bCubeX[l];
                    y[0] = tr_bCubeY[k]; y[1] = tr_bCubeY[l];
                    plot->line(x, y, 2);
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: TesT(TR2 *t, Vector<Real>& rTrx, Vector<Real>& rTry)
{
    Real    xt[3], yt[3], jac;
    
    xt[0] = rTrx[t->p1->node];  yt[0] = rTry[t->p1->node];
    xt[1] = rTrx[t->p2->node];  yt[1] = rTry[t->p2->node];
    xt[2] = rTrx[t->p3->node];  yt[2] = rTry[t->p3->node];

    jac = GrJacobi(xt[0],yt[0],xt[1],yt[1],xt[2],yt[2]);
    jac = Abs(jac);
    Minjac = Min(Minjac,jac);
    Maxjac = Max(Maxjac,jac);
    Totjac += jac;
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: SortTri(TR2 *t, Real currFarZ)
{
    L_Tr_E *le, *len;
    Real midFarZ;
    int lwb, upb, mid, InsPos, Poff, i;
    int NotFound;
    
    if (No_L_Tr_E == 0)
    {
        le = &TrL[No_L_Tr_E++];
        (le->t) = t;
        (le->FarZ) = currFarZ;
    }
    else
    {
        lwb = 0; upb = No_L_Tr_E - 1;
        NotFound = True;
        while (NotFound)
	{
            mid = (upb + lwb) / 2;
            le = &TrL[mid];
            midFarZ = (le->FarZ);
            if (currFarZ > midFarZ)
	    { lwb = mid; Poff = 1; }
            else 
	    { upb = mid; Poff = 0; }
            if (upb == lwb)
	    {
                NotFound = False;
	    }
            else if (upb == lwb+1)
	    {
                if (Poff > 0)
		{
                    le = &TrL[upb];
                    if (currFarZ >= (le->FarZ))  Poff += 1;
		}
                else
		{
                    le = &TrL[lwb];
                    if (currFarZ <= (le->FarZ))  Poff -= 1;
		}
                NotFound = False;
	    }
	} /* while */
        if (No_L_Tr_E < len_TrL)
	{
            InsPos = mid + Poff;
            for (i = No_L_Tr_E; i > InsPos; i--)
	    {
                le = &TrL[i-1]; len = &TrL[i];
                (len->t) = (le->t);
                (len->FarZ) = (le->FarZ);
	    }
            le = &TrL[InsPos];
            (le->t) = t;
            (le->FarZ) = currFarZ;
            No_L_Tr_E += 1;
	}
    }
}
//-------------------------------------------------------------------------

void FEPlotMESH2:: QuartSort(TR2 *t,Real p1x,Real p1y,Real p1z,Real p2x,Real p2y,
			Real p2z,Real p3x,Real p3y,Real p3z,Real jac)
{
    Real p12x, p12y, p12z, p23x, p23y, p23z, p31x, p31y, p31z, njac;
    Real currFarZ;
    
    p12x = (p1x + p2x) / 2.0; p12y = (p1y + p2y) / 2.0; p12z = (p1z + p2z) / 2.0;
    p23x = (p2x + p3x) / 2.0; p23y = (p2y + p3y) / 2.0; p23z = (p2z + p3z) / 2.0;
    p31x = (p3x + p1x) / 2.0; p31y = (p3y + p1y) / 2.0; p31z = (p3z + p1z) / 2.0;
    njac = jac / 4.0;
    if (njac > Meanjac)
    {
        QuartSort(t,p1x,p1y,p1z,p12x,p12y,p12z,p31x,p31y,p31z,njac);
        QuartSort(t,p2x,p2y,p2z,p23x,p23y,p23z,p12x,p12y,p12z,njac);
        QuartSort(t,p3x,p3y,p3z,p31x,p31y,p31z,p23x,p23y,p23z,njac);
        QuartSort(t,p12x,p12y,p12z,p23x,p23y,p23z,p31x,p31y,p31z,njac);
    }
    else
    {
        currFarZ = (p1z + p12z + p31z) / 3.0;
        (t->mark) += 1;
        SortTri(t,currFarZ);
        currFarZ = (p2z + p23z + p12z) / 3.0;
        (t->mark) += 1;
        SortTri(t,currFarZ);
        currFarZ = (p3z + p31z + p23z) / 3.0;
        (t->mark) += 1;
        SortTri(t,currFarZ);
        currFarZ = (p12z + p23z + p31z) / 3.0;
        (t->mark) += 1;
        SortTri(t,currFarZ);
    }
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: DivEtSort(TR2 *t, Vector<Real>& rTrx, Vector<Real>& rTry,
			Vector<Real>& rTrz)
{
    Real xt[3], yt[3], jac, currFarZ;
    
    t->mark = 0;

    xt[0] = rTrx[t->p1->node];  yt[0] = rTry[t->p1->node];
    xt[1] = rTrx[t->p2->node];  yt[1] = rTry[t->p2->node];
    xt[2] = rTrx[t->p3->node];  yt[2] = rTry[t->p3->node];

    jac = GrJacobi(xt[0],yt[0],xt[1],yt[1],xt[2],yt[2]);
    jac = Abs(jac);
    if (jac > Meanjac)
    {
        QuartSort(t, 
		  xt[0], yt[0], rTrz[t->p1->node], 
		  xt[1], yt[1], rTrz[t->p2->node], 
		  xt[2], yt[2], rTrz[t->p3->node], jac);
    }
    else
    {
        currFarZ = (rTrz[t->p1->node]+rTrz[t->p2->node]+rTrz[t->p3->node]) / 3.0;
	(t->mark) += 1;
	SortTri(t,currFarZ);
    }
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: Dr3D(TR2 *t, Vector<Real>& rTrx, Vector<Real>& rTry)
{
    Real xt[4], yt[4], jac;
    int col, drawOK;
    
    drawOK = (DrawatF && ((t->mark) > 0)) || (!DrawatF && (--(t->mark) == 0));
    if (drawOK)
    {
        (t->mark) = 0; 

	xt[0] = rTrx[t->p1->node];  yt[0] = rTry[t->p1->node];
	xt[1] = rTrx[t->p2->node];  yt[1] = rTry[t->p2->node];
	xt[2] = rTrx[t->p3->node];  yt[2] = rTry[t->p3->node];
        xt[3] = xt[0], yt[3] = yt[0];

        jac = GrJacobi(xt[0],yt[0],xt[1],yt[1],xt[2],yt[2]);

        col = (jac<=0.0) ? InnerCol : WHITE;

        plot->set(FILLCOL, col);
        plot->fill(xt, yt, 4);
        plot->line(xt, yt, 4);
    }
}
//-------------------------------------------------------------------------


void FEPlotMESH2:: plot3D(Vector<Real>& x)
{
    PT2 *p;
    TR2* t;
    Real MaxZ, MinZ, sFact;
    Real SinTheta, CosTheta, SinPhi, CosPhi;
    Real factor = REALPI / 180.;
    double Theta, Phi;
    int i, node;
    
    Vector<Real> rTrx(x.size()); 
    Vector<Real> rTry(x.size()); 
    Vector<Real> rTrz(x.size()); 
    
    InnerCol = GREEN;

    Real rotX = 22, rotY = 33, rotZ = 0.0;
    Cmd.get("rotX", &rotX);
    Cmd.get("rotY", &rotY);
    Cmd.get("rotZ", &rotZ);


    Theta = rotX * factor;
    Phi   = rotY * factor;
    DrawatF = int(rotZ + 0.5);

    SinTheta = sin(Theta); CosTheta = cos(Theta);
    SinPhi   = sin(Phi);   CosPhi   = cos(Phi);
    V11 =  CosPhi * CosTheta; V12 =  SinPhi * CosTheta; V13 = -SinTheta;
    V21 = -SinPhi           ; V22 =  CosPhi;            V23 = 0;
    V31 =  CosPhi * SinTheta; V32 =  SinPhi * SinTheta; V33 =  CosTheta;
    V34 = 0.0; // (curGraph->rotZ);

    MaxZ = x[1];
    MinZ = x[1];
    
    FORALL(x,i) 
    {
	if (x[i] > MaxZ) MaxZ = x[i];
	if (x[i] < MinZ) MinZ = x[i];
    }

    Xoff = xMin;
    Yoff = yMin;
    Zoff = MinZ;

    scalX = (xMax - xMin) / 2.0;
    scalY = (yMax - yMin) / 2.0;
    scalZ = (MaxZ - MinZ) / 2.0;

    if (Cmd.get("scaleX",&sFact))  scalX /= sFact;
    if (Cmd.get("scaleY",&sFact))  scalY /= sFact;
    if (Cmd.get("scaleZ",&sFact))  scalZ /= sFact;

    //xMin = -1.7321; yMin = -1.7321;  // wozu hatte ich das geschrieben?
    //xMax =  1.7321; yMax =  1.7321;

    plot->set (MINX, -1.7321);
    plot->set (MAXX,  1.7321);
    plot->set (MINY, -1.7321);
    plot->set (MAXY,  1.7321);
    plot->set (PENCOL,  BLACK);
    plot->set (PENSIZE, SMALL);

    VTr_bCube();

    DListIter<PT2> iterPT(mesh->points());
    while((p = iterPT.all()))
    {
	node = p->node;
	ViewTransPoint(p, x[node], rTrx[node], rTry[node], rTrz[node]);
    }
    
    len_TrL = 8 * mesh->noOfTriangles;
    TrL = new L_Tr_E[len_TrL];
    if (TrL==0) { 
	cout << "\n*** plot3D: allocation failure for TrL\n";  cout.flush(); abort(); 
    }

    No_L_Tr_E = 0;
    Minjac = 9999999.; 
    Maxjac = 0.0; 
    Totjac = 0.0;


    DListIter<TR2> iterTR(mesh->triangles());
    while((t = iterTR.all())) t->mark = 0;

    iterTR.reset();
    while((t = iterTR.all()))  TesT(t, rTrx, rTry);

    Meanjac = (Totjac/mesh->noOfTriangles) / 2.0;

    iterTR.reset();
    while((t = iterTR.all()))  DivEtSort(t, rTrx, rTry, rTrz);

    for (i=0; i<No_L_Tr_E; i++)  Dr3D((TrL[i]).t, rTrx, rTry); 

    // Dr_bCube(True);
    // Dr_bCube(False);

    plot->flush();
    delete TrL;    
}