BGGridMesh.cpp
#include "BGMesh.h"
#include "BGVertex.h"
#include <math.h>
#include <algorithm>
#include "coreGeometry.h"
void BGGridMesh::
initialize(std::vector< GeometryNode* >& nodes, std::vector< GeometryEdge* >& edges)
{
initialized = true;
std::cout << "Calculating background mesh" << std::endl;
int i;
int len = nodes.size();
double mean = MAP(nodes[0]->delta);
double mindelta = nodes[0]->delta;
double minx = nodes[0]->x, maxx = minx, miny = nodes[0]->y, maxy = miny;
for( i = 1; i < len; ++i )
{
mean += MAP(nodes[i]->delta);
if (nodes[i]->delta < mindelta)
mindelta = nodes[i]->delta;
if (nodes[i]->x < minx)
minx = nodes[i]->x;
if (nodes[i]->x > maxx)
maxx = nodes[i]->x;
if (nodes[i]->y < miny)
miny = nodes[i]->y;
if (nodes[i]->y > maxy)
maxy = nodes[i]->y;
}
mean /= len;
width = maxx - minx;
height = maxy - miny;
ox = minx;
oy = miny;
cellsize = mindelta * cellsize;
nh = (int)ceil(width / cellsize) + 1;
width = (nh - 1) * cellsize;
nv = (int)ceil(height / cellsize) + 1;
height = (nv - 1) * cellsize;
nh += 2;
nv += 2;
ox -= cellsize;
oy -= cellsize;
width += 2.0 * cellsize;
height += 2.0 * cellsize;
grid = new double[nh * nv];
bool *mask = new bool[nh * nv];
for (i = 0; i < nh * nv; i++)
{
grid[i] = mean;
mask[i] = false;
}
for (i = 0; i < len; i++)
{
int h = (int)((nodes[i]->x - ox) / cellsize + 0.5);
int v = (int)((nodes[i]->y - oy) / cellsize + 0.5);
grid[v*nh+h] = MAP(nodes[i]->delta);
mask[v*nh+h] = true;
}
len = edges.size();
for (i = 0; i < len; i++)
{
std::vector<int> x, y;
std::vector<double> delta;
edges[i]->getGridPoints(ox, oy, cellsize, x, y, delta);
for (int j = 0; j < x.size(); j++)
{
grid[y[j]*nh+x[j]] = delta[j];
mask[y[j]*nh+x[j]] = true;
}
}
std::vector<std::pair<int, int> > stack;
stack.push_back(std::make_pair(0, 0));
while (!stack.empty())
{
std::pair<int, int> p = stack.back();
stack.pop_back();
int h = p.first, v = p.second;
mask[v*nh+h] = true;
if (h > 0 && !mask[v*nh+h-1])
stack.push_back(std::make_pair(h-1, v));
if (h < nh - 1 && !mask[v*nh+h+1])
stack.push_back(std::make_pair(h+1, v));
if (v > 0 && !mask[(v-1)*nh+h])
stack.push_back(std::make_pair(h, v-1));
if (v < nv - 1 && !mask[(v+1)*nh+h])
stack.push_back(std::make_pair(h, v+1));
}
int niters = 0;
double error = 1e100, last;
do
{
niters++;
last = error;
error = 0.0;
int i, j;
for (int c = 0; c < 200; c++)
{
for (j = 0; j < nv; j++)
{
for (i = 0; i < nh; i++)
{
if (mask[j*nh+i]) continue;
double v = 0.0;
int n = 0;
if (i > 0)
{
v += grid[j*nh+i-1];
n++;
}
if (i < nh - 1)
{
v += grid[j*nh+i+1];
n++;
}
if (j > 0)
{
v += grid[(j-1)*nh+i];
n++;
}
if (j < nv - 1)
{
v += grid[(j+1)*nh+i];
n++;
}
grid[j*nh+i] = v / n;
}
}
}
for (j = 0; j < nv; j++)
{
for (i = 0; i < nh; i++)
{
if (mask[j*nh+i]) continue;
double residual = 0.0;
double v = grid[j*nh+i];
if (i > 0)
residual += grid[j*nh+i-1] - v;
if (i < nh - 1)
residual += grid[j*nh+i+1] - v;
if (j > 0)
residual += grid[(j-1)*nh+i] - v;
if (j < nv - 1)
residual += grid[(j+1)*nh+i] - v;
residual /= cellsize;
double e = fabs(residual);
if (e > error)
error = e;
}
}
} while (error > 0.00001);
}
double BGGridMesh::
interpolate( const double ix, const double iy )
{
int h = (int)((ix - ox) / cellsize);
int v = (int)((iy - oy) / cellsize);
if (h < 0) h = 0;
if (h > nh - 1) h = nh - 1;
if (v < 0) v = 0;
if (v > nv - 1) v = nv - 1;
double xoff = (ix - (ox + h * cellsize)) / cellsize,
yoff = (iy - (oy + v * cellsize)) / cellsize;
return UNMAP((1.0 - xoff) * (1.0 - yoff) * grid[v*nh+h] +
xoff * (1.0 - yoff) * grid[v*nh+h+1] +
(1.0 - xoff) * yoff * grid[(v+1)*nh+h] +
xoff * yoff * grid[(v+1)*nh+h+1]);
}