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
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software Foundation
18 59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
19 */
20
21 #include "tochnog.h"
22
23 #define EPS_ATAN 1.e-10
24
bounda()25 void bounda( )
26
27 {
28 long int in=0, ready=0, iuknwn=0, inod=0, iboun=0, found=0,
29 max_bounda=0, max_bounda_unknown=0, max_bounda_force=0,
30 bounda_time_user=0, ind1=0, ipuknwn=0, iu=0, iu_start=0, iu_end=0,
31 inc=0, ninc=0, length=0, range_length=0, unknown=0, force=0,
32 idim=0, time=0, sine=0, user=0, rotate=0, rotate_axis=0,
33 use_geom=0, swit=0, swit_node=0, bounda_length=0,
34 use_range=0, use_all=0, use_node_set=0, ifreq=0, nfreq=0, indx=0, max_node=0,
35 bounda_time_file=0, length_bounda_time=0, ldum=0, idum[1],
36 *val=NULL, *integer_range=NULL,
37 *dof_label=NULL, *dof_type=NULL, *node_bounded=NULL;
38 double load=0., time0=0., time1=0., load0=0., load1=0., factor=0.,
39 time_current=0., time_total=0., dtime=0., amplitude=0., frequency=0.,
40 time_start=0., radius=0., angle_start=0.,
41 angle_total=0., rdum=0., ddum[MDIM], coord_start[MDIM], coord_total[MDIM],
42 *bounda_time=NULL, *new_node_dof=NULL,
43 *node_dof=NULL, *bounda_sine=NULL, *node_rhside=NULL;
44
45 swit = set_swit(-1,-1,"bounda");
46 if ( swit ) pri( "In routine BOUNDA" );
47
48 val = get_new_int(MBOUNDA);
49 integer_range = get_new_int(MRANGE);
50 dof_label = get_new_int(MUKNWN);
51 dof_type = get_new_int(MUKNWN);
52 bounda_time = get_new_dbl(DATA_ITEM_SIZE);
53
54 db_set_int( NODE_BOUNDED, VERSION_NORMAL );
55
56 groundflow_phreatic_apply();
57
58 db_max_index( BOUNDA_UNKNOWN, max_bounda_unknown, VERSION_NORMAL, GET );
59 db_max_index( BOUNDA_FORCE, max_bounda_force, VERSION_NORMAL, GET );
60 if ( max_bounda_unknown<0 && max_bounda_force<0 ) goto end_of_bounda;
61
62 db_max_index( NODE, max_node, VERSION_NORMAL, GET );
63 if ( max_node<0 ) goto end_of_bounda;
64
65 nodes_in_geometry = get_new_int( 1+max_node );
66
67 db( DOF_LABEL, 0, dof_label, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
68 db( DOF_TYPE, 0, dof_type, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
69 db( DTIME, 0, idum, &dtime, ldum, VERSION_NEW, GET_IF_EXISTS );
70 db( TIME_CURRENT, 0, idum, &time_current, ldum,
71 VERSION_NORMAL, GET_IF_EXISTS );
72 time_total = time_current + dtime;
73 if ( swit ) pri( "time_total", time_total );
74
75 if ( max_bounda_unknown>max_bounda ) max_bounda = max_bounda_unknown;
76 if ( max_bounda_force>max_bounda ) max_bounda = max_bounda_force;
77 for ( iboun=0; iboun<=max_bounda; iboun++ ) {
78 unknown = db_active_index( BOUNDA_UNKNOWN, iboun, VERSION_NORMAL );
79 force = db_active_index( BOUNDA_FORCE, iboun, VERSION_NORMAL );
80 bounda_time_user = -NO; db( BOUNDA_TIME_USER, iboun, &bounda_time_user, ddum,
81 ldum, VERSION_NORMAL, GET_IF_EXISTS );
82 if ( unknown || force ) {
83 if ( swit ) pri( "iboun", iboun );
84 time = sine = user = 0; ninc = 2;
85 if ( db_active_index( BOUNDA_SINE, iboun, VERSION_NORMAL ) ) {
86 bounda_sine = db_dbl( BOUNDA_SINE, iboun, VERSION_NORMAL );
87 time_start = bounda_sine[0];
88 nfreq = ( db_len( BOUNDA_SINE, iboun, VERSION_NORMAL ) - 1 ) / 2;
89 ninc = 2;
90 sine = 1;
91 }
92 else if ( bounda_time_user==-YES ) {
93 ninc = 2;
94 user = 1;
95 }
96 else if ( db_active_index( BOUNDA_TIME_FILE, iboun, VERSION_NORMAL ) ) {
97 db( BOUNDA_TIME_FILE, iboun, &bounda_time_file, ddum, ldum,
98 VERSION_NORMAL, GET_IF_EXISTS );
99 if ( bounda_time_file==-YES ) {
100 bounda_time_file_apply( iboun, time_total, bounda_time, ninc );
101 time = 1;
102 }
103 else
104 db_error( BOUNDA_TIME_FILE, iboun );
105 }
106 else if ( db_active_index( BOUNDA_TIME, iboun, VERSION_NORMAL ) ) {
107 db( BOUNDA_TIME, iboun, idum, bounda_time, length_bounda_time,
108 VERSION_NORMAL, GET );
109 time = 1;
110 if ( length_bounda_time==1 )
111 ninc = 2;
112 else {
113 if ( length_bounda_time<4 ) db_error( BOUNDA_TIME, iboun );
114 ninc = length_bounda_time / 2;
115 }
116 }
117 else {
118 ninc = 2;
119 time = 1;
120 length_bounda_time = 0;
121 }
122
123 if ( unknown ) {
124 db( BOUNDA_UNKNOWN, iboun, val, ddum, bounda_length,
125 VERSION_NORMAL, GET );
126 if ( bounda_length<2 ) db_error( BOUNDA_UNKNOWN, iboun );
127 rotate = 0; rotate_axis = val[bounda_length-1];
128 if ( rotate_axis==-ROTATION_X_AXIS ) {
129 rotate = 1;
130 val[bounda_length-1] = dof_label[vel_indx+1*nder];
131 val[bounda_length-0] = dof_label[vel_indx+2*nder];
132 bounda_length++;
133 }
134 else if ( rotate_axis==-ROTATION_Y_AXIS ) {
135 rotate = 1;
136 val[bounda_length-1] = dof_label[vel_indx+0*nder];
137 val[bounda_length-0] = dof_label[vel_indx+2*nder];
138 bounda_length++;
139 }
140 else if ( rotate_axis==-ROTATION_Z_AXIS ) {
141 rotate = 1;
142 val[bounda_length-1] = dof_label[vel_indx+0*nder];
143 val[bounda_length-0] = dof_label[vel_indx+1*nder];
144 bounda_length++;
145 }
146 }
147 else {
148 assert( force );
149 db( BOUNDA_FORCE, iboun, val, ddum, bounda_length,
150 VERSION_NORMAL, GET );
151 if ( bounda_length<2 ) db_error( BOUNDA_FORCE, iboun );
152 }
153 if ( val[0]<0 && db_data_class(val[0])==GEOMETRY ) {
154 array_move( val, geometry_ent, 2 );
155 array_set( nodes_in_geometry, 0, 1+max_node );
156 parallel_sys_routine( ¶llel_geometry );
157 }
158
159 found = 0;
160 for ( inc=0; !found && inc<ninc-1; inc++ ) {
161
162 if ( time ) {
163 if ( length_bounda_time==0 ) {
164 load = 0.;
165 found = 1;
166 }
167 else if ( length_bounda_time==1 ) {
168 load = bounda_time[0];
169 found = 1;
170 }
171 else {
172 time0 = bounda_time[inc*2+0];
173 load0 = bounda_time[inc*2+1];
174 time1 = bounda_time[inc*2+2];
175 load1 = bounda_time[inc*2+3];
176 if ( time0>=time1 ) db_error( BOUNDA_TIME, iboun );
177 if ( swit ) {
178 pri( "time_total", time_total );
179 pri( "time0", time0 );
180 pri( "time1", time1 );
181 pri( "load0", load0 );
182 pri( "load1", load1 );
183 }
184 if ( time_total>=(time0-1.e-10) && time_total<=time1 ) {
185 found = 1;
186 if ( time0==time1 ) load = load0;
187 else load = load0 + (load1-load0)*(time_total-time0)/(time1-time0);
188 }
189 }
190 }
191 else if ( user ) {
192 load = 0.;
193 found = 1;
194 user_bounda_time( iboun, time_total, load );
195 }
196 else if ( sine )
197 found = time_current>=time_start;
198 else {
199 load = 0.;
200 found = 1;
201 }
202
203 if ( found ) {
204
205 iu_end = bounda_length - 1;
206 use_geom = use_range = use_all = use_node_set = 0;
207 if ( val[0]==-RA ) {
208 use_range = 1;
209 range_expand( val, integer_range, length, range_length );
210 iu_start = length;
211 }
212 else if ( val[0]==-ALL ) {
213 use_all = 1;
214 iu_start = 1;
215 }
216 else if ( val[0]==-NODE_SET ) {
217 use_node_set = 1;
218 iu_start = 1;
219 }
220 else if ( val[0]<0 && db_data_class(val[0])==GEOMETRY ) {
221 use_geom = 1;
222 iu_start = 2;
223 }
224 else
225 iu_start = 1;
226
227 for ( in=0, ready=0; !ready; in++ ) {
228 found = 0;
229 if ( use_range ) {
230 inod = integer_range[in];
231 factor = 1.;
232 found = 1;
233 }
234 else if ( use_all ) {
235 if ( db_active_index( NODE, in, VERSION_NORMAL ) ) {
236 inod = in;
237 factor = 1.;
238 found = 1;
239 }
240 }
241 else if ( use_geom ) {
242 if ( nodes_in_geometry[in] ) {
243 found = 1;
244 inod = in;
245 geometry( in, ddum, val, found, factor, ddum, rdum,
246 ddum, NODE_START_REFINED, PROJECT_EXACT, VERSION_NORMAL );
247 }
248 }
249 else if ( use_node_set ) {
250 inod = in;
251 if ( db_active_index( NODE_SET, inod, VERSION_NORMAL ) ) {
252 found = 1;
253 factor = 1.;
254 }
255 }
256 else {
257 inod = val[0];
258 factor = 1.;
259 found = 1;
260 }
261 if ( found ) {
262 if ( db_active_index( NODE, inod, VERSION_NORMAL ) ) {
263 swit_node = swit;
264 swit = set_swit(-1,inod,"bounda");
265 if ( swit ) pri( "inod", inod );
266 for ( iu=iu_start; iu<=iu_end; iu++ ) {
267 array_member( dof_label, val[iu], nuknwn, iuknwn );
268 if ( iuknwn<0 ) {
269 if ( unknown )
270 db_error( BOUNDA_UNKNOWN, iboun );
271 else
272 db_error( BOUNDA_FORCE, iboun );
273 }
274 ipuknwn = iuknwn / nder;
275 if ( unknown ) {
276 node_bounded = db_int( NODE_BOUNDED, inod, VERSION_NORMAL );
277 new_node_dof = db_dbl( NODE_DOF, inod, VERSION_NEW );
278 node_dof = db_dbl( NODE_DOF, inod, VERSION_NORMAL );
279 if ( swit ) pri( "node_dof", node_dof, nuknwn );
280 ind1 = ipuknwn*nder + nder - 1;
281 node_bounded[ipuknwn] = 1;
282 if ( sine ) {
283 new_node_dof[iuknwn] = new_node_dof[ind1] = 0.;
284 for ( ifreq=0; ifreq<nfreq; ifreq++ ) {
285 frequency = bounda_sine[1+ifreq*2+0];
286 amplitude = bounda_sine[1+ifreq*2+1];
287 new_node_dof[iuknwn] += factor * amplitude *
288 sin( 2. * PIRAD * frequency * time_total );
289 if ( derivatives ) new_node_dof[ind1] +=
290 factor * amplitude * 2. * PIRAD * frequency *
291 cos( 2. * PIRAD * frequency * time_total );
292 }
293 }
294 else if ( rotate ) {
295 db( NODE, inod, idum, coord_start, ldum,
296 VERSION_NORMAL, GET );
297 if ( materi_displacement ) {
298 for ( idim=0; idim<ndim; idim++ )
299 coord_start[idim] += node_dof[dis_indx+idim*nder];
300 }
301 if ( rotate_axis==-ROTATION_X_AXIS ) {
302 radius = sqrt( coord_start[1]*coord_start[1] +
303 coord_start[2]*coord_start[2] );
304 if ( scalar_dabs(coord_start[1])<EPS_ATAN ) {
305 if ( coord_start[2]>=0. )
306 angle_start = 1.*PIRAD/2.;
307 else
308 angle_start = 3.*PIRAD/2.;
309 }
310 else {
311 angle_start =
312 atan(scalar_dabs(coord_start[2]/coord_start[1]));
313 if ( coord_start[1]<0. && coord_start[2]>=0. )
314 angle_start = 1.*PIRAD - angle_start;
315 else if ( coord_start[1]<0. && coord_start[2]<0. )
316 angle_start = 1.*PIRAD + angle_start;
317 else if ( coord_start[1]>0. && coord_start[2]<0. )
318 angle_start = 2.*PIRAD - angle_start;
319 }
320 angle_total = angle_start +
321 factor * load * dtime * PIRAD / 180.;
322 coord_total[0] = coord_start[0];
323 coord_total[1] = radius * cos(angle_total);
324 coord_total[2] = radius * sin(angle_total);
325 }
326 else if ( rotate_axis==-ROTATION_Y_AXIS ) {
327 radius = sqrt( coord_start[0]*coord_start[0] +
328 coord_start[2]*coord_start[2] );
329 if ( scalar_dabs(coord_start[0])<EPS_ATAN ) {
330 if ( coord_start[2]>=0. )
331 angle_start = 1.*PIRAD/2.;
332 else
333 angle_start = 3.*PIRAD/2.;
334 }
335 else {
336 angle_start =
337 atan(scalar_dabs(coord_start[2]/coord_start[0]));
338 if ( coord_start[0]<0. && coord_start[2]>=0. )
339 angle_start = 1.*PIRAD - angle_start;
340 else if ( coord_start[0]<0. && coord_start[2]<0. )
341 angle_start = 1.*PIRAD + angle_start;
342 else if ( coord_start[0]>0. && coord_start[2]<0. )
343 angle_start = 2.*PIRAD - angle_start;
344 }
345 angle_total = angle_start -
346 factor * load * dtime * PIRAD / 180.;
347 coord_total[1] = coord_start[1];
348 coord_total[0] = radius * cos(angle_total);
349 coord_total[2] = radius * sin(angle_total);
350 }
351 else {
352 assert( rotate_axis==-ROTATION_Z_AXIS );
353 radius = sqrt( coord_start[0]*coord_start[0] +
354 coord_start[1]*coord_start[1] );
355 if ( scalar_dabs(coord_start[0])<EPS_ATAN ) {
356 if ( coord_start[1]>=0. )
357 angle_start = 1.*PIRAD/2.;
358 else
359 angle_start = 3.*PIRAD/2.;
360 }
361 else {
362 angle_start =
363 atan(scalar_dabs(coord_start[1]/coord_start[0]));
364 if ( coord_start[0]<0. && coord_start[1]>=0. )
365 angle_start = 1.*PIRAD - angle_start;
366 else if ( coord_start[0]<0. && coord_start[1]<0. )
367 angle_start = 1.*PIRAD + angle_start;
368 else if ( coord_start[0]>0. && coord_start[1]<0. )
369 angle_start = 2.*PIRAD - angle_start;
370 }
371 angle_total = angle_start +
372 factor * load * dtime * PIRAD / 180.;
373 coord_total[2] = coord_start[2];
374 coord_total[0] = radius * cos(angle_total);
375 coord_total[1] = radius * sin(angle_total);
376 }
377 idim = ( iuknwn - vel_indx ) / nder;
378 new_node_dof[iuknwn] =
379 ( coord_total[idim] - coord_start[idim] ) / dtime;
380 if ( derivatives ) new_node_dof[ind1] =
381 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
382 }
383 else {
384 new_node_dof[iuknwn] = factor * load;
385 if ( derivatives ) new_node_dof[ind1] =
386 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
387 }
388 if ( dof_type[iuknwn]==-MATERI_DISPLACEMENT ) {
389 indx = vel_indx + ( iuknwn - dis_indx );
390 new_node_dof[indx] =
391 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
392 node_bounded[indx/nder] = 1;
393 }
394 if ( dof_type[iuknwn]==-MATERI_VELOCITY_INTEGRATED ) {
395 indx = vel_indx + ( iuknwn - veli_indx );
396 new_node_dof[indx] =
397 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
398 node_bounded[indx/nder] = 1;
399 }
400 if ( dof_type[iuknwn]==-MAXWELL_E ) {
401 indx = maxfe_indx + ( iuknwn - maxe_indx );
402 new_node_dof[indx] =
403 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
404 node_bounded[indx/nder] = 1;
405 }
406 if ( dof_type[iuknwn]==-WAVE_SCALAR ) {
407 new_node_dof[fscal_indx] =
408 ( new_node_dof[iuknwn] - node_dof[iuknwn] ) / dtime;
409 node_bounded[fscal_indx/nder] = 1;
410 }
411 if ( swit ) {
412 pri( "node_bounded", node_bounded, npuknwn );
413 pri( "new_node_dof", new_node_dof, nuknwn );
414 }
415 }
416 else if ( force ) {
417 if ( dof_type[iuknwn]==-MATERI_DISPLACEMENT )
418 indx = vel_indx + ( iuknwn - dis_indx );
419 else if ( dof_type[iuknwn]==-MATERI_VELOCITY_INTEGRATED )
420 indx = vel_indx + ( iuknwn - veli_indx );
421 else if ( dof_type[iuknwn]==-MATERI_DISPLACEMENT )
422 indx = fscal_indx + ( iuknwn - scal_indx );
423 else
424 indx = iuknwn;
425 indx /= nder;
426 node_rhside = db_dbl( NODE_RHSIDE, inod, VERSION_NORMAL );
427 if ( sine ) {
428 node_rhside[indx] = 0.;
429 for ( ifreq=0; ifreq<nfreq; ifreq++ ) {
430 frequency = bounda_sine[1+ifreq*2+0];
431 amplitude = bounda_sine[1+ifreq*2+1];
432 node_rhside[indx] += factor*amplitude*
433 sin(2.*PIRAD*frequency*time_total);
434 }
435 }
436 else {
437 node_rhside[indx] = factor * load;
438 }
439 }
440 }
441 swit = swit_node;
442 }
443 }
444 if ( use_range )
445 ready = ( (in+1)==range_length );
446 else if ( use_all )
447 ready = ( in==max_node );
448 else if ( use_geom )
449 ready = ( in==max_node );
450 else if ( use_node_set )
451 ready = ( in==max_node );
452 else
453 ready = 1;
454 }
455 }
456 }
457 }
458 }
459
460 delete[] nodes_in_geometry;
461
462 end_of_bounda:
463
464 delete[] val;
465 delete[] integer_range;
466 delete[] dof_label;
467 delete[] dof_type;
468 delete[] bounda_time;
469
470 if ( swit ) pri( "Out routine BOUNDA" );
471
472 }
473
bounda_time_file_apply(long int iboun,double time_total,double bounda_time[],long int & ninc)474 void bounda_time_file_apply( long int iboun, double time_total,
475 double bounda_time[], long int &ninc )
476
477 {
478 long int n=0, found=0;
479 char str[MCHAR], filename[MCHAR];
480 double old_load=0., old_time=0., new_load=0., new_time=0.;
481
482 strcpy( filename, long_to_a(iboun,str) );
483 strcat( filename, ".bounda" );
484 ifstream in( filename );
485 if ( !in ) db_error( BOUNDA_TIME_FILE, iboun );
486
487 ninc = 0;
488 for ( ; !in.eof() && !found ; ) {
489 n++;
490 old_load = new_load;
491 old_time = new_time;
492 if ( !(in>>new_time) || !(in>>new_load) )
493 db_error( BOUNDA_TIME_FILE, iboun );
494 if ( n>1 && time_total>=old_time && time_total<=new_time ) {
495 found = 1;
496 ninc = 2;
497 bounda_time[0] = old_time;
498 bounda_time[1] = old_load;
499 bounda_time[2] = new_time;
500 bounda_time[3] = new_load;
501 }
502 }
503 if ( ninc<2 ) db_error( BOUNDA_TIME_FILE, iboun );
504
505 in.close();
506 }
507