QuadLayer.cpp
#include "QuadLayer.h"
#include "QuadElement.h"
#include "coreGeometry.h"
#include <iostream>
void getGridline(NodeMap& allNodes, std::vector<Node*> &strip, MeshNode** grid, int offset, double *places)
{
double prevx = strip[0]->x, prevy = strip[0]->y, totLength = 0.0;
int len = strip.size(), i;
for (i = 0; i < len; i++)
{
totLength += distance(prevx, prevy, strip[i]->x, strip[i]->y);
places[i] = totLength;
prevx = strip[i]->x;
prevy = strip[i]->y;
grid[i * offset] = static_cast< MeshNode * >(strip[i]);
allNodes[strip[i]->tag] = strip[i];
}
for (i = 0; i < len; i++)
places[i] /= totLength;
}
// grid generation changed for curved bodies by Reino Ruusu 24.1.2000
void QuadLayer::
makeGrid( NodeMap& allNodes, const int m, const int n)
{
int i,j;
grid = new MeshNode*[m*n];
// relative place of nodes on edges (0..1)
double *places[4];
places[0] = new double[n];
places[1] = new double[m];
places[2] = new double[n];
places[3] = new double[m];
std::vector<Node*> strip;
// edge 0 is the column m
edges[0]->exportNodes(strip, directions[0]);
getGridline(allNodes, strip, &grid[m-1], m, places[0]);
strip.clear();
// edge 1 is the row n (reversed)
edges[1]->exportNodes(strip, -directions[1]);
getGridline(allNodes, strip, &grid[(n-1)*m], 1, places[1]);
strip.clear();
// edge 2 is the column 1 (reversed)
edges[2]->exportNodes(strip, -directions[2]);
getGridline(allNodes, strip, &grid[0], m, places[2]);
strip.clear();
// edge 3 is the row 1
edges[3]->exportNodes(strip, directions[3]);
getGridline(allNodes, strip, &grid[0], 1, places[3]);
strip.clear();
double *_offh = new double[8*m],
*_offv = new double[8*n];
double *xhll = _offh+0*m, *xhlr = _offh+1*m, *xhul = _offh+2*m, *xhur = _offh+3*m;
double *yhll = _offh+4*m, *yhlr = _offh+5*m, *yhul = _offh+6*m, *yhur = _offh+7*m;
double *xvll = _offv+0*n, *xvlr = _offv+1*n, *xvul = _offv+2*n, *xvur = _offv+3*n;
double *yvll = _offv+4*n, *yvlr = _offv+5*n, *yvul = _offv+6*n, *yvur = _offv+7*n;
// Precalculate offsets from corner points
// rows
for (i = 0; i < m; i++)
{
// LL
xhll[i] = grid[i]->x - grid[0]->x;
yhll[i] = grid[i]->y - grid[0]->y;
// LR
xhlr[i] = grid[i]->x - grid[m-1]->x;
yhlr[i] = grid[i]->y - grid[m-1]->y;
// UL
xhul[i] = grid[(n-1)*m+i]->x - grid[(n-1)*m]->x;
yhul[i] = grid[(n-1)*m+i]->y - grid[(n-1)*m]->y;
// UR
xhur[i] = grid[(n-1)*m+i]->x - grid[(n-1)*m+m-1]->x;
yhur[i] = grid[(n-1)*m+i]->y - grid[(n-1)*m+m-1]->y;
}
// columns
for (j = 0; j < n; j++)
{
// LL
xvll[j] = grid[j*m]->x - grid[0]->x;
yvll[j] = grid[j*m]->y - grid[0]->y;
// UL
xvul[j] = grid[j*m]->x - grid[(n-1)*m]->x;
yvul[j] = grid[j*m]->y - grid[(n-1)*m]->y;
// LR
xvlr[j] = grid[j*m+m-1]->x - grid[m-1]->x;
yvlr[j] = grid[j*m+m-1]->y - grid[m-1]->y;
// UR
xvur[j] = grid[j*m+m-1]->x - grid[(n-1)*m+m-1]->x;
yvur[j] = grid[j*m+m-1]->y - grid[(n-1)*m+m-1]->y;
}
for( j = 1; j < (n - 1); ++j)
{
for( i = 1; i < (m - 1); ++i)
{
double pl = places[2][j],
pr = places[0][j],
pb = places[3][i],
pt = places[1][i];
/* Calculate weights for interpolation
wt = wl * pl + wr * pr
wb = 1 - wt
wr = wb * pb + wt * pt
wl = 1 - wr
*/
double wt = ((pr - pl) * pb + pl) / (1.0 - (pr - pl) * (pt - pb));
double wb = 1.0 - wt;
double wr = wb * pb + wt * pt;
double wl = 1.0 - wr;
// Interpolate coordinates from the four gridline endpoints
Node *left = grid[j*m], *right = grid[j*m+m-1],
*bottom = grid[i], *top = grid[(n-1)*m+i];
double x, y;
x = ((wl * (left->x + wb * xhll[i] + wt * xhul[i]) +
wr * (right->x + wb * xhlr[i] + wt * xhur[i])) +
(wb * (bottom->x + wl * xvll[j] + wr * xvlr[j]) +
wt * (top->x + wl * xvul[j] + wr * xvur[j]))) / 2.0;
y = ((wl * (left->y + wb * yhll[i] + wt * yhul[i]) +
wr * (right->y + wb * yhlr[i] + wt * yhur[i])) +
(wb * (bottom->y + wl * yvll[j] + wr * yvlr[j]) +
wt * (top->y + wl * yvul[j] + wr * yvur[j]))) / 2.0;
grid[ j*m+i ] = new MeshNode( x, y );
allNodes[grid[ j*m+i ]->tag] = grid[ j*m+i ];
}
}
delete[] _offh;
delete[] _offv;
for (i = 0; i < 4; i++)
delete[] places[i];
}
void QuadLayer::
discretize(NodeMap& fixedNodes, NodeMap& allNodes, std::list< Element* >& allElements)
{
int i, j;
int m = edges[1]->size();
int n = edges[0]->size();
QuadElement **elements = new QuadElement*[(n-1)*(m-1)];
makeGrid( allNodes, m, n );
for( j = 0; j < (n - 1); ++j)
for( i = 0; i < (m - 1); ++i)
{
QuadElement *el = new QuadElement( grid[j*m+i+1], grid[(j+1)*m+i+1], grid[(j+1)*m+i], grid[j*m+i]);
elements[j*(m-1)+i] = el;
allElements.push_back(el);
}
// set boundary element relations
std::vector<BoundaryElement*> bels;
// edge 0 is the column m
edges[0]->elements(bels, directions[0]);
for (j = 0; j < n - 1; j++)
{
bels[j]->setLeft(elements[j*(m-1)+m-2]->elementId());
}
bels.clear();
// edge 1 is the row n (reversed)
edges[1]->elements(bels, -directions[1]);
for (i = 0; i < m - 1; i++)
{
bels[i]->setRight(elements[(n-2)*(m-1)+i]->elementId());
}
bels.clear();
// edge 2 is the column 1 (reversed)
edges[2]->elements(bels, -directions[2]);
for (j = 0; j < n - 1; j++)
{
bels[j]->setRight(elements[j*(m-1)]->elementId());
}
bels.clear();
// edge 3 is the row 1
edges[3]->elements(bels, directions[3]);
for (i = 0; i < m - 1; i++)
{
bels[i]->setLeft(elements[i]->elementId());
}
delete [] grid;
delete [] elements;
}