/* $Log: post.c,v $
* Revision 1.1 1995/02/21 07:53:52 martenl
* Initial revision
*
* Revision 1.5 1994/10/26 13:01:55 martenl
* Added memory allocation for vecX and vecY
*
* Revision 1.4 1994/10/26 12:58:59 martenl
* Added plotVectors
*
* Revision 1.3 1994/10/25 09:25:56 martenl
* Major cleanup
*
* Revision 1.2 1994/10/14 13:24:23 martenl
* (readSolution): Added a filename argument
* */
#define min(x, y) (((x) < (y)) ? (x) : (y))
#define max(x, y) (((x) > (y)) ? (x) : (y))
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include "libsx.h"
#include "globals.h"
#include "geodata.h"
#include "kraftwerk.h"
#include "grid.h"
#include "crosshair.h"
#include "graph.h"
int setUpMycolorMap(int ncolor);
void lingraph();
void dolingraph();
float xlin1, xlin2, ylin1, ylin2;
float xold[1000], fold[1000];
int xplin1, xplin2, yplin1, yplin2;
int nlin, nold;
int linitem;
void drawText(char *text, float x, float y);
int clickedInElement(int x, int y, int *thisEl);
void updateLevelButtons()
{
int i;
for(i = 0; i < MAXLEVELS; i++){
SetWidgetState(solButt[i], FALSE);
}
if(state == HAVE_SOLUTION){
for(i = 0; i < NoOfMeshes; i++)
SetWidgetState(solButt[i], TRUE);
}
}
/*
*
* +-------------------------------------------+
* | Allocate memory for the mesh and solution |
* +-------------------------------------------+
*
*/
int allocateMemory(int imesh, int elnum, int nnode)
{
int ierr;
ierr = 0;
/* allocate new mwmory */
mySol[imesh].nodpek = malloc(NOBF*elnum * sizeof(int));
if(mySol[imesh].nodpek == NULL){
errorHandler("Unable to allocate memory");
ierr = 1;
}
mySol[imesh].coord = malloc(2*nnode * sizeof(float));
if(mySol[imesh].coord == NULL){
errorHandler("Unable to allocate memory");
ierr = 1;
}
mySol[imesh].u = malloc(nnode * sizeof(float));
if(mySol[imesh].u == NULL){
errorHandler("Unable to allocate memory");
ierr = 1;
}
mySol[imesh].vecX = malloc(nnode * sizeof(float));
if(mySol[imesh].vecX == NULL){
errorHandler("Unable to allocate memory");
ierr = 1;
}
mySol[imesh].vecY = malloc(nnode * sizeof(float));
if(mySol[imesh].vecY == NULL){
errorHandler("Unable to allocate memory");
ierr = 1;
}
return ierr;
}
/*
*
* +---------------+
* | Plot the mesh |
* +---------------+
*
*/
void plotMesh()
{
float x[4], y[4];
int i, el, nod, ncells;
float cfl, xmin, ymin, umax, vmax;
getDrawingScale();
cfl = -1000.;
mySetColor(GREY);
for( el = 0; el < mySol[currentSolution].elnum; el++){
xmin = -1000.;
ymin = -1000.;
umax = -1000.;
vmax = -1000.;
for( i = 0; i < NOBF; i++){
nod = mySol[currentSolution].nodpek[el][i];
x[i] = mySol[currentSolution].coord[ nod ][0];
y[i] = mySol[currentSolution].coord[ nod ][1];
}
x[NOBF] = x[0];
y[NOBF] = y[0];
/*
for( i = 0; i < NOBF; i++){
nod = mySol[currentSolution].nodpek[el][i];
xmin = max(xmin, (float) fabs( (double) (x[i]-x[i+1])));
umax = max(umax, (float) fabs( (double) mySol[currentSolution].vecX[ nod ]));
ymin = max(ymin, (float) fabs( (double) (y[i]-y[i+1])));
vmax = max(vmax, (float) fabs( (double) mySol[currentSolution].vecY[ nod ]));
x[i] += mySol[currentSolution].vecX[ nod ]*arrwscl;
y[i] += mySol[currentSolution].vecY[ nod ]*arrwscl;
}
cfl = max(cfl, arrwscl * umax / xmin);
cfl = max(cfl, arrwscl * vmax / ymin);
x[NOBF] = x[0];
y[NOBF] = y[0];
*/
line(x, y, NOBF+1);
}
mySetColor(WHITE);
drawTitle("Mesh");
/* printf("dt = %f, cfl = %f\n", arrwscl, cfl); */
}
void plotElementNumbers()
{
char number[3];
int el, i, node;
float xx, yy;
int x, y;
for( el = 0; el < mySol[currentSolution].elnum; el++){
xx = 0.;
yy = 0.;
for(i = 0; i < NOBF; i++){
node = mySol[currentSolution].nodpek[el][i];
xx += mySol[currentSolution].coord[node][0];
yy += mySol[currentSolution].coord[node][1];
}
xx = xx / NOBF;
yy = yy / NOBF;
sprintf(number, "%i", el + 1);
phys2pix(&x, &y, xx, yy);
myDrawText(number, x, y);
}
}
void plotNodeNumbers()
{
char number[3];
int el, i, node;
float xx, yy;
int x, y;
for( el = 0; el < mySol[currentSolution].elnum; el++){
for(i = 0; i < NOBF; i++){
node = mySol[currentSolution].nodpek[el][i];
xx = mySol[currentSolution].coord[node][0];
yy = mySol[currentSolution].coord[node][1];
sprintf(number, "%i", node + 1);
phys2pix(&x, &y, xx, yy);
myDrawText(number, x, y);
}
}
}
void plotPhysBoundary()
{
int i,j;
int n, nb;
float xx[2], yy[2];
GEO_NODE *g;
BND_NODE *s;
FOR_ALL_BOUNDARIES(g){
FOR_ALL_SEGMENTS(g, s){
xx[0] = ((s->seg)->pts1)->physicalX;
yy[0] = ((s->seg)->pts1)->physicalY;
xx[1] = ((s->seg)->pts2)->physicalX;
yy[1] = ((s->seg)->pts2)->physicalY;
line(xx, yy, 2);
}
}
}
void plotVectors(int this)
{
int i, ptr, ncells;
float xx, xy, xu, xv;
float umax, uu;
ncells = setUpMycolorMap(ncolor);
getDrawingScale();
ptr = mySol[currentSolution].ptr[this];
umax = 0.0;
for( i = 0; i < mySol[currentSolution].nnode; i++){
xu = mySol[currentSolution].f[i + ptr];
xv = mySol[currentSolution].f[i + ptr + mySol[currentSolution].nnode];
umax = max(umax, xu*xu+xv*xv);
}
umax = (float)sqrt((double) umax);
for( i = 0; i < mySol[currentSolution].nnode; i++){
xu = mySol[currentSolution].f[i + ptr];
xv = mySol[currentSolution].f[i + ptr + mySol[currentSolution].nnode];
uu = (float) sqrt((double)xu*xu+xv*xv) ;
if(uu <= umax*maxscalarX && uu >= umax*minscalarX){
xx = mySol[currentSolution].coord[i][0];
xy = mySol[currentSolution].coord[i][1];
xu = xx + maxarrwscl*arrwscl*xu;
xv = xy + maxarrwscl*arrwscl*xv;
mySetMyColor( (uu-umax*minscalarX)/
(umax*maxscalarX-umax*minscalarX));
arrow(xx, xy, xu, xv);
}
}
mySetColor(WHITE);
plotPhysBoundary();
drawTitle(mySol[currentSolution].name[this]);
#ifdef NAVIERSTOKES
dolingraph();
#endif
}
void plotIsolines(int this)
{
double maxscalar, minscalar, scl;
double xx[NOBF+1], yy[NOBF+1], sc[NOBF+1];
float xp[10], yp[10];
double level;
int el, i, ilev, ii, node, ptr, ntext, ncells;
char text[100];
float xf, yf, dx, dy;
int last_level, last_levely[40];
int ix, iy,x1, y1, nx, ny;
ntext = 0;
linitem = this;
ncells = setUpMycolorMap(ncolor);
getDrawingScale();
ptr = mySol[currentSolution].ptr[this];
maxscalar = -100000.;
minscalar = 100000.;
for(i = 0; i < mySol[currentSolution].nnode; i++){
maxscalar = max(maxscalar, (double) mySol[currentSolution].f[i+ptr]);
minscalar = min(minscalar, (double) mySol[currentSolution].f[i+ptr]);
}
scl = maxscalar - minscalar;
maxscalar = minscalar +(double) maxscalarX * scl;
minscalar = minscalar +(double) minscalarX * scl;
scl = maxscalar - minscalar;
sprintf(text, "Min scalar: %f , Max scalar: %f", minscalar, maxscalar);
infoHandler(text);
for( el = 0; el < mySol[currentSolution].elnum; el++){
for(i = 0; i < NOBF; i++){
node = mySol[currentSolution].nodpek[el][i];
xx[i] = (double)mySol[currentSolution].coord[node][0];
yy[i] = (double)mySol[currentSolution].coord[node][1];
sc[i] = (double)mySol[currentSolution].f[node + ptr];
}
xx[NOBF] = xx[0];
yy[NOBF] = yy[0];
sc[NOBF] = sc[0];
for(ilev = 0; ilev < nlev; ilev++){
level = -scl + (double) ilev*2*scl/nlev;
level = 1.*ilev*scl / nlev + minscalar;
level = minscalar + scl * (double) ilev/nlev;
ii = 0;
for(i = 0; i < NOBF; i++){
if( (sc[i]-level)*(sc[i+1]-level) <= 0.0){
if(sc[i+1]-sc[i] != 0.){
xp[ii] = (float) xx[i] +(level-sc[i])/(sc[i+1]-sc[i])*(xx[i+1]-xx[i]);
yp[ii] = (float) yy[i] +(level-sc[i])/(sc[i+1]-sc[i])*(yy[i+1]-yy[i]);
ii += 1;
}
}
}
if(ii > 0){
mySetMyColor( (float) ilev / (float) nlev);
line(xp, yp, ii);
}
}
} /* next element */
mySetColor(WHITE);
plotPhysBoundary();
drawTitle(mySol[currentSolution].name[this]);
}
void plotIsolineLabels(int this)
{
double maxscalar, minscalar, scl;
double xx[NOBF+1], yy[NOBF+1], sc[NOBF+1];
float xp[10], yp[10];
double level;
int el, i, ilev, ii, node, ptr, ntext, ncells;
char text[100];
float xf, yf, dx, dy;
int last_level, last_levely[40];
int ix, iy,x1, y1, nx, ny;
ntext = 0;
linitem = this;
ncells = setUpMycolorMap(ncolor);
getDrawingScale();
ptr = mySol[currentSolution].ptr[this];
maxscalar = -100000.;
minscalar = 100000.;
for(i = 0; i < mySol[currentSolution].nnode; i++){
maxscalar = max(maxscalar, (double) mySol[currentSolution].f[i+ptr]);
minscalar = min(minscalar, (double) mySol[currentSolution].f[i+ptr]);
}
scl = maxscalar - minscalar;
maxscalar = minscalar +(double) maxscalarX * scl;
minscalar = minscalar +(double) minscalarX * scl;
scl = maxscalar - minscalar;
/* Put up labels on the isolines */
nx = 15;
ny = 15;
dx = (float)(xmax - xmin ) / nx;
dy = (float)(ymax - ymin ) / ny;
ntext = 0;
last_level = 1233;
for(ix = 0; ix <=nx; ix++)
last_levely[ix] = 121222;
for( iy = 0; iy <= ny; iy=iy+2){
for( ix = 0; ix <= nx; ix=ix+2){
xf = xmin + (float)ix*dx;
yf = ymin + (float)iy*dy;
phys2pix(&x1, &y1, xf, yf);
if(clickedInElement(x1, y1, &el) == TRUE){
for(i = 0; i < NOBF; i++){
node = mySol[currentSolution].nodpek[el][i];
xx[i] = (double)mySol[currentSolution].coord[node][0];
yy[i] = (double)mySol[currentSolution].coord[node][1];
sc[i] = (double)mySol[currentSolution].f[node + ptr];
}
xx[NOBF] = xx[0];
yy[NOBF] = yy[0];
sc[NOBF] = sc[0];
for(ilev = 0; ilev <= nlev; ilev=ilev+1){
level = minscalar + scl * (double) ((ilev+ntext)%nlev)/nlev;
ii = 0;
for(i = 0; i < NOBF; i++){
if( (sc[i]-level)*(sc[i+1]-level) <= 0.0){
if(sc[i+1]-sc[i] != 0.){
xp[ii] = (float) xx[i] +(level-sc[i])/(sc[i+1]-sc[i])*(xx[i+1]-xx[i]);
yp[ii] = (float) yy[i] +(level-sc[i])/(sc[i+1]-sc[i])*(yy[i+1]-yy[i]);
ii += 1;
}
}
}
if(ii > 0){
mySetMyColor( (float)((ilev+ntext)%nlev/(float)nlev));
if((ilev+ntext)%nlev != last_level){
if((ilev+ntext)%nlev != last_levely[ix]){
phys2pix(&x1, &y1, xp[0], yp[0]);
sprintf(text, "%.3g", level);
myDrawText(text, x1+2, y1-2);
myDrawCross(x1, y1, 4);
last_level = (ilev+ntext)%nlev;
last_levely[ix] = (ilev+ntext)%nlev;
}
}
ntext++;
if(ntext > nlev)
ntext = 0;
break;
}
}
}
}
}
mySetColor(WHITE);
}
void plotIsoColorLines(int this)
{
float maxscalar, minscalar, scl;
int xx[NOBF+1], yy[NOBF+1], sc[NOBF+1];
int ix, iy;
float x, y;
int level;
int lmin, lmax;
int start, end;
int el, i, ilev, ii, node, ncells, ptr;
XPoint *pp, *pts;
char text[100];
pts = (XPoint *) malloc( 10 * sizeof(XPoint) );
pp = pts;
ncells = setUpMycolorMap(ncolor);
getDrawingScale();
ptr = mySol[currentSolution].ptr[this];
maxscalar = -100000.;
minscalar = 100000.;
for(i = 0; i < mySol[currentSolution].nnode; i++){
maxscalar = max(maxscalar, mySol[currentSolution].f[i+ptr]);
minscalar = min(minscalar, mySol[currentSolution].f[i+ptr]);
}
scl = maxscalar - minscalar;
sprintf(text, "Min scalar: %f , Max scalar: %f", minscalar, maxscalar);
infoHandler(text);
for( el = 0; el < mySol[currentSolution].elnum; el++){
lmax = -10000;
lmin = 10000;
for(i = 0; i < NOBF; i++){
node = mySol[currentSolution].nodpek[el][i];
sc[i] =(int) floor((double)(ncells-1)* (mySol[currentSolution].f[ptr+node]-minscalar)
/ scl);
x = mySol[currentSolution].coord[node][0];
y = mySol[currentSolution].coord[node][1];
phys2pix(&ix, &iy, x, y);
xx[i] = ix;
yy[i] = iy;
lmax = max(lmax, sc[i]);
lmin = min(lmin, sc[i]);
}
xx[NOBF] = xx[0];
yy[NOBF] = yy[0];
sc[NOBF] = sc[0];
start = max(0, lmin);
end = min(ncells - 1, lmax);
for(ilev = start; ilev < end + 1; ilev++){
level = ilev;
ii = 0;
for(i = 0; i < NOBF; i++){
if(sc[i] >= level){
pp[ii].x = xx[i];
pp[ii].y = yy[i];
ii += 1;
}
if( (sc[i]-level)*(sc[i+1]-level) <= 0){
if(sc[i+1]-sc[i] != 0){
pp[ii].x = xx[i] +(level-sc[i])*(xx[i+1]-xx[i])/(sc[i+1]-sc[i]);
pp[ii].y = yy[i] +(level-sc[i])*(yy[i+1]-yy[i])/(sc[i+1]-sc[i]);
ii++;
}
else{
pp[ii].x = xx[i];
pp[ii].y = yy[i];
ii++;
}
}
}
mySetMyColor((float)ilev/ (float)ncells);
if(ii > 0){
pp[ii].x = pp[0].x;
pp[ii].y = pp[0].y;
ii += 1;
myDrawFilledPolygon(pp, ii);
}
}
} /* next element */
plotPhysBoundary();
free(pts);
}
void
solfcn(Widget w, void *data)
{
int hh = 10, ww = 10;
int i;
for( i = 0; i < NoOfMeshes; i++){
if( GetToggleState( solButt[i])){
currentSolution = i;
redisplay(w, hh, ww, data);
}
}
}
void putUpLevelButtons()
{
int i;
char text[2];
currentSolution = 0;
label = MakeLabel(" Refinment Level:");
sprintf(text,"%i",0);
solButt[0] = MakeToggle(text, TRUE, NULL, solfcn, NULL);
for(i = 1; i < MAXLEVELS; i++){
sprintf(text,"%i",i);
solButt[i] = MakeToggle(text, FALSE, solButt[0], solfcn, NULL);
}
SetWidgetPos(label, PLACE_RIGHT, plotMenu.name, NO_CARE, NULL);
SetWidgetPos(solButt[0], PLACE_RIGHT, label, NO_CARE, NULL);
for(i = 1; i < MAXLEVELS; i++)
SetWidgetPos(solButt[i], PLACE_RIGHT, solButt[i-1], NO_CARE, NULL);
/* ShowDisplay(); */
}
void linMouse(Widget w, int x, int y, void *data)
{
int width, height;
GetDrawAreaSize(&width, &height);
crossHair(x, y);
if(nlin == 1){
SetDrawMode(SANE_XOR);
DrawLine(xplin2, yplin2, xplin1, yplin1);
DrawLine(xplin2, yplin2, x, y);
xplin1 = x;
yplin1 = y;
SetDrawMode(GXcopy);
}
/*
if(nlin == 0){
nold = 0;
pix2phys(x, 0, &xlin1, &ylin1);
pix2phys(x+1, height, &xlin2, &ylin2);
SetDrawMode(SANE_XOR);
if(have_colors == TRUE)
SetColor(YELLOW);
if(have_colors == TRUE)
SetColor(WHITE);
SetDrawMode(GXcopy);
nlin = 1;
}
else{
SetDrawMode(SANE_XOR);
pix2phys(x, 0, &xlin1, &ylin1);
pix2phys(x+1, height, &xlin2, &ylin2);
if(have_colors == TRUE)
SetColor(YELLOW);
if(have_colors == TRUE)
SetColor(WHITE);
SetDrawMode(GXcopy);
nlin = 1;
}
*/
}
void linDown(Widget w, int button, int x, int y, void *data)
{
if(nlin == 0){
pix2phys(x, y, &xlin2, &ylin2);
xplin2 = x;
yplin2 = y;
xplin1 = x;
yplin1 = y;
nlin = 1;
}
else{
if(x == xplin2)
x = x + 1;
pix2phys(x, y, &xlin1, &ylin1);
deleteCrossHair();
SetButtonDownCB(drawWindow, NULL);
SetMouseMotionCB(drawWindow, NULL);
lingraph();
}
}
void findFunction(float *xx, float *yy, int *n)
{
double x1, y1; /* Starting point of line */
double x2, y2; /* End point of line */
double k0, y0, k, m;
double xp[NOBF+1], yp[NOBF+1];
double ux[NOBF+1];
double xcross, ycross, r, f, li, xi;
int indx, el, i, node;
int ptr;
ptr = mySol[currentSolution].ptr[linitem];
x1 = (double) xlin1;
x2 = (double) xlin2;
y1 = (double) ylin1;
y2 = (double) ylin2;
k0 = (y2 - y1) / (x2 - x1);
y0 = y1 - k0 * x1 ;
indx = 0;
for( el = 0; el < mySol[currentSolution].elnum; el++){
for(i = 0; i < NOBF; i++){
node = (double) mySol[currentSolution].nodpek[el][i];
xp[i] = (double) mySol[currentSolution].coord[node][0];
yp[i] = (double) mySol[currentSolution].coord[node][1];
ux[i] = (double) mySol[currentSolution].f[node + ptr];
}
xp[NOBF] = xp[0];
yp[NOBF] = yp[0];
ux[NOBF] = ux[0];
for(i = 0; i < NOBF; i++){
if(xp[i] != xp[i+1]){
k = (yp[i] - yp[i+1]) / (xp[i] - xp[i+1]);
m = yp[i] - k * xp[i];
if( k != k0){
xcross = (y0 - m) / (k - k0);
ycross = k*xcross + m;
if( (xcross - xp[i])*(xcross - xp[i+1]) <= 0.0){
if( (ycross - yp[i])*(ycross - yp[i+1]) <= 0.0){
if( (xcross-x1)*(xcross-x2) <=0.0 && (ycross-y1)*(ycross-y2) <=0.0){
r = (xcross - x1)*(xcross - x1) + (ycross - y1)*(ycross - y1);
li = (xp[i+1] - xp[i])*(xp[i+1] - xp[i]) +
(yp[i+1] - yp[i])*(yp[i+1] - yp[i]);
xi = (xcross - xp[i])*(xcross - xp[i]) +
(ycross - yp[i])*(ycross - yp[i]);
xi = sqrt(xi/li);
f = ux[i] + xi * (ux[i+1] - ux[i]);
xx[indx] = (float) sqrt(r);
yy[indx] = (float) f;
indx++;
}
}
}
}
}
}
}
*n = indx;
}
void dolingraph()
{
nlin = 0;
SetButtonDownCB(drawWindow, linDown);
SetMouseMotionCB(drawWindow, linMouse);
}
void lingraph()
{
float x[1000], f[1000];
float tst, tmp;
/* float xmin, xmax, ymin, ymax;*/
int n;
int width, height;
int i, j;
float xx[3], yy[3];
findFunction(x, f, &n);
/* printf("Plotting linear graph with %i points\n", n); */
xmax = -1000.;
ymax = -1000.;
xmin = 1000.;
ymin = 1000.;
for(i = 0; i < n; i++){
ymax = max(ymax, f[i]);
ymin = min(ymin, f[i]);
xmax = max(xmax, x[i]);
xmin = min(xmin, x[i]);
}
/*
printf("xmin, xmax : %f, %f\n", xmin, xmax);
printf("ymin, ymax : %f, %f\n", ymin, ymax);
*/
for(i = 0; i < n; i++){
tst = x[i];
for(j = i + 1; j < n; j++){
if(x[j] < tst){
tst = x[j];
tmp = x[i];
x[i] = x[j];
x[j] = tmp;
tmp = f[i];
f[i] = f[j];
f[j] = tmp;
}
}
}
ClearDrawArea();
GetDrawAreaSize(&width, &height);
getDrawingScale();
for(i = 0; i < n - 1; i++){
xx[0] = x[i];
yy[0] = f[i];
xx[1] = x[i+1];
yy[1] = f[i+1];
line(xx, yy, 2);
}
xx[0] = xmin;
yy[0] = ymin;
xx[1] = xmax;
yy[1] = ymin;
line(xx, yy, 2);
xx[0] = xmax;
yy[0] = ymin;
xx[1] = xmax;
yy[1] = ymax;
line(xx, yy, 2);
xx[0] = xmax;
yy[0] = ymax;
xx[1] = xmin;
yy[1] = ymax;
line(xx, yy, 2);
xx[0] = xmin;
yy[0] = ymax;
xx[1] = xmin;
yy[1] = ymin;
line(xx, yy, 2);
/*
DrawLine((int)(0.1*width),(int)( 0.1*height),(int)( 0.9*width),(int)( 0.1*height));
DrawLine((int)(0.1*width),(int)( 0.1*height),(int)( 0.1*width),(int)( 0.9*height));
DrawLine((int)(0.9*width),(int)( 0.1*height),(int)( 0.9*width),(int)( 0.9*height));
DrawLine((int)(0.1*width),(int)( 0.9*height),(int)( 0.9*width),(int)( 0.9*height));
if(xmin*xmax <= 0.0){
xp = (int) (width *(0.8 * (0.0 - xmin)/(xmax - xmin) + 0.1));
DrawLine(xp,(int)( 0.1*height), xp, (int)( 0.9*height));
}
*/
}
syntax highlighted by Code2HTML, v. 0.9.1