-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathpoly3d.c
3738 lines (3237 loc) · 111 KB
/
poly3d.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
//#define DEBUG
/*==================================================
SET TABSTOPS AT EVERY FOUR SPACES FOR PROPER DISPLAY
====================================================*/
/*****************************************************************************
* FILE: poly3d.c, version 0.0 (beta)
* DATE: June, 1993
* BY: Andrew L. Thomas
*
* A 3-D polygonal boundary element code based on superposition of angular
* dislocations (Yoffe, 1956; Comninou & Dunders, 1966). Translated and
* expanded from a triangular element FORTRAN code written by M. Jeyakumaran
* at Northwestern University (Jeyakumaran et al., 1992).
*
* See the "Poly3D Users Manual" for more information.
*****************************************************************************/
/**************************** Revision History *******************************
*
* If you change this code in ANY way, describe and initial those changes
* below. Put a comment line containing the word CHANGE (all uppercase)
* where the changes are made, so they will be easy to find in the source code.
*
* DATE FILE NAME INTLS
* -------- -------------- ---------------------- -----
* Jun-93 poly3d.c Andrew L. Thomas ALT
* Version 0.0 (beta) of Poly3D completed
* Nov-97 poly3d.c Yann Lagalaye
* Newpoly3d (execut) Correction of "shadow effect"
*****************************************************************************/
/******************************** Includes **********************************/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "matrix.h"
#include "safetan.h"
#include "getoptPoly3D.h"
#include "getwords.h"
#include "infcoeff.h"
#include "elastic.h"
#include "pi.h"
#include "nrutil.h"
#include "nr.h"
/******************************** Constants *********************************/
/*-------------
MISCELLANEOUS
--------------*/
#define MAXFILE 256 /* Max length of file names */
#define MAX_ERROR_MSG 256 /* Max length of error messages */
#define MAXLINE 256 /* Max line length for getwords() */
#define MAXWORDS 50 /* Max # words on getwords() line */
#define GLOBAL_NAME "global" /* Global coord system name */
#define ELT_CSYS_NAME "elocal" /* Element coord system name */
#define END_STMT "end" /* End flag for input file sections */
#define COMMENT_CHAR '*' /* Input file comment character */
#define CONTINUE_CHAR '\\' /* Input file line continuation char*/
#define ERROR -1 /* Return value for function errors */
#define FALSE 0 /* False flag */
#define TRUE 1 /* True flag */
#define BVECTOR_BC 0 /* Burgers vector BC component flag */
#define TRACTION_BC 1 /* Traction BC component flag */
/*-----------
PROGRAM INFO
-------------*/
#define PROGRAM "poly3d.c"
#define VERSION "Beta-Release"
#ifdef __DATE__
#define COMPILE_DATE __DATE__
#else
#define COMPILE_DATE "Date Unavailable"
#endif
/*--------------------------------------
NUMERICAL LIMITS USED WHEN CHECKING...
---------------------------------------*/
#define SWAP_TINY 1.0e-10 /* ...if vertices must be swapped */
#define TINY_ANGLE 0.5*PI/180. /* ...if elt coord sys can be calc */
#define BVERT_TINY 1.0e-10 /* ...if point lies below a vertex */
#define COPLANAR_LIMIT 30. /* ...if elt vertices are co-planar */
/*-------------------------------------
PRINT OPTION ARRAY SIZE AND POSITIONS
--------------------------------------*/
#define NUM_PR_OPTS 5
#define DISPL 0
#define STRAIN 1
#define PSTRAIN 2
#define STRESS 3
#define PSTRESS 4
/*--------------------------------------------------------------
CHARS USED IN INPUT FILE PRINT STRINGS TO ENABLE PRINT OPTIONS
---------------------------------------------------------------*/
#define DISPL_CHAR 'd'
#define STRESS_CHAR 's'
#define STRAIN_CHAR 'e'
#define TRACTION_CHAR 't'
#define BVECTOR_CHAR 'b'
#define PRINCIPAL_CHAR 'p'
/*----------------------------------------
INPUT FILE FORMAT FOR DEFINING CONSTANTS
-----------------------------------------*/
#define CONST_NAME_POS 0
#define CONST_VALUE_POS 2
#define CONST_NUM_PARAMS 3
/*------------------------------------------------------
INPUT FILE FORMAT FOR DEFINING USER COORDINATE SYSTEMS
-------------------------------------------------------*/
#define CS_NAME_POS 0
#define CS_PARENT_POS 1
#define CS_ORIGIN_POS 2
#define CS_ROT_POS 5
#define CS_ROT_ORDER_POS 8
#define CS_NUM_PARAMS 9
/*------------------------------------------------
INPUT FILE FORMAT FOR DEFINING OBSERVATION GRIDS
-------------------------------------------------*/
#define OG_NAME_POS 0
#define OG_DIMEN_POS 1
#define OG_PRINT_OPS_POS 2
#define OG_INPUT_CSYS_POS 3
#define OG_OBSPT_CSYS_POS 4
#define OG_DATA_CSYS_POS 5
#define OG_BEGIN_POS 6
#define OG_END_POS 9
#define OG_NUMPTS_POS 12
#define OG_MIN_NUM_PARAMS 9
/*---------------------------------------
INPUT FILE FORMAT FOR DEFINING VERTICES
----------------------------------------*/
#define V_CHAR 'v'
#define V_CHAR_POS 0
#define V_NAME_POS 1
#define V_CSYS_POS 2
#define V_X_POS 3
#define V_NUM_PARAMS 6
/*--------------------------------------
INPUT FILE FORMAT FOR DEFINING OBJECTS
---------------------------------------*/
#define OBJ_CHAR 'o'
#define OBJ_CHAR_POS 0
#define OBJ_NAME_POS 1
#define OBJ_PRINT_OPS_POS 2
#define OBJ_POS_CSYS_POS 3
#define OBJ_MIN_NUM_PARAMS 2
/*---------------------------------------
INPUT FILE FORMAT FOR DEFINING ELEMENTS
----------------------------------------*/
#define E_CHAR 'e'
#define E_CHAR_POS 0
#define E_NUM_VERT_POS 1
#define E_BC_CSYS_POS 2
#define E_BC_TYPE_POS 3
#define E_BC_POS 4
#define E_VERTEX_POS 7
#define E_MIN_NUM_PARAMS 10
/*------------------------------------
FORMAT FOR PRINTING ELEMENT GEOMETRY
-------------------------------------*/
#define ELT_GEOM_LABELS \
" ELT Vertex Name X1 X2 X3\n"
#define ELT_GEOM_UNDLNS \
"---- -------------------- ---------- ---------- ----------\n"
#define ELT_GEOM_FMT \
"%4d %-20s %10.4f %10.4f %10.4f\n"
/*-----------------------------------------
FORMAT FOR PRINTING OBSERVATION GRID DATA
------------------------------------------*/
#define OG_LOC_LABELS \
" X1 X2 X3 "
#define OG_LOC_UNDLNS \
"-------- -------- -------- "
#define OG_LOC_FMT \
"%8.3f %8.3f %8.3f "
#define OG_DISPL_TITLE \
"\nDISPLACEMENTS:\n\n"
#define OG_DISPL_LABELS \
" U1 U2 U3\n"
#define OG_DISPL_UNDLNS \
"---------- ---------- ----------\n"
#define OG_DISPL_FMT \
"%10.3e %10.3e %10.3e\n"
#define OG_STRAIN_TITLE \
"\nSTRAINS:\n\n"
#define OG_STRAIN_LABELS \
" E11 E22 E33 E12 E23 E13\n"
#define OG_STRAIN_UNDLNS \
"---------- ---------- ---------- ---------- ---------- ----------\n"
#define OG_STRAIN_FMT \
"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n"
#define OG_STRESS_TITLE \
"\nSTRESSES:\n\n"
#define OG_STRESS_LABELS \
" SIG11 SIG22 SIG33 SIG12 SIG23 SIG13\n"
#define OG_STRESS_UNDLNS \
"---------- ---------- ---------- ---------- ---------- ----------\n"
#define OG_STRESS_FMT \
"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n"
#define OG_PSTRAIN_TITLE \
"\nPRINCIPAL STRAINS:\n\n"
#define OG_PSTRAIN_LABELS \
" N1 N2 N3 E1 N1 N2 N3 E2 N1 N2 N3 E3\n"
#define OG_PSTRAIN_UNDLNS \
"------ ------ ------ ---------- ------ ------ ------ ---------- ------ ------ ------ ----------\n"
#define OG_PSTRAIN_FMT \
"%6.3f %6.3f %6.3f %10.3e %6.3f %6.3f %6.3f %10.3e %6.3f %6.3f %6.3f %10.3e\n"
#define OG_PSTRESS_TITLE \
"\nPRINCIPAL STRESSES:\n\n"
#define OG_PSTRESS_LABELS \
" N1 N2 N3 SIG1 N1 N2 N3 SIG2 N1 N2 N3 SIG3\n"
#define OG_PSTRESS_UNDLNS \
"------ ------ ------ ---------- ------ ------ ------ ---------- ------ ------ ------ ----------\n"
#define OG_PSTRESS_FMT \
"%6.3f %6.3f %6.3f %10.3e %6.3f %6.3f %6.3f %10.3e %6.3f %6.3f %6.3f %10.3e\n"
/*-------------------------------
FORMAT FOR PRINTING OBJECT DATA
--------------------------------*/
#define OBJ_LOC_LABELS \
" ELT X1C X2C X3C "
#define OBJ_LOC_UNDLNS \
"---- -------- -------- -------- "
#define OBJ_LOC_FMT \
"%4d %8.3f %8.3f %8.3f "
#define OBJ_DISPL_TITLE \
"\nDISPLACEMENTS:\n\n"
#define OBJ_DISPL_LABELS \
" B1 U1(+) U1(-) B2 U2(+) U2(-) B3 U3(+) U3(-) "
#define OBJ_DISPL_UNDLNS \
"---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- "
#define OBJ_DISPL_FMT \
"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e "
#define OBJ_STRESS_TITLE \
"\nSTRESSES (TRACTIONS):\n\n"
#define OBJ_STRESS_LABELS \
" T1 T2 T3 "
#define OBJ_STRESS_UNDLNS \
"---------- ---------- ---------- "
#define OBJ_STRESS_FMT \
"%10.3e %10.3e %10.3e "
#define OBJ_BC_CSYS_LABELS \
"Coord Sys\n"
#define OBJ_BC_CSYS_UNDLNS \
"---------\n"
#define OBJ_BC_CSYS_FMT \
"%s\n"
/********************************** Macros **********************************/
#define RADIANS(A) ((A)*PI/180.) /* Convert degrees to radians */
#define MAX(A,B) (((A) > (B)) ? (A):(B))/* MAX macro */
/******************************** Structures ********************************/
struct csys_s { /* -- COORDINATE SYSTEM STRUCTURE - */
char *name; /* Coordinate system name */
double origin[3]; /* Coord sys origin (global) */
double local_rot[3][3]; /* (To) global rotation matrix */
struct csys_s *next; /* Ptr to next c.s. in linked list */
};
typedef struct csys_s csys_t;
struct obs_grid_s { /* -- OBSERVATION GRID STRUCTURE -- */
char *name; /* Observation grid name */
int dimension; /* Dimension of observation grid */
double begin[3]; /* Obs grid beginning coords */
double end[3]; /* Obs grid ending coords */
int numpts[3]; /* No of obs pts along x1,x2,x3 */
int print[NUM_PR_OPTS]; /* Print options array */
csys_t *endpt_csys; /* Input coordinate system */
csys_t *obspt_csys; /* Observation grid coord system */
csys_t *outp_csys; /* Output coord sys for obs grid */
struct obs_grid_s *next; /* Ptr to next o.l. in linked list */
};
typedef struct obs_grid_s obs_grid_t;
struct vert_s { /* ------- VERTEX STRUCTURE ------- */
char *name; /* Vertex name */
double x[3]; /* Vertex coordinates (global) */
csys_t *csys; /* Coordinate system for vertex */
struct vert_s *next; /* Ptr to next vertex in linked list*/
};
typedef struct vert_s vert_t;
struct disloc_seg_s { /* --- DISLOC SEGMENT STRUCTURE --- */
double elt_b[3][3]; /* Proj of element b to segment b */
double trend; /* Strike of plunging leg of d.s. */
double plunge; /* Plunge of plunging leg of d.s. */
double local_rot[3][3]; /* Local-to-global rotation matrix */
vert_t *vert[2]; /* Dislocation segment vertices */
};
typedef struct disloc_seg_s disloc_seg_t;
struct elt_s { /* ------ ELEMENT STRUCTURE ------- */
int num_vertices; /* Number of vertices */
int bc_type[3]; /* Boundary condition type array */
double bc[3]; /* Boundary condition magnitudes */
csys_t elt_csys; /* Element-local coordinate system */
double *b[3]; /* Burgers vector array */
disloc_seg_t *disloc_seg; /* Dislocation segment array */
csys_t *bc_csys; /* Ptr to coord sys for element BCs */
struct elt_s *next; /* Ptr to next elt in linked list */
};
typedef struct elt_s elt_t;
struct obj_s { /* ------ OBJECT STRUCTURE -------- */
char *name; /* Object type */
int print[NUM_PR_OPTS]; /* Print options */
csys_t *pos_csys; /* Position coordinate system */
elt_t *first_elt; /* Pointer to first element */
elt_t *last_elt; /* Pointer to last element */
struct obj_s *next; /* Ptr to next obj in linked list */
};
typedef struct obj_s obj_t;
/*************************** External Variables *****************************/
char *title1_E = NULL; /* Problem title */
char *title2_E = NULL; /* Problem subtitle */
int half_space_E = TRUE; /* Half/whole space flag */
int check_cond_num_E = TRUE; /* Check matrix condition num flag */
double cond_num_E = -1.0; /* Matrix condition number */
char infile_E[MAXFILE]; /* Input file name */
char outfile_E[MAXFILE]; /* Output file name */
int linenum_E = 0; /* Current line # in input file */
int num_elts_E = 0; /* Number of elements */
int below_vertex_E = FALSE; /* Below vertex flag */
double null_value_E = -999.0; /* Null output value */
/*********************************************************************************************************/
/************************* NEW: added 98-12-09 to reflect the problem of observation point near vertices */
double coef_exclu_E = 0.0; /* coef exclusion value */
int near_vertex_E = FALSE; /* nearnest vertex flag for observation point */
/* EXPLANATIONS:
Let's obs_pt be an observation point.
Put flag near_vertex_E to FALSE.
For each vertex in current project:
Let's d be the mean length of segments containing v
Let's l be the distance from v to a choosen observation point
Then if l<d*coef_exclu => near_vertex_E = TRUE
Then, before printing computed values for this observation point:
if near_vertex_E=TRUE => print null_value_E
otherwise print computed values
So coef_exclu_E = 1.0, means 100% of the mean length of segments containing v.
*/
/*********************************************************************************************************/
FILE *ifp_E = stdin; /* Input file ptr (default stdin) */
FILE *ofp_E = stdout; /* Output file ptr (default stdout) */
FILE *tempfp_E[NUM_PR_OPTS]; /* Temporary file ptrs */
double shear_mod_E = -1.0; /* Shear modulus */
double psn_ratio_E = -1.0; /* Poisson's ratio */
double youngs_mod_E = -1.0; /* Young's modulus */
double bulk_mod_E = -1.0; /* Bulk modulus */
double lame_lambda_E = -1.0; /* Lame's lambda */
int print_elt_geom_E = FALSE; /* Print element geometry flag */
char *elt_geom_csys_name_E =NULL;/* Element geometry coord sys name */
int rem_stress_bc_E = TRUE; /* Remote stress vs strain bc flag */
double rem_stress_E[3][3]; /* Remote stress tensor */
double rem_strain_E[3][3]; /* Remote strain tensor */
double *b_vector_E; /* Burger's vector array */
double **ic_matrix_E; /* Influence coeff matrix */
csys_t *first_csys_E = NULL; /* 1st memb of csys linked list */
obs_grid_t *first_obs_grid_E = NULL; /* 1st memb of obs grid linked list */
obj_t *first_obj_E = NULL; /* 1st memb of obj linked list */
elt_t *first_elt_E = NULL; /* 1st memb of elt linked list */
vert_t *first_vert_E = NULL; /* 1st memb of vert linked list */
/************************* ANSI Function Declarations ***********************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
double array_max_norm(double **a, int start_row, int end_row, int start_col,
int end_col);
int calc_elt_parameters(elt_t *current_elt);
void close_temp_files(void);
void copy_temp_files(void);
void determine_burgers_vectors(void);
void displ_strain(int calc_displ, int calc_strain, double x[3],
double displ[3], double strain[3][3], elt_t *omit_elt);
void displ_strain_poly_elt(int calc_displ, int calc_strain,
elt_t *current_elt, double x[3], double displ[3],
double strain[3][3], elt_t *omit_elt, int under);
void print_obj_data(void);
csys_t *find_csys(char *name);
vert_t *find_vert(char *name);
int get_double_var(double *var, char *var_name, char *word[],
int numwords);
int get_boolean_var(int *var, char *var_name, char *true_string,
char *false_string, char *word[], int numwords, char *line);
int get_text_var(char **var, char *var_name, char *word[], int numwords,
char *line);
void get_elt_info(elt_t **current_elt, obj_t *current_obj, int numwords,
char *word[], char *line);
void get_obj_info(obj_t **current_obj, int numwords, char *word[],
char *line);
void get_program_args(void);
void get_vert_info(vert_t **current_vert, int numwords, char *word[],
char *line);
void displ_strain_ics_poly_elt(int calc_displ, int calc_strain,
elt_t *current_elt, double x[3],
double displ_ic[3][3], double strain_ic[3][3][3],
elt_t *omit_elt);
void print_obs_grid_data(void);
int open_files();
void open_temp_files(int print[]);
int parse_command_line_args(int argc, char *argv[]);
void p_error(char *error_msg, char *line);
void display_msg(char *_msg);
void print_elt_data(obj_t *current_obj,
elt_t *current_elt, int elt_num, double displ[3],
double stress[3][3]);
void print_elt_geometry(void);
void print_obj_titles(obj_t *current_obj);
void print_obs_grid_titles(obs_grid_t *current_obs_grid);
void print_obs_pt_data(obs_grid_t *current_obs_grid, double x[3],
double displ[3], double strain[3][3]);
void print_problem_info(void);
int read_csystems();
int read_objs_elts_verts();
void read_infile(void);
int read_line(char *line, char *word[]);
int read_observation_grids();
int read_constants();
void setup_global_coords(void);
/************************** K&R Function Declarations ***********************/
#else
double array_max_norm();
int calc_elt_parameters();
void close_temp_files();
void copy_temp_files();
void determine_burgers_vectors();
void displ_strain();
void displ_strain_poly_elt();
void print_obj_data();
csys_t *find_csys();
vert_t *find_vert();
int get_double_var();
int get_boolean_var();
int get_text_var();
void get_elt_info();
void get_obj_info();
void get_program_args();
void get_vert_info();
void displ_strain_ics_poly_elt();
void print_obs_grid_data();
int open_files();
void open_temp_files();
int parse_command_line_args();
void p_error();
void display_msg();
void print_elt_data();
void print_elt_geometry();
void print_obj_titles();
void print_obs_grid_titles();
void print_obs_pt_data();
void print_problem_info();
int read_csystems();
int read_objs_elts_verts();
void read_infile();
int read_line();
int read_observation_grids();
int read_constants();
void setup_global_coords();
#endif
/******************************* Function: main ******************************
* In: argc - number of command line arguments
* argv - array of command line arguments
*****************************************************************************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
main(int argc, char *argv[])
#else
main(argc, argv)
int argc;
char *argv[];
#endif
{
time_t start, t_reading, t_matrix, t_problem, t_object, t_obs,finish;
double elapsed_time;
/* Use get_program_args() or parse_command_line_args() to get
file names and program options.
-------------------------------------------------------------*/
#if FPROMPT
get_program_args();
#else
parse_command_line_args(argc,argv);
#endif
/* Start counter */
time( &start );
/* Open the input and output files
----------------------------------*/
open_files();
/* Read input file and set up the problem
-----------------------------------------*/
read_infile();
time( &t_reading );
elapsed_time = difftime( t_reading, start );
printf( "\nReading file : %6.0f seconds.", elapsed_time );
/* Solve for burger's vector for each element
---------------------------------------------*/
determine_burgers_vectors();
time( &t_matrix );
elapsed_time = difftime( t_matrix,t_reading );
printf( "\nEquations/matrix : %6.0f seconds.", elapsed_time );
/* Print the problem info
-------------------------*/
print_problem_info();
time( &t_problem );
elapsed_time = difftime( t_problem,t_matrix );
printf( "\nPrint pb info : %6.0f seconds.", elapsed_time );
/* Calculate displs and tractions on elements
---------------------------------------------*/
print_obj_data();
time( &t_object );
elapsed_time = difftime( t_object,t_problem );
printf( "\nD and T on elements: %6.0f seconds.", elapsed_time );
/* Calculate displacements and stresses along observation grids
---------------------------------------------------------------*/
print_obs_grid_data();
time( &t_obs );
elapsed_time = difftime( t_obs,t_object );
printf( "\nD and S on obs : %6.0f seconds.", elapsed_time );
/* End counter, and print elapsed time */
time( &finish );
elapsed_time = difftime( finish, start );
printf( "\nProgram takes %6.0f seconds.\n", elapsed_time );
return(0);
}
/*********************** Function: print_problem_info ***********************
* Prints general problem information to the output file.
****************************************************************************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
void print_problem_info(void)
#else
void print_problem_info()
#endif
{
/* Print program name, version, and date
----------------------------------------*/
fprintf(ofp_E,"OUTPUT FROM: %s, version %s\n",PROGRAM, VERSION);
fprintf(ofp_E," COMPILED: %s\n",COMPILE_DATE);
/* Print input file name and problem titles
-------------------------------------------*/
fprintf(ofp_E,"\n INPUT FILE: %s\n",infile_E);
fprintf(ofp_E, " TITLE1: %s\n",title1_E);
fprintf(ofp_E, " TITLE2: %s\n",title2_E);
/* Print elastic constant values
--------------------------------*/
fprintf(ofp_E,"\nELASTIC CONSTANTS:\n");
fprintf(ofp_E, " Shear Modulus = %f\n",shear_mod_E);
fprintf(ofp_E, " Poisson's Ratio = %f\n",psn_ratio_E);
fprintf(ofp_E, " Young's Modulus = %f\n",youngs_mod_E);
fprintf(ofp_E, " Bulk Modulus = %f\n",bulk_mod_E);
fprintf(ofp_E, " Lame's Lambda = %f\n",lame_lambda_E);
/* Print the null output value
------------------------------*/
fprintf(ofp_E,"\nNULL OUPUT VALUE = %f\n",null_value_E);
/* Print the coef exclusion value
------------------------------*/
fprintf(ofp_E,"\nCOEF EXCLUSION VALUE = %f\n",coef_exclu_E);
/* Print condition number of the influence coefficient matrix
-------------------------------------------------------------*/
fprintf(ofp_E,"\nCONDITION NUMBER = ");
if (cond_num_E < 0.0) {
fprintf(ofp_E,"(no traction bc's -> no matrix needed)\n");
} else {
if (check_cond_num_E)
fprintf(ofp_E,"%f\n",cond_num_E);
else
fprintf(ofp_E,"(not requested)\n");
}
/* Print element geometries (if requested)
------------------------------------------*/
if (print_elt_geom_E)
print_elt_geometry();
}
/************************** Function: print_obj_data ************************
* Calculates and prints object data to the output file.
*****************************************************************************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
void print_obj_data(void)
#else
void print_obj_data()
#endif
{
int elt_num;
int calc_displ;
int calc_strain;
obj_t *current_obj;
double displ[3];
double strain[3][3];
double stress[3][3];
elt_t *current_elt;
current_obj = first_obj_E;
current_elt = first_elt_E;
while (current_elt != NULL) {
/* If this element starts new object......
-----------------------------------------*/
if (current_elt == current_obj->first_elt) {
/* Reset element number
----------------------*/
elt_num = 1;
if (current_obj->print[DISPL] || current_obj->print[STRESS]) {
/* Open temp files
------------------*/
open_temp_files(current_obj->print);
/* Print object titles
----------------------*/
print_obj_titles(current_obj);
} else {
/* Skip this object if no output requested
------------------------------------------*/
current_elt = current_obj->last_elt;
}
}
if (current_obj->print[DISPL] || current_obj->print[STRESS]) {
calc_displ = current_obj->print[DISPL];
calc_strain = current_obj->print[STRESS];
displ_strain(calc_displ,calc_strain,
current_elt->elt_csys.origin,displ,strain,
current_elt);
strain_to_stress(strain,shear_mod_E,lame_lambda_E,
stress);
// Print the data for this elt
print_elt_data(current_obj, current_elt, elt_num,
displ, stress);
elt_num++;
}
if (current_elt == current_obj->last_elt) {
if (current_obj->print[DISPL] || current_obj->print[STRESS]) {
/* Copy temp files to main output file, then close them.
--------------------------------------------------------*/
copy_temp_files();
close_temp_files();
}
current_obj = current_obj->next;
}
current_elt = current_elt->next;
}
}
/********************** Function: print_obs_grid_data ******************
* Loops through linked list of observation grids, caculating and printing
* the requested displacement, strain, and stress data for each to the
* output file.
************************************************************************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
void print_obs_grid_data(void)
#else
void print_obs_grid_data()
#endif
{
obs_grid_t *current_obs_grid;
int i, j, k;
double x[3];
double dx[3];
double displ[3];
double strain[3][3];
double *begin;
double *end;
int *numpts;
char error_msg[MAX_ERROR_MSG];
int *print;
int calc_displ;
int calc_strain;
/* Declaration for the bug correction */
double x_copy[3];
/* Loop over each observation grid
----------------------------------*/
current_obs_grid = first_obs_grid_E;
while (current_obs_grid != NULL) {
/* Determine data required by obs grid print options
----------------------------------------------------*/
print = current_obs_grid->print;
calc_displ = print[DISPL];
calc_strain = (print[STRAIN] || print[PSTRAIN] || print[STRESS]
|| print[PSTRESS]);
/* Remenber that begin_pt & end_pt are in GLOBAL coordinate system.*/
begin = current_obs_grid->begin;
end = current_obs_grid->end;
numpts = current_obs_grid->numpts;
subtract_vectors(end,begin,dx);
/* Open temp files
------------------*/
open_temp_files(current_obs_grid->print);
/* Print observation grid titles
--------------------------------*/
print_obs_grid_titles(current_obs_grid);
/* Process the observation grid
-------------------------------*/
switch (current_obs_grid->dimension) {
case 0:
copy_vector(begin,x);
/* Transform current observation point into global CSys
----------------------------------------------------*/
transform_position_vector(INVERSE_TRANS,
current_obs_grid->endpt_csys->origin,
current_obs_grid->endpt_csys->local_rot,
x);
displ_strain(calc_displ,calc_strain,x,displ,strain,
NULL);
print_obs_pt_data(current_obs_grid,x,displ,strain);
break;
case 1:
dx[2] /= (numpts[0]-1);
dx[1] /= (numpts[0]-1);
dx[0] /= (numpts[0]-1);
for (i=0; i < numpts[0]; i++) {
x[2] = begin[2] + dx[2]*i;
x[1] = begin[1] + dx[1]*i;
x[0] = begin[0] + dx[0]*i;
/* Transform current observation point into global CSys
----------------------------------------------------*/
transform_position_vector(INVERSE_TRANS,
current_obs_grid->endpt_csys->origin,
current_obs_grid->endpt_csys->local_rot,
x);
displ_strain(calc_displ,calc_strain,x,displ,
strain,NULL);
print_obs_pt_data(current_obs_grid,x,displ,strain);
}
break;
case 2:
case 3:
for (i=0; i < 3; i++) {
dx[i] /= (numpts[i] == 1) ? 1 : (numpts[i]-1);
}
for (i=0; i < numpts[2]; i++) {
x[2] = begin[2] + dx[2]*i;
for (j=0; j < numpts[1]; j++) {
x[1] = begin[1] + dx[1]*j;
for (k=0; k < numpts[0]; k++) {
x[0] = begin[0] + dx[0]*k;
copy_vector(x,x_copy);
/* Transform current observation point into global CSys
----------------------------------------------------*/
transform_position_vector(INVERSE_TRANS,
current_obs_grid->endpt_csys->origin,
current_obs_grid->endpt_csys->local_rot,
x_copy);
displ_strain(calc_displ,calc_strain,/****/x_copy/****/,displ,strain,NULL);
print_obs_pt_data(current_obs_grid, /****/x_copy/****/,displ,strain);
}
}
}
break;
default:
sprintf(error_msg,
"Invalid dimension (%d) for observation grid",
current_obs_grid->dimension);
p_error(error_msg,NULL);
}
/* Copy temp files to main output file, then close them.
--------------------------------------------------------*/
copy_temp_files();
close_temp_files();
current_obs_grid = current_obs_grid->next;
} /*while*/
}
/************************ Function: displ_strain ****************************
* Calculates the total displacement and/or strain at a point due to ALL
* elements.
*
* In: calc_displ - calculate displacements flag
* calc_strain - calculate strains flag
* x - coords (global) of pt at which to calc displ & strain
* omit_elt - element to omit when calculating displs (NULL = none)
*
* Out: displ - displacement vector (global coords)
* strain - strain tensor (global coords)
*
*****************************************************************************/
#if defined(__STDC__) || defined(ANSI) /* ANSI */
void displ_strain(int calc_displ, int calc_strain, double x[3],double displ[3], double strain[3][3], elt_t *omit_elt)
#else
void displ_strain(calc_displ, calc_strain, x, displ, strain, omit_elt)
int calc_displ;
int calc_strain;
double x[3];
double displ[3];
double strain[3][3];
elt_t *omit_elt;
#endif
{
elt_t *current_elt;
double elt_displ[3];
double elt_strain[3][3];
/*Declarations for correction of the "shadow effect" */
int i;
double orient;
double under_plane;
int under;
double normal [3];
double data1 [3];
double data [3];
vert_t *verta;
vert_t *vertb;
vert_t *vertc;
double seg1 [3];
double seg2 [3];
double inside;
double inside_vector [3];
int inside_test;
int inside_test_fin;
double x3_global [3];
disloc_seg_t *disloc_seg;
/* end of the declaration */
/*Declarations for coef_exclu */
double dist;
/* end of the declaration */
initialize_vector(displ,0.0);
initialize_matrix(strain,0.0);
/* Loop over each element
-------------------------*/
current_elt = first_elt_E;
near_vertex_E = FALSE;
while (current_elt != NULL)
{
/* Test UNDER to determine if the data point is inside the "shadow zone"
---------------------------------------------------------------------*/
under = 0;
inside_test_fin = 1;
/* Determine if the element has a positive side up or down (orient)
and determine if the data is under the plane
defined by the element (under_plane)
-----------------------------------------------------------------*/
disloc_seg = current_elt->disloc_seg;
x3_global [0] = 0;
x3_global [1] = 0;
x3_global [2] = 1;
verta = disloc_seg [0].vert[0];
vertb = disloc_seg [0].vert[1];
vertc = disloc_seg [1].vert[1];
subtract_vectors (vertb->x, verta->x, seg1);
subtract_vectors (vertc->x, verta->x, seg2);
cross_product (seg1, seg2, normal);
orient = dot_product (normal, x3_global);
subtract_vectors (x, verta->x, data1);
under_plane = dot_product (normal, data1);
/* Determine if obs_pt is nearnest vertex verta, vertb or vertc
1) compute d = mean length of the 3 disloc_seg
2) compute distance from obs_pt to verta and see if this diance is < d*coef_exclu
3) compute distance from obs_pt to vertb and see if this diance is < d*coef_exclu
4) compute distance from obs_pt to vertc and see if this diance is < d*coef_exclu
if condition is TRUE, break the loop over elements and set the flag near_vertex_E=TRUE
and return.
We use a new function "distance" define in matrix.h & matrix.c
*/
dist = 0.0;
dist = distance(verta->x,vertb->x);
dist += distance(vertb->x,vertc->x);
dist += distance(vertc->x,verta->x);
dist *= coef_exclu_E/3.0;
if (distance(x,verta->x)<=dist || distance(x,vertb->x)<=dist || distance(x,vertc->x)<=dist)
{
near_vertex_E = TRUE;
break;
}
/* Determine if the data point is in the "rigid body"
(cf. explanation of the bug)
---------------------------------------------------*/
for (i = 0; i < current_elt->num_vertices; i++)
{
inside = 0;
verta = disloc_seg [i].vert[0];
vertb = disloc_seg [i].vert[1];
subtract_vectors (vertb->x, verta->x, seg1);
subtract_vectors (x, verta->x, data);
cross_product (seg1, data, inside_vector);
inside = dot_product (x3_global, inside_vector);
if (orient > 0)
{
if (inside > 0 && under_plane < 0)
inside_test = 1;
else
inside_test = 0;
}
if (orient < 0)
{
if (inside < 0 && under_plane > 0)
inside_test = 1;
else
inside_test = 0;
}
if (orient == 0)
under = 0;
inside_test_fin *= inside_test;
}
/*Gives a value 1 to under for data points under the element
and positive side up
Gives a value 2 to under for data points under the element
and positive side down
Gives a value 0 to under for data points not under the element
---------------------------------------------------------------*/
if (inside_test_fin > 0)