GeometryEdge.cpp
#include "GeometryEdge.h"
#include "BGMesh.h"
#include <math.h>
#include <algorithm>
#include <iostream>
#include "coreGeometry.h"
static int nextTag = 1;
GeometryEdge::GeometryEdge( int nSegments )
{
tag = nextTag;
segments = nSegments;
boundaryTag = 0;
type = OUTER;
++nextTag;
}
GeometryEdge::GeometryEdge( const int t, int nSegments )
{
tag = t;
segments = nSegments;
boundaryTag = 0;
type = OUTER;
++nextTag;
}
void GeometryEdge::
exportNodes(std::vector<Node*>& strip, int direction)
{
if( direction < 0 )
{
std::reverse_copy( nodes.begin(), nodes.end(), std::back_inserter( strip ) );
}
else
{
std::copy( nodes.begin(), nodes.end(), std::back_inserter( strip ) );
}
}
void GeometryEdge::
elements(std::vector< BoundaryElement * >& strip, int direction)
{
int i, len = bels.size();
if( len == 0 ) // Not constructed yet
{
int nlen = nodes.size();
for( i = 0; i < (nlen - 1); ++i )
{
bels.push_back( new BoundaryElement( boundaryTag, nodes[i], nodes[i+1] ) );
}
len = bels.size();
}
if( direction < 0 )
{
for( i = 0; i < len; ++i ) bels[i]->flip( true );
std::reverse_copy( bels.begin(), bels.end(), std::back_inserter( strip ) );
}
else
{
for( i = 0; i < len; ++i ) bels[i]->flip( false );
std::copy( bels.begin(), bels.end(), std::back_inserter( strip ) );
}
}
void GeometryEdge::
midNodes( Node*& a, Node*& b, int dir )
{
int len = nodes.size();
int ta, tb;
ta = len / 2 - 1;
tb = ta + 1;
if( dir == 1 )
{
a = nodes[ta];
b = nodes[tb];
}
else
{
a = nodes[tb];
b = nodes[ta];
}
}
double GeometryEdge::interpolate(double x, double y)
{
double value = 0.0;
for (std::vector<BGMesh*>::iterator it = bgMeshes.begin(); it != bgMeshes.end(); it++)
{
double h = (*it)->interpolate(x, y);
if (value == 0.0 || h < value)
value = h;
}
return value;
}
void GeometryEdge::
discretize( NodeMap &allNodes )
{
if (segments == 0)
{
int len = dots.size() - 1;
for (int i = 0; i < len; i++)
{
discretizeGradedSegment(allNodes, dots[i], dots[i+1]);
if (i < len - 1)
nodes.pop_back();
}
}
else
{
int len = dots.size() - 1, i;
if (segments < len) segments = len;
double totLength = 0.0;
double *lengths = new double[len];
for (i = 0; i < len; i++)
{
lengths[i] = distance(dots[i]->x, dots[i]->y, dots[i+1]->x, dots[i+1]->y);
totLength += lengths[i];
}
int *segs = new int[len];
for (i = 0; i < len; i++)
segs[i] = 1;
int remaining = segments - len;
while (remaining > 0)
{
int longest = 0;
double length = lengths[0] / segs[0];
for (i = 1; i < len; i++)
{
double l = lengths[i] / segs[i];
if (l > length)
{
length = l;
longest = i;
}
}
segs[longest]++;
remaining--;
}
delete [] lengths;
for (i = 0; i < len; i++)
{
discretizeConstantSegment(segs[i], allNodes, dots[i], dots[i+1]);
if (i < len - 1)
nodes.pop_back();
}
delete [] segs;
}
}
void GeometryEdge::
discretizeConstantSegment( int nSeg, NodeMap& allNodes, GeometryNode *from, GeometryNode *to )
{
Node *first = allNodes[from->tag];
Node *last = allNodes[to->tag];
double dx = last->x - first->x;
double dy = last->y - first->y;
double length = sqrt(dx*dx + dy*dy);
double sx = dx / length;
double sy = dy / length;
double step = length / nSeg;
double bx = first->x;
double by = first->y;
from->setDelta( step );
to->setDelta( step );
nodes.push_back( first );
for( int i = 1; i < nSeg; ++i )
{
MeshNode *nd = new MeshNode(bx+i*step*sx, by+i*step*sy);
nd->boundarynode = true;
nodes.push_back( nd );
}
nodes.push_back( last );
}
void GeometryEdge::
discretizeGradedSegment( NodeMap& allNodes, GeometryNode *from, GeometryNode *to )
{
int i;
double lineLength = distance( from->x, from->y, to->x, to->y );
double grading[2];
double sx = (to->x - from->x) / lineLength;
double sy = (to->y - from->y) / lineLength;
double limit = lineLength;
double at = 0.0;
double predicted = interpolate(from->x, from->y);
std::vector< double > lengths;
int step = 0;
while( 1 )
{
double prev;
do
{
prev = predicted;
double center = at + predicted / 2.0;
// Avoid interpolating outside the domain
if (center > lineLength)
center = lineLength;
predicted = 0.5 * (predicted + interpolate(from->x + center * sx, from->y + center * sy));
}
while( fabs((prev - predicted)/lineLength) > 0.0000001 );
at += predicted;
if (at > limit - 0.01 * predicted)
break;
lengths.push_back( at );
++step;
}
double mult = lineLength / at;
for( i = 0; i < step; ++i)
{
lengths[i] *= mult;
}
Node *first = allNodes[from->tag];
Node *last = allNodes[to->tag];
nodes.push_back( first );
for( i = 0; i < step; ++i )
{
MeshNode *nd = new MeshNode( lengths[i]*sx + from->x, lengths[i]*sy + from->y );
nd->boundarynode = true;
nodes.push_back( nd );
}
nodes.push_back( last );
}
std::ostream& operator<< (std::ostream& o, const GeometryEdge& A)
{
const int len = A.nodes.size();
o << A.tag << ' ' << len << ' ';
for(int i = 0; i < len; ++i)
{
o << A.nodes[i]->tag << ' ';
}
o << std::endl;
return o;
}
void GeometryEdge::
getGridPoints(double ox, double oy, double cellsize, std::vector<int> &x, std::vector<int> &y, std::vector<double> &delta)
{
int i, len = dots.size();
for (i = 0; i < len - 1; i++)
{
double ax = dots[i]->x, ay = dots[i]->y;
double bx = dots[i+1]->x, by = dots[i+1]->y;
double ad = MAP(dots[i]->delta), bd = MAP(dots[i+1]->delta);
int h = (int)((ax - ox) / cellsize + 0.5);
int v = (int)((ay - oy) / cellsize + 0.5);
x.push_back(h);
y.push_back(v);
delta.push_back(ad);
double dx = bx - ax;
double dy = by - ay;
double d = sqrt(dx * dx + dy * dy);
dx /= d;
dy /= d;
double dd = bd - ad;
int endh = (int)((bx - ox) / cellsize + 0.5);
int endv = (int)((by - oy) / cellsize + 0.5);
if (endh == h && endv == v)
continue;
int hstep = endh > h ? 1 : -1;
int vstep = endv > v ? 1 : -1;
if (abs(endh - h) > abs(endv - v))
{
int pos = 0, num = abs(endv - v), div = abs(endh - h);
while ((h += hstep) != endh)
{
pos += num;
if (pos > div / 2)
{
v += vstep;
pos -= div;
}
double nx = ox + h * cellsize;
double ny = oy + v * cellsize;
double nd = ad + ((nx - ax) * dx + (ny - ay) * dy) / d * dd;
x.push_back(h);
y.push_back(v);
delta.push_back(nd);
}
}
else
{
int pos = 0, num = abs(endh - h), div = abs(endv - v);
while ((v += vstep) != endv)
{
pos += num;
if (pos > div / 2)
{
h += hstep;
pos -= div;
}
double nx = ox + h * cellsize;
double ny = oy + v * cellsize;
double nd = ad + ((nx - ax) * dx + (ny - ay) * dy) / d * dd;
x.push_back(h);
y.push_back(v);
delta.push_back(nd);
}
}
if (i == len - 2)
{
x.push_back(endh);
y.push_back(endv);
delta.push_back(bd);
}
}
}