/*
$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;
}
syntax highlighted by Code2HTML, v. 0.9.1