Attached is cyl1.inp I happened to have a 3D duct with
a cylinder. Change to suit you needs.
This is a boundary file, not a refinement file.
No DOF included. It has both Dirichlet and Neumann
boundary conditions included.
I also attached a plot UCD source code to plot UCD .inp
files. Oh, again for 3D with roll, pitch and heading rotation.
It does need OpenGL and glut to compile and link.
plot_ucd.c Just a plain C program.
Jon [email protected] www.csee.umbc.edu/~squire
On Thu, 7 Oct 2010, Keith Brauss wrote:
Hello all,
I am interested in solving a incompressible fluid flow
(in a pipe or rectangular duct domain) around a cylinder
contained within the pipe or duct domain.
I was wondering if anybody could point me in the right
direction concerning making the mesh and implementing the
boundary conditions on the cylinder. Thanks.
Dan Brauss
# cyl1.inp very coarse boundary, UCD format
# cells are boundary faces
# cylinder radius=1.0 in duct 10 by 10 input and output, 20 long
# X coordinates of duct -5.0 to +5.0
# Y coordinates of duct -5.0 to +5.0
# Z length coordinates of duct -10.0 to 10.0
# cylinder axis -5.0, 0.0, 0.0 to +5.0, 0.0, 0.0
# Dirichlet pressure at nodes
# Neumann velocity on faces, positive outward normal
24 14 1 1 0
1 -5.0 +5.0 -10.0
2 +5.0 +5.0 -10.0
3 +5.0 -5.0 -10.0
4 -5.0 -5.0 -10.0
5 -5.0 +5.0 +10.0
6 +5.0 +5.0 +10.0
7 +5.0 -5.0 +10.0
8 -5.0 -5.0 +10.0
9 -5.0 +1.0000 0.0000
10 -5.0 +0.7071 +0.7071
11 -5.0 0.0000 +1.0000
12 -5.0 -0.7071 +0.7071
13 -5.0 -1.0000 0.0000
14 -5.0 -0.7071 -0.7071
15 -5.0 0.0000 -1.0000
16 -5.0 +0.7071 -0.7071
17 +5.0 +1.0000 0.0000
18 +5.0 +0.7071 +0.7071
19 +5.0 0.0000 +1.0000
20 +5.0 -0.7071 +0.7071
21 +5.0 -1.0000 0.0000
22 +5.0 -0.7071 -0.7071
23 +5.0 0.0000 -1.0000
24 +5.0 +0.7071 -0.7071
1 0 quad 1 2 3 4
2 0 quad 1 5 6 2
3 0 quad 2 6 7 3
4 0 quad 4 8 7 3
5 0 quad 1 5 8 4
6 0 quad 5 6 7 8
7 0 quad 9 10 18 17
8 0 quad 10 11 19 18
9 0 quad 11 12 20 19
10 0 quad 12 13 21 20
11 0 quad 13 14 22 21
12 0 quad 14 15 23 22
13 0 quad 15 16 24 23
14 0 quad 16 9 17 24
1 1
Dirichlet Boundary, pressure PSI
1 14.7
2 14.7
3 14.7
4 14.7
5 14.7
6 14.7
7 14.7
8 14.7
9 14.7
10 14.7
11 14.7
12 14.7
13 14.7
14 14.7
15 14.7
16 14.7
17 14.7
18 14.7
19 14.7
20 14.7
21 14.7
22 14.7
23 14.7
24 14.7
1 1
Neumann Boundary, velocity outward normal MPH
1 -89.0
2 0.0
3 0.0
4 0.0
5 0.0
6 +89.0
7 0.0
8 0.0
9 0.0
10 0.0
11 0.0
12 0.0
13 0.0
14 0.0
/* plot_ucd.c Rotating, wire frame or with lighting */
/* UCD format 3D boundary file *.inp */
/* link with opengl(mesa) and glut */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <GL/glut.h>
static FILE * infile; /* input file name on command line */
static FILE * outfile; /* optional, second, output file on command line */
static int num_pts;
static int num_pys;
typedef struct {GLfloat x; GLfloat y; GLfloat z;} dpts;
static dpts * data_points;
static int * data_verts;
static int debug = 0;
static GLfloat size = 1.0; /* auto scale */
static GLfloat roll = 0.0; /* rotate about Z axis 'r' */
static GLfloat pitch = 0.0; /* rotate about X axis 'p' */
static GLfloat heading = 0.0; /* rotate about Y axis 'h' */
static GLfloat pos[3] = {0.0, 0.0, 0.0};
static int background = 1; /* background color 'b' */
static int wire = 1; /* draw wireframe 'w' */
static int vert = 0; /* draw vertices 'v' */
static int colorz = 0; /* color based on z 'c' */
static GLfloat zavg = 0.0;
static GLfloat zmin = 0.0;
static GLfloat zmax = 0.0;
static GLfloat amb[4] = { 0.329412, 0.223529, 0.027451, 1.0}; /* brass */
static GLfloat dif[4] = { 0.780392, 0.568627, 0.113725, 1.0}; /* material */
static GLfloat spe[4] = { 0.992157, 0.941176, 0.807843, 1.0};
static GLfloat shiney = 0.21794872;
/* light source */
static GLfloat ambient[] = {0.0, 0.0, 0.0, 1.0}; /* no ambient */
static GLfloat diffuse[] = {0.4, 0.4, 0.4, 1.0}; /* mid diffuse */
static GLfloat specular[] = {0.4, 0.4, 0.4, 1.0}; /* mid specular */
static GLfloat position[] = {0.0, 0.0, 3.0, 0.0};
static GLdouble m[16]; /* 4 by 4 transformation matrix */
static int win_w, win_h;
static void draw_solid_dat(void)
{
int i, j, k, pts, pt;
GLfloat v[3][3];
int kk[3];
GLfloat ax, ay, az, bx, by, bz, nx, ny, nz, s;
glLightfv(GL_LIGHT0, GL_AMBIENT, ambient);
glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse);
glLightfv(GL_LIGHT0, GL_POSITION, position);
glFrontFace(GL_CCW);
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
glEnable(GL_NORMALIZE);
glEnable(GL_AUTO_NORMAL);
glEnable(GL_DEPTH_TEST);
glLightfv(GL_LIGHT0, GL_POSITION, position);
glPushMatrix();
glLoadIdentity();
glRotatef(pitch, 1.0, 0.0, 0.0);
glRotatef(heading, 0.0, 1.0, 0.0);
glRotatef(roll, 0.0, 0.0, 1.0);
glTranslatef(pos[0], pos[1], pos[2]);
glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, amb);
glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, dif);
glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, spe);
glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, shiney * 64.0);
/* use pys to do surfaces */
k = 0;
for(i=0; i<num_pys; i++)
{
pts = data_verts[k++];
if(debug>2) printf("pts=%d\n", pts);
glBegin(GL_POLYGON);
for(j=0; j<3; j++)
{
pt = data_verts[k++];
if(debug>2) printf("j=%d, pt=%d\n", j, pt);
kk[j] = pt;
v[j][0] = data_points[pt-1].x;
v[j][1] = data_points[pt-1].y;
v[j][2] = data_points[pt-1].z;
}
/* compute normals, assuming planar, normalize */
ax = v[2][0] - v[1][0];
ay = v[2][1] - v[1][1];
az = v[2][2] - v[1][2];
bx = v[1][0] - v[0][0];
by = v[1][1] - v[0][1];
bz = v[1][2] - v[0][2];
nx = ay*bz-az*by;
ny = az*bx-ax*bz;
nz = ax*by-ay*bx;
s = sqrt(nx*nx+ny*ny+nz*nz);
nx = nx / s;
ny = ny / s;
nz = nz / s;
if(debug>1)
{
printf("\nv0x=%f, v0y=%f, v0z=%f\n", v[0][0], v[0][1], v[0][2]);
printf("v1x=%f, v1y=%f, v1z=%f\n", v[1][0], v[1][1], v[1][2]);
printf("v2x=%f, v2y=%f, v2z=%f\n", v[2][0], v[2][1], v[2][2]);
printf("nx=%f, ny=%f, nz=%f\n", nx, ny, nz);
}
for(j=0; j<3; j++)
{
glNormal3f(nx, ny, nz);
glVertex3f(data_points[kk[j]-1].x, data_points[kk[j]-1].y,
data_points[kk[j]-1].z);
}
for(j=3; j<pts; j++)
{
pt = data_verts[k++];
if(debug>2) printf("j=%d, pt=%d\n", j, pt);
glNormal3f(nx, ny, nz);
glVertex3f(data_points[pt-1].x, data_points[pt-1].y,
data_points[pt-1].z);
}
glEnd();
} /* end loop on num_pys */
glPopMatrix();
} /* end draw_solid_dat */
static void draw_wire_dat(void)
{
int i, j, k, pts, pt, pt0;
glDisable(GL_LIGHTING);
glDisable(GL_LIGHT0);
glDisable(GL_NORMALIZE);
glDisable(GL_AUTO_NORMAL);
if(background) glColor3f(0.0, 0.0, 0.0);
if(!background) glColor3f(1.0, 1.0, 1.0);
glPointSize(3.0);
glPushMatrix();
glLoadIdentity();
glRotatef(pitch, 1.0, 0.0, 0.0);
glRotatef(heading, 0.0, 1.0, 0.0);
glRotatef(roll, 0.0, 0.0, 1.0);
glTranslatef(pos[0], pos[1], pos[2]);
glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */
/* use pys to do edges */
k = 0;
for(i=0; i<num_pys; i++)
{
pts = data_verts[k++];
if(debug>2) printf("pts=%d\n", pts);
glBegin(GL_LINE_STRIP);
for(j=0; j<pts; j++)
{
pt = data_verts[k++]; /* packed data */
if(j==0) pt0 = pt; /* close polygon */
if(debug>2) printf("j=%d, pt=%d\n", j, pt);
glVertex3f(data_points[pt-1].x, data_points[pt-1].y,
data_points[pt-1].z);
}
glVertex3f(data_points[pt0-1].x, data_points[pt0-1].y,
data_points[pt0-1].z);
glEnd();
} /* end loop on num_pys */
glPopMatrix();
} /* end draw_wire_dat */
static void draw_vert_dat(void)
{
int i, j, k, pts, pt;
glDisable(GL_LIGHTING);
glDisable(GL_LIGHT0);
glDisable(GL_NORMALIZE);
glDisable(GL_AUTO_NORMAL);
if(background) glColor3f(0.0, 0.0, 0.0);
if(!background) glColor3f(1.0, 1.0, 1.0);
glPointSize(3.0);
glPushMatrix();
glLoadIdentity();
glRotatef(pitch, 1.0, 0.0, 0.0);
glRotatef(heading, 0.0, 1.0, 0.0);
glRotatef(roll, 0.0, 0.0, 1.0);
glTranslatef(pos[0], pos[1], pos[2]);
glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */
/* use pys to do surfaces */
k = 0;
for(i=0; i<num_pys; i++)
{
pts = data_verts[k++];
if(debug>2) printf("pts=%d\n", pts);
glBegin(GL_POINTS);
for(j=0; j<pts; j++)
{
pt = data_verts[k++];
if(debug>2) printf("j=%d, pt=%d\n", j, pt);
glVertex3f(data_points[pt-1].x, data_points[pt-1].y,
data_points[pt-1].z);
}
glEnd();
} /* end loop on num_pys */
glPopMatrix();
} /* end draw_vert_dat */
static writeText(GLfloat x, GLfloat y, GLfloat z, GLfloat size, char * msg)
{ /* uses background and wire frame */
char * p;
GLfloat wht[] = {1.0, 1.0, 1.0, 1.0};
GLfloat blk[] = {0.0, 0.0, 0.0, 0.0};
glPushMatrix();
glLoadIdentity ();
if(background)
{
if(!wire)
{
glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, blk);
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, blk);
}
else
{
glColor3f(0.0, 0.0, 0.0);
}
}
else
{
if(!wire)
{
glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, wht);
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, wht);
}
else
{
glColor3f(1.0, 1.0, 1.0);
}
}
glEnable(GL_LINE_SMOOTH);
glEnable(GL_BLEND);
glTranslatef(x, y, z);
glScalef(size/1000.0, size/1000.0, 0.0);
for(p=msg; *p; p++)
glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
glPopMatrix();
}
void display(void)
{
if(background) glClearColor(1.0, 1.0, 1.0, 1.0);
if(!background) glClearColor(0.0, 0.0, 0.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
if(vert) draw_vert_dat();
else if(wire) draw_wire_dat();
else draw_solid_dat();
writeText(-size*0.95, -size*0.95+size/16.0, 2.0*size-1.0, size/2.5,
"Key r,R roll, p,P pitch, h,H heading, b toggles background");
writeText(-size*0.95, -size*0.95, 2.0*size-1.0, size/2.5,
"Key x,X y,Y z,Z to position, w for wireframe, v for vertices");
glFlush();
glutSwapBuffers();
} /* end display */
static GLdouble dist(GLdouble a1x, GLdouble a1y, GLdouble a1z, /* line p1 */
GLdouble a2x, GLdouble a2y, GLdouble a2z, /* line p2 */
GLdouble a3x, GLdouble a3y, GLdouble a3z) /* point */
{
GLdouble d, ax, ay, az, bx, by, bz, cx, cy, cz;
ax = a2x-a1x; ay = a2y-a1y; az = a2z-a1z;
bx = a1x-a3x; by = a1y-a3y; bz = a1z-a3z;
cx = ay*bz-az*by; cy = ax*bz-az*bx; cz = ax*by-ay*bx;
d = sqrt(cx*cx+cy*cy+cz*cz)/sqrt(ax*ax+ay*ay+az*az);
return d;
}
void mouse(int button, int state, int x, int y)
{
GLint viewport[4];
GLdouble mvmatrix[16], projmatrix[16];
GLint realy; /* OpenGL y coordinate position */
GLdouble wx0, wy0, wz0; /* returned world x, y, z coords at near */
GLdouble wx1, wy1, wz1; /* returned world x, y, z coords at far */
GLdouble d, dmin;
int i, imin;
if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
{
glGetIntegerv (GL_VIEWPORT, viewport);
glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);
/* note viewport[3] is height of window in pixels */
realy = viewport[3] - (GLint) y - 1;
gluUnProject((GLdouble) x, (GLdouble) realy, 0.0, /* near */
m, projmatrix, viewport, &wx0, &wy0, &wz0);
gluUnProject((GLdouble) x, (GLdouble) realy, 1.0, /* far */
m, projmatrix, viewport, &wx1, &wy1, &wz1);
imin = 0;
for(i=0; i<num_pts; i++)
{
d = dist(wx0, wy0, wz0, wx1, wy1, wz1,
data_points[i].x, data_points[i].y, data_points[i].z);
if(i==0) dmin = d;
if(debug) printf("distance, d to point %d is %g\n", i, d);
if(d < dmin) { dmin = d; imin = i; }
}
printf("%g close to data_points[%d]\n", dmin, imin);
}
if(button==GLUT_MIDDLE_BUTTON && state == GLUT_DOWN) ; /* commit */
if(button==GLUT_RIGHT_BUTTON && state == GLUT_DOWN) ; /* drag */
}
void myReshape(int w, int h)
{
win_w = w;
win_h = h;
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
if(w <= h)
glOrtho(-size, size, -size * (GLfloat)h / (GLfloat)w,
size * (GLfloat)h / (GLfloat)w, -5.0*size, 5.0*size);
else
glOrtho(-size * (GLfloat)w / (GLfloat)h,
size * (GLfloat)w / (GLfloat)h, -size, size, -5.0*size, 5.0*size);
glMatrixMode(GL_MODELVIEW);
}
/* x,X,y,Y,z and Z keys for moving position
r,R,p,P,h and H keys for rool, pitch and heading
w toggle wire frame, v toggle vertices
b toggle background, c toggle color in z */
void keyboard (unsigned char key, int x, int y)
{
if(key == 'x') pos[0] -= 0.125;
if(key == 'X') pos[0] += 0.125;
if(key == 'y') pos[1] -= 0.125;
if(key == 'Y') pos[1] += 0.125;
if(key == 'z') pos[2] -= 0.125;
if(key == 'Z') pos[2] += 0.125;
if(key == 'r') roll -= 5.0;
if(key == 'R') roll += 5.0;
if(key == 'p') pitch -= 5.0;
if(key == 'P') pitch += 5.0;
if(key == 'h') heading -= 5.0;
if(key == 'H') heading += 5.0;
if(key == 'b' || key == 'B') background = 1-background;
if(key == 'w' || key == 'W') wire = 1-wire;
if(key == 'v' || key == 'V') vert = 1-vert;
if(key == 'c' || key == 'C') colorz = 1-colorz;
glutPostRedisplay();
} /* end keyboard */
#undef abs
#define abs(x) (((x)<0.0)?(-(x)):(x))
static void ucd_read(int argc, char *argv[])
{
int i, j, k, pts, pt;
int kk, len;
char line[255];
char *stat = NULL;
float x,y,z;
int icell, matl;
char typ[8];
/* cell type name, integer number, number of vertices */
int cvert[8];
char *ctyps[9] = {"none", "pt", "line", "tri", "quad",
"tet", "pyr", "prism", "hex"};
int ctypi[9] = { 0, 1, 2, 3, 4,
5, 6, 7, 8};
int ctypn[9] = { 0, 1, 2, 3, 4,
4, 5, 6, 8};
/* get file name of data to render */
if(argc<2)
{
printf("supply a file name of a UCD .inp file\n");
exit(0);
}
printf("open %s for reading\n", argv[1]);
infile = fopen(argv[1], "r");
if(infile == NULL)
{
printf("could not open %s for reading\n", argv[1]);
exit(0);
}
/* skip comments */
while(1) /* read comments from UCD file */
{
stat = fgets(line, 254, infile);
if(stat==NULL)
{
printf("plot_ucd, premature EOF in %s\n", argv[1]);
exit(1);
}
printf("input %s",line); /* optional */
if(line[0]!='#') break;
} /* end while to get past comments */
sscanf(line, "%d %d", &num_pts, &num_pys);
printf("nodes=%d, cells=%d \n", num_pts, num_pys);
data_points = (dpts *)malloc(sizeof(dpts)*(num_pts+2));
data_verts = (int *)malloc(sizeof(int)*8*(num_pys+2)); /* only room pts<8 */
if(data_points == NULL || data_verts == NULL)
{
printf("can not allocate enough space for data\n");
exit(0);
}
for(i=0; i<num_pts; i++)
{
fscanf(infile, "%d %f %f %f\n", &k, &data_points[i].x,
&data_points[i].y, &data_points[i].z);
if(debug) printf("k=%d, x=%f, y=%f, z=%f\n", k, data_points[i].x,
data_points[i].y, data_points[i].z);
if(abs(data_points[i].x)>size) size = abs(data_points[i].x);
if(abs(data_points[i].y)>size) size = abs(data_points[i].y);
if(abs(data_points[i].z)>size) size = abs(data_points[i].z);
}
printf("scaling to size=%g\n", size);
size = size * 1.2; /* for spinning */
position[2] = 1.2*size;
/* use pys to read in cells = boundary surfaces */
kk = 0;
for(i=0; i<num_pys; i++)
{
stat = fgets(line, 254, infile);
if(debug) printf("%s", line);
sscanf(line,"%d %d %s", &icell, &matl, typ);
/* ignored, could test matl==0 for boundary */
k = 0;
for(j=1; j<9; j++)
{
if(strcmp(ctyps[j],typ)==0)
{
k = ctypn[j];
break;
}
}
switch(k)
{
case 1:
data_verts[kk++] = 1;
sscanf(line,"%d %d %s %d", &icell, &matl, typ, &data_verts[kk++]);
break;
case 2:
data_verts[kk++] = 2;
sscanf(line,"%d %d %s %d %d", &icell, &matl, typ, &data_verts[kk++],
&data_verts[kk++]);
break;
case 3:
data_verts[kk++] = 3;
sscanf(line,"%d %d %s %d %d %d", &icell, &matl, typ, &data_verts[kk++],
&data_verts[kk++], &data_verts[kk++]);
break;
case 4:
data_verts[kk++] = 4;
sscanf(line,"%d %d %s %d %d %d %d", &icell, &matl, typ, &data_verts[kk++],
&data_verts[kk++], &data_verts[kk++], &data_verts[kk++]);
break;
default:
printf("bad cell type %s\n", line);
} /* end switch */
} /* end cells */
/* just ignore extra stuff in file */
fclose(infile);
} /* end ucdread */
static void ucd_write(int argc, char *argv[])
{
/* get file name to write */
if(argc<3)
{
printf("no output of UCD .inp file\n");
return;
}
printf("open %s for writing\n", argv[2]);
outfile = fopen(argv[2], "w");
if(outfile == NULL)
{
printf("could not open %s for writing\n", argv[2]);
return;
}
/* fprintf(outfile, "%d %d \n", num_pts, num_pys); */
printf("not implemented yet.\n");
/* fflush(outfile); */
/* fclose(outfile); */
} /* end ucd_write */
int main(int argc, char *argv[])
{
glutInit(&argc, argv);
/* need both double buffering and z buffer */
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
glutInitWindowSize(600, 600);
glutCreateWindow(argv[0]);
glutReshapeFunc(myReshape);
glutDisplayFunc(display);
glutMouseFunc(mouse); /* "mouse" setup */
glutKeyboardFunc(keyboard); /* "keyboard" setup */
glEnable(GL_DEPTH_TEST); /* Enable hidden--surface--removal */
ucd_read(argc, argv); /* read a UCD .inp file */
/* ucd_write(argc, argv);*/ /* write a .inp file when second file name */
glutMainLoop();
return 0;
} /* end plot_ucd.c */
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
