/*
This file is part of the FElt finite element analysis package.
Copyright (C) 1993-2000 Jason I. Gobat and Darren C. Atkinson
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
/****************************************************************************
*
* File: contour.c
*
***************************************************************************/
# include <stdio.h>
# include <string.h>
# include <ctype.h>
# include "fe.h"
# include "allocate.h"
# include "error.h"
# include "region.h"
# include "colormap.h"
# define HORIZ 1
# define VERT 2
# define WIDTH 400
# define HEIGHT 400
# define NUMCELLS 240
# define MaxNodesPerElement 24
static unsigned char red [256];
static unsigned char green [256];
static unsigned char blue [256];
static unsigned char white, black;
# define Sign(x) ((x) >= 0 ? 1 : -1)
static void BresenhamLine (xs, ys, xe, ye, image, nrows, ncols)
int xs, ys, xe, ye;
unsigned char **image;
int nrows;
int ncols;
{
unsigned i;
int x,y,
dx,dy,
e;
int signx, signy,
temp;
unsigned change;
if (xe == xs && ye == ys)
return;
signx = Sign (xe - xs);
signy = Sign (ye - ys);
dx = signx*(xe - xs);
dy = signy*(ye - ys);
if (dy <= dx)
change = 0;
else {
change = 1;
temp = dx;
dx = dy;
dy = temp;
}
e = 2*dy - dx;
x = xs;
y = ys;
for (i = 0 ; i <= dx ; i++) {
if (y > 0 && y < nrows && x > 0 && x < ncols)
image [y][x] = black;
while (e >= 0) {
if (!change)
y += signy;
else
x += signx;
e -= 2*dx;
}
if (!change)
x += signx;
else
y += signy;
e += 2*dy;
}
return;
}
static void InitializeColormapTable ( )
{
int i;
for (i = 0 ; i < NUMCELLS ; i++) {
red [i] = color_values [i].r;
green [i] = color_values [i].g;
blue [i] = color_values [i].b;
}
white = NUMCELLS;
black = white + 1;
red [white] = 255;
blue [white] = 255;
green [white] = 255;
red [black] = 0;
blue [black] = 0;
green [black] = 0;
for (i = black + 1 ; i < 256 ; i++) {
red [i] = 0;
blue [i] = 0;
green [i] = 0;
}
return;
}
static void HequalImageData (image, nrows, ncols)
unsigned char **image;
int nrows;
int ncols;
{
unsigned i,j;
unsigned char point;
int histogram [NUMCELLS];
float p [NUMCELLS],
v [NUMCELLS],
vstar [NUMCELLS],
vmin, area, sum;
for (i = 0 ; i < NUMCELLS ; i++)
histogram [i] = 0;
area = 0;
for (i = 0 ; i < nrows ; i++) {
for (j = 0 ; j < ncols ; j++) {
point = image [i][j];
if (point != white) {
area++;
histogram [point]++;
}
}
}
for (i = 0 ; i < NUMCELLS ; i++)
p [i] = histogram [i] / area;
sum = 0.0;
vmin = p[0];
for (i = 0 ; i < NUMCELLS ; i++) {
sum += p[i];
v[i] = sum;
if (v[i] < vmin)
vmin = v[i];
}
for (i = 0 ; i < NUMCELLS ; i++)
vstar[i] = ((((v[i] - vmin)/(1.0 - vmin))*(NUMCELLS - 1.0)) + 0.5);
for (i = 0 ; i < nrows ; i++) {
for (j = 0 ; j < ncols ; j++) {
point = image [i][j];
if (point != white) {
point = (unsigned char) vstar [point];
image [i][j] = point;
}
}
}
}
static unsigned char **LoadImage (z, width, height, mask, equalize, wedge_orient)
float **z;
unsigned char **mask;
int width;
int height;
int equalize;
unsigned wedge_orient;
{
unsigned x_max,y_max,x_min,y_min;
float max,min;
unsigned i,j;
float slope;
int point;
unsigned flag;
int max_width;
int min_width;
int sx1,sy1,sx2,sy2;
unsigned char **image;
x_max = x_min = y_max = y_min = 0;
max = min = max_width = min_width = 0;
flag = 0;
for (i = 0 ; i < height ; i++) {
for (j = 0 ; j < width ; j++) {
if (mask[i][j]) {
if (!flag) {
x_max = x_min = j;
y_max = y_min = i;
max = min = z[i][j];
flag = 1;
}
else {
if (z[i][j] > max)
max = z[i][j];
else if (z[i][j] < min)
min = z[i][j];
if (i > y_max)
y_max = i;
else if (i < y_min)
y_min = i;
if (j > x_max)
x_max = j;
else if (j < x_min)
x_min = j;
}
}
}
}
image = Allocate (unsigned char *, height);
for (i = 0 ; i < height ; i++)
image [i] = Allocate (unsigned char, width);
slope = (float) NUMCELLS / (max - min);
for (i = 0 ; i < height ; i++) {
for (j = 0 ; j < width ; j++) {
if (mask[i][j]) {
if (z[i][j] != max)
point = (int) (slope*(z[i][j] - min));
else
point = NUMCELLS - 1;
}
else
point = white;
image [height - i - 1][j] = point;
}
}
if (wedge_orient == HORIZ) {
for (j = x_min ; j <= x_max ; j++) {
point = (int) ((float) (j-x_min)/
(float) (x_max - x_min)*(NUMCELLS - 1.0));
for (i = (int) (0.88*height) ; i < height - 12 ; i++)
image [i][j] = point;
}
sx1 = x_min;
sx2 = x_max - max_width;
sy1 = sy2 = height - 2;
}
else {
for (i = y_min ; i <= y_max ; i++) {
point = (int) ((float) (i-y_min)/
(float) (y_max - y_min)*(NUMCELLS - 1.0));
for (j = (int) (0.9*width) ; j < (int) (0.98*width) ; j++)
image [height - i - 1][j] = point;
}
sy1 = y_max + 10;
sy2 = y_min - 2;
sx1 = width - min_width - 1;
sx2 = width - max_width - 1;
}
if (equalize)
HequalImageData (image, height, width);
return image;
}
static void PlotContourField (filename, node, numnodes, element,
numelts, nd, component, equalize, plot_elt,
width, height)
char *filename;
Node *node;
unsigned numnodes;
Element *element;
unsigned numelts,
nd;
int component;
int equalize;
int plot_elt;
int width;
int height;
{
unsigned i,j;
unsigned n;
float x_min, y_min,
x_max, y_max;
float maxX,maxY,
minX,minY;
float Xscale, Yscale;
float *xd,*yd,*zd;
float *xi,*yi,**zi;
unsigned result;
unsigned char **mask;
unsigned wedge_orient;
Region current, global;
XPoint rgpnt [MaxNodesPerElement];
unsigned char **image;
InitializeColormapTable ();
/*
* some preliminaries to set-up a properly
* proportioned window
*/
x_min = x_max = element[1] -> node[1] -> x;
y_min = y_max = element[1] -> node[1] -> y;
for (i = 1 ; i <= numelts ; i++) {
for (j = 1 ; j <= element[i] -> definition -> shapenodes ; j++) {
if (element[i] -> node[j] -> x > x_max)
x_max = element[i] -> node[j] -> x;
else if (element[i] -> node[j] -> x < x_min)
x_min = element[i] -> node[j] -> x;
if (element[i] -> node[j] -> y > y_max)
y_max = element[i] -> node[j] -> y;
else if (element[i] -> node[j] -> y < y_min)
y_min = element[i] -> node[j] -> y;
}
}
minX = x_min - 0.05*(x_max - x_min);
maxX = x_max + 0.05*(x_max - x_min);
minY = y_min - 0.05*(y_max - y_min);
maxY = y_max + 0.05*(y_max - y_min);
if (x_max - x_min < y_max - y_min) {
maxX += 0.25*(x_max - x_min);
wedge_orient = VERT;
}
else {
minY -= 0.25*(y_max - y_min);
wedge_orient = HORIZ;
}
if (width && height) {
Xscale = (float) width / (maxX - minX);
Yscale = (float) height / (maxY - minY);
}
else {
if (maxX - minX > maxY - minY) {
Xscale = (float) WIDTH / (maxX - minX);
width = WIDTH;
Yscale = Xscale;
height = Yscale*(maxY - minY);
}
else {
Yscale = (float) HEIGHT / (maxY - minY);
height = HEIGHT;
Xscale = Yscale;
width = Xscale*(maxX - minX);
}
}
/*
* set-up to do the bivariate interpolation so
* we can fill in the entire stress field
*/
zi = Allocate (float *, height);
mask = Allocate (unsigned char *, height);
for (i = 0 ; i < height ; i++) {
zi [i] = Allocate (float, width);
mask [i] = Allocate (unsigned char, width);
}
xi = Allocate (float, width);
UnitOffset (xi);
yi = Allocate (float, height);
UnitOffset (yi);
if (numnodes != 0 && node != NULL)
nd = numnodes;
xd = Allocate (float, nd);
yd = Allocate (float, nd);
zd = Allocate (float, nd);
if (xd == NULL || yd == NULL || zd == NULL)
Fatal ("allocation error plotting contour field");
UnitOffset (xd);
UnitOffset (yd);
UnitOffset (zd);
global = XCreateRegion ();
n = 0;
for (i = 1 ; i <= numelts ; i++) {
for (j = 1; j <= element[i] -> definition -> shapenodes ; j++) {
rgpnt [j - 1].x = (int) ((element[i] -> node[j] -> x - minX)*Xscale);
rgpnt [j - 1].y = (int) ((element[i] -> node[j] -> y - minY)*Yscale);
}
current=XPolygonRegion (rgpnt, element[i] -> definition -> shapenodes,
EvenOddRule);
XUnionRegion (global, current, global);
XDestroyRegion (current);
if (numnodes == 0) {
for (j = 1; j <= element[i] -> ninteg ; j++) {
xd[n+1] = (element[i] -> stress[j] -> x - minX)*Xscale;
yd[n+1] = (element[i] -> stress[j] -> y - minY)*Yscale;
zd[n+1] = element[i] -> stress[j] -> values[component];
n++;
}
}
}
if (numnodes > 0 && node != NULL) {
for (i = 1 ; i <= numnodes ; i++) {
xd[i] = (node[i] -> x - minX)*Xscale;
yd[i] = (node[i] -> y - minY)*Yscale;
zd[i] = node[i] -> dx [component];
}
}
for (i = 0 ; i < height ; i++) {
for (j = 0 ; j < width ; j++) {
if (XPointInRegion(global,j,i))
mask[i][j] = 1;
else
mask[i][j] = 0;
}
}
XDestroyRegion (global);
for (i = 1 ; i <= height; i++)
yi [i] = (float) i - 1.0;
for (j = 1 ; j <= width ; j++)
xi [j] = (float) j - 1.0;
result = BivariateInterp (nd, xd, yd, zd, width, height,xi,yi,zi,mask);
if (result) {
error ("could not interpolate to form stress field");
image = NULL;
}
else { /* interpolation ok - draw everything */
image = LoadImage (zi, width, height, mask, equalize, wedge_orient);
if (plot_elt) {
for (i = 1 ; i <= numelts ; i++) {
for (j = 1 ; j < element[i] -> definition -> shapenodes ; j++)
BresenhamLine((int) ((element[i] -> node[j] -> x - minX)*Xscale),
height - (int) ((element[i] -> node[j] -> y - minY)*Yscale),
(int) ((element[i] -> node[j+1] -> x - minX)*Xscale),
height - (int) ((element[i] -> node[j+1] -> y - minY)*Yscale),
image, height, width);
BresenhamLine((int) ((element[i] -> node[j] -> x - minX)*Xscale),
height - (int) ((element[i] -> node[j] -> y - minY)*Yscale),
(int) ((element[i] -> node[1] -> x - minX)*Xscale),
height - (int) ((element[i] -> node[1] -> y - minY)*Yscale),
image, height, width);
}
}
}
/*
* cleanup
*/
ZeroOffset (xd); Deallocate (xd);
ZeroOffset (yd); Deallocate (yd);
ZeroOffset (zd); Deallocate (zd);
ZeroOffset (xi); Deallocate (xi);
ZeroOffset (yi); Deallocate (yi);
for (i = 0 ; i < height ; i++) {
Deallocate (zi [i]);
Deallocate (mask [i]);
}
Deallocate (mask);
Deallocate (zi);
if (image)
ImageDataToBMP(filename, image, height, width, red, green, blue);
return;
}
void PlotStressField (out, element, numelts, comp, equalize, plot_elt,
width, height)
char *out;
Element *element;
unsigned numelts;
int comp;
int equalize;
int plot_elt;
int width;
int height;
{
int nd;
int i, j;
int flag;
nd = 0;
flag = 1;
for (i = 1; i <= numelts ; i++) {
if (element [i] -> definition -> shape != Planar) {
error ("cannot plot stresses for non-planar elements");
return;
}
if (element [i] -> stress == NULL) {
error ("could not get stresses for all elements");
return;
}
if (comp > element [i] -> definition -> numstresses) {
error ("invalid stress component for element %d",i);
return;
}
if (flag) {
for (j = 1 ; j <= element [i] -> ninteg ; j++) {
if (element [i] -> stress [j] -> values [comp] != 0)
flag = 0;
}
}
nd += element [i] -> ninteg;
}
if (flag) {
error ("all stresses are zero for specified component");
return;
}
PlotContourField (out, NULL, 0, element, numelts, nd, comp,
equalize, plot_elt, width, height);
}
void PlotDisplacementField (out, node, numnodes, element, numelts, comp,
equalize, plot_elt, width, height)
char *out;
Node *node;
unsigned numnodes;
Element *element;
unsigned numelts;
int comp;
int equalize;
int plot_elt;
int width;
int height;
{
int i;
int flag;
for (i = 1; i <= numelts ; i++) {
if (element [i] -> definition -> shape != Planar) {
error ("cannot plot stresses for non-planar elements");
return;
}
}
flag = 1;
for (i = 1 ; i <= numnodes ; i++) {
if (node [i] -> dx [comp] != 0) {
flag = 0;
break;
}
}
if (flag) {
error ("all displacements are zero for specified component");
return;
}
PlotContourField (out, node, numnodes, element, numelts, 0, comp,
equalize, plot_elt, width, height);
}
syntax highlighted by Code2HTML, v. 0.9.1