1 /*
2 Copyright (c) 1994 - 2010, Lawrence Livermore National Security, LLC.
3 LLNL-CODE-425250.
4 All rights reserved.
5
6 This file is part of Silo. For details, see silo.llnl.gov.
7
8 Redistribution and use in source and binary forms, with or without
9 modification, are permitted provided that the following conditions
10 are met:
11
12 * Redistributions of source code must retain the above copyright
13 notice, this list of conditions and the disclaimer below.
14 * Redistributions in binary form must reproduce the above copyright
15 notice, this list of conditions and the disclaimer (as noted
16 below) in the documentation and/or other materials provided with
17 the distribution.
18 * Neither the name of the LLNS/LLNL nor the names of its
19 contributors may be used to endorse or promote products derived
20 from this software without specific prior written permission.
21
22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
25 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE
26 LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
27 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
29 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
30 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
31 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
32 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34
35 This work was produced at Lawrence Livermore National Laboratory under
36 Contract No. DE-AC52-07NA27344 with the DOE.
37
38 Neither the United States Government nor Lawrence Livermore National
39 Security, LLC nor any of their employees, makes any warranty, express
40 or implied, or assumes any liability or responsibility for the
41 accuracy, completeness, or usefulness of any information, apparatus,
42 product, or process disclosed, or represents that its use would not
43 infringe privately-owned rights.
44
45 Any reference herein to any specific commercial products, process, or
46 services by trade name, trademark, manufacturer or otherwise does not
47 necessarily constitute or imply its endorsement, recommendation, or
48 favoring by the United States Government or Lawrence Livermore
49 National Security, LLC. The views and opinions of authors expressed
50 herein do not necessarily state or reflect those of the United States
51 Government or Lawrence Livermore National Security, LLC, and shall not
52 be used for advertising or product endorsement purposes.
53 */
54
55 /*
56 * Robb Matzke, Fri Dec 9 13:05:38 EST 1994
57 * Changed for device independence. This is a support file and is not
58 * strictly part of SILO or the SILO-Taurus device driver. However,
59 * we are including the `silo_taurus_private.h' header file to get the
60 * memory management macros, `taurus.h' header file, and function
61 * prototype for `db_taur_extface'.
62 */
63
64 #define SILO_NO_CALLBACKS
65 #include "config.h"
66 #include "silo_taurus_private.h"
67
68 #if HAVE_FCNTL_H
69 #include <fcntl.h> /*open */
70 #endif
71 #if HAVE_SYS_FCNTL_H
72 #include <sys/fcntl.h> /*open */
73 #endif
74 #include <math.h> /*sqrt, acos */
75 #if HAVE_SYS_STAT_H
76 #include <sys/stat.h> /*stat */
77 #endif
78 #include <ctype.h> /*isspace */
79
80 #define MAXBUF 100000
81
82 #ifndef FALSE
83 #define FALSE 0
84 #endif
85 #ifndef TRUE
86 #define TRUE 1
87 #endif
88
89 /*
90 * ftpi = (4/3)*pi, ttpi = (2/3)*pi
91 */
92 #define ftpi 4.188790205
93 #define ttpi 2.094395102
94
95 /*
96 * The list of taurus variables.
97 */
98 var_list_s taur_var_list[] =
99 {
100 {"disp_x", "mesh1", 3, NODAL_VAR,
101 VAL_COORDX, VAR_DISPX},
102 {"disp_y", "mesh1", 3, NODAL_VAR,
103 VAL_COORDY, VAR_DISPY},
104 {"disp_z", "mesh1", 3, NODAL_VAR,
105 VAL_COORDZ, VAR_DISPZ},
106 {"disp_mag", "mesh1", 3, NODAL_VAR,
107 VAL_COORDX, VAR_DISP_MAG},
108 {"vel_x", "mesh1", 3, NODAL_VAR,
109 VAL_VELX, VAR_NORMAL},
110 {"vel_y", "mesh1", 3, NODAL_VAR,
111 VAL_VELY, VAR_NORMAL},
112 {"vel_z", "mesh1", 3, NODAL_VAR,
113 VAL_VELZ, VAR_NORMAL},
114 {"vel_mag", "mesh1", 3, NODAL_VAR,
115 VAL_VELX, VAR_VEL_MAG},
116 {"acc_x", "mesh1", 3, NODAL_VAR,
117 VAL_ACCX, VAR_NORMAL},
118 {"acc_y", "mesh1", 3, NODAL_VAR,
119 VAL_ACCY, VAR_NORMAL},
120 {"acc_z", "mesh1", 3, NODAL_VAR,
121 VAL_ACCZ, VAR_NORMAL},
122 {"acc_mag", "mesh1", 3, NODAL_VAR,
123 VAL_ACCX, VAR_ACC_MAG},
124 {"temp_x", "mesh1", 3, NODAL_VAR,
125 VAL_TEMPX, VAR_NORMAL},
126 {"temp_y", "mesh1", 3, NODAL_VAR,
127 VAL_TEMPY, VAR_NORMAL},
128 {"temp_z", "mesh1", 3, NODAL_VAR,
129 VAL_TEMPZ, VAR_NORMAL},
130 {"m_xx_bending", "shell_mesh", 4, ZONAL_VAR,
131 VAL_SHELL_RES1, VAR_NORMAL},
132 {"m_yy_bending", "shell_mesh", 4, ZONAL_VAR,
133 VAL_SHELL_RES2, VAR_NORMAL},
134 {"m_xy_bending", "shell_mesh", 4, ZONAL_VAR,
135 VAL_SHELL_RES3, VAR_NORMAL},
136 {"q_xx_shear", "shell_mesh", 4, ZONAL_VAR,
137 VAL_SHELL_RES4, VAR_NORMAL},
138 {"q_yy_shear", "shell_mesh", 4, ZONAL_VAR,
139 VAL_SHELL_RES5, VAR_NORMAL},
140 {"n_xx_normal", "shell_mesh", 4, ZONAL_VAR,
141 VAL_SHELL_RES6, VAR_NORMAL},
142 {"n_yy_normal", "shell_mesh", 4, ZONAL_VAR,
143 VAL_SHELL_RES7, VAR_NORMAL},
144 {"n_xy_normal", "shell_mesh", 4, ZONAL_VAR,
145 VAL_SHELL_RES8, VAR_NORMAL},
146 {"thickness", "shell_mesh", 4, ZONAL_VAR,
147 VAL_SHELL_THICKNESS, VAR_NORMAL},
148 {"int_energy", "shell_mesh", 4, ZONAL_VAR,
149 VAL_SHELL_INT_ENG, VAR_NORMAL},
150 {"surf_stress_1", "shell_mesh", 4, ZONAL_VAR,
151 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_1},
152 {"surf_stress_2", "shell_mesh", 4, ZONAL_VAR,
153 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_2},
154 {"surf_stress_3", "shell_mesh", 4, ZONAL_VAR,
155 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_3},
156 {"surf_stress_4", "shell_mesh", 4, ZONAL_VAR,
157 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_4},
158 {"surf_stress_5", "shell_mesh", 4, ZONAL_VAR,
159 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_5},
160 {"surf_stress_6", "shell_mesh", 4, ZONAL_VAR,
161 VAL_SHELL_MID_SIGX, VAR_SURF_STRESS_6},
162 {"eff_upp_stress", "shell_mesh", 4, ZONAL_VAR,
163 VAL_SHELL_MID_SIGX, VAR_UP_STRESS},
164 {"eff_low_stress", "shell_mesh", 4, ZONAL_VAR,
165 VAL_SHELL_MID_SIGX, VAR_LOW_STRESS},
166 {"eff_max_stress", "shell_mesh", 4, ZONAL_VAR,
167 VAL_SHELL_MID_SIGX, VAR_MAX_STRESS},
168 {"upp_surf_eps", "shell_mesh", 4, ZONAL_VAR,
169 VAL_SHELL_OUT_EPS_EFF, VAR_NORMAL},
170 {"low_surf_eps", "shell_mesh", 4, ZONAL_VAR,
171 VAL_SHELL_IN_EPS_EFF, VAR_NORMAL},
172 {"low_xx_strain", "shell_mesh", 4, ZONAL_VAR,
173 VAL_SHELL_IN_SIGX, VAR_NORMAL},
174 {"low_yy_strain", "shell_mesh", 4, ZONAL_VAR,
175 VAL_SHELL_IN_SIGY, VAR_NORMAL},
176 {"low_zz_strain", "shell_mesh", 4, ZONAL_VAR,
177 VAL_SHELL_IN_SIGZ, VAR_NORMAL},
178 {"low_xy_strain", "shell_mesh", 4, ZONAL_VAR,
179 VAL_SHELL_IN_SIGXY, VAR_NORMAL},
180 {"low_yz_strain", "shell_mesh", 4, ZONAL_VAR,
181 VAL_SHELL_IN_SIGYZ, VAR_NORMAL},
182 {"low_zx_strain", "shell_mesh", 4, ZONAL_VAR,
183 VAL_SHELL_IN_SIGZX, VAR_NORMAL},
184 {"upp_xx_strain", "shell_mesh", 4, ZONAL_VAR,
185 VAL_SHELL_OUT_SIGX, VAR_NORMAL},
186 {"upp_yy_strain", "shell_mesh", 4, ZONAL_VAR,
187 VAL_SHELL_OUT_SIGY, VAR_NORMAL},
188 {"upp_zz_strain", "shell_mesh", 4, ZONAL_VAR,
189 VAL_SHELL_OUT_SIGZ, VAR_NORMAL},
190 {"upp_xy_strain", "shell_mesh", 4, ZONAL_VAR,
191 VAL_SHELL_OUT_SIGXY, VAR_NORMAL},
192 {"upp_yz_strain", "shell_mesh", 4, ZONAL_VAR,
193 VAL_SHELL_OUT_SIGYZ, VAR_NORMAL},
194 {"upp_zx_strain", "shell_mesh", 4, ZONAL_VAR,
195 VAL_SHELL_OUT_SIGZX, VAR_NORMAL},
196 {"mid_xx_strain", "shell_mesh", 4, ZONAL_VAR,
197 VAL_SHELL_MID_SIGX, VAR_NORMAL},
198 {"mid_yy_strain", "shell_mesh", 4, ZONAL_VAR,
199 VAL_SHELL_MID_SIGY, VAR_NORMAL},
200 {"mid_zz_strain", "shell_mesh", 4, ZONAL_VAR,
201 VAL_SHELL_MID_SIGZ, VAR_NORMAL},
202 {"mid_xy_strain", "shell_mesh", 4, ZONAL_VAR,
203 VAL_SHELL_MID_SIGXY, VAR_NORMAL},
204 {"mid_yz_strain", "shell_mesh", 4, ZONAL_VAR,
205 VAL_SHELL_MID_SIGYZ, VAR_NORMAL},
206 {"mid_zx_strain", "shell_mesh", 4, ZONAL_VAR,
207 VAL_SHELL_MID_SIGZX, VAR_NORMAL},
208 {"stress_xx", "hs_mesh", 5, ZONAL_VAR,
209 VAL_HEX_SIGX, VAR_SIGX},
210 {"stress_yy", "hs_mesh", 5, ZONAL_VAR,
211 VAL_HEX_SIGY, VAR_SIGY},
212 {"stress_zz", "hs_mesh", 5, ZONAL_VAR,
213 VAL_HEX_SIGZ, VAR_SIGZ},
214 {"stress_xy", "hs_mesh", 5, ZONAL_VAR,
215 VAL_HEX_SIGXY, VAR_SIGXY},
216 {"stress_yz", "hs_mesh", 5, ZONAL_VAR,
217 VAL_HEX_SIGYZ, VAR_SIGYZ},
218 {"stress_zx", "hs_mesh", 5, ZONAL_VAR,
219 VAL_HEX_SIGZX, VAR_SIGZX},
220 {"stress_eps", "hs_mesh", 5, ZONAL_VAR,
221 VAL_HEX_EPS_EFF, VAR_EPS},
222 {"pressure", "hs_mesh", 5, ZONAL_VAR,
223 VAL_HEX_SIGX, VAR_PRESSURE},
224 {"stress_eff", "hs_mesh", 5, ZONAL_VAR,
225 VAL_HEX_SIGX, VAR_SIG_EFF},
226 {"princ_dev_stress_1", "hs_mesh", 5, ZONAL_VAR,
227 VAL_HEX_SIGX, VAR_DEV_STRESS_1},
228 {"princ_dev_stress_2", "hs_mesh", 5, ZONAL_VAR,
229 VAL_HEX_SIGX, VAR_DEV_STRESS_2},
230 {"princ_dev_stress_3", "hs_mesh", 5, ZONAL_VAR,
231 VAL_HEX_SIGX, VAR_DEV_STRESS_3},
232 {"max_shear_stress", "hs_mesh", 5, ZONAL_VAR,
233 VAL_HEX_SIGX, VAR_MAX_SHEAR_STR},
234 {"princ_stress_1", "hs_mesh", 5, ZONAL_VAR,
235 VAL_HEX_SIGX, VAR_PRINC_STRESS_1},
236 {"princ_stress_2", "hs_mesh", 5, ZONAL_VAR,
237 VAL_HEX_SIGX, VAR_PRINC_STRESS_2},
238 {"princ_stress_3", "hs_mesh", 5, ZONAL_VAR,
239 VAL_HEX_SIGX, VAR_PRINC_STRESS_3},
240 {"temperature", "mesh1", 8, NODAL_VAR,
241 VAL_TEMP, VAR_NORMAL},
242 {"flux_x", "mesh1", 8, NODAL_VAR,
243 VAL_FLUXX, VAR_NORMAL},
244 {"flux_y", "mesh1", 8, NODAL_VAR,
245 VAL_FLUXY, VAR_NORMAL},
246 {"flux_z", "mesh1", 8, NODAL_VAR,
247 VAL_FLUXZ, VAR_NORMAL},
248 {"vel_x", "mesh1", 9, NODAL_VAR,
249 VAL_VELX, VAR_NORMAL},
250 {"vel_y", "mesh1", 9, NODAL_VAR,
251 VAL_VELY, VAR_NORMAL},
252 {"vel_z", "mesh1", 9, NODAL_VAR,
253 VAL_VELZ, VAR_NORMAL},
254 {"vel_mag", "mesh1", 9, NODAL_VAR,
255 VAL_VELX, VAR_VEL_MAG},
256 {"vort_x", "mesh1", 9, NODAL_VAR,
257 VAL_VORTX, VAR_NORMAL},
258 {"vort_y", "mesh1", 9, NODAL_VAR,
259 VAL_VORTY, VAR_NORMAL},
260 {"vort_z", "mesh1", 9, NODAL_VAR,
261 VAL_VORTZ, VAR_NORMAL},
262 {"vort_mag", "mesh1", 9, NODAL_VAR,
263 VAL_VORTX, VAR_VORT_MAG},
264 {"pressure", "hs_mesh", 9, ZONAL_VAR,
265 VAL_PRESSURE, VAR_PRESSURE},
266 {"dummy", "dummy", 10, NODAL_VAR,
267 0, VAR_NORMAL}
268 };
269
270 /*-------------------------------------------------------------------------
271 * Function: fam_name
272 *
273 * Purpose:
274 *
275 * Return: Success: void
276 *
277 * Failure:
278 *
279 * Programmer:
280 *
281 * Modifications:
282 *
283 * Jim Reus, 23 Apr 97
284 * Changed to proto type form.
285 *
286 *-------------------------------------------------------------------------
287 */
288 static void
fam_name(char * basename,int filenumber,char * filename)289 fam_name (char *basename, int filenumber, char *filename)
290 {
291
292 if (filenumber == 0)
293 strcpy(filename, basename);
294 else if (filenumber < 100)
295 sprintf(filename, "%s%02d", basename, filenumber);
296 else
297 sprintf(filename, "%s%03d", basename, filenumber);
298 }
299
300 /*-------------------------------------------------------------------------
301 * Function: fix_title
302 *
303 * Purpose: Fixes the screwy title. For some reason, the title
304 * contains extra padding at character positions
305 * n*8+6 and n*8+7. The title may also contain leading
306 * and trailing whitespace. This routine will remove
307 * all of this. This routine replaces a buggy version
308 * from the old Taurus stuff (see modification section
309 * below).
310 *
311 * Return: void
312 *
313 * Programmer: Robb Matzke, Wed Jan 25 15:46:03 PST 1995
314 *
315 * Modifications:
316 *
317 * Robb Matzke, Wed Jan 25 15:20:03 PST 1995
318 * Removed final `for' loop since the title should have
319 * only 41 characters counting the null terminator.
320 * Fixed final while loop for cases where the title contains
321 * no printable characters or all whitespce. After removing
322 * the junk part of the title and shortening the title, we
323 * zero-fill everything to the right. This is because the
324 * title size is initialized before we know what the title
325 * is and the browser output routines will print all bytes
326 * of the title, including the stuff after the null terminator.
327 *
328 * Eric Brugger, Fri Apr 28 10:09:40 PDT 1995
329 * I modified the routine to correct screwy titles and leave normal
330 * ones alone.
331 *
332 * Jim Reus, 23 Apr 97
333 * Changed to prototype form.
334 *
335 *-------------------------------------------------------------------------
336 */
337 PRIVATE void
fix_title(char * title)338 fix_title (char *title)
339 {
340 int i, j;
341 int fixit;
342
343 /*
344 * Determine if the title needs fixing.
345 */
346 fixit = TRUE;
347 for (i = 0; i < 40; i += 8) {
348 if (title[i + 6] != ' ' || title[i + 7] != ' ')
349 fixit = FALSE;
350 }
351
352 /*
353 * Fix it if necessary.
354 */
355 if (fixit == TRUE) {
356 j = 0;
357 for (i = 0; i < 6; i++)
358 title[j++] = title[i];
359 for (i = 8; i < 14; i++)
360 title[j++] = title[i];
361 for (i = 16; i < 22; i++)
362 title[j++] = title[i];
363 for (i = 24; i < 30; i++)
364 title[j++] = title[i];
365 for (i = 32; i < 38; i++)
366 title[j++] = title[i];
367 for (i = 0; i < 10; i++)
368 title[j++] = ' ';
369 }
370
371 /*
372 * Remove trailing blanks.
373 */
374 j = 39;
375 while (j > 0 && isspace(title[j]))
376 j--;
377 title[j + 1] = '\0';
378 }
379
380 /*-------------------------------------------------------------------------
381 * Function: init_file_info
382 *
383 * Purpose:
384 *
385 * Return: Success: void
386 *
387 * Failure:
388 *
389 * Programmer:
390 *
391 * Modifications:
392 *
393 * Jim Reus, 23 Apr 97
394 * Changed to prototype form.
395 *
396 *-------------------------------------------------------------------------
397 */
398 static void
init_file_info(TAURUSfile * taurus)399 init_file_info (TAURUSfile *taurus)
400 {
401 int i;
402 int nfiles;
403 struct stat statbuf;
404
405 /*
406 * Determine the number of files.
407 */
408 nfiles = 0;
409 fam_name(taurus->basename, nfiles, taurus->filename);
410 while (stat(taurus->filename, &statbuf) != -1) {
411 nfiles++;
412 fam_name(taurus->basename, nfiles, taurus->filename);
413 }
414 taurus->nfiles = nfiles;
415
416 /*
417 * Determine the size of each file in the family.
418 */
419 taurus->filesize = ALLOC_N(int, nfiles);
420
421 for (i = 0; i < nfiles; i++) {
422 fam_name(taurus->basename, i, taurus->filename);
423 stat(taurus->filename, &statbuf);
424 taurus->filesize[i] = statbuf.st_size;
425 }
426 }
427
428 /*-------------------------------------------------------------------------
429 * Function: taurus_read
430 *
431 * Purpose:
432 *
433 * Return: Success: 0
434 *
435 * Failure: -1
436 *
437 * Programmer:
438 *
439 * Modifications:
440 *
441 * Eric Brugger, Fri Apr 28 10:41:17 PDT 1995
442 * I modified the routine to handle the case where the address is
443 * not in the specified file.
444 *
445 * Jim Reus, 23 Apr 97
446 * Changed to prototype form.
447 *
448 *-------------------------------------------------------------------------
449 */
450 static int
taurus_read(TAURUSfile * taurus,int ifile,int iadd,int length,char * buffer)451 taurus_read (TAURUSfile *taurus, int ifile, int iadd, int length, char *buffer)
452 {
453 int n;
454 int ibuf;
455 int idisk;
456
457 /*
458 * Skip to the correct file if the address is not in the
459 * specified file.
460 */
461 while (iadd > taurus->filesize[ifile]) {
462 iadd -= taurus->filesize[ifile];
463 ifile++;
464 }
465
466 /*
467 * Read the file.
468 */
469 ibuf = 0;
470 idisk = iadd;
471 while (length > 0) {
472 /*
473 * If the desired file is not open, close the current file
474 * and open the desired file.
475 */
476 if (taurus->ifile != ifile) {
477 if (taurus->fd != -1)
478 close(taurus->fd);
479 fam_name(taurus->basename, ifile, taurus->filename);
480 if ((taurus->fd = open(taurus->filename, O_RDONLY)) < 0) {
481 return (-1);
482 }
483 taurus->ifile = ifile;
484 }
485
486 /*
487 * Read the maximum amount from the current file.
488 */
489 n = MIN(taurus->filesize[ifile] - idisk, length);
490 lseek(taurus->fd, idisk, SEEK_SET);
491
492 if (read(taurus->fd, &buffer[ibuf], n) != n) {
493 return (-2);
494 }
495
496 ibuf += n;
497 length -= n;
498 ifile++;
499 idisk = 0;
500 }
501
502 return (0);
503 }
504
505 /*-------------------------------------------------------------------------
506 * Function: init_state_info
507 *
508 * Purpose:
509 *
510 * Return: Success: void
511 *
512 * Failure:
513 *
514 * Programmer:
515 *
516 * Modifications:
517 * Eric Brugger, Tue Mar 28 15:03:37 PST 1995
518 * I modified the routines to read a topaz3d data file.
519 *
520 * Eric Brugger, Wed Apr 26 13:52:00 PDT 1995
521 * I modified the routine to set the directory to none.
522 *
523 * Eric Brugger, Fri Apr 28 10:12:02 PDT 1995
524 * I modified the routine to handle states that overlapped several
525 * files.
526 *
527 * Eric Brugger, Thu Jul 27 12:49:40 PDT 1995
528 * I modified the routine to handle files generated by hydra.
529 *
530 * Jim Reus, 23 Apr 97
531 * Changed to prototype form.
532 *
533 *-------------------------------------------------------------------------
534 */
535 static void
init_state_info(TAURUSfile * taurus)536 init_state_info (TAURUSfile *taurus)
537 {
538 int i;
539 int geomsize;
540 int statesize;
541 int totsize;
542 int maxstates;
543 int nstates;
544 int loc;
545 int ifile;
546 int nfiles;
547 int nv1dact, nv2dact, nv3dact;
548
549 nfiles = taurus->nfiles;
550
551 nv1dact = taurus->nv1d;
552 nv2dact = taurus->nv2d;
553 nv3dact = taurus->nv3d;
554 if (taurus->activ >= 1000 && taurus->activ <= 1005) {
555 if (taurus->nel2 > 0)
556 nv1dact++;
557 if (taurus->nel4 > 0)
558 nv2dact++;
559 if (taurus->nel8 > 0)
560 nv3dact++;
561 }
562
563 /*
564 * Determine the file and disk address for the start of each
565 * state in the database.
566 */
567 geomsize = (taurus->ndim * taurus->numnp +
568 9 * taurus->nel8 +
569 5 * taurus->nel4 +
570 6 * taurus->nel2) * sizeof(int);
571
572 switch (taurus->icode) {
573 /*
574 * topaz3d.
575 */
576 case 1:
577 /*
578 * This is an extension to the taurus data base. If it
579 * is four then fluxes are present.
580 */
581 if (taurus->it == 4)
582 statesize = (4 * taurus->numnp + 1) * sizeof(int);
583
584 else
585 statesize = (1 * taurus->numnp + 1) * sizeof(int);
586
587 break;
588 /*
589 * dyna3d or nike3d.
590 */
591 case 2:
592 case 6:
593 case 200:
594 statesize = (taurus->it * taurus->numnp +
595 taurus->ndim * taurus->numnp *
596 (taurus->iu + taurus->iv + taurus->ia) +
597 taurus->nel8 * nv3dact +
598 taurus->nel4 * nv2dact +
599 taurus->nel2 * nv1dact +
600 taurus->nglbv + 1) * sizeof(int);
601
602 break;
603 }
604
605 totsize = 0;
606 for (i = 0; i < nfiles; i++)
607 totsize += taurus->filesize[i];
608 maxstates = (totsize / statesize) + 1;
609 taurus->state_file = ALLOC_N(int, maxstates);
610 taurus->state_loc = ALLOC_N(int, maxstates);
611 taurus->state_time = ALLOC_N(float, maxstates);
612
613 loc = 64 * sizeof(int) + geomsize;
614
615 ifile = 0;
616 while (loc >= taurus->filesize[ifile]) {
617 loc -= taurus->filesize[ifile];
618 ifile++;
619 }
620
621 for (nstates = 0;; nstates++) {
622 if (nfiles == 1) {
623 if ((loc + statesize) > taurus->filesize[ifile])
624 break;
625 }
626 else if (statesize <= taurus->filesize[1]) {
627 if ((loc + statesize) > taurus->filesize[ifile]) {
628 ifile++;
629 loc = 0;
630 if (ifile >= nfiles)
631 break;
632 }
633 }
634 else {
635 while (loc > 0) {
636 loc -= taurus->filesize[ifile];
637 ifile++;
638 }
639 loc = 0;
640 if (ifile >= nfiles)
641 break;
642 }
643
644 taurus->state_file[nstates] = ifile;
645 taurus->state_loc[nstates] = loc;
646
647 loc += statesize;
648 }
649
650 taurus->nstates = nstates;
651 taurus->state = -1;
652 taurus->idir = -1;
653
654 /*
655 * Read in the time for each state.
656 */
657 for (i = 0; i < nstates; i++)
658 taurus_read(taurus, taurus->state_file[i], taurus->state_loc[i],
659 sizeof(float), (char*)&taurus->state_time[i]);
660 }
661
662 /*-------------------------------------------------------------------------
663 * Function: init_var_info
664 *
665 * Purpose:
666 *
667 * Return: Success: void
668 *
669 * Failure:
670 *
671 * Programmer:
672 *
673 * Modifications:
674 * Eric Brugger, Tue Mar 28 15:03:37 PST 1995
675 * I modified the routines to read a topaz3d data file.
676 *
677 * Eric Brugger, Thu Apr 27 08:47:01 PDT 1995
678 * I modified the code to work properly with state directories.
679 *
680 * Eric Brugger, Thu Jul 27 12:49:40 PDT 1995
681 * I modified the routine to handle files generated by hydra.
682 *
683 * Jim Reus, 23 Apr 97
684 * Changed to prototype form.
685 *
686 *-------------------------------------------------------------------------
687 */
688 static void
init_var_info(TAURUSfile * taurus)689 init_var_info (TAURUSfile *taurus)
690 {
691 int i;
692 int loc;
693
694 /*
695 * Set up the default values.
696 */
697 for (i = 0; i < MAX_VAL; i++)
698 taurus->var_start[i] = -1;
699
700 loc = (1 + taurus->nglbv) * sizeof(float);
701
702 /*
703 * Topaz3d data.
704 */
705 if (taurus->icode == 1) {
706 taurus->var_start[VAL_TEMP] = loc;
707 taurus->var_ncomps[VAL_TEMP] = 1;
708 taurus->var_len[VAL_TEMP] = taurus->numnp;
709 taurus->var_offset[VAL_TEMP] = 0;
710 loc += taurus->numnp * sizeof(float);
711
712 /*
713 * This is an extension to the taurus data base. If it
714 * is four then fluxes are present.
715 */
716 if (taurus->it == 4) {
717 for (i = VAL_FLUXX; i <= VAL_FLUXZ; i++) {
718 taurus->var_start[i] = loc;
719 taurus->var_ncomps[i] = 3;
720 taurus->var_len[i] = taurus->numnp;
721 }
722 taurus->var_offset[VAL_FLUXX] = 0;
723 taurus->var_offset[VAL_FLUXY] = 1;
724 taurus->var_offset[VAL_FLUXZ] = 2;
725
726 loc += 3 * taurus->numnp * sizeof(float);
727 }
728
729 return;
730 }
731
732 /*
733 * Hydra data.
734 */
735 if (taurus->icode == 200) {
736 if (taurus->iv != 0) {
737 for (i = VAL_VELX; i <= VAL_VELZ; i++) {
738 taurus->var_start[i] = loc;
739 taurus->var_ncomps[i] = 3;
740 taurus->var_len[i] = taurus->numnp;
741 }
742 taurus->var_offset[VAL_VELX] = 0;
743 taurus->var_offset[VAL_VELY] = 1;
744 taurus->var_offset[VAL_VELZ] = 2;
745 loc += taurus->ndim * taurus->numnp * sizeof(float);
746 }
747
748 if (taurus->ia != 0) {
749 for (i = VAL_VORTX; i <= VAL_VORTZ; i++) {
750 taurus->var_start[i] = loc;
751 taurus->var_ncomps[i] = 3;
752 taurus->var_len[i] = taurus->numnp;
753 }
754 taurus->var_offset[VAL_VORTX] = 0;
755 taurus->var_offset[VAL_VORTY] = 1;
756 taurus->var_offset[VAL_VORTZ] = 2;
757 loc += taurus->ndim * taurus->numnp * sizeof(float);
758 }
759
760 if (taurus->nel8 > 0 && taurus->nv3d == 7) {
761 taurus->var_start[VAL_PRESSURE] = loc;
762 taurus->var_ncomps[VAL_PRESSURE] = 7;
763 taurus->var_len[VAL_PRESSURE] = taurus->nel8;
764 taurus->var_offset[VAL_PRESSURE] = 0;
765 loc += taurus->nv3d * taurus->nel8 * sizeof(float);
766 }
767
768 return;
769 }
770
771 /*
772 * Initial nodal coordinates (for displacement calculations).
773 */
774 if (taurus->iu != 0) {
775 for (i = VAL_COORDX; i <= VAL_COORDZ; i++) {
776 taurus->var_start[i] = 64 * sizeof(int);
777
778 taurus->var_ncomps[i] = 3;
779 taurus->var_len[i] = taurus->numnp;
780 }
781 taurus->var_offset[VAL_COORDX] = 0;
782 taurus->var_offset[VAL_COORDY] = 1;
783 taurus->var_offset[VAL_COORDZ] = 2;
784 loc += taurus->ndim * taurus->numnp * sizeof(float);
785 }
786
787 /*
788 * Nodal velocities.
789 */
790 if (taurus->iv != 0) {
791 for (i = VAL_VELX; i <= VAL_VELZ; i++) {
792 taurus->var_start[i] = loc;
793 taurus->var_ncomps[i] = 3;
794 taurus->var_len[i] = taurus->numnp;
795 }
796 taurus->var_offset[VAL_VELX] = 0;
797 taurus->var_offset[VAL_VELY] = 1;
798 taurus->var_offset[VAL_VELZ] = 2;
799 loc += taurus->ndim * taurus->numnp * sizeof(float);
800 }
801
802 /*
803 * Nodal accelerations.
804 */
805 if (taurus->ia != 0) {
806 for (i = VAL_ACCX; i <= VAL_ACCZ; i++) {
807 taurus->var_start[i] = loc;
808 taurus->var_ncomps[i] = 3;
809 taurus->var_len[i] = taurus->numnp;
810 }
811 taurus->var_offset[VAL_ACCX] = 0;
812 taurus->var_offset[VAL_ACCY] = 1;
813 taurus->var_offset[VAL_ACCZ] = 2;
814 loc += taurus->ndim * taurus->numnp * sizeof(float);
815 }
816
817 /*
818 * Nodal temperatures.
819 */
820 if (taurus->it != 0) {
821 for (i = VAL_TEMPX; i <= VAL_TEMPZ; i++) {
822 taurus->var_start[i] = loc;
823 taurus->var_ncomps[i] = 3;
824 taurus->var_len[i] = taurus->numnp;
825 }
826 taurus->var_offset[VAL_TEMPX] = 0;
827 taurus->var_offset[VAL_TEMPY] = 1;
828 taurus->var_offset[VAL_TEMPZ] = 2;
829 loc += taurus->ndim * taurus->numnp * sizeof(float);
830 }
831
832 /*
833 * Brick data.
834 */
835 if (taurus->nel8 > 0 && taurus->nv3d == 7) {
836 for (i = VAL_HEX_SIGX; i <= VAL_HEX_EPS_EFF; i++) {
837 taurus->var_start[i] = loc;
838 taurus->var_ncomps[i] = 7;
839 taurus->var_len[i] = taurus->nel8;
840 }
841 taurus->var_offset[VAL_HEX_SIGX] = 0;
842 taurus->var_offset[VAL_HEX_SIGY] = 1;
843 taurus->var_offset[VAL_HEX_SIGZ] = 2;
844 taurus->var_offset[VAL_HEX_SIGXY] = 3;
845 taurus->var_offset[VAL_HEX_SIGYZ] = 4;
846 taurus->var_offset[VAL_HEX_SIGZX] = 5;
847 taurus->var_offset[VAL_HEX_EPS_EFF] = 6;
848 loc += taurus->nv3d * taurus->nel8 * sizeof(float);
849 }
850
851 /*
852 * Beam data.
853 */
854 if (taurus->nel2 > 0 && taurus->nv1d == 6) {
855 loc += taurus->nv1d * taurus->nel2 * sizeof(float);
856 }
857
858 /*
859 * Shell data.
860 */
861 if (taurus->nel4 > 0 && taurus->nv2d >= 33) {
862 for (i = VAL_SHELL_MID_SIGX; i <= VAL_SHELL_INT_ENG; i++) {
863 taurus->var_start[i] = loc;
864 taurus->var_ncomps[i] = taurus->nv2d;
865 taurus->var_len[i] = taurus->nel4;
866 }
867 taurus->var_offset[VAL_SHELL_MID_SIGX] = 0;
868 taurus->var_offset[VAL_SHELL_MID_SIGY] = 1;
869 taurus->var_offset[VAL_SHELL_MID_SIGZ] = 2;
870 taurus->var_offset[VAL_SHELL_MID_SIGXY] = 3;
871 taurus->var_offset[VAL_SHELL_MID_SIGYZ] = 4;
872 taurus->var_offset[VAL_SHELL_MID_SIGZX] = 5;
873 taurus->var_offset[VAL_SHELL_MID_EPS_EFF] = 6; /* not used */
874 taurus->var_offset[VAL_SHELL_IN_SIGX] = 7;
875 taurus->var_offset[VAL_SHELL_IN_SIGY] = 8;
876 taurus->var_offset[VAL_SHELL_IN_SIGZ] = 9;
877 taurus->var_offset[VAL_SHELL_IN_SIGXY] = 10;
878 taurus->var_offset[VAL_SHELL_IN_SIGYZ] = 11;
879 taurus->var_offset[VAL_SHELL_IN_SIGZX] = 12;
880 taurus->var_offset[VAL_SHELL_IN_EPS_EFF] = 13;
881 taurus->var_offset[VAL_SHELL_OUT_SIGX] = 14;
882 taurus->var_offset[VAL_SHELL_OUT_SIGY] = 15;
883 taurus->var_offset[VAL_SHELL_OUT_SIGZ] = 16;
884 taurus->var_offset[VAL_SHELL_OUT_SIGXY] = 17;
885 taurus->var_offset[VAL_SHELL_OUT_SIGYZ] = 18;
886 taurus->var_offset[VAL_SHELL_OUT_SIGZX] = 19;
887 taurus->var_offset[VAL_SHELL_OUT_EPS_EFF] = 20;
888 taurus->var_offset[VAL_SHELL_RES1] = 21;
889 taurus->var_offset[VAL_SHELL_RES2] = 22;
890 taurus->var_offset[VAL_SHELL_RES3] = 23;
891 taurus->var_offset[VAL_SHELL_RES4] = 24;
892 taurus->var_offset[VAL_SHELL_RES5] = 25;
893 taurus->var_offset[VAL_SHELL_RES6] = 26;
894 taurus->var_offset[VAL_SHELL_RES7] = 27;
895 taurus->var_offset[VAL_SHELL_RES8] = 28;
896 taurus->var_offset[VAL_SHELL_THICKNESS] = 29;
897 taurus->var_offset[VAL_SHELL_ELDEP1] = 30; /* not used */
898 taurus->var_offset[VAL_SHELL_ELDEP2] = 31; /* not used */
899 taurus->var_offset[VAL_SHELL_INT_ENG] = 32;
900 /*
901 * The variables in the following if block are not used.
902 */
903 if (taurus->nv2d == 45 || taurus->nv2d == 46) {
904 for (i = VAL_SHELL_EPSX_IN; i <= VAL_SHELL_EPSZX_OUT; i++) {
905 taurus->var_start[i] = loc;
906 taurus->var_ncomps[i] = taurus->nv2d;
907 taurus->var_len[i] = taurus->nel4;
908 }
909 taurus->var_offset[VAL_SHELL_EPSX_IN] = 33;
910 taurus->var_offset[VAL_SHELL_EPSY_IN] = 34;
911 taurus->var_offset[VAL_SHELL_EPSZ_IN] = 35;
912 taurus->var_offset[VAL_SHELL_EPSXY_IN] = 36;
913 taurus->var_offset[VAL_SHELL_EPSYZ_IN] = 37;
914 taurus->var_offset[VAL_SHELL_EPSZX_IN] = 38;
915 taurus->var_offset[VAL_SHELL_EPSX_OUT] = 39;
916 taurus->var_offset[VAL_SHELL_EPSY_OUT] = 40;
917 taurus->var_offset[VAL_SHELL_EPSZ_OUT] = 41;
918 taurus->var_offset[VAL_SHELL_EPSXY_OUT] = 42;
919 taurus->var_offset[VAL_SHELL_EPSYZ_OUT] = 43;
920 taurus->var_offset[VAL_SHELL_EPSZX_OUT] = 44;
921 }
922 loc += taurus->nv2d * taurus->nel4 * sizeof(float);
923 }
924 }
925
926 /*-------------------------------------------------------------------------
927 * Function: init_mat_info
928 *
929 * Purpose:
930 *
931 * Return: Success: void
932 *
933 * Failure:
934 *
935 * Programmer:
936 *
937 * Modifications:
938 * Eric Brugger, Mon Aug 28 13:37:09 PDT 1995
939 * I modified the routine to find the actual materials referenced
940 * rather than the maximum material number and assume that they were
941 * all used.
942 *
943 * Jim Reus, 23 Apr 97
944 * Changed to prototype form.
945 *
946 *-------------------------------------------------------------------------
947 */
948 static void
init_mat_info(TAURUSfile * taurus)949 init_mat_info (TAURUSfile *taurus)
950 {
951 int i;
952 int iadd, ibuf, imat;
953 int len;
954 int lbuf;
955 int *buf, *buf2;
956 int maxmat;
957 int nmat;
958 int *matnos;
959
960 maxmat = 0;
961
962 /*
963 * Set up for reading the nodelists and material data.
964 */
965 iadd = 64 * sizeof(int) + taurus->numnp * taurus->ndim * sizeof(float);
966
967 ibuf = 0;
968 lbuf = (9 * taurus->nel8) + (5 * taurus->nel4) + (6 * taurus->nel2);
969 buf = ALLOC_N(int, lbuf);
970
971 /*
972 * Read the hexahedron elements.
973 */
974 if (taurus->nel8 > 0) {
975 len = 9 * taurus->nel8 * sizeof(int);
976
977 taurus_read(taurus, 0, iadd, len, (char*)buf);
978
979 for (i = 0; i < taurus->nel8; i++)
980 if (buf[i * 9 + 8] > maxmat)
981 maxmat = buf[i * 9 + 8];
982
983 iadd += len;
984 ibuf += 9 * taurus->nel8;
985 }
986
987 /*
988 * Read the beam elements.
989 */
990 if (taurus->nel2 > 0) {
991 len = 6 * taurus->nel2 * sizeof(int);
992
993 taurus_read(taurus, 0, iadd, len, (char*)(&buf[ibuf]));
994
995 for (i = 0; i < taurus->nel2; i++)
996 if (buf[ibuf + i * 6 + 5] > maxmat)
997 maxmat = buf[ibuf + i * 6 + 5];
998
999 iadd += len;
1000 ibuf += 6 * taurus->nel2;
1001 }
1002
1003 /*
1004 * Read the shell elements.
1005 */
1006 if (taurus->nel4 > 0) {
1007 len = 5 * taurus->nel4 * sizeof(int);
1008
1009 taurus_read(taurus, 0, iadd, len, (char*)(&buf[ibuf]));
1010
1011 for (i = 0; i < taurus->nel4; i++)
1012 if (buf[ibuf + i * 5 + 4] > maxmat)
1013 maxmat = buf[ibuf + i * 5 + 4];
1014 }
1015
1016 /*
1017 * Find all the materials referenced.
1018 */
1019 buf2 = ALLOC_N(int, maxmat);
1020
1021 for (i = 0; i < maxmat; i++)
1022 buf2[i] = 0;
1023
1024 ibuf = 0;
1025 if (taurus->nel8 > 0) {
1026 for (i = 0; i < taurus->nel8; i++)
1027 buf2[buf[i * 9 + 8] - 1] = 1;
1028 ibuf += taurus->nel8 * 9;
1029 }
1030 if (taurus->nel2 > 0) {
1031 for (i = 0; i < taurus->nel2; i++)
1032 buf2[buf[ibuf + i * 6 + 5] - 1] = 1;
1033 ibuf += taurus->nel2 * 6;
1034 }
1035 if (taurus->nel4 > 0) {
1036 for (i = 0; i < taurus->nel4; i++)
1037 buf2[buf[ibuf + i * 5 + 4] - 1] = 1;
1038 }
1039
1040 /*
1041 * Find nmat, and set the matnos array.
1042 */
1043 nmat = 0;
1044 for (i = 0; i < maxmat; i++)
1045 nmat += buf2[i];
1046 matnos = ALLOC_N(int, nmat);
1047
1048 imat = 0;
1049 for (i = 0; i < maxmat; i++)
1050 if (buf2[i] == 1)
1051 matnos[imat++] = i + 1;
1052
1053 FREE(buf);
1054 FREE(buf2);
1055
1056 taurus->nmat = nmat;
1057 taurus->matnos = matnos;
1058 }
1059
1060 /*-------------------------------------------------------------------------
1061 * Function: taurus_calc
1062 *
1063 * Purpose:
1064 *
1065 * Return: Success: void
1066 *
1067 * Failure:
1068 *
1069 * Programmer:
1070 *
1071 * Modifications:
1072 * Eric Brugger, Thu Apr 27 08:47:01 PDT 1995
1073 * I modified the code to work properly with state directories.
1074 *
1075 * Eric Brugger, Fri Jul 28 08:41:03 PDT 1995
1076 * I added a calculation for the vorticity magnitude.
1077 *
1078 * Jim Reus, 23 Apr 97
1079 * Changed to prototype form.
1080 *
1081 *-------------------------------------------------------------------------
1082 */
1083 static void
taurus_calc(TAURUSfile * taurus,float * buf,int lbuf,int ncomps,int offset,int var_id,float * var,int ivar)1084 taurus_calc (TAURUSfile *taurus, float *buf, int lbuf, int ncomps, int offset,
1085 int var_id, float *var, int ivar)
1086 {
1087 int ibuf;
1088 double t1, t2, t3, t4, t5;
1089 double pr, aa, bb, cc, dd, angp;
1090 float *buf2, *buf3, *buf4;
1091
1092 switch (var_id) {
1093 case VAR_NORMAL:
1094 case VAR_SIGX:
1095 case VAR_SIGY:
1096 case VAR_SIGZ:
1097 case VAR_SIGXY:
1098 case VAR_SIGYZ:
1099 case VAR_SIGZX:
1100 case VAR_EPS:
1101 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1102 var[ivar] = buf[ibuf];
1103 ivar++;
1104 }
1105 break;
1106 case VAR_PRESSURE:
1107 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1108 var[ivar] = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1109 ivar++;
1110 }
1111 break;
1112 case VAR_SIG_EFF:
1113 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1114 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1115 buf[ibuf] += pr;
1116 buf[ibuf + 1] += pr;
1117 buf[ibuf + 2] += pr;
1118 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1119 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1120 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1121 var[ivar] = sqrt(3.0 * fabs(aa));
1122 ivar++;
1123 }
1124 break;
1125 case VAR_DEV_STRESS_1:
1126 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1127 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1128 buf[ibuf] += pr;
1129 buf[ibuf + 1] += pr;
1130 buf[ibuf + 2] += pr;
1131 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1132 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1133 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1134 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1135 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1136 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1137 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1138 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1139 buf[ibuf] = aa;
1140 buf[ibuf + 1] = bb;
1141 }
1142 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1143 if (buf[ibuf] < 1.0e-25)
1144 var[ivar] = 0.;
1145 else {
1146 aa = buf[ibuf];
1147 bb = buf[ibuf + 1];
1148 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1149 cc = MAX(MIN(cc, 1.0), -1.0);
1150 angp = acos(cc) / 3.0;
1151 dd = 2.0 * sqrt(aa / 3.0);
1152 var[ivar] = dd * cos(angp);
1153 }
1154 ivar++;
1155 }
1156 break;
1157 case VAR_DEV_STRESS_2:
1158 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1159 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1160 buf[ibuf] += pr;
1161 buf[ibuf + 1] += pr;
1162 buf[ibuf + 2] += pr;
1163 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1164 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1165 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1166 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1167 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1168 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1169 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1170 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1171 buf[ibuf] = aa;
1172 buf[ibuf + 1] = bb;
1173 }
1174 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1175 if (buf[ibuf] < 1.0e-25)
1176 var[ivar] = 0.;
1177 else {
1178 aa = buf[ibuf];
1179 bb = buf[ibuf + 1];
1180 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1181 cc = MAX(MIN(cc, 1.0), -1.0);
1182 angp = acos(cc) / 3.0;
1183 dd = 2.0 * sqrt(aa / 3.0);
1184 var[ivar] = dd * cos(angp + ftpi);
1185 }
1186 ivar++;
1187 }
1188 break;
1189 case VAR_DEV_STRESS_3:
1190 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1191 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1192 buf[ibuf] += pr;
1193 buf[ibuf + 1] += pr;
1194 buf[ibuf + 2] += pr;
1195 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1196 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1197 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1198 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1199 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1200 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1201 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1202 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1203 buf[ibuf] = aa;
1204 buf[ibuf + 1] = bb;
1205 }
1206 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1207 if (buf[ibuf] < 1.0e-25)
1208 var[ivar] = 0.;
1209 else {
1210 aa = buf[ibuf];
1211 bb = buf[ibuf + 1];
1212 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1213 cc = MAX(MIN(cc, 1.0), -1.0);
1214 angp = acos(cc) / 3.0;
1215 dd = 2.0 * sqrt(aa / 3.0);
1216 var[ivar] = dd * cos(angp + ttpi);
1217 }
1218 ivar++;
1219 }
1220 break;
1221 case VAR_MAX_SHEAR_STR:
1222 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1223 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1224 buf[ibuf] += pr;
1225 buf[ibuf + 1] += pr;
1226 buf[ibuf + 2] += pr;
1227 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1228 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1229 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1230 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1231 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1232 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1233 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1234 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1235 buf[ibuf] = aa;
1236 buf[ibuf + 1] = bb;
1237 }
1238 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1239 if (buf[ibuf] < 1.0e-25)
1240 var[ivar] = 0.;
1241 else {
1242 aa = buf[ibuf];
1243 bb = buf[ibuf + 1];
1244 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1245 cc = MAX(MIN(cc, 1.0), -1.0);
1246 angp = acos(cc) / 3.0;
1247 dd = sqrt(aa / 3.0);
1248 var[ivar] = dd * (cos(angp) - cos(angp + ttpi));
1249 }
1250 ivar++;
1251 }
1252 break;
1253 case VAR_PRINC_STRESS_1:
1254 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1255 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1256 buf[ibuf] += pr;
1257 buf[ibuf + 1] += pr;
1258 buf[ibuf + 2] += pr;
1259 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1260 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1261 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1262 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1263 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1264 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1265 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1266 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1267 buf[ibuf] = aa;
1268 buf[ibuf + 1] = bb;
1269 buf[ibuf + 2] = pr;
1270 }
1271 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1272 if (buf[ibuf] < 1.0e-25)
1273 var[ivar] = -buf[ibuf + 2];
1274 else {
1275 aa = buf[ibuf];
1276 bb = buf[ibuf + 1];
1277 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1278 cc = MAX(MIN(cc, 1.0), -1.0);
1279 angp = acos(cc) / 3.0;
1280 dd = 2.0 * sqrt(aa / 3.0);
1281 var[ivar] = dd * cos(angp) - buf[ibuf + 2];
1282 }
1283 ivar++;
1284 }
1285 break;
1286 case VAR_PRINC_STRESS_2:
1287 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1288 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1289 buf[ibuf] += pr;
1290 buf[ibuf + 1] += pr;
1291 buf[ibuf + 2] += pr;
1292 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1293 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1294 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1295 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1296 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1297 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1298 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1299 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1300 buf[ibuf] = aa;
1301 buf[ibuf + 1] = bb;
1302 buf[ibuf + 2] = pr;
1303 }
1304 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1305 if (buf[ibuf] < 1.0e-25)
1306 var[ivar] = -buf[ibuf + 2];
1307 else {
1308 aa = buf[ibuf];
1309 bb = buf[ibuf + 1];
1310 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1311 cc = MAX(MIN(cc, 1.0), -1.0);
1312 angp = acos(cc) / 3.0;
1313 dd = 2.0 * sqrt(aa / 3.0);
1314 var[ivar] = dd * cos(angp + ftpi) - buf[ibuf + 2];
1315 }
1316 ivar++;
1317 }
1318 break;
1319 case VAR_PRINC_STRESS_3:
1320 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1321 pr = -(buf[ibuf] + buf[ibuf + 1] + buf[ibuf + 2]) / 3.0;
1322 buf[ibuf] += pr;
1323 buf[ibuf + 1] += pr;
1324 buf[ibuf + 2] += pr;
1325 aa = buf[ibuf + 3] * buf[ibuf + 3] + buf[ibuf + 4] * buf[ibuf + 4] +
1326 buf[ibuf + 5] * buf[ibuf + 5] - buf[ibuf] * buf[ibuf + 1] -
1327 buf[ibuf + 1] * buf[ibuf + 2] - buf[ibuf] * buf[ibuf + 2];
1328 bb = buf[ibuf] * buf[ibuf + 4] * buf[ibuf + 4] +
1329 buf[ibuf + 1] * buf[ibuf + 5] * buf[ibuf + 5] +
1330 buf[ibuf + 2] * buf[ibuf + 3] * buf[ibuf + 3] -
1331 buf[ibuf] * buf[ibuf + 1] * buf[ibuf + 2] - 2.0 *
1332 buf[ibuf + 3] * buf[ibuf + 4] * buf[ibuf + 5];
1333 buf[ibuf] = aa;
1334 buf[ibuf + 1] = bb;
1335 buf[ibuf + 2] = pr;
1336 }
1337 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1338 if (buf[ibuf] < 1.0e-25)
1339 var[ivar] = -buf[ibuf + 2];
1340 else {
1341 aa = buf[ibuf];
1342 bb = buf[ibuf + 1];
1343 cc = -sqrt(27.0 / aa) * bb * 0.5 / aa;
1344 cc = MAX(MIN(cc, 1.0), -1.0);
1345 angp = acos(cc) / 3.0;
1346 dd = 2.0 * sqrt(aa / 3.0);
1347 var[ivar] = dd * cos(angp + ttpi) - buf[ibuf + 2];
1348 }
1349 ivar++;
1350 }
1351 break;
1352 case VAR_DISPX:
1353 buf2 = taurus->coords[0];
1354 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1355 var[ivar] = buf2[ivar] - buf[ibuf];
1356 ivar++;
1357 }
1358 break;
1359 case VAR_DISPY:
1360 buf2 = taurus->coords[1];
1361 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1362 var[ivar] = buf2[ivar] - buf[ibuf];
1363 ivar++;
1364 }
1365 break;
1366 case VAR_DISPZ:
1367 buf2 = taurus->coords[2];
1368 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1369 var[ivar] = buf2[ivar] - buf[ibuf];
1370 ivar++;
1371 }
1372 break;
1373 case VAR_DISP_MAG:
1374 buf2 = taurus->coords[0];
1375 buf3 = taurus->coords[1];
1376 buf4 = taurus->coords[2];
1377 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1378 var[ivar] = sqrt((buf2[ivar] - buf[ibuf]) *
1379 (buf2[ivar] - buf[ibuf]) +
1380 (buf3[ivar] - buf[ibuf + 1]) *
1381 (buf3[ivar] - buf[ibuf + 1]) +
1382 (buf4[ivar] - buf[ibuf + 2]) *
1383 (buf4[ivar] - buf[ibuf + 2]));
1384 ivar++;
1385 }
1386 break;
1387 case VAR_VEL_MAG:
1388 case VAR_ACC_MAG:
1389 case VAR_VORT_MAG:
1390 for (ibuf = offset; ibuf < lbuf; ibuf += ncomps) {
1391 var[ivar] = sqrt(buf[ibuf] * buf[ibuf] +
1392 buf[ibuf + 1] * buf[ibuf + 1] +
1393 buf[ibuf + 2] * buf[ibuf + 2]);
1394 ivar++;
1395 }
1396 break;
1397 case VAR_SURF_STRESS_1:
1398 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1399 var[ivar] = (buf[ibuf + 26] / buf[ibuf + 29]) + 6.0 *
1400 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1401 ivar++;
1402 }
1403 break;
1404 case VAR_SURF_STRESS_2:
1405 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1406 var[ivar] = (buf[ibuf + 26] / buf[ibuf + 29]) - 6.0 *
1407 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1408 ivar++;
1409 }
1410 break;
1411 case VAR_SURF_STRESS_3:
1412 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1413 var[ivar] = (buf[ibuf + 27] / buf[ibuf + 29]) + 6.0 *
1414 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1415 ivar++;
1416 }
1417 break;
1418 case VAR_SURF_STRESS_4:
1419 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1420 var[ivar] = (buf[ibuf + 27] / buf[ibuf + 29]) - 6.0 *
1421 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1422 ivar++;
1423 }
1424 break;
1425 case VAR_SURF_STRESS_5:
1426 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1427 var[ivar] = (buf[ibuf + 28] / buf[ibuf + 29]) + 6.0 *
1428 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1429 ivar++;
1430 }
1431 break;
1432 case VAR_SURF_STRESS_6:
1433 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1434 var[ivar] = (buf[ibuf + 28] / buf[ibuf + 29]) - 6.0 *
1435 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1436 ivar++;
1437 }
1438 break;
1439 case VAR_UP_STRESS:
1440 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1441 t1 = (buf[ibuf + 26] / buf[ibuf + 29]) + 6.0 *
1442 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1443 t2 = (buf[ibuf + 27] / buf[ibuf + 29]) + 6.0 *
1444 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1445 t3 = (buf[ibuf + 28] / buf[ibuf + 29]) + 6.0 *
1446 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1447 var[ivar] = sqrt(t1 * t1 - t1 * t2 + t2 * t2 + 3.0 * t3 * t3);
1448 ivar++;
1449 }
1450 break;
1451 case VAR_LOW_STRESS:
1452 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1453 t1 = (buf[ibuf + 26] / buf[ibuf + 29]) - 6.0 *
1454 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1455 t2 = (buf[ibuf + 27] / buf[ibuf + 29]) - 6.0 *
1456 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1457 t3 = (buf[ibuf + 28] / buf[ibuf + 29]) - 6.0 *
1458 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1459 var[ivar] = sqrt(t1 * t1 - t1 * t2 + t2 * t2 + 3.0 * t3 * t3);
1460 ivar++;
1461 }
1462 break;
1463 case VAR_MAX_STRESS:
1464 for (ibuf = 0; ibuf < lbuf; ibuf += ncomps) {
1465 t1 = (buf[ibuf + 26] / buf[ibuf + 29]) + 6.0 *
1466 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1467 t2 = (buf[ibuf + 27] / buf[ibuf + 29]) + 6.0 *
1468 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1469 t3 = (buf[ibuf + 28] / buf[ibuf + 29]) + 6.0 *
1470 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1471 t4 = sqrt(t1 * t1 - t1 * t2 + t2 * t2 + 3.0 * t3 * t3);
1472 t1 = (buf[ibuf + 26] / buf[ibuf + 29]) - 6.0 *
1473 (buf[ibuf + 21] / (buf[ibuf + 29] * buf[ibuf + 29]));
1474 t2 = (buf[ibuf + 27] / buf[ibuf + 29]) - 6.0 *
1475 (buf[ibuf + 22] / (buf[ibuf + 29] * buf[ibuf + 29]));
1476 t3 = (buf[ibuf + 28] / buf[ibuf + 29]) - 6.0 *
1477 (buf[ibuf + 23] / (buf[ibuf + 29] * buf[ibuf + 29]));
1478 t5 = sqrt(t1 * t1 - t1 * t2 + t2 * t2 + 3.0 * t3 * t3);
1479 var[ivar] = MAX(t4, t5);
1480 ivar++;
1481 }
1482 break;
1483 }
1484 }
1485
1486 /*-------------------------------------------------------------------------
1487 * Function: taurus_readblockvar
1488 *
1489 * Purpose:
1490 *
1491 * Return: Success: 0
1492 *
1493 * Failure: -1
1494 *
1495 * Programmer:
1496 *
1497 * Modifications:
1498 * Eric Brugger, Fri Apr 28 09:13:43 PDT 1995
1499 * I removed some debugging code that I accidently left in previously.
1500 *
1501 * Jim Reus, 23 Apr 97
1502 * Changed to prototype form.
1503 *
1504 *-------------------------------------------------------------------------
1505 */
1506 static int
taurus_readblockvar(TAURUSfile * taurus,int var_id,int val_id,float * var)1507 taurus_readblockvar (TAURUSfile *taurus, int var_id, int val_id, float *var)
1508 {
1509 int size;
1510 int lbuf;
1511 float *buf;
1512 int ivar;
1513 int n;
1514 int ifile, iadd, nel, offset, ncomps;
1515
1516 /*
1517 * When reading displacements we are reading the coordinates in
1518 * the first file. This isn't very pretty but it fits the model.
1519 */
1520 if (var_id >= VAR_DISPX && var_id <= VAR_DISP_MAG) {
1521 ifile = 0;
1522 iadd = taurus->var_start[val_id];
1523 }
1524 else {
1525 ifile = taurus->state_file[taurus->state];
1526 iadd = taurus->state_loc[taurus->state] +
1527 taurus->var_start[val_id];
1528 }
1529
1530 nel = taurus->var_len[val_id];
1531 offset = taurus->var_offset[val_id];
1532 ncomps = taurus->var_ncomps[val_id];
1533
1534 /*
1535 * Allocate space for the buffer.
1536 */
1537 size = nel * ncomps;
1538
1539 if (size < MAXBUF)
1540 lbuf = size;
1541 else
1542 lbuf = (MAXBUF / ncomps) * ncomps;
1543
1544 buf = ALLOC_N(float, lbuf);
1545
1546 /*
1547 * Read a buffer full of data at a time, and transfer it to
1548 * the real variable using the appropriate stride.
1549 */
1550 ivar = 0;
1551 while (size > 0) {
1552 n = MIN(size, lbuf);
1553 taurus_read(taurus, ifile, iadd, n * sizeof(int), (char*)buf);
1554
1555 taurus_calc(taurus, buf, n, ncomps, offset, var_id, var, ivar);
1556 ivar += n / ncomps;
1557 iadd += n * sizeof(int);
1558
1559 size -= n;
1560 }
1561
1562 FREE(buf);
1563
1564 return (0);
1565 }
1566
1567 /*-------------------------------------------------------------------------
1568 * Function: db_taur_open
1569 *
1570 * Purpose: Open a taurus file.
1571 *
1572 * Return: Success: pointer to new file.
1573 *
1574 * Failure: NULL
1575 *
1576 * Programmer:
1577 *
1578 * Modifications:
1579 * Eric Brugger, Tue Mar 28 15:03:37 PST 1995
1580 * I modified the routines to read a topaz3d data file.
1581 *
1582 * Eric Brugger, Wed Apr 26 13:14:42 PDT 1995
1583 * I modified the routine to read the title in the taurus file
1584 * and put it in the taurus structure.
1585 *
1586 * Eric Brugger, Thu Jul 27 12:49:40 PDT 1995
1587 * I modified the routine to handle files generated by hydra.
1588 *
1589 * Jim Reus, 23 Apr 97
1590 * Changed to prototype form.
1591 *
1592 *-------------------------------------------------------------------------
1593 */
1594 TAURUSfile *
db_taur_open(char const * basename)1595 db_taur_open (char const *basename)
1596 {
1597 int fd;
1598 int loc, size;
1599 int ctl[40];
1600 char title[48];
1601 TAURUSfile *taurus;
1602
1603 /*
1604 * Create the structure to hold the open file information.
1605 */
1606 taurus = ALLOC_N(TAURUSfile, 1);
1607
1608 /*
1609 * Open the file and read the header.
1610 */
1611 if ((fd = open(basename, O_RDONLY)) < 0) {
1612 FREE(taurus);
1613 return (NULL);
1614 }
1615
1616 taurus->ifile = 0;
1617 taurus->fd = fd;
1618 taurus->basename = ALLOC_N(char, strlen(basename) + 1);
1619 strcpy(taurus->basename, basename);
1620 taurus->filename = ALLOC_N(char, strlen(basename) + 4);
1621
1622 taurus->mesh_read = 0;
1623
1624 loc = 15 * sizeof(int);
1625
1626 lseek(fd, loc, SEEK_SET);
1627 size = 40 * sizeof(int);
1628
1629 if (read(fd, ctl, size) != size) {
1630 FREE(taurus->basename);
1631 FREE(taurus->filename);
1632 FREE(taurus);
1633 close(fd);
1634 return (NULL);
1635 }
1636
1637 /*
1638 * Do a simple check to see that this is indeed a taurus file.
1639 * ctl [0] should really be a 4, which indicates a 3d mesh with
1640 * an unpacked node list. 3 indicates a 3d mesh with a packed
1641 * node list, but this can still have an unpacked node list, so
1642 * we will assume that it is unpacked regardless of the flag.
1643 */
1644 if (!((ctl[0] == 3 || ctl[0] == 4) &&
1645 (ctl[2] == 1 || ctl[2] == 2 || ctl[2] == 6 || ctl[2] == 200))) {
1646 FREE(taurus->basename);
1647 FREE(taurus->filename);
1648 FREE(taurus);
1649 close(fd);
1650 return (NULL);
1651 }
1652
1653 taurus->ndim = ctl[0];
1654 taurus->numnp = ctl[1];
1655 taurus->icode = ctl[2];
1656 taurus->nglbv = ctl[3];
1657 taurus->it = ctl[4];
1658 taurus->iu = ctl[5];
1659 taurus->iv = ctl[6];
1660 taurus->ia = ctl[7];
1661 taurus->nel8 = ctl[8];
1662 taurus->nummat8 = ctl[9];
1663 taurus->nv3d = ctl[12];
1664 taurus->nel2 = ctl[13];
1665 taurus->nummat2 = ctl[14];
1666 taurus->nv1d = ctl[15];
1667 taurus->nel4 = ctl[16];
1668 taurus->nummat4 = ctl[17];
1669 taurus->nv2d = ctl[18];
1670 taurus->activ = ctl[20];
1671
1672 /*
1673 * if ndim is 4 then the nodelist are unpacked. We cannot handle
1674 * packed values so let us assume no files exist that are packed.
1675 * If the values were packed they would presumably be three values
1676 * per integer, which means a very small number of nodes.
1677 */
1678 if (taurus->ndim == 4)
1679 taurus->ndim = 3;
1680
1681 /*
1682 * Read the title.
1683 */
1684 loc = 0;
1685 lseek(fd, loc, SEEK_SET);
1686 size = 40 * sizeof(char);
1687
1688 if (read(fd, title, size) != size) {
1689 FREE(taurus->basename);
1690 FREE(taurus->filename);
1691 FREE(taurus);
1692 close(fd);
1693 return (NULL);
1694 }
1695 title[40] = '\0';
1696 fix_title(title);
1697 strcpy(taurus->title, title);
1698
1699 /*
1700 * Initialize the file information.
1701 */
1702 init_file_info(taurus);
1703
1704 /*
1705 * Initialize the state information.
1706 */
1707 init_state_info(taurus);
1708
1709 /*
1710 * Initialize the variable information.
1711 */
1712 init_var_info(taurus);
1713
1714 /*
1715 * Initialize the material information.
1716 */
1717 init_mat_info(taurus);
1718
1719 return (taurus);
1720 }
1721
1722 /*-------------------------------------------------------------------------
1723 * Function: db_taur_close
1724 *
1725 * Purpose: Close a taurus file pointer and free the memory that
1726 * belongs to the taurus driver.
1727 *
1728 * Return: Success: 0
1729 *
1730 * Failure: -1
1731 *
1732 * Programmer: robb@cloud
1733 * Fri Dec 9 12:56:43 EST 1994
1734 *
1735 * Modifications:
1736 * Eric Brugger, Wed Dec 20 12:04:03 PST 1995
1737 * I modified the code to handle the activity data.
1738 *
1739 * Jim Reus, 23 Apr 97
1740 * Changed to prototype form.
1741 *
1742 *-------------------------------------------------------------------------
1743 */
1744 int
db_taur_close(TAURUSfile * taurus)1745 db_taur_close (TAURUSfile *taurus)
1746 {
1747
1748 close(taurus->fd);
1749 FREE(taurus->basename);
1750 FREE(taurus->filename);
1751 FREE(taurus->filesize);
1752
1753 FREE(taurus->state_file);
1754 FREE(taurus->state_loc);
1755 FREE(taurus->state_time);
1756
1757 FREE(taurus->hex_nodelist);
1758 FREE(taurus->shell_nodelist);
1759 FREE(taurus->beam_nodelist);
1760 FREE(taurus->hex_facelist);
1761 FREE(taurus->hex_zoneno);
1762 FREE(taurus->hex_matlist);
1763 FREE(taurus->shell_matlist);
1764 FREE(taurus->beam_matlist);
1765 FREE(taurus->hex_activ);
1766 FREE(taurus->shell_activ);
1767 FREE(taurus->beam_activ);
1768 if (taurus->coords != NULL) {
1769 FREE(taurus->coords[0]);
1770 FREE(taurus->coords[1]);
1771 if (taurus->ndim > 2)
1772 FREE(taurus->coords[2]);
1773 }
1774 FREE(taurus->coords);
1775 FREE(taurus);
1776 return (0);
1777 }
1778
1779 /*-------------------------------------------------------------------------
1780 * Function: init_mesh_info
1781 *
1782 * Purpose:
1783 *
1784 * Return: Success: void
1785 *
1786 * Failure:
1787 *
1788 * Programmer:
1789 *
1790 * Modifications:
1791 * Eric Brugger, Mon Aug 28 12:00:10 PDT 1995
1792 * I modified the routine to not use unreferenced nodes in the
1793 * mesh extent calculations.
1794 *
1795 * Eric Brugger, Wed Dec 20 12:04:03 PST 1995
1796 * I modified the code to handle the activity data.
1797 *
1798 * Jim Reus, 23 Apr 97
1799 * Changed to prototype form.
1800 *
1801 *-------------------------------------------------------------------------
1802 */
1803 void
init_mesh_info(TAURUSfile * taurus)1804 init_mesh_info (TAURUSfile *taurus)
1805 {
1806 int i;
1807 int iadd;
1808 int len;
1809 int lbuf;
1810 int *buf;
1811 float *rbuf;
1812 int *zones, *faces, *mats;
1813 int *zoneno;
1814 int nzones, nfaces;
1815 float minval, maxval;
1816 int ndim, numnp;
1817 float *coords;
1818 float xval, yval, zval;
1819
1820 if (taurus->mesh_read == 1)
1821 return;
1822
1823 taurus->nhex = 0;
1824 taurus->nshell = 0;
1825 taurus->nbeam = 0;
1826 taurus->coords = NULL;
1827 taurus->coord_state = -1;
1828
1829 /*
1830 * Read the coordinate information.
1831 */
1832 ndim = taurus->ndim;
1833 numnp = taurus->numnp;
1834
1835 /*
1836 * Allocate storage if not already allocated.
1837 */
1838 if (taurus->coords == NULL) {
1839 taurus->coords = ALLOC_N(float *, ndim);
1840 taurus->coords[0] = ALLOC_N(float, numnp);
1841 taurus->coords[1] = ALLOC_N(float, numnp);
1842
1843 if (ndim > 2)
1844 taurus->coords[2] = ALLOC_N(float, numnp);
1845 }
1846
1847 iadd = 64 * sizeof(int);
1848
1849 lbuf = numnp * ndim;
1850 rbuf = ALLOC_N(float, lbuf);
1851 len = lbuf * sizeof(float);
1852
1853 taurus_read(taurus, 0, iadd, len, (char*)rbuf);
1854
1855 coords = taurus->coords[0];
1856 for (i = 0; i < numnp; i++)
1857 coords[i] = rbuf[i * ndim];
1858
1859 coords = taurus->coords[1];
1860 for (i = 0; i < numnp; i++)
1861 coords[i] = rbuf[i * ndim + 1];
1862
1863 if (ndim > 2) {
1864 coords = taurus->coords[2];
1865 for (i = 0; i < numnp; i++)
1866 coords[i] = rbuf[i * ndim + 2];
1867 }
1868
1869 FREE(rbuf);
1870
1871 /*
1872 * Set up for reading the nodelists and material data.
1873 */
1874 iadd = 64 * sizeof(int) + taurus->numnp * taurus->ndim * sizeof(float);
1875
1876 lbuf = MAX(9 * taurus->nel8, MAX(5 * taurus->nel4, 6 * taurus->nel2));
1877 buf = ALLOC_N(int, lbuf);
1878
1879 /*
1880 * Read the hexahedron elements.
1881 */
1882 if (taurus->nel8 > 0) {
1883 zones = ALLOC_N(int, 8 * taurus->nel8);
1884 mats = ALLOC_N(int, taurus->nel8);
1885
1886 len = 9 * taurus->nel8 * sizeof(int);
1887
1888 taurus_read(taurus, 0, iadd, len, (char*)buf);
1889 for (i = 0; i < taurus->nel8; i++) {
1890 zones[i * 8] = buf[i * 9] - 1;
1891 zones[i * 8 + 1] = buf[i * 9 + 1] - 1;
1892 zones[i * 8 + 2] = buf[i * 9 + 2] - 1;
1893 zones[i * 8 + 3] = buf[i * 9 + 3] - 1;
1894 zones[i * 8 + 4] = buf[i * 9 + 4] - 1;
1895 zones[i * 8 + 5] = buf[i * 9 + 5] - 1;
1896 zones[i * 8 + 6] = buf[i * 9 + 6] - 1;
1897 zones[i * 8 + 7] = buf[i * 9 + 7] - 1;
1898 }
1899 for (i = 0; i < taurus->nel8; i++) {
1900 mats[i] = buf[i * 9 + 8];
1901 }
1902
1903 taurus->nhex = taurus->nel8;
1904 taurus->hex_nodelist = zones;
1905 taurus->hex_matlist = mats;
1906 taurus->hex_activ = NULL;
1907
1908 if (taurus->activ >= 1000 || taurus->activ <= 1005) {
1909 taurus->nhex_faces = 0;
1910 taurus->hex_facelist = NULL;
1911 taurus->hex_zoneno = NULL;
1912 }
1913 else {
1914 db_taur_extface(zones, taurus->numnp, taurus->nel8,
1915 mats, &faces, &nfaces, &zoneno);
1916
1917 taurus->nhex_faces = nfaces;
1918 taurus->hex_facelist = faces;
1919 taurus->hex_zoneno = zoneno;
1920 }
1921
1922 iadd += len;
1923 }
1924
1925 /*
1926 * Read the beam elements.
1927 */
1928 if (taurus->nel2 > 0) {
1929 zones = ALLOC_N(int, 2 * taurus->nel2);
1930 mats = ALLOC_N(int, taurus->nel2);
1931
1932 len = 6 * taurus->nel2 * sizeof(int);
1933
1934 taurus_read(taurus, 0, iadd, len, (char*)buf);
1935 for (i = 0; i < taurus->nel2; i++) {
1936 zones[i * 2] = buf[i * 6] - 1;
1937 zones[i * 2 + 1] = buf[i * 6 + 1] - 1;
1938 }
1939 for (i = 0; i < taurus->nel2; i++) {
1940 mats[i] = buf[i * 6 + 5];
1941 }
1942
1943 taurus->nbeam = taurus->nel2;
1944 taurus->beam_nodelist = zones;
1945 taurus->beam_matlist = mats;
1946 taurus->beam_activ = NULL;
1947
1948 iadd += len;
1949 }
1950
1951 /*
1952 * Read the shell elements.
1953 */
1954 if (taurus->nel4 > 0) {
1955 zones = ALLOC_N(int, 4 * taurus->nel4);
1956 mats = ALLOC_N(int, taurus->nel4);
1957
1958 len = 5 * taurus->nel4 * sizeof(int);
1959
1960 taurus_read(taurus, 0, iadd, len, (char*)buf);
1961 for (i = 0; i < taurus->nel4; i++) {
1962 zones[i * 4] = buf[i * 5] - 1;
1963 zones[i * 4 + 1] = buf[i * 5 + 1] - 1;
1964 zones[i * 4 + 2] = buf[i * 5 + 2] - 1;
1965 zones[i * 4 + 3] = buf[i * 5 + 3] - 1;
1966 }
1967 for (i = 0; i < taurus->nel4; i++) {
1968 mats[i] = buf[i * 5 + 4];
1969 }
1970
1971 taurus->nshell = taurus->nel4;
1972 taurus->shell_nodelist = zones;
1973 taurus->shell_matlist = mats;
1974 taurus->shell_activ = NULL;
1975
1976 iadd += len;
1977 }
1978
1979 /*
1980 * Find the unreferenced coordinates.
1981 */
1982 for (i = 0; i < numnp; i++)
1983 buf[i] = 0;
1984
1985 zones = taurus->hex_nodelist;
1986 nzones = taurus->nel8;
1987 for (i = 0; i < nzones * 8; i++)
1988 buf[zones[i]] = 1;
1989
1990 zones = taurus->beam_nodelist;
1991 nzones = taurus->nel2;
1992 for (i = 0; i < nzones * 2; i++)
1993 buf[zones[i]] = 1;
1994
1995 zones = taurus->shell_nodelist;
1996 nzones = taurus->nel4;
1997 for (i = 0; i < nzones * 4; i++)
1998 buf[zones[i]] = 1;
1999
2000 for (i = 0; i < numnp && buf[i] == 0; i++)
2001 /* do nothing */ ;
2002 if (i < numnp) {
2003 xval = taurus->coords[0][i];
2004 yval = taurus->coords[1][i];
2005 if (ndim > 2)
2006 zval = taurus->coords[2][i];
2007 }
2008 else {
2009 xval = 0.;
2010 yval = 0.;
2011 zval = 0.;
2012 }
2013
2014 coords = taurus->coords[0];
2015 for (i = 0; i < numnp; i++)
2016 if (buf[i] == 0)
2017 coords[i] = xval;
2018 coords = taurus->coords[1];
2019 for (i = 0; i < numnp; i++)
2020 if (buf[i] == 0)
2021 coords[i] = yval;
2022 if (ndim > 2) {
2023 coords = taurus->coords[2];
2024 for (i = 0; i < numnp; i++)
2025 if (buf[i] == 0)
2026 coords[i] = zval;
2027 }
2028
2029 /*
2030 * Determine the extents of the data. The extents will change
2031 * over the course of the problem but the extents at the beginning
2032 * will be used for all the states.
2033 */
2034 coords = taurus->coords[0];
2035 minval = coords[0];
2036 maxval = coords[0];
2037 for (i = 0; i < numnp; i++) {
2038 minval = MIN(minval, coords[i]);
2039 maxval = MAX(maxval, coords[i]);
2040 }
2041 taurus->min_extents[0] = minval;
2042 taurus->max_extents[0] = maxval;
2043
2044 coords = taurus->coords[1];
2045 minval = coords[0];
2046 maxval = coords[0];
2047 for (i = 0; i < numnp; i++) {
2048 minval = MIN(minval, coords[i]);
2049 maxval = MAX(maxval, coords[i]);
2050 }
2051 taurus->min_extents[1] = minval;
2052 taurus->max_extents[1] = maxval;
2053
2054 if (ndim > 2) {
2055 coords = taurus->coords[2];
2056 minval = coords[0];
2057 maxval = coords[0];
2058 for (i = 0; i < numnp; i++) {
2059 minval = MIN(minval, coords[i]);
2060 maxval = MAX(maxval, coords[i]);
2061 }
2062 taurus->min_extents[2] = minval;
2063 taurus->max_extents[2] = maxval;
2064 }
2065 else {
2066 taurus->min_extents[2] = 0.;
2067 taurus->max_extents[2] = 0.;
2068 }
2069
2070 FREE(buf);
2071
2072 taurus->mesh_read = 1;
2073 }
2074
2075 /*-------------------------------------------------------------------------
2076 * Function: init_coord_info
2077 *
2078 * Purpose:
2079 *
2080 * Return: Success:
2081 *
2082 * Failure:
2083 *
2084 * Programmer:
2085 *
2086 * Modifications:
2087 *
2088 * Jim Reus, 23 Apr 97
2089 * Changed to prototype form.
2090 *
2091 *-------------------------------------------------------------------------
2092 */
2093 void
init_coord_info(TAURUSfile * taurus)2094 init_coord_info (TAURUSfile *taurus)
2095 {
2096 int i;
2097 int ndim, numnp;
2098 int state_loc, state_file;
2099 int loc, len, lbuf;
2100 float *buf;
2101 float *coords;
2102
2103 ndim = taurus->ndim;
2104 numnp = taurus->numnp;
2105
2106 /*
2107 * Allocate storage if not already allocated.
2108 */
2109 if (taurus->coords == NULL) {
2110 taurus->coords = ALLOC_N(float *, ndim);
2111 taurus->coords[0] = ALLOC_N(float, numnp);
2112 taurus->coords[1] = ALLOC_N(float, numnp);
2113
2114 if (ndim > 2)
2115 taurus->coords[2] = ALLOC_N(float, numnp);
2116 }
2117
2118 /*
2119 * Set up the addresses of where to read the data from.
2120 */
2121 if (taurus->iu == 1) {
2122 state_loc = taurus->state_loc[taurus->state];
2123 state_file = taurus->state_file[taurus->state];
2124 loc = state_loc + (1 + taurus->nglbv) * sizeof(float);
2125 }
2126 else {
2127 state_loc = 64 * sizeof(int);
2128
2129 state_file = 0;
2130 loc = state_loc;
2131 }
2132
2133 lbuf = taurus->numnp * taurus->ndim;
2134 buf = ALLOC_N(float, lbuf);
2135
2136 len = lbuf * sizeof(float);
2137
2138 taurus_read(taurus, state_file, loc, len, (char*)buf);
2139
2140 coords = taurus->coords[0];
2141 for (i = 0; i < numnp; i++)
2142 coords[i] = buf[i * ndim];
2143
2144 coords = taurus->coords[1];
2145 for (i = 0; i < numnp; i++)
2146 coords[i] = buf[i * ndim + 1];
2147
2148 if (taurus->ndim > 2) {
2149 coords = taurus->coords[2];
2150 for (i = 0; i < numnp; i++)
2151 coords[i] = buf[i * ndim + 2];
2152 }
2153
2154 FREE(buf);
2155
2156 taurus->coord_state = taurus->state;
2157 }
2158
2159 /*-------------------------------------------------------------------------
2160 * Function: init_zone_info
2161 *
2162 * Purpose:
2163 *
2164 * Return: Success:
2165 *
2166 * Failure:
2167 *
2168 * Programmer: Eric Brugger
2169 * Date: December 20, 1995
2170 *
2171 * Modifications:
2172 *
2173 * Jim Reus, 23 Apr 97
2174 * Changed to prototype form.
2175 *
2176 *-------------------------------------------------------------------------
2177 */
2178 void
init_zone_info(TAURUSfile * taurus)2179 init_zone_info (TAURUSfile *taurus)
2180 {
2181 int i, j;
2182 int ifile, loc;
2183 int *zones, *mats;
2184 int *faces, nfaces, *zoneno;
2185
2186 /*
2187 * Check if any work needs to be done.
2188 */
2189 if (taurus->activ < 1000 && taurus->activ > 1005)
2190 return;
2191
2192 if (taurus->state < 0 && taurus->state >= taurus->nstates)
2193 return;
2194
2195 /*
2196 * Allocate storage if not already allocated.
2197 */
2198 if (taurus->hex_activ == NULL && taurus->nel8 > 0)
2199 taurus->hex_activ = ALLOC_N (int, taurus->nel8);
2200
2201 if (taurus->beam_activ == NULL && taurus->nel2 > 0)
2202 taurus->beam_activ = ALLOC_N (int, taurus->nel2);
2203
2204 if (taurus->shell_activ == NULL && taurus->nel4 > 0)
2205 taurus->shell_activ = ALLOC_N (int, taurus->nel4);
2206
2207 /*
2208 * Read the activity data from the file if it is present.
2209 */
2210 ifile = taurus->state_file [taurus->state];
2211 loc = taurus->state_loc [taurus->state];
2212 loc += (taurus->it * taurus->numnp +
2213 taurus->ndim * taurus->numnp *
2214 (taurus->iu + taurus->iv + taurus->ia) +
2215 taurus->nel8 * taurus->nv3d +
2216 taurus->nel4 * taurus->nv2d +
2217 taurus->nel2 * taurus->nv1d +
2218 taurus->nglbv + 1) * sizeof(int);
2219
2220 taurus_read (taurus, ifile, loc, sizeof (int)*taurus->nel8,
2221 (char*)(taurus->hex_activ));
2222 loc += sizeof (int) * taurus->nel8;
2223
2224 taurus_read (taurus, ifile, loc, sizeof (int)*taurus->nel2,
2225 (char*)(taurus->beam_activ));
2226 loc += sizeof (int) * taurus->nel2;
2227
2228 taurus_read (taurus, ifile, loc, sizeof (int)*taurus->nel4,
2229 (char*)(taurus->shell_activ));
2230 loc += sizeof (int) * taurus->nel4;
2231
2232 /*
2233 * Create the face list for the hex elements.
2234 */
2235 FREE (taurus->hex_facelist);
2236 FREE (taurus->hex_zoneno);
2237 zones = ALLOC_N (int, 8 * taurus->nel8);
2238 mats = ALLOC_N (int, taurus->nel8);
2239 j = 0;
2240 for (i = 0; i < taurus->nel8; i++) {
2241 if (taurus->hex_activ [i] != 0) {
2242 zones[j * 8] = taurus->hex_nodelist[i * 8];
2243 zones[j * 8 + 1] = taurus->hex_nodelist[i * 8 + 1];
2244 zones[j * 8 + 2] = taurus->hex_nodelist[i * 8 + 2];
2245 zones[j * 8 + 3] = taurus->hex_nodelist[i * 8 + 3];
2246 zones[j * 8 + 4] = taurus->hex_nodelist[i * 8 + 4];
2247 zones[j * 8 + 5] = taurus->hex_nodelist[i * 8 + 5];
2248 zones[j * 8 + 6] = taurus->hex_nodelist[i * 8 + 6];
2249 zones[j * 8 + 7] = taurus->hex_nodelist[i * 8 + 7];
2250 mats[j] = taurus->hex_matlist [i];
2251 j++;
2252 }
2253 }
2254 db_taur_extface(zones, taurus->numnp, j,
2255 mats, &faces, &nfaces, &zoneno);
2256 taurus->nhex_faces = nfaces;
2257 taurus->hex_facelist = faces;
2258 taurus->hex_zoneno = zoneno;
2259
2260 FREE (zones);
2261 FREE (mats);
2262 }
2263
2264 /*-------------------------------------------------------------------------
2265 * Function: taurus_readvar
2266 *
2267 * Purpose:
2268 *
2269 * Return: Success: 0
2270 *
2271 * Failure: -1
2272 *
2273 * Programmer:
2274 *
2275 * Modifications:
2276 * Eric Brugger, Wed Apr 26 15:21:29 PDT 1995
2277 * I modified the routine to return the title in the file
2278 * as _fileinfo.
2279 *
2280 * Eric Brugger, Thu Jul 27 12:49:40 PDT 1995
2281 * I modified the routine to handle files generated by hydra.
2282 *
2283 * Eric Brugger, Thu Dec 21 09:57:09 PST 1995
2284 * I modified the routine to return the meshname of the variable.
2285 *
2286 * Jim Reus, 23 Apr 97
2287 * Changed to prototype form.
2288 *
2289 *-------------------------------------------------------------------------
2290 */
2291 int
taurus_readvar(TAURUSfile * taurus,char const * varname,float ** var,int * length,int * center,char * meshname)2292 taurus_readvar (TAURUSfile *taurus, char const *varname, float **var, int *length,
2293 int *center, char *meshname)
2294 {
2295 int i;
2296 int idir;
2297 int ivar;
2298 int var_id, val_id;
2299
2300 if (taurus->icode == 1)
2301 idir = 8;
2302 else if (taurus->icode == 200)
2303 idir = 9;
2304 else
2305 idir = taurus->idir;
2306
2307 if (idir == -1)
2308 return (-1);
2309
2310 /*
2311 * Find the variable name in the variable list.
2312 */
2313 for (i = 0; taur_var_list[i].idir < idir; i++)
2314 /* do nothing */ ;
2315
2316 for (i = i; taur_var_list[i].idir == idir &&
2317 strcmp(taur_var_list[i].name, varname) != 0; i++)
2318 /* do nothing */ ;
2319
2320 if (taur_var_list[i].idir != idir)
2321 return (-1);
2322
2323 ivar = i;
2324 var_id = taur_var_list[ivar].ivar;
2325 val_id = taur_var_list[ivar].ival;
2326
2327 if (taurus->var_start[val_id] == -1)
2328 return (-1);
2329
2330 /*
2331 * Set the return values.
2332 */
2333 *center = taur_var_list[ivar].centering;
2334 if (var_id >= VAR_SIGX && var_id <= VAR_PRINC_STRESS_3) {
2335 *length = taurus->nel8 + taurus->nel4;
2336 }
2337 else {
2338 *length = taurus->var_len[val_id];
2339 }
2340 strcpy (meshname, taur_var_list[ivar].mesh);
2341
2342 /*
2343 * Allocate space for the variable.
2344 */
2345 *var = ALLOC_N(float, *length);
2346
2347 /*
2348 * Read the variable.
2349 */
2350 taurus_readblockvar(taurus, var_id, val_id, *var);
2351 if (var_id >= VAR_SIGX && var_id <= VAR_PRINC_STRESS_3) {
2352 val_id += (VAL_SHELL_MID_SIGX - VAL_HEX_SIGX);
2353 taurus_readblockvar(taurus, var_id, val_id,
2354 &(var[0][taurus->nel8]));
2355 }
2356
2357 return (0);
2358 }
2359