/*
$Id: triang2.cc,v 1.9 1996/11/07 10:56:04 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "triang2.h"
#if NOKASKPLOT == 0
#include "feplot2.h"
#endif
#include "utils.h"
#include "numerics.h"
#include "mzibutil.h"
#include "alloc.h"
#include "triangtempl.h"
#include "cmdpars.h"
extern CmdPars Cmd;
Vector<PT*> EDG2::p(2);
Vector<PT*> TR2::p(3);
Vector<EDG*> TR2::e(3);
Vector<PATCH*> TR2::f(3);
//-------------------------------------------------------------------------
MESH2:: MESH2 (const char* inFileName, Bool readMesh)
: MESH(inFileName), ptList(0,25), edgList(0,25), trList(0,25),
ptAlloc(ElementsInBlock), trAlloc(2*ElementsInBlock),
edgAlloc(3*ElementsInBlock), oldGreenEdges(25), oldGreenTriangles(25)
{
noOfPoints = 0;
noOfEdges = 0;
noOfTriangles = 0;
// if (Cmd.isTrue("Circle")) circle = True;
// else circle = False;
if (readMesh)
{
if (Cmd.isTrue("ReadOldFormat")) readTriangulationOldFormat(inFileName);
else readTriangulation(inFileName);
}
}
//-------------------------------------------------------------------------
MESH2:: ~MESH2()
{
int i;
if (Cmd.isTrue("ReadOldFormat"))
{
delete [] trList[0]->first;
delete [] ptList[0]->first;
delete [] edgList[0]->first;
//RemoveGreen();
//FORALL( trList,i) trList[i]->deleteAll();
//FORALL( ptList,i) ptList[i]->deleteAll();
//FORALL(edgList,i) edgList[i]->deleteAll();
}
FORALL( trList,i) delete trList[i];
FORALL( ptList,i) delete ptList[i];
FORALL(edgList,i) delete edgList[i];
}
//-------------------------------------------------------------------------
FEPlot* MESH2:: newFEPlot(int plotType, char* caption, float size)
{
#if NOKASKPLOT == 0
return new FEPlotMESH2(this, plotType, caption, size);
#else
return 0;
#endif
}
//-------------------------------------------------------------------------
void MESH2:: readTriangulation(const char* fileName0)
{
int i, maxIndex;
TR2 *t;
DList<PT2>* ptL = new DList<PT2>; ptList.push(ptL);
DList<TR2>* trL = new DList<TR2>; trList.push(trL);
DList<EDG2>* edgL = new DList<EDG2>; edgList.push(edgL);
BufferedParser parser(fileName0, "end", CommentFlag);
if (parser.findWord("Poin") != parser.ok) missingItem("points",fileName);
if (parser.findWord("max") != parser.ok) missingItem("maxIndex",fileName);
parser.getValue(&maxIndex);
checkMaxIndex(maxIndex);
Vector<PT2*> ptAdr(0,maxIndex); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
readPoints(parser, ptAdr);
if (parser.findWord("Elem") != parser.ok) missingItem("triangles",fileName);
readElements(parser, ptAdr);
// -- compute edges:
Vector<Stack<EDG2*>*> edsOfPts(noOfPoints);
FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG2*>(10);
DListIter<TR2> iter(triangles());
while ((t=iter.all())) CompleteTriangle(t, edsOfPts);
if (parser.findWord("Boun") != parser.ok) missingItem("boundary",fileName);
readBoundary(parser, ptAdr, edsOfPts);
FORALL(edsOfPts,i) delete edsOfPts[i];
MakeGreen();
iter.reset();
while ((t=iter.all())) t->setBoundary();
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH2:: readPoints(BufferedParser& parser, Vector<PT2*>& ptAdr)
{
int ptNo, count=0;
PT2* p;
while (parser.getValue(&ptNo) == parser.ok)
{
if (ptNo > ptAdr.high()) pointIndexError(ptNo);
p = ptAlloc.Get(); // new PT2;
ptAdr[ptNo] = p;
parser.getValue(&p->x);
parser.getValue(&p->y);
p->node = ++count; // used in Find...
ptList[0]->add(p);
++noOfPoints;
}
}
//-------------------------------------------------------------------------
void MESH2:: readElements(BufferedParser& parser, Vector<PT2*>& ptAdr)
{
int pNo1, pNo2, pNo3, classA;
TR2 *t;
while (parser.getValue(&pNo1) == parser.ok)
{
parser.getValue(&pNo2);
parser.getValue(&pNo3);
t = getTR2(); // new TR2;
t->p1 = ptAdr[pNo1];
t->p2 = ptAdr[pNo2];
t->p3 = ptAdr[pNo3];
parser.getValue(&classA);
t->classA = classA; // int -> char
trList[0]->add(t);
++noOfTriangles;
}
}
//-------------------------------------------------------------------------
void MESH2:: readBoundary(BufferedParser& parser, Vector<PT2*>& ptAdr,
Vector<Stack<EDG2*>*>& edsOfPts)
{
int i, pNo1, pNo2, pNo;
PT2 *pt1, *pt2, *p;
EDG2 *e=0;
while (parser.getValue(&pNo1) == parser.ok) // edges on boundaries
{
parser.getValue(&pNo2);
pt1 = ptAdr[pNo1];
pt2 = ptAdr[pNo2];
Stack<EDG2*>& edges = *edsOfPts[pt1->node];
FORALL(edges,i)
{
e = edges[i];
if (pointOnEdge(pt2,e)) break;
}
if (!e->onBoundary()) boundaryWarning(pNo1,pNo2);
e->readBC(parser);
}
if (parser.findWord("Poin") == parser.ok) // additional BCs on points
{
while (parser.getValue(&pNo) == parser.ok)
{
p = ptAdr[pNo];
p->readBC(parser);
}
}
}
//-------------------------------------------------------------------------
void MESH2:: CompleteTriangle(TR2* t, Vector<Stack<EDG2*>*>& edsOfPts)
{
setRightHand(t);
PT2 *p1 = t->p1, *p2 = t->p2, *p3 = t->p3;
t->e1 = FindEdgByPts(p2,p3,edsOfPts);
t->e2 = FindEdgByPts(p3,p1,edsOfPts);
t->e3 = FindEdgByPts(p1,p2,edsOfPts);
SetTE(t);
}
//-------------------------------------------------------------------------
EDG2* MESH2:: FindEdgByPts(PT2* p1, PT2* p2, Vector<Stack<EDG2*>*>& edsOfPts)
{
int i;
EDG2* ed;
Stack<EDG2*>& edges = *edsOfPts[p1->node];
FORALL(edges,i)
{
ed = edges[i];
if (pointOnEdge(p2,ed)) return ed;
}
ed = edgAlloc.Get(); // new EDG2;
ed->p1 = p1;
ed->p2 = p2;
edsOfPts[p1->node]->push(ed);
edsOfPts[p2->node]->push(ed);
edgList[0]->add(ed);
++noOfEdges;
return ed;
}
//-------------------------------------------------------------------------
void MESH2:: setRightHand(TR2* t)
{
Jacobian J;
t->compJ(J,False);
if (J.detJ < 0.0) { PT2* p=t->p2; t->p2=t->p3; t->p3=p; }
}
//-------------------------------------------------------------------------
Bool MESH2:: pointOnEdge(PT2* p, EDG2* e)
{
return (e->p1 == p || e->p2 == p);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
#define GETLINE() line=tp;lineNo++;while(((*tp)!='\n')&&((*tp)!=NUL))tp++;if ((*tp)!=NUL){*tp=NUL;tp++;}
Bool MESH2:: readTriangulationOldFormat(const char* inFileName)
{
int ptdim, edgdim, trdim, bounddim;
int countP, countE, countT;
double x,y;
int rc, i, k1, k2, k3,
lineNo = 0, classA, partner;
char ch, *line, *tp;
char globBuf[256];
PT2 *p;
EDG2 *ed;
TR2 *t;
char* inFileText = 0;
ZIBReadFile(inFileName, &inFileText);
tp = inFileText;
GETLINE();
GETLINE();
rc = sscanf(line,"Dimension:(%d,%d,%d,%d)", &ptdim, &edgdim,
&trdim, &bounddim);
if (rc<=0)
{
cout << "\n*** ReadTriangulation: Error at top of file\n";
return False;
}
// allocate lists for depth 0, reset and push them onto the stacks:
DList<PT2>* ptL = new DList<PT2>;
ptList.push(ptL);
PT2* PArray = new PT2[ptdim];
for (i=0; i<ptdim; ++i)
{
PArray[i].reset();
ptList[0]->add(&PArray[i]);
}
DList<TR2>* trL = new DList<TR2>;
trList.push(trL);
TR2* TArray = new TR2[trdim];
for (i=0; i<trdim; ++i)
{
TArray[i].reset();
trList[0]->add(&TArray[i]);
}
DList<EDG2>* edgL = new DList<EDG2>;
edgList.push(edgL);
EDG2* EArray = new EDG2[edgdim];
for (i=0; i<edgdim; ++i)
{
EArray[i].reset();
edgList[0]->add(&EArray[i]);
}
if (PArray==0 || EArray==0 || TArray==0)
cout << "\n*** Read Triangulation: error in new for Arrays\n";
countP = 0; countE = 0; countT = 0;
while (True)
{
GETLINE();
for (i=0; i<256; ++i) globBuf[i] = '\0';
rc = sscanf(line, "%d:(%le,%le),%c%d", &i, &x, &y, &ch, &classA);
if (rc==0) break;
p = &PArray[i];
p->x = x;
p->y = y;
p->node = countP + 1;
p->boundP = ((ch=='D')||(ch=='d'))?DIRICHLET:
(((ch=='I')||(ch=='i'))?INTERIOR:
(((ch=='N')||(ch=='n'))?NEUMANN:CAUCHY));
if (rc==5) p->classA = classA;
++countP;
}
if (strncmp(line,"END",3) != 0)
{
printf("\n*** ReadTri: syntax error reading points (line %d)\n", lineNo);
return False;
}
noOfPoints = countP;
/* edges */
while (True)
{
GETLINE();
for (i=0; i<256; ++i) globBuf[i] = '\0';
rc = sscanf(line,"%d:(%d,%d),%c%s", &i,&k1,&k2,&ch,globBuf);
if (rc==0) break;
if ((k1>=ptdim)||(k2>=ptdim))
{
printf("\n *** ReadTri: bound check error reading edge (line %d)\n",
lineNo);
return False;
}
countE++;
ed = &EArray[i];
ed->p1 = &PArray[k1];
ed->p2 = &PArray[k2];
ed->boundP = ((ch=='D')||(ch=='d'))?DIRICHLET:
(((ch=='I')||(ch=='i'))?INTERIOR:
(((ch=='N')||(ch=='n'))?NEUMANN:CAUCHY));
if ((rc==5)&&((*globBuf)!=NUL))
{
rc = sscanf(globBuf,"%d,%c%d", &classA, &ch, &partner);
if (rc==0)
{
rc = sscanf(globBuf,",%c%d", &ch, &partner);
if (rc==0)
printf("\n Not a number after boundary type (line %d)\n",
lineNo);
}
else ed->classA = classA;
if (rc>1)
{
if ((ch=='S')||(ch=='s'))
ed->firstSon = &EArray[partner];
else if ((ch=='F')||(ch=='f'))
{
EArray[i+1].father = ed->father = &EArray[partner];
EArray[partner].pm = &PArray[k2];
}
else
{
printf("\n *** Unknown Character '%c' (line %d)\n",
ch, lineNo);
}
}
}
}
if (strncmp(line,"END",3)!=0)
{
printf("\n*** ReadTri: syntax error reading edges (line %d)\n", lineNo);
return False;
}
noOfEdges=countE;
/* triangles */
while (True)
{
GETLINE();
for (i=0; i<256; ++i) globBuf[i] = '\0';
rc = sscanf(line,"%d:(%d,%d,%d)%s",&i,&k1,&k2,&k3,globBuf);
if (rc==0) break;
if ((k1>=edgdim)||(k2>=edgdim)||(k3>=edgdim))
{
printf(
"\n *** ReadTri: bound check error reading triangle (line %d)\n",
lineNo);
cout.flush(); abort();
}
countT++;
t = &TArray[i];
t->e1 = &EArray[k1];
t->e2 = &EArray[k2];
t->e3 = &EArray[k3];
// t->next = &TArray[i+1];
// t->prev = &TArray[i-1];
// CallProcs(NewTriangle,(ptr)t,sizeof(TR2),0);
SetTP(t, 0);
SetTE(t);
if ((rc==5)&&((*globBuf)!=NUL))
{
rc = sscanf(globBuf,",%d,%d", &classA, &partner);
if (rc==0)
{
printf("\n *** Not a number after triangle type (line %d)\n",
lineNo);
}
else t->classA = classA;
/*
if (rc>1)
{
if ((partner>=0)&&(partner<trdim))
t->partner = &TArray[partner];
else
{
sprintf(globBuf, "ReadTri: bound check error reading
triangle (line %d)\n", lineNo);
ZIBStdOut(globBuf);
CleanUp(inFile);
return False;
}
}
*/
}
}
if (strncmp(line,"END",3)!=0)
{
printf("\n *** ReadTri: syntax error reading triangles (line %d)\n",
lineNo);
return False;
}
noOfTriangles = countT;
/* boundary */
/*
if (bounddim!=0)
{
for (k = 0;k<bounddim;k++)
{
GETLINE();
rc = sscanf(line,"(%le,%le):%d\n",&x,&y,&k1);
midPoint = NewArcVec();
if (midPoint==0) break;
midPoint[0] = x;
midPoint[1] = y;
if (rc==0) break;
for (i=0;i<k1;i++)
{
GETLINE();
rc = sscanf(line,"%d",&k2);
AE(&EArray[k2],arcData) = (ptr)midPoint;
}
}
GETLINE();
if (strncmp(line,"END",3)!=0)
{
printf("\n *** ReadTri: syntax error reading bounds (line %d)\n",
lineNo);
return False;
}
}
*/
MakeGreen();
DListIter<TR2> iter(triangles());
while ((t=iter.all())) t->setBoundary();
if (infoRefinement) Info();
delete inFileText;
return True;
}
//-------------------------------------------------------------------------
void MESH2:: writeTriangulation(const char* /*outFileText*/, const Vector<Real>& u)
const
{
if (Cmd.isTrue("writeGrapeFormat")) writeGrapeFormat (fileName, u);
if (Cmd.isTrue("writeZIBFormat")) writeZIBFormat (fileName, u);
}
//-------------------------------------------------------------------------
void MESH2:: writeZIBFormat(const char* outFileText, const Vector<Real>& u)
const
{
int i;
int node;
TR2 *t;
EDG2 *ed;
PT2 *p;
FILE *outFile;
Vector<int> PNodes(noOfPoints);
// saving the node numbers for later reuse
i = 0;
DListIter<PT2> PIter(points());
while ((p=PIter.all())) PNodes[++i] = p->node;
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".zib.geo");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** writeTriangulation: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
fprintf(outFile, "# %s\n\n", fileName);
fprintf(outFile, "Points: maxIndex %d\n", noOfPoints);
node = 0;
DListIter<PT2> iterPts(points());
while ((p=iterPts.all())) {
p->node = node+1;
fprintf(outFile, " %d %e %e\n", p->node, p->x, p->y);
node++;
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nElements:\n");
DListIter<TR2> iterTrs(triangles());
while ((t=iterTrs.all())) {
fprintf(outFile, "%d %d %d %d\n",
(t->p1)->node, (t->p2)->node,
(t->p3)->node, t->classA );
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nBoundary: \n");
DListIter<EDG2> iterEds(edges());
while ((ed=iterEds.all())) {
if (!ed->onBoundary()) continue;
if (ed->isDirichlet())
fprintf(outFile, "%d %d d %d\n",
(ed->p1)->node, (ed->p2)->node,ed->classA );
else if (ed->isNeumann())
fprintf(outFile, "%d %d n %d\n",
(ed->p1)->node, (ed->p2)->node,ed->classA );
else if (ed->isCauchy())
fprintf(outFile, "%d %d c %d\n",
(ed->p1)->node, (ed->p2)->node, ed->classA );
}
fprintf(outFile, "END\n");
// output of values at the nodes located in the points
fprintf(outFile, "\nValues in nodes: \n");
i = 0;
iterPts.reset();
while ((p=iterPts.all())) fprintf(outFile, "%e\n", u[PNodes[++i]]);
fprintf(outFile, "END\n");
cout << "\n WriteZIB-Data ok.; written to " << outText << "\n";
delete outText;
fclose(outFile);
// restore of node numbers
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
}
//-------------------------------------------------------------------------
void MESH2:: writeGrapeFormat(const char* outFileText, const Vector<Real>& u)
const
{
FILE *outFile;
int i;
int node;
int k[3];
TR2 *t;
PT2 *p;
Vector<int> PNodes(noOfPoints), TNodes(noOfTriangles);
i = 0;
DListIter<PT2> PIter(points());
while ((p=PIter.all())) PNodes[++i] = p->node;
i = 0;
DListIter<TR2> TIter(triangles());
while ((t=TIter.all())) TNodes[++i] = t->node;
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".grape");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** grapeTriangulation: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
fprintf(outFile, "%d %d\n", noOfPoints, noOfTriangles);
node = 0;
DListIter<PT2> iterPts(points());
while ((p=iterPts.all())) {
p->node = node;
fprintf(outFile, "%e %e\n", p->x, p->y);
node++;
}
node = 0;
DListIter<TR2> iterT(triangles());
while ((t=iterT.all())) {
t->node = node;
fprintf(outFile, "%d %d %d \n",
(t->p1)->node, (t->p2)->node,
(t->p3)->node);
node++;
}
iterT.reset();
while ((t=iterT.all()))
{
EDG2* ed;
ed = t->e1;
if ( (ed->t1 == nil) || (ed->t2 == nil) )
k[0] = -1;
else
if ( (ed->t1 == t) )
k[0] = ed->t2->node;
else k[0] = ed->t1->node;
ed = t->e2;
if ( (ed->t1 == nil) || (ed->t2 == nil) )
k[1] = -1;
else
if ( (ed->t1 == t) ) k[1] = ed->t2->node;
else k[1] = ed->t1->node;
ed = t->e3;
if ( (ed->t1 == nil) || (ed->t2 == nil) )
k[2] = -1;
else
if ( (ed->t1 == t) ) k[2] = ed->t2->node;
else k[2] = ed->t1->node;
fprintf(outFile, "%d %d %d\n",
k[0], k[1], k[2]);
}
// fprintf(outFile, "1\n");
i = 0;
iterPts.reset();
while ((p=iterPts.all())) fprintf(outFile, "%e\n", u[PNodes[++i]]);
cout << "\n GRAPE-Data ok.; written to " << outText << "\n";
delete outText;
fclose(outFile);
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
i = 0;
TIter.reset();
while ((t=TIter.all())) t->node = TNodes[++i];
return;
}
//-------------------------------------------------------------------------
void MESH2:: Refine()
{
TR2* t;
int k, noPt[3], noEdg[3], noTr[3];
Timer timer, accTimer;
noPt[0] = noOfPoints;
noTr[0] = noOfTriangles;
noEdg[0] = noOfEdges;
RemoveGreen(); // delete all green triangless
noPt[1] = noOfPoints;
noTr[1] = noOfTriangles;
noEdg[1] = noOfEdges;
DList<TR2>* trL = new DList<TR2>;
trList.push(trL);
FORALL(oldGreenTriangles,k) RefineTriangle(oldGreenTriangles[k]);
DListIter<TR2> iter(triangles());
while ((t=iter.all())) RefineTriangle(t);
noPt[2] = noOfPoints;
noTr[2] = noOfTriangles;
noEdg[2] = noOfEdges;
MakeGreen();
++refStep;
if (infoRefinement) refinementInfo (noPt, noEdg, noTr);
timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------
void MESH2:: Flatten()
{
int noReleasedEdges=0, noReleasedTriangles=0;
int noPoints=0, noEdges=0, noTriangles=0;
int k, greenTriangles, greenEdges;
const int impossibleDepth=255;
TR2* t;
EDG2* ed;
PT2* p;
DList<TR2>* trL = new DList<TR2>;
DList<EDG2>* edgL = new DList<EDG2>;
DList<PT2>* ptL = new DList<PT2>;
Stack<TR2*> trToRetain(25), trToRelease(25);
Stack<EDG2*> edgToRetain(25), edgToRelease(25);
Stack<PT2*> ptToRetain(25);
DListIter<TR2> iterTr(triangles());
DListIter<EDG2> iterEdg(edges());
DListIter<PT2> iterPt(points());
greenEdges = noOfEdges;
greenTriangles = noOfTriangles;
RemoveGreen(); // delete all green triangless
greenEdges -= noOfEdges;
greenTriangles -= noOfTriangles;
while ((ed=iterEdg.allOfHistory())) ed->mark=False;
iterEdg.reset();
for (k=oldGreenEdges.l; k<=oldGreenEdges.h; k++)
{
ed = oldGreenEdges[k];
ed->mark=False;
}
while (!oldGreenTriangles.empty())
{
t = oldGreenTriangles.pop();
if (t->refined())
{
trToRelease.push(t);
noReleasedTriangles++;
}
else
{
trToRetain.push(t);
noTriangles++;
(t->e1)->mark = True;
(t->e2)->mark = True;
(t->e3)->mark = True;
}
}
while ((t=iterTr.allOfHistory()))
{
if (t->refined())
{
trToRelease.push(t);
noReleasedTriangles++;
}
else
{
trToRetain.push(t);
noTriangles++;
(t->e1)->mark = True;
(t->e2)->mark = True;
(t->e3)->mark = True;
}
}
for (k=oldGreenEdges.l; k<=oldGreenEdges.h; k++)
{
ed = oldGreenEdges[k];
if (ed->mark)
{
edgToRetain.push(ed);
noEdges++;
}
else
{
edgToRelease.push(ed);
noReleasedEdges++;
}
}
while ((ed=iterEdg.allOfHistory()))
{
if (ed->mark)
{
edgToRetain.push(ed);
noEdges++;
}
else
{
edgToRelease.push(ed);
noReleasedEdges++;
}
}
while (!oldGreenEdges.empty()) oldGreenEdges.pop();
while ((p=iterPt.allOfHistory()))
{
ptToRetain.push(p);
noPoints++;
}
while (!trList.empty()) { delete trList.pop(); }
trList.push(trL);
while (!edgList.empty()) { delete edgList.pop(); }
edgList.push(edgL);
while (!ptList.empty()) { delete ptList.pop(); }
ptList.push(ptL);
for (k=ptToRetain.l; k<=ptToRetain.h; k++)
{
p = ptToRetain[k];
ptL->add(p);
p->depth =0;
}
for (k=edgToRetain.l; k<=edgToRetain.h; k++)
{
ed = edgToRetain[k];
if (ed->firstSon)
{
ed->depth = impossibleDepth;
oldGreenEdges.push(ed);
noEdges--;
ed->father = nil;
continue;
}
if ((ed->father) && (!((ed->father)->mark))) ed->father = nil;
ed->depth = 0;
edgL->add(ed);
}
for (k=trToRetain.l; k<=trToRetain.h; k++)
{
t = trToRetain[k];
t->father = nil;
if (((t->e1)->firstSon)||((t->e2)->firstSon)||((t->e3)->firstSon))
{
t->depth = impossibleDepth;
t->mark = False;
oldGreenTriangles.push(t);
}
else
{
t->depth = 0;
trL->add(t);
}
}
while (!trToRelease.empty())
{
t = trToRelease.pop();
DelTE(t,3);
returnTR2(t);
}
while (!edgToRelease.empty())
{
ed = edgToRelease.pop();
edgAlloc.Return(ed);
}
refStep = 0;
maxDepth = 0;
ptIterLevel = 0;
iterLevel = 0;
noOfPoints = noPoints;
noOfEdges = noEdges;
noOfTriangles = noTriangles;
greenEdges = noOfEdges;
greenTriangles = noOfTriangles;
MakeGreen();
greenEdges -= noOfEdges;
greenTriangles -= noOfTriangles;
if (infoRefinement)
{
cout << "\n Flatten: " << noReleasedEdges << " edges, " <<
noReleasedTriangles << " triangles released\n";
cout << " " << noPoints << " points, " <<
noEdges+greenEdges << " edges, " <<
noTriangles+greenTriangles << " triangles retained\n";
}
return;
}
//-------------------------------------------------------------------------
#define NEIGH_TR2(ed,t) (((ed->t1)==t)?(ed->t2):(ed->t1))
/* ***************************************************************
RefineTriangle (DoRefTr) internal routine to refine a triangle, loTrues for
neighbors do be refined too.
--------
TR2 *t Adress of Triangle
return int True or False
*************************************************************** */
int MESH2:: RefineTriangle(TR2 *t)
{
TR2 *neighbor, *bigNeighbor;
EDG2 *ed, *edFather;
EDG2* eds[1+3];
int k;
if (t->mark) t->mark = False; // return if not marked
else return True;
if (t->firstSon != 0) return True; // return if already refined
if (RefineTriangleRed(t) != True) // refine it
{
cout << "\n*** RefineTriangle: triangle not refined \n"; cout.flush(); abort();
}
/*
Algorithm: For all three edges:
(1) look if neighbor triangle has 2 refined neighbors
and refine it this case,
(2) if no neighbor triangle exists, "father neighbor"
triangle needs to be refined (green refined at
the last level),
This routine uses recursion.
*/
t->getEdges(eds);
for (k=1; k<=3; k++)
{
ed=eds[k];
if ((ed->boundP)==INTERIOR)
{
neighbor = NEIGH_TR2(ed,t);
if (neighbor!=0)
{
if (Mark2(neighbor)==True) // (1)
{
neighbor->mark=True;
if (!RefineTriangle(neighbor)) return False;
}
}
else // (2)
{
edFather = ed->father;
if (edFather==0) continue;
if ((edFather->t1)==0) continue;
bigNeighbor = (edFather->t1->firstSon == 0) ?
edFather->t1 : edFather->t2;
if (bigNeighbor == 0) continue; // return True; ???
bigNeighbor->mark=True;
if (!RefineTriangle(bigNeighbor))
{ printf("\nHuch\n"); return False; }
}
}
}
return True;
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
Mark2(t) tests if triangle has two refined neighbors
----------
TR2 *t Adress of the triangle
return int True or False
****************************************************************** */
int MESH2:: Mark2(TR2 *t)
{
int k, anz=0;
EDG2* eds[1+3];
t->getEdges(eds);
if (t->mark) return True;
for (k=1; k<=3; k++) if ((eds[k]->firstSon)!=0) anz++;
return (anz > 1) ? True : False;
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
RefineTriangleRed (RefineTr) does the actual refinement
----------
TR2 *t Adress of the triangle
return int True or False
****************************************************************** */
int MESH2:: RefineTriangleRed(TR2 *t)
{
TR2 *t1, *t2, *t3, *t4;
EDG2 *e1 = t->e1, *e2 = t->e2, *e3 = t->e3, *ereg,
*ie1,*ie2,*ie3, *e11, *e12, *e21, *e22, *e31, *e32;
// Get 4 new triangles, 3 new edges:
t1 = getTR2();
t2 = getTR2();
t3 = getTR2();
t4 = getTR2();
noOfTriangles += 3;
ie1 = edgAlloc.Get();
ie2 = edgAlloc.Get();
ie3 = edgAlloc.Get();
noOfEdges += 3;
// Set the "son" edges eik, i=1,2,3; k=1,2;
// if neeccessary refine the edges
if ((e1->firstSon)==0) /* Edge not refined */
{
e11 = RefineEdge(e1);
if (e11==0) return False; /* No space for Refinement */
}
else e11 = e1->firstSon; /* Take the existing "son" edge */
e12 = e11->next;
if ((e11->p1)==(t->p3))
{ ereg = e11; e11 = e12; e12 = ereg; }
if ((e2->firstSon)==0)
{
e21 = RefineEdge(e2);
if (e21==0) return False;
}
else e21 = e2->firstSon;
e22 = e21->next;
if ((e21->p1)==(t->p1))
{ ereg = e21; e21 = e22; e22 = ereg; }
if ((e3->firstSon)==0)
{
e31 = RefineEdge(e3);
if (e31==0) return False;
}
else e31 = e3->firstSon;
e32 = e31->next;
if ((e31->p1)==(t->p2))
{ ereg = e31; e31 = e32; e32 = ereg; }
/* Build the new triangles */
InnerEdges(t,e1,e2,e3,ie1,ie2,ie3);
t1->e1 = e32; t1->e2 = e11; t1->e3 = ie1;
t2->e1 = e12; t2->e2 = e21; t2->e3 = ie2;
t3->e1 = e22; t3->e2 = e31; t3->e3 = ie3;
t4->e3 = ie1; t4->e2 = ie3; t4->e1 = ie2;
SetTP(t1,t); SetTE(t1);
SetTP(t2,t); SetTE(t2);
SetTP(t3,t); SetTE(t3);
SetTP(t4,t); SetTE(t4);
t->firstSon = t1;
t->type = T_RED;
ie1->depth = ie2->depth = ie3->depth = t1->depth;
trList[trList.h]->add(t1);
trList[trList.h]->add(t2);
trList[trList.h]->add(t3);
trList[trList.h]->add(t4);
if (ie1->depth > edgList.h)
{
DList<EDG2>* edgL = new DList<EDG2>;
edgList.push(edgL);
}
edgList[ie1->depth]->add(ie1);
edgList[ie1->depth]->add(ie2);
edgList[ie1->depth]->add(ie3);
return True;
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
RefineEdge(ed) does the actual refinement of an edge
----------
EDG2 *ed Adress of the edge
return int True or False
****************************************************************** */
EDG2* MESH2:: RefineEdge(EDG2 *ed)
{
EDG2 *e1, *e2;
PT2 *pm;
e1 = edgAlloc.Get();
e2 = edgAlloc.Get();
++noOfEdges;
pm = ptAlloc.Get();
++noOfPoints;
pm->classA = ed->classA;
pm->boundP = ed->boundP;
pm->depth = ed->depth + 1;
if ( (ed->isoP == CIRCULAR) && (ed->boundP != INTERIOR) )
setXYpm(pm, ed);
else
{
pm->x = 0.5*(ed->p1->x + ed->p2->x);
pm->y = 0.5*(ed->p1->y + ed->p2->y);
}
e1->p1 = ed->p1; e1->p2 = pm;
e1->boundP = ed->boundP; e1->father = ed; e1->type = ed->type;
e2->p1 = pm; e2->p2 = ed->p2;
e2->boundP = ed->boundP; e2->father = ed; e2->type = ed->type;
e1->classA = e2->classA = ed->classA;
e1->isoP = e2->isoP = ed->isoP;
e1->isBlueDiagonal = e2->isBlueDiagonal = ed->isBlueDiagonal;
ed->firstSon = e1; ed->pm =pm;
e1->depth = e2->depth = ed->depth+1;
if (e1->depth > edgList.h) // create new list at depth if necessary
{
DList<EDG2>* edgL = new DList<EDG2>;
edgList.push(edgL);
}
edgList[e1->depth]->add(e1);
edgList[e1->depth]->add(e2);
if (pm->depth > ptList.h)
{
DList<PT2>* ptL = new DList<PT2>;
ptList.push(ptL);
}
ptList[pm->depth]->add(pm);
return e1;
}
/*-------------------------------------------------------------------------*/
void MESH2:: setXYpm(PT2* pm, const EDG2* e)
{
Real x1, x2, y1, y2, phi1, phi2, phi, radius;
Real xx, yy;
int quadrant;
Real realPi=3.1415926535897932384, realPi2=1.57079632679489661923;
// dieser Teil ist von einem frueheren Programm uebernommen.
// fuer Kreisboegen geht das sicher einfacher.
x1 = e->p1->x; y1 = e->p1->y;
x2 = e->p2->x; y2 = e->p2->y;
radius = sqrt(x1*x1 + y1*y1);
if (x1>=0.0)
quadrant = (y1>0)?0:3;
else
quadrant = (y1>=0)?1:2;
xx = fabs((quadrant==0)||(quadrant==2)?x1:y1);
yy = fabs((quadrant==0)||(quadrant==2)?y1:x1);
phi1 = atan2(yy,xx)+quadrant*realPi2;
if (x2>=0.0)
quadrant = (y2>0)?0:3;
else
quadrant = (y2>=0)?1:2;
xx = fabs((quadrant==0)||(quadrant==2)?x2:y2);
yy = fabs((quadrant==0)||(quadrant==2)?y2:x2);
phi2 = atan2(yy,xx)+quadrant*realPi2;
phi = 0.5*(phi1 + phi2);
if (fabs(phi1 - phi2) > realPi) phi = phi - realPi;
pm->x = radius * cos(phi);
pm->y = radius * sin(phi);
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
InnerEdges(..) sets the points of the inner edges
----------
EDG2 *e1,*e2,*e3 Edges of the triangle to be refined
EDG2 *ie1,*ie2,*ie3 Inner edges
****************************************************************** */
void MESH2:: InnerEdges (TR2 *t,EDG2 *e1,EDG2 *e2,EDG2 *e3,EDG2 *ie1,
EDG2 *ie2,EDG2 *ie3)
{
ie1->p1 = e1->pm; ie1->p2 = e3->pm;
ie1->type = T_RED; ie1->classA = t->classA;
ie2->p1 = e2->pm; ie2->p2 = e1->pm;
ie2->type = T_RED; ie2->classA = t->classA;
ie3->p1 = e3->pm; ie3->p2 = e2->pm;
ie3->type = T_RED; ie3->classA = t->classA;
ie1->isBlueDiagonal = e2->isBlueDiagonal;
ie2->isBlueDiagonal = e3->isBlueDiagonal;
ie3->isBlueDiagonal = e1->isBlueDiagonal;
return;
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
MakeTriangleGreen(t) generates two green triangles
----------
TR2 *t Adress of the triangle
return int True or False
****************************************************************** */
int MESH2:: MakeGreen()
{
TR2* t;
int k;
DListIter<TR2> iter(triangles());
while ((t=iter.all())) MakeTriangleGreen(t);
FORALL(oldGreenTriangles,k)
if (!(oldGreenTriangles[k]->firstSon))
MakeTriangleGreen(oldGreenTriangles[k]);
return True;
}
//-------------------------------------------------------------------------
int MESH2:: MakeTriangleGreen(TR2 *t)
{
PT2 *p1, *p2;
EDG2 *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew;
TR2 *t1New,*t2New;
// Testing if one neighbor triangle is refined and setting pi, ei
// corresponding to this neighbor:
if ((e1->firstSon)!=0)
{
p1 = t->p1; p2 = t->p2;
e1 = t->e1; e2 = t->e2; e3 = t->e3;
}
else
{
if ((e2->firstSon)!=0)
{
p1 = t->p2; p2 = t->p3;
e1 = t->e2; e2 = t->e3; e3 = t->e1;
}
else
{
if ((e3->firstSon)!=0)
{
p1 = t->p3; p2 = t->p1;
e1 = t->e3; e2 = t->e1; e3 = t->e2;
}
else return True; // nothing to refine
}
}
// e1->pm->mark = T_GREEN; // "green point"
// Get one new edge and two new triangles
eNew = edgAlloc.Get();
noOfEdges++;
t1New = getTR2();
t2New = getTR2();
noOfTriangles += 1;
t1New->next = t2New;
t2New->prev = t1New;
// Set the fields for the new objects
eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;
eNew->type = T_GREEN;
eNew->classA = t->classA;
eNew->depth = t->depth+1;
t1New->e1 = e3; t1New->e3 = eNew; t1New->type = T_GREEN; // !!! T_WHITE;
t2New->e1 = e2; t2New->e2 = eNew; t2New->type = T_GREEN; // !!! T_WHITE;
if ((e1->p1)==p2) // test for orientation
{
t1New->e2 = e1->firstSon; t2New->e3 = (e1->firstSon)->next;
}
else
{
t1New->e2 = (e1->firstSon)->next; t2New->e3 = e1->firstSon;
}
DelTE(t,3);
SetTP(t1New,t); SetTP(t2New,t);
SetTE(t1New); SetTE(t2New);
t->firstSon = t1New;
// if (t->type == T_GREEN) cout << "*** double green triangle! \n"; // !!!
t->type = T_GREEN;
// Append two new ones to their lists
trList[trList.h]->add(t1New);
trList[trList.h]->add(t2New);
if (eNew->depth > edgList.h)
{
DList<EDG2>* edgL = new DList<EDG2>;
edgList.push(edgL);
}
edgList[eNew->depth]->add(eNew);
return True;
}
//-------------------------------------------------------------------------
/* *******************************************************************
RemoveGreen() remove all green triangles
----------
****************************************************************** */
void MESH2:: RemoveGreen()
{
TR2 *t1old,*t2old,*tFather;
EDG2 *eMid = 0;
int k;
// Loop for all triangles of last refinement step (green triangles only here)
TR2* t = trList[refStep]->first;
while (t != 0)
{
k = 0;
tFather = t->father;
if (tFather == 0) { t = t->next; continue; }
if (tFather->type != T_GREEN) { t = t->next; continue; }
// Remove two green triangles, print internal error message
// if the green edge is not found
t1old = t; t2old = t->next;
if (t->e1->type==T_GREEN) { eMid = t->e1; k++;}
if (t->e2->type==T_GREEN) { eMid = t->e2; k++;}
if (t->e3->type==T_GREEN) { eMid = t->e3; k++;}
if (k != 1)
cout << "\n*** RemoveGreen: Triangle with impossible state k="
<< k << "\n";
if (t1old->mark) tFather->mark = True; // refine father
else if (t2old->mark) tFather->mark = True; // refine father
else tFather->mark = False;
t = t2old->next;
// -- Remove the two green triangles; return the space
trList[refStep]->remove(t1old);
trList[refStep]->remove(t2old);
DelTE(t1old,3);
DelTE(t2old,3);
returnTR2(t1old);
returnTR2(t2old);
tFather->firstSon = 0;
tFather->type = T_WHITE;
SetTE(tFather);
edgList[eMid->depth]->remove(eMid);
edgAlloc.Return(eMid);
--noOfTriangles;
--noOfEdges;
}
return;
}
/*-------------------------------------------------------------------------*/
/* *******************************************************************
SetTP(t,f) set the points of a triangle, computed from the
edges
----------
TR2 *t Adress of the triangle
TR2 *f Father of t
****************************************************************** */
void MESH2:: SetTP(TR2 *t, TR2 *f)
{
PT2 *P11 = (t->e1)->p1, *P12 = (t->e1)->p2,
*P21 = (t->e2)->p1, *P22 = (t->e2)->p2,
*P31 = (t->e3)->p1, *P32 = (t->e3)->p2;
t->p1 = ((P31==P11)||(P31==P12)) ? P32:P31;
t->p2 = ((P11==P21)||(P11==P22)) ? P12:P11;
t->p3 = ((P21==P31)||(P21==P32)) ? P22:P21;
t->father = f;
if (f != 0)
{
t->depth = (f->depth)+1;
t->classA = f->classA;
}
if (t->depth > maxDepth && t->type != T_GREEN) maxDepth = t->depth;
return;
}
//-------------------------------------------------------------------------
/* *******************************************************************
SetTE(t) set the triangles of an edge, computed from a new
triangle, makes checks for programming errors
----------
TR2 *t Adress of the triangle
****************************************************************** */
void MESH2:: SetTE(TR2 *t)
{
EDG2 *ed;
int k = 0;
ed = t->e1;
if ((ed->t1)==0) { ed->t1 = t; k++; }
else if ((ed->t2)==0) { ed->t2 = t; k++; }
ed = t->e2;
if ((ed->t1)==0) { ed->t1 = t; k++; }
else if ((ed->t2)==0) { ed->t2 = t; k++; }
ed = t->e3;
if ((ed->t1)==0) { ed->t1 = t; k++; }
else if ((ed->t2)==0) { ed->t2 = t; k++; }
if (k!=3)
{
cout << "\n *** Tri: impossible state (SetTE), k = " << k;
cout << " (type: " << (int) t->type << ")\n";
printf("\n adress t: %p (%d,%d,%d)",t,
t->p1->node,t->p2->node,t->p3->node);
printf("\n adress e1: %p (%d-%d)",t->e1,
t->e1->p1->node,t->e1->p2->node);
printf("\n adress e2: %p (%d-%d)",t->e2,
t->e2->p1->node,t->e2->p2->node);
printf("\n adress e3: %p (%d-%d)",t->e3,
t->e3->p1->node,t->e3->p2->node);
printf("\n t->e1->t1: %p",(t->e1)->t1);
printf("\n t->e1->t2: %p",(t->e1)->t2);
printf("\n t->e2->t1: %p",(t->e2)->t1);
printf("\n t->e2->t2: %p",(t->e2)->t2);
printf("\n t->e3->t1: %p",(t->e3)->t1);
printf("\n t->e3->t2: %p",(t->e3)->t2);
#if NOKASKPLOT == 0
FEPlotMESH2* fePlot = new FEPlotMESH2(this, SCREEN);
fePlot->plotElements();
fePlot->plotTriangleNodes();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e1);
fePlot->plotEdge(t->e2);
fePlot->plotEdge(t->e3);
fePlot->plot->flush();
Pause();
#endif
}
if ( t->e1->onBoundary() || t->e2->onBoundary() || t->e3->onBoundary() )
t->boundP = BOUNDARY;
return;
}
/* *******************************************************************
DelTE(t, check) delete a triangle on its edges
----------
TR2 *t Adress of the triangle
****************************************************************** */
void MESH2:: DelTE(TR2 *t, int check)
{
EDG2 *ed;
int k = 0;
/*
if triangle on edge, remove it, if not generate internal
error message
*/
ed = t->e1;
if ((ed->t1)==t) { ed->t1 = 0; k++; }
else if ((ed->t2)==t) { ed->t2 = 0; k++; }
ed = t->e2;
if ((ed->t1)==t) { ed->t1 = 0; k++; }
else if ((ed->t2)==t) { ed->t2 = 0; k++; }
ed = t->e3;
if ((ed->t1)==t) { ed->t1 = 0; k++; }
else if ((ed->t2)==t) { ed->t2 = 0; k++; }
if (k!=check)
{
cout << "\n *** Tri: impossible state (DelTE), k = " << k;
cout << " (type: " << (int) t->type << ")\n";
printf("\n t (%d,%d,%d): %p",t->p1->node,t->p2->node,t->p3->node,t);
printf("\n e1 (%d-%d): %p",t->e1->p1->node,t->e1->p2->node,t->e1);
printf("\n e2 (%d-%d): %p",t->e2->p1->node,t->e2->p2->node,t->e2);
printf("\n e3 (%d-%d): %p",t->e3->p1->node,t->e3->p2->node,t->e3);
printf("\n t->e1->t1: %p",(t->e1)->t1);
printf("\n t->e1->t2: %p",(t->e1)->t2);
printf("\n t->e2->t1: %p",(t->e2)->t1);
printf("\n t->e2->t2: %p",(t->e2)->t2);
printf("\n t->e3->t1: %p",(t->e3)->t1);
printf("\n t->e3->t2: %p",(t->e3)->t2);
#if NOKASKPLOT == 0
FEPlotMESH2* fePlot = new FEPlotMESH2(this, SCREEN);
fePlot->plotElements();
fePlot->plotTriangleNodes();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e1);
fePlot->plotEdge(t->e2);
fePlot->plotEdge(t->e3);
fePlot->plot->flush();
Pause();
#endif
}
return;
}
//-------------------------------------------------------------------------
ostream& operator << (ostream& os, MESH2& t)
{
PT2* p;
EDG2* e;
TR2* tr;
os << "\nTriangulation";
if (t.fileName) os << " "<<t.fileName;
os <<":\n";
DListIter<PT2> iterPT (t.points());
DListIter<TR2> iterTR (t.triangles());
DListIter<EDG2> iterEDG(t.edges());
while ((p =iterPT.all())) p->print(cout);
while ((e =iterEDG.all())) e->print(cout);
while ((tr=iterTR.all())) tr->print(cout);
os << "\nRefinement step " << t.refStep;
os << " -- points: " << t.noOfPoints;
os << " edges: " << t.noOfEdges;
os << " triangles: " << t.noOfTriangles;
os << "\n";
return os;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void MESH2:: Info()
{
cout << "\n Triangulation: points edges triangles "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%10d%10d%10d", " ", noOfPoints, noOfEdges, noOfTriangles)
<< "\n";
}
//-------------------------------------------------------------------------
int MESH2:: MemSpace()
{
int space = varAllocator.MemSpace();
space += edgAlloc.MemSpace();
space += ptAlloc.MemSpace();
space += trAlloc.MemSpace();
return space;
}
//-------------------------------------------------------------------------
void MESH2:: refinementInfo(int* noPt, int* noEdg, int* noTr)
{
int space = MemSpace();
int blocks = varAllocator.NoOfBlocks();
blocks += edgAlloc.NoOfBlocks();
blocks += ptAlloc.NoOfBlocks();
blocks += trAlloc.NoOfBlocks();
cout << "\n Refine: points edges triangles "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%10d%10d%10d", "previous", noPt[0], noEdg[0], noTr[0]);
if (infoRefinement > 1)
{
cout << Form("\n%15s%10d%10d%10d", "removed",
noPt[1]-noPt[0], noEdg[1]-noEdg[0], noTr[1]-noTr[0])
<< Form("\n%15s%10d%10d%10d", "generated",
noPt[2]-noPt[1], noEdg[2]-noEdg[1], noTr[2]-noTr[1])
<< Form("\n%15s%10d%10d%10d", "added",
noOfPoints-noPt[2], noOfEdges-noEdg[2], noOfTriangles-noTr[2])
<< " for green closure";
}
cout << Form("\n%15s%10d%10d%10d Space : %1.6f MB","total",
noOfPoints, noOfEdges, noOfTriangles, float(space)/1e6);
cout << Form("\n%15s%10.2f%10.2f%10.2f Blocks: %1d\n","ratio",
float(noOfPoints)/noPt[0], float(noOfEdges)/noEdg[0],
float(noOfTriangles)/noTr[0], blocks);
if (infoRefinement > 2)
{
int npDepth, neDepth, step, np=0, ne=0, nt=0;
cout << "\n Grid Statistics:"
<< Form("\n%15s%10s%10s%10s%10s", "depth/step", "points", "edges",
"triangles","(total)");
FORALL (trList, step)
{
npDepth = (step<= ptList.h) ? ptList[step]->noOfElements : 0;
neDepth = (step<=edgList.h) ? edgList[step]->noOfElements : 0;
cout << Form("\n%15d%10d%10d%10d", step,
npDepth, neDepth, trList[step]->noOfElements);
np += npDepth;
ne += neDepth;
nt += trList[step]->noOfElements;
}
cout << Form("\n%15s%10d%10d%10d", "total sum", np, ne, nt);
cout << "\n";
}
}
//-------------------------------------------------------------------------
void MESH2:: consistencyCheck()
{
/*
TR2* t;
EDG2* e;
PT2* p;
DListIter<PT2> iterPT (points());
DListIter<TR2> iterTR (triangles());
DListIter<EDG2> iterEDG(edges());
while (p=iterPT.all())
{
if (p->depth != depth)
cout << "\n*** Check Pt: Point with wrong depth\n";
}
while (e=iterEDG.all())
{
if (e->depth != depth)
cout << "\n*** Check Edg: Edge with wrong depth\n";
}
while (t=iterTR.all())
{
if (t->depth != depth)
cout << "\n*** Check Tri: Tri with wrong depth\n";
if (t->type != T_GREEN)
{
if (t->e1->depth != depth)
{
cout << "\n*** Check Tri: edge 1 with wrong depth "
<< t->e1->depth << " " << depth << "\n";
FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
fePlot->plotElements();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e1);
fePlot->plot->flush();
Pause();
}
if (t->e2->depth != depth)
{
cout << "\n*** Check Tri: edge 2 with wrong depth "
<< t->e2->depth << " " << depth << "\n";
FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
fePlot->plotElements();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e2);
fePlot->plot->flush();
Pause();
}
if (t->e3->depth != depth)
{
cout << "\n*** Check Tri: edge 3 with wrong depth "
<< t->e3->depth << " " << depth << "\n";
FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
fePlot->plotElements();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e3);
fePlot->plot->flush();
Pause();
}
}
int k = 0;
ed = t->e1;
if ((ed->t1)==t) ++k;
if ((ed->t2)==t) ++k;
ed = t->e2;
if ((ed->t1)==t) ++k;
if ((ed->t2)==t) ++k;
ed = t->e3;
if ((ed->t1)==t) ++k;
if ((ed->t2)==t) ++k;
if (k != 3)
{
cout << "\n *** Check Tri: impossible state, k = " << k << "\n";
cout << "type: " << t->type << "\n";
printf("\n adress t: %p",t);
printf("\n t->e1->t1: %p",(t->e1)->t1);
printf("\n t->e1->t2: %p",(t->e1)->t2);
printf("\n t->e2->t1: %p",(t->e2)->t1);
printf("\n t->e2->t2: %p",(t->e2)->t2);
printf("\n t->e3->t1: %p",(t->e3)->t1);
printf("\n t->e3->t2: %p",(t->e3)->t2);
FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
fePlot->plotElements();
fePlot->plot->set(PENSIZE, BIG);
fePlot->plotEdge(t->e1);
fePlot->plotEdge(t->e2);
fePlot->plotEdge(t->e3);
fePlot->plot->flush();
Pause();
}
*/
}
//-------------------------------------------------------------------------
void MESH2:: Resolve(Real reqhMin, int nSteps)
{
Real hMin, h;
int depth;
PT2* pt;
TR2* tr;
EDG2* e;
Timer timer;
if (nSteps==0) // -> reqhMin will be used
{
// -- get minimal h: --
DListIter<EDG2> iterEDG(edges());
hMin = machMax(Real(0.0));
while ((e=iterEDG.all()))
{
h = e->lengthSqr();
if (h < hMin) hMin = h;
}
hMin = sqrt(hMin);
nSteps = Round(hMin/reqhMin);
}
if (nSteps > 10)
{
cout << "\n*** MESH2::Resolve: nSteps = " << nSteps << " too large\n";
cout.flush(); abort();
}
// resolve mesh:
RemoveGreen();
for (depth=0; depth < nSteps; ++depth)
{
DListIter<TR2> iter(triangles(), depth);
while ((tr=iter.all())) tr->mark = 1;
iter.reset();
while ((tr=iter.all())) RefineTriangle(tr);
}
/*
// delete lower levels:
for (depth=0; depth < nSteps; ++depth)
{
edgList[depth]->deleteAll();
trList [depth]->deleteAll();
}
edgList[0] = edgList[nSteps];
trList [0] = trList [nSteps];
edgList.h = 0;
trList.h = 0;
for (TrListProc tl(trList[0], all,&tr); tr; tr=tl.All())
{
tr->father = 0;
tr->depth = 0;
}
for (EdgListProc el(edgList[0],all,&e); e; e=el.All())
{
e->father = 0;
e->depth = 0;
}
for (depth=0; depth < nSteps; ++depth) ptList[0]->concat(*ptList[depth+1]);
ptList.h = 0;
*/
// set node numbers:
DListIter<PT2> iterPT(points());
int node = 0;
while ((pt=iterPT.all())) { pt->node = ++node; } // pt->depth = 0; }
/*
noOfPoints = node;
noOfEdges = edgList[0]->noOfElements;
noOfTriangles = trList [0]->noOfElements;
maxDepth = 0;
*/
MakeGreen();
if (Cmd.isTrue("timeResolve")) { cout << "\tResolve: "; timer.cpu(); }
if (infoRefinement) Info();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void PT2:: print(ostream& os) const
{
os << "Point " << node << " BC " << (int) boundP << " class "
<< (int) classA << " mark " << (int) mark
<< " x: " << x << " y: " << y << "\n";
}
void PT2:: reset()
{
boundP = INTERIOR;
mark = classA = depth = NUL;
node = 0;
}
void PT2:: readBC(BufferedParser& parser)
{
readBCValues(&boundP, &classA, parser);
}
//-------------------------------------------------------------------------
void PT2:: getNodeCoord(Vector<Real>& v) const { v[1]=x; v[2]=y; }
//-------------------------------------------------------------------------
void PT2:: getCoordinates(Vector<Real>& v) const { v[1]=x; v[2]=y; }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void EDG2:: print(ostream& os) const
{
os << "Edge with points " << p1->node << " " << p2->node
<< " type " << (int) type << " BC " << (int) boundP
<< " class " << (int) classA << " depth " << (int) depth << "\n";
}
void EDG2:: reset()
{
boundP = INTERIOR;
mark = classA = type = depth = isoP = NUL;
isBlueDiagonal = NUL;
p1 = p2 = pm = 0;
t1 = t2 = 0;
node = 0;
father = firstSon = 0;
}
//-------------------------------------------------------------------------
void EDG2:: compNormalVector(Vector<Real>& n, Real* length) const
{
Vector<Real> u(3),v(3),t(3);
p1->getCoordinates(u);
p2->getCoordinates(v);
t[1] = v[1] - u[1];
t[2] = v[2] - u[2];
n[1] = t[2]; // rotate 90 deg. clockwise
n[2] = -t[1];
*length = sqrt(n[1]*n[1]+n[2]*n[2]);
n[1] /= (*length);
n[2] /= (*length);
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void EDG2:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
Real L1 = 1.0 - u[1];
Real L2 = u[1];
x[1] = L1*p1->x + L2*p2->x;
x[2] = L1*p1->y + L2*p2->y;
}
//-------------------------------------------------------------------------
void EDG2:: getNodeCoord(Vector<Real>& v) const
{
v[1] = 0.5*(p1->x + p2->x);
v[2] = 0.5*(p1->y + p2->y);
}
//-------------------------------------------------------------------------
Real EDG2:: hMax() const
{
return sqrt(sqr(p2->x - p1->x) + sqr(p2->y - p1->y));
}
//-------------------------------------------------------------------------
void EDG2:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
}
//-------------------------------------------------------------------------
void EDG2:: readBC(BufferedParser& parser)
{
readBCValues(&boundP, &classA, &isoP, parser);
if (boundP == DIRICHLET)
{
p1->boundP = p2->boundP = boundP;
p1->classA = p2->classA = classA;
}
else if (boundP == CAUCHY || boundP == NEUMANN || boundP == INTERIOR)
{
if (p1->boundP != DIRICHLET)
{
p1->boundP = boundP;
p1->classA = classA;
}
if (p2->boundP != DIRICHLET)
{
p2->boundP = boundP;
p2->classA = classA;
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void TR2:: compJ(Jacobian& j, Bool checkDetJ) const
{
j.J(1,1) = p2->x - p1->x; // Jacobian j.J(i,k) = dx(k)/dr(i)
j.J(2,1) = p3->x - p1->x;
j.J(1,2) = p2->y - p1->y;
j.J(2,2) = p3->y - p1->y;
j.detJ = j.J(1,1)*j.J(2,2) - j.J(1,2)*j.J(2,1);
if (checkDetJ)
if (j.detJ <= 0.0)
{
cout << "\n*** Jacobi-Det.<= 0 " << j.detJ << " encountered\n";
cout.flush(); abort();
}
}
void TR2:: compJinv(Jacobian& j) const
{
int i,k;
compJ(j);
j.Jinv(1,1) = j.J(2,2); j.Jinv(1,2) = -j.J(1,2);
j.Jinv(2,1) = -j.J(2,1); j.Jinv(2,2) = j.J(1,1);
for (i=1; i<=2; ++i)
for (k=1; k<=2; ++k) j.Jinv(i,k) /= j.detJ;
}
//-------------------------------------------------------------------------
void TR2:: print(ostream& os) const
{
TR::print(os);
os << " type " << (int) type << " class " << (int)classA
<< " mark "<< (int) mark << " depth " << (int) depth << "\n";
}
void TR2:: reset()
{
p1 = p2 = p3 = 0;
e1 = e2 = e3 = 0;
mark = type = classA = depth = NUL;
boundP = INTERIOR;
father = firstSon = 0;
node = 0;
error = 0.0;
}
void TR2:: setBoundary()
{
if ( e1->onBoundary() || e2->onBoundary() || e3->onBoundary() )
boundP = BOUNDARY;
}
//-------------------------------------------------------------------------
int TR2:: NoOfSons() const
{
if (type == MESH::T_RED) return 4;
else if (type == MESH::T_GREEN) return 2;
else
{
cout << "\n*** TR2::NoOfSons(): unknown type " << type << "\n";
cout.flush(); abort();
return 0;
}
}
//-------------------------------------------------------------------------
void TR2:: getEdges(EDG2* eds[]) const
{
eds[1] = e1; eds[2] = e2; eds[3] = e3;
}
//-------------------------------------------------------------------------
void TR2:: getEdgeOrientation(Vector<Bool>& orient) const
{
int i;
FORALL(orient,i) orient[i] = False;
if (p1->node == 0 || p2->node == 0 || p3->node == 0) // !!!
{
cout <<"\n*** TR2:: getEdgeOrientation: node == 0 !\n";
cout.flush(); abort();
}
if (p2->node < p3->node) orient[1] = True;
if (p3->node < p1->node) orient[2] = True;
if (p1->node < p2->node) orient[3] = True;
}
//-------------------------------------------------------------------------
void TR2:: getAllFaces(Vector<PATCH*>& faces) const
{
faces[1] = e1; faces[2] = e2; faces[3] = e3;
}
//-------------------------------------------------------------------------
void TR2:: getNeighbours(Vector<PATCH*>& neighbours) const
{
if (e1->t1 != this) neighbours[1] = e1->t1;
else neighbours[1] = e1->t2;
if (e2->t1 != this) neighbours[2] = e2->t1;
else neighbours[2] = e2->t2;
if (e3->t1 != this) neighbours[3] = e3->t1;
else neighbours[3] = e3->t2;
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void TR2:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
Real L1 = 1.0 - u[1] - u[2];
Real L2 = u[1];
Real L3 = u[2];
x[1] = L1*p1->x + L2*p2->x + L3*p3->x;
x[2] = L1*p1->y + L2*p2->y + L3*p3->y;
}
//-------------------------------------------------------------------------
// -- point at x -> unit coordinates
void TR2:: unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const
{
int i, k;
Jacobian J;
Vector<Real> xn(2);
compJinv(J);
xn[1] = x[1] - p1->x;
xn[2] = x[2] - p1->y;
for (i=1; i<=2; ++i) // jacobian is transposed!
{
u[i] = 0.0;
for (k=1; k<=2; ++k) u[i] += J.Jinv(k,i) * xn[k];
}
}
//-------------------------------------------------------------------------
void TR2:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2,
Vector<Real>& cP3) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
p3->getCoordinates(cP3);
}
//-------------------------------------------------------------------------
Real TR2:: volume() const
{
Matrix<Real> J(2,2);
J(1,1) = p2->x - p1->x; // Jacobian J(i,k) = dx(k)/dr(i)
J(2,1) = p3->x - p1->x;
J(1,2) = p2->y - p1->y;
J(2,2) = p3->y - p1->y;
return 0.5 * Abs(J(1,1)*J(2,2) - J(1,2)*J(2,1));
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void MESH2:: resetElemIter(Bool lastStep)
{
if (elemIter != 0)
{
cout << "\n*** resetElemIter: elemIter != 0"
<< " \t (previous operation not completed?)\n";
cout.flush(); abort();
}
if (lastStep) iterLevel = trList.high();
else iterLevel = 0;
}
//-------------------------------------------------------------------------
PATCH* MESH2:: elemIterAll()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > trList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=trList[iterLevel]->first; elemIter;
elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > trList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH2:: elemIterAllOfHistory()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
return elemIter;
if (++iterLevel > trList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=trList[iterLevel]->first; elemIter;
elemIter=elemIter->Next()) return elemIter;
if (++iterLevel > trList.high()) return 0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
PATCH* MESH2:: edgIterAll()
{
if (edgIter != 0) // iterate rest of list
{
for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next())
if (!edgIter->refined()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (edgIter=edgList[edgIterLevel]->first; edgIter;
edgIter=edgIter->Next())
if (!edgIter->refined()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH2:: edgIterAllOfHistory()
{
if (edgIter != 0) // iterate rest of list
{
for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next())
return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (edgIter=edgList[edgIterLevel]->first; edgIter;
edgIter=edgIter->Next()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH2:: ptIterAll()
{
if (ptIter != 0) // iterate rest of list
{
for (ptIter=ptIter->Next(); ptIter; ptIter=ptIter->Next())
if (!ptIter->refined()) return ptIter;
if (++ptIterLevel > ptList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (ptIter=ptList[ptIterLevel]->first; ptIter;
ptIter=ptIter->Next())
if (!ptIter->refined()) return ptIter;
if (++ptIterLevel > ptList.high()) return 0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
MESH2BISECT:: MESH2BISECT (const char* inFileName, Bool /*readMesh*/)
: MESH2(inFileName), newTriangles(25)
{ }
//-------------------------------------------------------------------------
void MESH2BISECT:: Refine()
{
TR2 *t, *tNext;
int noPt[3], noEdg[3], noTr[3];
Timer timer, accTimer;
while (!newTriangles.empty()) newTriangles.pop();
noPt[0] = noOfPoints;
noTr[0] = noOfTriangles;
noEdg[0] = noOfEdges;
noPt[1] = noOfPoints;
noTr[1] = noOfTriangles;
noEdg[1] = noOfEdges;
DListIter<TR2> iter(triangles());
while ((t=iter.all())) EdgBisection(t);
int depth;
FORALL(trList,depth)
{
for (t=trList[depth]->first; t; t=tNext)
{
tNext = t->next; // a tri may be removed!
RefineTriangle(t);
}
}
int moreTrianglesToRefine;
while (1)
{
moreTrianglesToRefine = 0;
iter.reset();
while ((t=iter.all())) moreTrianglesToRefine += TestBisection(t);
if (moreTrianglesToRefine == 0) break;
if (infoRefinement>2)
cout << "\n moreTrianglesToRefine = " << moreTrianglesToRefine << "\n";
FORALL(trList,depth)
{
for (t=trList[depth]->first; t; t=tNext)
{
tNext = t->next; // a tri may be removed!
RefineTriangle(t);
}
}
}
noPt[2] = noOfPoints;
noTr[2] = noOfTriangles;
noEdg[2] = noOfEdges;
++refStep;
if (infoRefinement) refinementInfo (noPt, noEdg, noTr);
timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------
int MESH2BISECT:: RefineTriangle(TR2 *t)
{
if (t->mark) t->mark = False; // return if not marked
else return True;
if (t->firstSon) return True; // return if already refined
if (!RefineTriangleRed(t)) // refine it
{
cout << "\n*** RefineTriangle: triangle not bisected \n"; cout.flush(); abort();
}
return True;
}
//-------------------------------------------------------------------------
int MESH2BISECT:: TestBisection(TR2 *t)
{
EDG2* eds[1+3];
int k;
t->getEdges(eds);
for (k=1; k<=3; k++)
if (eds[k]->firstSon != 0) { t->mark=True; return 1; }
return 0;
}
//-------------------------------------------------------------------------
void MESH2BISECT:: EdgBisection(TR2 *t)
{
EDG2* eds[1+3];
int k;
if (t->mark)
{
t->getEdges(eds);
for (k=1; k<=3; k++) if (eds[k]->firstSon == 0) RefineEdge(eds[k]);
}
}
//-------------------------------------------------------------------------
/* *******************************************************************
MESH2BISECT::RefineTriangleRed does the actual bisection of a triangle
----------
TR2 *t Adress of the triangle
return int True or False
****************************************************************** */
int MESH2BISECT:: RefineTriangleRed(TR2 *t)
{
PT2 *p1=0, *p2=0;
EDG2 *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew, *e11;
TR2 *t1New, *t2New;
TR2* tFather = t->father;
int longestEdge;
if (!t->isVolatile()) // compute longest edge
{
Real lengthEdg1 = e1->hMax();
Real lengthEdg2 = e2->hMax();
Real lengthEdg3 = e3->hMax();
longestEdge = 1;
Real maxLength = lengthEdg1;
if (lengthEdg2 > maxLength) { longestEdge=2; maxLength=lengthEdg2; }
if (lengthEdg3 > maxLength) { longestEdge=3; }
}
else // a volatile son: the edge with lower depth MUST be refined
{
if (e1->depth < t->depth) longestEdge = 1;
else if (e2->depth < t->depth) longestEdge = 2;
else longestEdge = 3;
}
if (longestEdge==1)
{
p1 = t->p1; p2 = t->p2;
e1 = t->e1; e2 = t->e2; e3 = t->e3;
}
else if (longestEdge==2)
{
p1 = t->p2; p2 = t->p3;
e1 = t->e2; e2 = t->e3; e3 = t->e1;
}
else
{
p1 = t->p3; p2 = t->p1;
e1 = t->e3; e2 = t->e1; e3 = t->e2;
}
// Get one new edge and two new triangles
eNew = edgAlloc.Get();
++noOfEdges;
t1New = getTR2();
t2New = getTR2();
++noOfTriangles;
// Set the fields for the new objects
if (e1->firstSon == 0) e11 = RefineEdge(e1);
else e11 = e1->firstSon;
eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;
eNew->type = T_RED;
eNew->classA = t->classA;
eNew->depth = e11->depth;
t1New->e1 = e3; t1New->e3 = eNew;
t2New->e1 = e2; t2New->e2 = eNew;
if ((e1->p1)==p2) // test for orientation
{
t1New->e2 = e1->firstSon;
t2New->e3 = (e1->firstSon)->next;
}
else
{
t1New->e2 = (e1->firstSon)->next;
t2New->e3 = e1->firstSon;
}
DelTE(t,3);
if (!t->isVolatile()) // a non-volatile son creates two volatile sons
{
t1New->type = t2New->type = T_VOLATILE;
SetTP(t1New,t);
SetTP(t2New,t);
t->firstSon = t1New;
t->type = T_RED1;
if (t1New->depth > trList.h) trList.push(new DList<TR2>);
trList[t1New->depth]->add(t1New);
trList[t2New->depth]->add(t2New);
}
else // a volatile son creates two non-volatile sons and is removed
{
t1New->type = t2New->type = T_WHITE;
SetTP(t1New,tFather);
SetTP(t2New,tFather);
if (tFather->type == T_RED1) tFather->type = T_RED2;
else if (tFather->type == T_RED2) tFather->type = T_RED3;
else { cout << "\n*** error in tFather->type\n"; cout.flush(); abort(); }
if (tFather->firstSon == t) tFather->firstSon = t1New;
trList[t->depth]->substitute(t,t1New,t2New);
returnTR2(t);
}
SetTE(t1New);
SetTE(t2New);
newTriangles.push(t1New);
newTriangles.push(t2New);
if (eNew->depth > edgList.h) edgList.push(new DList<EDG2>);
edgList[eNew->depth]->add(eNew);
return True;
}
syntax highlighted by Code2HTML, v. 0.9.1