1 /*
2 Copyright (C) 1998 Dennis Roddeman
3 email: dennis.roddeman@feat.nl
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software Foundation
17 59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
18 */
19
20 #include "tochnog.h"
21
element_middle_radius_set()22 void element_middle_radius_set( )
23
24 {
25 long int element=0, max_element=0, idim=0, inol=0, inod=0,
26 nnol=0, length=0, ldum=0, idum[1], el[1+MNOL], nodes[MNOL];
27 double distance=0., radius=0., middle[MDIM], coord[MDIM],
28 work[MDIM], ddum[1], *node_dof=NULL;
29
30 db_max_index( ELEMENT, max_element, VERSION_NORMAL, GET );
31 for ( element=0; element<=max_element; element++ ) {
32 if ( db_active_index( ELEMENT, element, VERSION_NORMAL ) ) {
33 db( ELEMENT, element, el, ddum, length, VERSION_NORMAL, GET );
34 nnol = length - 1; array_move( &el[1], nodes, nnol );
35 array_set( middle, 0., MDIM );
36 for ( inol=0; inol<nnol; inol++ ) {
37 inod = nodes[inol];
38 db( NODE, inod, idum, coord, ldum, VERSION_NORMAL, GET );
39 if ( materi_displacement ) {
40 node_dof = db_dbl( NODE_DOF, inod, VERSION_NORMAL );
41 for ( idim=0; idim<ndim; idim++ ) {
42 coord[idim] += node_dof[dis_indx+idim*nder];
43 }
44 }
45 array_add( coord, middle, middle, ndim );
46 }
47 array_multiply( middle, middle, (double) 1/nnol, ndim );
48 db( ELEMENT_MIDDLE, element, idum, middle, ndim, VERSION_NORMAL, PUT );
49 radius = 0.;
50 for ( inol=0; inol<nnol; inol++ ) {
51 inod = nodes[inol];
52 db( NODE, inod, idum, coord, ldum, VERSION_NORMAL, GET );
53 distance = array_distance( coord, middle, work, ndim );
54 if ( distance>radius ) radius = distance;
55 }
56 length = 1;
57 db( ELEMENT_RADIUS, element, idum, &radius, length, VERSION_NORMAL, PUT );
58 }
59 }
60 }
61
element_volume_set(long int name,long int nodes[],long int version,double & element_volume)62 void element_volume_set( long int name, long int nodes[], long int version,
63 double &element_volume )
64
65 {
66
67 long int inod=0, jnod=0, knod=0, lnod=0, ldum=0, idum[1];
68 double coord0[MDIM], coord1[MDIM], coord2[MDIM], coord3[MDIM], work[MDIM];
69
70 if ( ndim==1 ) {
71 inod = nodes[0];
72 jnod = nodes[1];
73 db( NODE, inod, idum, coord0, ldum, version, GET );
74 db( NODE, jnod, idum, coord1, ldum, version, GET );
75 element_volume = array_distance( coord0, coord1, work, ndim );
76 }
77 else if ( ndim==2 ) {
78 if ( name==-TRIA3 ) {
79 inod = nodes[0];
80 jnod = nodes[1];
81 knod = nodes[2];
82 }
83 else {
84 assert( name==-TRIA6 );
85 inod = nodes[0];
86 jnod = nodes[2];
87 knod = nodes[5];
88 }
89 db( NODE, inod, idum, coord0, ldum, version, GET );
90 db( NODE, jnod, idum, coord1, ldum, version, GET );
91 db( NODE, knod, idum, coord2, ldum, version, GET );
92 element_volume = triangle_area( coord0, coord1, coord2 );
93 }
94 else {
95 assert( ndim==3 );
96 if ( name==-TET4 ) {
97 inod = nodes[0];
98 jnod = nodes[1];
99 knod = nodes[2];
100 lnod = nodes[3];
101 }
102 else {
103 assert( name==-TET10 );
104 inod = nodes[0];
105 jnod = nodes[2];
106 knod = nodes[5];
107 lnod = nodes[9];
108 }
109 db( NODE, inod, idum, coord0, ldum, version, GET );
110 db( NODE, jnod, idum, coord1, ldum, version, GET );
111 db( NODE, knod, idum, coord2, ldum, version, GET );
112 db( NODE, lnod, idum, coord3, ldum, version, GET );
113 element_volume = tetrahedron_volume( coord0, coord1, coord2, coord3 );
114 }
115
116 }
117
get_element_matrix_unknowns(long int element,long int element_matrix_unknowns[])118 void get_element_matrix_unknowns( long int element,
119 long int element_matrix_unknowns[] )
120
121 {
122 long int i=0, n=0, inol=0, jnol=0, inod=0, jnod=0,
123 ipuknwn=0, jpuknwn=0, length=0, nnol=0, element_group=0, ldum=0,
124 el[1+MNOL], nodes[MNOL], *group_matrix_unknowns=NULL;
125 double ddum[1];
126
127 if ( db_active_index( ELEMENT_MATRIX_UNKNOWNS, element, VERSION_NORMAL ) )
128 db( ELEMENT_MATRIX_UNKNOWNS, element, element_matrix_unknowns, ddum, ldum,
129 VERSION_NORMAL, GET );
130 else {
131 db( ELEMENT, element, el, ddum, length, VERSION_NORMAL, GET );
132 nnol = length - 1; array_move( &el[1], nodes, nnol );
133 db( ELEMENT_GROUP, element, &element_group, ddum, ldum,
134 VERSION_NORMAL, GET_IF_EXISTS );
135 group_matrix_unknowns =
136 db_int( GROUP_MATRIX_UNKNOWNS, element_group, VERSION_NORMAL );
137 length = db_len( GROUP_MATRIX_UNKNOWNS, element_group, VERSION_NORMAL );
138 n = length / 4;
139 for ( i=0; i<n; i++ ) {
140 inol = group_matrix_unknowns[4*i+0]; inod = nodes[inol];
141 ipuknwn = group_matrix_unknowns[4*i+1];
142 jnol = group_matrix_unknowns[4*i+2]; jnod = nodes[jnol];
143 jpuknwn = group_matrix_unknowns[4*i+3];
144 element_matrix_unknowns[2*i+0] = inod*npuknwn + ipuknwn;
145 element_matrix_unknowns[2*i+1] = jnod*npuknwn + jpuknwn;
146 }
147 }
148
149 }
150
get_new_char(long int n)151 char *get_new_char( long int n )
152
153 {
154 char *ptr=NULL;
155
156 if ( n<=0 ) n = 1;
157 if ( !(ptr = new char[n] ) ) {
158 pri( "Error: cannot allocate enough memory." );
159 exit(TN_EXIT_STATUS);
160 }
161 return ptr;
162 }
163
get_new_dbl(long int n)164 double *get_new_dbl( long int n )
165
166 {
167 double *ptr=NULL;
168
169 if ( n<=0 ) n = 1;
170 if ( !(ptr = new double[n] ) ) {
171 pri( "Error: cannot allocate enough memory." );
172 exit(TN_EXIT_STATUS);
173 }
174 return ptr;
175 }
176
get_new_int(long int n)177 long int *get_new_int( long int n )
178
179 {
180 long int *ptr=NULL;
181
182 if ( n<=0 ) n = 1;
183 if ( !(ptr = new long int[n] ) ) {
184 pri( "Error: cannot allocate enough memory." );
185 exit(TN_EXIT_STATUS);
186 }
187 return ptr;
188 }
189
get_new_int_short(long int n)190 int *get_new_int_short( long int n )
191
192 {
193 int *ptr=NULL;
194
195 if ( n<=0 ) n = 1;
196 if ( !(ptr = new int[n] ) ) {
197 pri( "Error: cannot allocate enough memory." );
198 exit(TN_EXIT_STATUS);
199 }
200 return ptr;
201 }
202
203
set_environment(void)204 void set_environment( void )
205
206 {
207 long int length=0, options_processors=1,
208 options_solver=-MATRIX_ITERATIVE_BICG;
209 double ddum[1];
210 char *str=NULL;
211
212 if ( !db_active_index( OPTIONS_PROCESSORS, 0, VERSION_NORMAL ) ) {
213 str = getenv("TOCHNOG_OPTIONS_PROCESSORS");
214 if ( str!=NULL ) {
215 options_processors = atoi( str );
216 length = 1;
217 db( OPTIONS_PROCESSORS, 0, &options_processors, ddum,
218 length, VERSION_NORMAL, PUT );
219 }
220 }
221 if ( !db_active_index( OPTIONS_SOLVER, 0, VERSION_NORMAL ) ) {
222 str = getenv("TOCHNOG_OPTIONS_SOLVER");
223 if ( str!=NULL ) {
224 options_solver = -db_number( str );
225 length = 1;
226 db( OPTIONS_SOLVER, 0, &options_solver, ddum,
227 length, VERSION_NORMAL, PUT );
228 }
229 }
230
231 }
232
set_swit(long int element,long int inod,const char * routine)233 long int set_swit( long int element, long int inod, const char* routine )
234
235 {
236 long int i=0, result=0, iteration=0, icontrol=0, ldum=0;
237 double ddum[1];
238 char *str=NULL;
239
240 if ( strcmp(routine,swit_routine_stack[0]) ) {
241 for ( i=MSTACK-2; i>=0; i-- ) {
242 strcpy( swit_routine_stack[i+1], swit_routine_stack[i] );
243 }
244 strcpy( swit_routine_stack[0], routine );
245 }
246 if ( element>=0 )
247 swit_element_stack = element;
248 else
249 swit_element_stack = -1;
250 if ( inod>=0 )
251 swit_node_stack = inod;
252 else
253 swit_node_stack = -1;
254
255 if ( getenv("TOCHNOG_REPORT_ROUTINE")!=NULL ) {
256 pri( "In routine ", routine );
257 }
258
259 if ( getenv("TOCHNOG_DEBUG")!=NULL ) {
260 str = getenv("TOCHNOG_DEBUG");
261 if ( !strcmp("yes",str ) ) result = 1;
262 }
263 if ( result==0 ) return 0;
264
265 if ( element>=0 ) {
266 if ( getenv("TOCHNOG_ELEMENT")!=NULL ) {
267 str = getenv("TOCHNOG_ELEMENT");
268 if ( element!=atol(str) ) result=0;
269 }
270 }
271
272 if ( db_active_index( ICONTROL, 0, VERSION_NORMAL ) ) {
273 if ( getenv("TOCHNOG_ICONTROL")!=NULL ) {
274 db( ICONTROL, 0, &icontrol, ddum, ldum, VERSION_NORMAL, GET );
275 str = getenv("TOCHNOG_ICONTROL");
276 if ( icontrol!=atol(str) ) result=0;
277 }
278 }
279
280 if ( db_active_index( NUMBER_ITERATIONS, 0, VERSION_NORMAL ) ) {
281 if ( getenv("TOCHNOG_ITERATION")!=NULL ) {
282 db( NUMBER_ITERATIONS, 0, &iteration, ddum, ldum, VERSION_NORMAL, GET );
283 str = getenv("TOCHNOG_ITERATION");
284 if ( iteration!=atol(str) ) result=0;
285 }
286 }
287
288 if ( inod>=0 ) {
289 if ( getenv("TOCHNOG_NODE")!=NULL ) {
290 str = getenv("TOCHNOG_NODE");
291 if ( inod!=atol(str) ) result=0;
292 }
293 }
294
295 if ( strcmp(routine,"") ) {
296 if ( getenv("TOCHNOG_ROUTINE")!=NULL ) {
297 str = getenv("TOCHNOG_ROUTINE");
298 if ( strcmp(routine,str) ) result=0;
299 }
300 }
301
302 return result;
303
304 }
305
stress_indx(long int idim,long int jdim)306 long int stress_indx( long int idim, long int jdim )
307
308 {
309 long int kdim=0, ldim=0, indx=0;
310
311 if ( idim<jdim ) {
312 kdim = idim;
313 ldim = jdim;
314 }
315 else {
316 kdim = jdim;
317 ldim = idim;
318 }
319
320 if ( kdim==0 ) indx = 0 + ldim;
321 else if ( kdim==1 ) indx = 2 + ldim;
322 else if ( kdim==2 ) indx = 3 + ldim;
323
324 return indx;
325 }
326
long_to_a(long int n,char s[])327 char *long_to_a( long int n, char s[] )
328
329 {
330
331 int i=0, sign=0, m=0;
332 char *ptr=NULL;
333
334 m = (int) n;
335
336 if ( ( sign = m ) < 0 )
337 m = -m;
338
339 i = 0;
340 do {
341 s[i++] = m % 10 + '0';
342 } while ( ( m/= 10 ) > 0 );
343
344 if ( sign<0 )
345 s[i++] = '-';
346 s[i] = '\0';
347 string_reverse( s );
348
349 ptr = s;
350 return ptr;
351
352 }
353
exit_tn(long int print_database_type)354 void exit_tn( long int print_database_type )
355
356 {
357
358 long int itarget=0, ntarget=0, data_item_name=0,
359 data_item_index=0, number=0, correct=0, length=0, ldum=0,
360 dof_label[MUKNWN], *target_item=NULL, *int_data=NULL;
361 double value=0., tolerance=0., ddum[1], *dbl_data=NULL, *target_value=NULL;
362
363 if ( any_runtime ) {
364 pri( "*** Runtime file used in calculation. ***" );
365 }
366
367 if ( print_database_type==-RESTART ) {
368 print_restart( -1 );
369 }
370 else {
371 assert( print_database_type==-YES );
372 print_database( -1, VERSION_NORMAL, -YES );
373 }
374
375 print_gid( -YES );
376 if ( db_partialname_any("control_print_dx") ) print_dx( -YES );
377
378 db( DOF_LABEL, 0, dof_label, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
379 db_max_index( TARGET_ITEM, ntarget, VERSION_NORMAL, GET );
380 if ( ntarget>=0 ) {
381 for ( itarget=0; itarget<=ntarget; itarget++ ) {
382 if ( db_active_index( TARGET_ITEM, itarget, VERSION_NORMAL ) &&
383 db_active_index( TARGET_VALUE, itarget, VERSION_NORMAL ) ) {
384 target_item = db_int( TARGET_ITEM, itarget, VERSION_NORMAL );
385 data_item_name = labs(target_item[0]);
386 data_item_index = target_item[1];
387 if ( !db_active_index( data_item_name, data_item_index, VERSION_NORMAL ) ) {
388 ofstream out( "tn.log", ios::app );
389 out << "\nError in calculation with data file " << data_file << ".";
390 out.close();
391 exit(TN_EXIT_STATUS);
392 }
393 else {
394 length = db_len( data_item_name, data_item_index, VERSION_NORMAL );
395 if ( target_item[2]<0 ) {
396 array_member(dof_label,target_item[2],nuknwn,number);
397 if ( db_len(data_item_name,data_item_index,VERSION_NORMAL)==npuknwn )
398 number /= nder;
399 }
400 else
401 number = target_item[2];
402 if ( number>length-1 ) db_error( TARGET_ITEM, itarget );
403 target_value = db_dbl( TARGET_VALUE, itarget, VERSION_NORMAL );
404 value = target_value[0];
405 tolerance = target_value[1];
406 correct = 1;
407 if ( db_type(data_item_name)==INTEGER ) {
408 int_data = db_int( data_item_name, data_item_index, VERSION_NORMAL );
409 if ( int_data[number]<(value-tolerance) ||
410 int_data[number]>(value+tolerance) ) correct = 0;
411 }
412 else {
413 dbl_data = db_dbl( data_item_name, data_item_index, VERSION_NORMAL );
414 if ( dbl_data[number]<(value-tolerance) ||
415 dbl_data[number]>(value+tolerance) ) correct = 0;
416 }
417 if ( !correct ) {
418 ofstream out( "tn.log", ios::app );
419 out << "\nError in calculation with data file " << data_file << ".";
420 out << "\nTarget value for: ";
421 out << "data item " << db_name(data_item_name) << " ";
422 out << "with index " << data_item_index << " ";
423 if ( target_item[2]<0 ) {
424 out << "for " << db_name(target_item[2]) << " ";
425 }
426 else
427 out << "and value number " << number << " ";
428 out << "is " << value << ".";
429 if ( db_type(data_item_name)==INTEGER )
430 out << "\nThe actual value is " << int_data[number] << ".";
431 else
432 out << "\nThe actual value is " << dbl_data[number] << ".";
433 out << "\n";
434 out.close();
435 exit(TN_EXIT_STATUS);
436 }
437 }
438 }
439 }
440 }
441
442 db_close();
443
444 ofstream out( "tn.log", ios::app );
445 out << "\nCalculation with data file " << data_file << " ready.\n";
446 out.close();
447 exit(0);
448
449 }
450
451
exit_tn_on_error(void)452 void exit_tn_on_error( void )
453
454 {
455 long int i=0, element_group=0, iarea=0, ldum=0;
456 double ddum[1];
457
458 if ( strcmp(swit_routine_stack[0],"") ) {
459 pri( "\n\nLast routines:" );
460 for ( i=0; i<MSTACK; i++ ) {
461 if ( strcmp( swit_routine_stack[i], "" ) ) pri( swit_routine_stack[i] );
462 }
463 }
464 if ( swit_node_stack>=0 ) pri( "\n\nLast node:" , swit_node_stack );
465 if ( swit_element_stack>=0 ) {
466 pri( "\n\nLast element:" , swit_element_stack );
467 element_group = 0;
468 db( ELEMENT_GROUP, swit_element_stack, &element_group,
469 ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
470 pri( "Last element has ELEMENT_GROUP", element_group );
471 if ( db_active_index( ELEMENT_GROUP_AREA_ELEMENT_GROUP, swit_element_stack, VERSION_NORMAL ) ) {
472 db( ELEMENT_GROUP_AREA_ELEMENT_GROUP, swit_element_stack, &iarea,
473 ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
474 pri( "Last element got ELEMENT_GROUP from AREA_ELEMENT_GROUP with index", iarea );
475 }
476 if ( db_active_index( ELEMENT_GROUP_AREA_ELEMENT_GROUP_SEQUENCE_ELEMENTGROUP,
477 swit_element_stack, VERSION_NORMAL ) ) {
478 db( ELEMENT_GROUP_AREA_ELEMENT_GROUP_SEQUENCE_ELEMENTGROUP, swit_element_stack, &iarea,
479 ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
480 pri( "Last element got ELEMENT_GROUP from AREA_ELEMENT_GROUP_SEQUENCE_ELEMENTGROUP with index", iarea );
481 }
482 }
483
484 parallel_sys_lock();
485 parallel_active = 0;
486 print_database( -1, VERSION_NORMAL, -EVERYTHING );
487 print_gid( -YES );
488 cout << flush;
489 exit(TN_EXIT_STATUS);
490 parallel_sys_unlock();
491 }
492