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 SLIP 1
24 #define STICK 2
25
parallel_contact(void)26 void parallel_contact( void )
27
28 {
29 long int icont=0, inum=0, itar=0, idim=0, inol=0, inod=0,
30 nnol=0, iside=0, nside=0, nnol_side=0, name=0,
31 max_node=0, max_element=0, swit=0, ipuknwn=0, ldum=0,
32 length=0, use_element_target=0, use_contact_geometry=0,
33 max_contact_geometry=0, max_total=0, inside_target=0,
34 contact_geometry_switch=0, contact_stick=-YES, status=0, iteration=0,
35 number_of_side_boundary_nodes=0, iloop=0,
36 nloop=0, idat=0, max=0, any_contact_data=0, ithread=0,
37 max_node_boundary=0, idum[1], contact_geometry[10],
38 *el=NULL, *nodes=NULL, *side_nodes=NULL, *next_of_loop=NULL;
39 double penetration=0.,
40 tmp=0., temp=0., pres=0., pressure_penalty=0., temperature_penalty=0.,
41 velocity_penalty=0., dtime=0., normal_force=0., friction_force=0.,
42 pressure_force=0., temperature_force=0., friction_energy=0.,
43 contact_heat_generation=0., slip_size=0., contact_friction=0.,
44 fac=0., contact_relaxation=1., rdum=0.,
45 ddum[MNOL], *tar_coord=NULL, *node_force=NULL,
46 *normal_dir=NULL, *vec1=NULL, *vec2=NULL, *average_tar_coord=NULL,
47 *average_tar_side_coord=NULL, *cont_coord=NULL, *tar_side_coord=NULL,
48 *cont_vel=NULL, *average_tar_vel=NULL, *tar_side_vel=NULL,
49 *slip_dir=NULL, *penetration_coord=NULL, *weight=NULL,
50 *node_lhside=NULL, *node_rhside=NULL, *new_node_dof=NULL;
51
52 for ( idat=0; idat<MDAT; idat++ ) {
53 db_highest_index( idat, max, VERSION_NORMAL );
54 if ( db_data_class(idat)==CONTACT && max>=0 ) any_contact_data = 1;
55 }
56
57 db_max_index( NODE, max_node, VERSION_NORMAL, GET );
58 if ( max_node >= 0 && any_contact_data ) {
59
60 db_max_index( ELEMENT, max_element, VERSION_NORMAL, GET );
61 db_max_index( NODE_BOUNDARY, max_node_boundary, VERSION_NORMAL, GET );
62 db_max_index( CONTACT_GEOMETRY, max_contact_geometry, VERSION_NORMAL, GET );
63 max_total = 1;
64 max_total += max_contact_geometry;
65 if ( max_node_boundary>=0 ) max_total += max_element;
66
67 el = get_new_int(MNOL+1);
68 nodes = get_new_int(MNOL);
69 side_nodes = get_new_int(MNOL);
70 tar_coord = get_new_dbl( MNOL*MDIM );
71 node_force = get_new_dbl( MDIM );
72 normal_dir = get_new_dbl( MDIM );
73 vec1 = get_new_dbl( MDIM );
74 vec2 = get_new_dbl( MDIM );
75 average_tar_coord = get_new_dbl( MDIM );
76 average_tar_side_coord = get_new_dbl( MDIM );
77 cont_coord = get_new_dbl( MDIM );
78 tar_side_coord = get_new_dbl( MDIM*MNOL );
79 cont_vel = get_new_dbl( MDIM );
80 average_tar_vel = get_new_dbl( MDIM );
81 tar_side_vel = get_new_dbl( MDIM );
82 slip_dir = get_new_dbl( MDIM );
83 penetration_coord = get_new_dbl( MDIM );
84 weight = get_new_dbl( MDIM );
85
86 area_node_dataitem();
87
88 swit = set_swit(-1,-1,"contact");
89 if ( swit ) pri( "In routine CONTACT" );
90
91 db( NUMBER_ITERATIONS, 0, &iteration, ddum, ldum, VERSION_NEW, GET );
92 db( CONTACT_HEATGENERATION, 0, idum,
93 &contact_heat_generation, ldum, VERSION_NORMAL, GET_IF_EXISTS );
94 db( CONTACT_PENALTY_PRESSURE, 0, idum,
95 &pressure_penalty, ldum, VERSION_NORMAL, GET_IF_EXISTS );
96 db( CONTACT_PENALTY_TEMPERATURE, 0, idum,
97 &temperature_penalty, ldum, VERSION_NORMAL, GET_IF_EXISTS );
98 db( CONTACT_PENALTY_VELOCITY, 0, idum,
99 &velocity_penalty, ldum, VERSION_NORMAL, GET_IF_EXISTS );
100 db( CONTACT_FRICTION, 0, idum, &contact_friction,
101 ldum, VERSION_NORMAL, GET_IF_EXISTS );
102 db( CONTACT_STICK, 0, &contact_stick, ddum,
103 ldum, VERSION_NORMAL, GET_IF_EXISTS );
104 db( CONTACT_RELAXATION, 0, idum, &contact_relaxation,
105 ldum, VERSION_NORMAL, GET_IF_EXISTS );
106 db( DTIME, 0, idum, &dtime, ldum, VERSION_NEW, GET );
107 if ( swit ) {
108 pri( "pressure_penalty", pressure_penalty );
109 pri( "temperature_penalty", temperature_penalty );
110 pri( "velocity_penalty", velocity_penalty );
111 pri( "contact_friction", contact_friction );
112 }
113
114 // loop over contacters
115 if ( max_node>=0 ) {
116 next_of_loop = get_new_int(1+max_node);
117 parallel_sys_next_of_loop( next_of_loop, max_node, nloop, ithread );
118 for ( iloop=0; iloop<nloop; iloop++ ) {
119 icont = next_of_loop[iloop];
120 if ( icont>max_node )
121 break;
122 else if ( db_active_index( NODE, icont, VERSION_NORMAL ) ) {
123 // contacter coordinates
124 db( NODE, icont, idum, cont_coord, ldum, VERSION_NEW, GET );
125 new_node_dof = db_dbl( NODE_DOF, icont, VERSION_NEW );
126 node_rhside = db_dbl( NODE_RHSIDE, icont, VERSION_NORMAL );
127 for ( idim=0; idim<ndim; idim++ ) {
128 if ( materi_displacement ) cont_coord[idim] +=
129 new_node_dof[dis_indx+idim*nder];
130 // explicit prediction
131 if ( iteration==1 ) cont_coord[idim] +=
132 new_node_dof[vel_indx+idim*nder]*dtime;
133 cont_vel[idim] = new_node_dof[vel_indx+idim*nder];
134 }
135 if ( swit ) {
136 pri( "node_contacter", icont );
137 pri( "cont_coord", cont_coord, ndim );
138 pri( "cont_vel", cont_vel, ndim );
139 }
140 // loop over target geometries and elements
141 for ( inum=0; inum<=max_total; inum++ ) {
142 inside_target = 0;
143 if ( inum<=max_contact_geometry ) {
144 use_contact_geometry = 1;
145 use_element_target = 0;
146 itar = inum;
147 }
148 else {
149 use_contact_geometry = 0;
150 use_element_target = 1;
151 itar = inum - max_contact_geometry - 1;
152 }
153 if ( (use_element_target&&db_active_index(ELEMENT,itar,
154 VERSION_NORMAL))||
155 (use_contact_geometry&&db_active_index(CONTACT_GEOMETRY,itar,
156 VERSION_NORMAL)) ) {
157 if ( use_element_target ) {
158 if ( swit ) pri( "itar", itar );
159 db( ELEMENT, itar, el, ddum, length, VERSION_NORMAL, GET );
160 name = el[0]; nnol = length - 1;
161 if ( name!=-BAR2 && name!=-TRIA3 && name!=-QUAD4 &&
162 name!=-TET4 && name!=-HEX8 ) {
163 cout << "\nError: " << db_name( name );
164 cout << " is not available for contact analysis.\n";
165 exit(TN_EXIT_STATUS);
166 }
167 array_move( &el[1], nodes, nnol );
168 // target coordinates and average
169 array_set( average_tar_coord, 0., ndim );
170 array_set( average_tar_vel, 0., ndim );
171 for ( inol=0; inol<nnol; inol++ ) {
172 inod = nodes[inol];
173 new_node_dof = db_dbl( NODE_DOF, inod, VERSION_NEW );
174 db( NODE, inod, idum, &tar_coord[inol*ndim],
175 ldum, VERSION_NEW, GET );
176 for ( idim=0; idim<ndim; idim++ ) {
177 if ( materi_displacement ) tar_coord[inol*ndim+idim] +=
178 new_node_dof[dis_indx+idim*nder];
179 if ( iteration==1 ) tar_coord[inol*ndim+idim] +=
180 new_node_dof[vel_indx+idim*nder]*dtime;
181 average_tar_coord[idim] += tar_coord[inol*ndim+idim]/nnol;
182 average_tar_vel[idim] += new_node_dof[vel_indx+idim*nder]/nnol;
183 }
184 }
185 // inside target element?
186 if ( !array_member(nodes,icont,nnol,ldum) &&
187 point_el( cont_coord, tar_coord, ddum, name, nnol ) ) {
188 if ( swit ) {
189 pri( "tar_coord", tar_coord, nnol, ndim );
190 pri( "average_tar_coord", average_tar_coord, ndim );
191 pri( "average_tar_vel", average_tar_vel, ndim );
192 }
193 // loop over target sides
194 if ( name==-BAR2 ) {
195 nside = 2;
196 nnol_side = 1;
197 }
198 else if ( name==-TRIA3 ) {
199 nside = 3;
200 nnol_side = 2;
201 }
202 else if ( name==-QUAD4 ) {
203 nside = 4;
204 nnol_side = 2;
205 }
206 else if ( name==-TET4 ) {
207 nside = 4;
208 nnol_side = 3;
209 }
210 else if ( name==-HEX8 ) {
211 nside = 6;
212 nnol_side = 4;
213 }
214 else {
215 nside = 0;
216 nnol_side = 0;
217 }
218 for ( iside=0; iside<nside && !inside_target; iside++ ) {
219 if ( swit ) pri( "iside", iside );
220 if ( name==-BAR2 ) {
221 if ( iside==0 )
222 side_nodes[0] = 0;
223 else {
224 assert( iside==1 );
225 side_nodes[0] = 1;
226 }
227 }
228 else if ( name==-TRIA3 ) {
229 if ( iside==0 ) {
230 side_nodes[0] = 0;
231 side_nodes[1] = 1;
232 }
233 else if ( iside==1 ) {
234 side_nodes[0] = 1;
235 side_nodes[1] = 2;
236 }
237 else {
238 assert ( iside==2 );
239 side_nodes[0] = 2;
240 side_nodes[1] = 0;
241 }
242 }
243 else if ( name==-QUAD4 ) {
244 if ( iside==0 ) {
245 side_nodes[0] = 0;
246 side_nodes[1] = 1;
247 }
248 else if ( iside==1 ) {
249 side_nodes[0] = 1;
250 side_nodes[1] = 3;
251 }
252 else if ( iside==2 ) {
253 side_nodes[0] = 3;
254 side_nodes[1] = 2;
255 }
256 else {
257 assert( iside==3 );
258 side_nodes[0] = 2;
259 side_nodes[1] = 0;
260 }
261 }
262 else if ( name==-TET4 ) {
263 if ( iside==0 ) {
264 side_nodes[0] = 0;
265 side_nodes[1] = 1;
266 side_nodes[2] = 2;
267 }
268 else if ( iside==1 ) {
269 side_nodes[0] = 0;
270 side_nodes[1] = 1;
271 side_nodes[2] = 3;
272 }
273 else if ( iside==2 ) {
274 side_nodes[0] = 2;
275 side_nodes[1] = 3;
276 side_nodes[2] = 0;
277 }
278 else {
279 assert ( iside==3 );
280 side_nodes[0] = 2;
281 side_nodes[1] = 3;
282 side_nodes[2] = 1;
283 }
284 }
285 else if ( name==-HEX8 ) {
286 if ( iside==0 ) {
287 side_nodes[0] = 0;
288 side_nodes[1] = 1;
289 side_nodes[2] = 2;
290 side_nodes[3] = 3;
291 }
292 else if ( iside==1 ) {
293 side_nodes[0] = 0;
294 side_nodes[1] = 1;
295 side_nodes[2] = 4;
296 side_nodes[3] = 5;
297 }
298 else if ( iside==2 ) {
299 side_nodes[0] = 1;
300 side_nodes[1] = 2;
301 side_nodes[2] = 5;
302 side_nodes[3] = 6;
303 }
304 else if ( iside==3 ) {
305 side_nodes[0] = 2;
306 side_nodes[1] = 3;
307 side_nodes[2] = 6;
308 side_nodes[3] = 7;
309 }
310 else if ( iside==4 ) {
311 side_nodes[0] = 3;
312 side_nodes[1] = 0;
313 side_nodes[2] = 7;
314 side_nodes[3] = 4;
315 }
316 else if ( iside==5 ) {
317 side_nodes[0] = 4;
318 side_nodes[1] = 5;
319 side_nodes[2] = 6;
320 side_nodes[3] = 7;
321 }
322 }
323 number_of_side_boundary_nodes = 0;
324 for ( inol=0; inol<nnol_side; inol++ ) {
325 inod = nodes[side_nodes[inol]];
326 if ( db_active_index( NODE_BOUNDARY, inod, VERSION_NORMAL ) )
327 number_of_side_boundary_nodes++;
328 }
329 if ( swit ) pri( "number_of_side_boundary_nodes",
330 number_of_side_boundary_nodes );
331 if ( number_of_side_boundary_nodes==ndim ) {
332 inside_target = 1;
333 // side properties
334 array_set( average_tar_side_coord, 0., ndim );
335 array_set( tar_side_vel, 0., ndim );
336 for ( inol=0; inol<nnol_side; inol++ ) {
337 inod = nodes[side_nodes[inol]];
338 new_node_dof = db_dbl( NODE_DOF, inod, VERSION_NEW );
339 db( NODE, inod, idum, &tar_side_coord[inol*ndim], ldum,
340 VERSION_NEW, GET );
341 for ( idim=0; idim<ndim; idim++ ) {
342 if ( materi_displacement )
343 tar_side_coord[inol*ndim+idim] +=
344 new_node_dof[dis_indx+idim*nder];
345 if ( iteration==1 )
346 tar_side_coord[inol*ndim+idim] +=
347 new_node_dof[vel_indx+idim*nder]*dtime;
348 average_tar_side_coord[idim] +=
349 tar_side_coord[inol*ndim+idim]/nnol_side;
350 tar_side_vel[idim] +=
351 new_node_dof[vel_indx+idim*nder] / nnol_side;
352 }
353 }
354 if ( swit ) {
355 pri( "average_tar_side_coord",
356 average_tar_side_coord, ndim );
357 pri( "tar_side_vel", tar_side_vel, ndim );
358 }
359 if ( nnol_side==1 ) {
360 weight[0] = 1.;
361 }
362 else if ( nnol_side==2 ) {
363 project_point_exactly_on_line( cont_coord,
364 &tar_side_coord[0*ndim],
365 &tar_side_coord[1*ndim], weight );
366 }
367 else if ( nnol_side==3 ) {
368 project_point_exactly_on_triangle( cont_coord,
369 &tar_side_coord[0*ndim],
370 &tar_side_coord[1*ndim],
371 &tar_side_coord[2*ndim],
372 weight );
373 }
374 else {
375 assert( nnol_side==4 );
376 project_point_exactly_on_quad( cont_coord,
377 &tar_side_coord[0*ndim],
378 &tar_side_coord[1*ndim],
379 &tar_side_coord[2*ndim],
380 &tar_side_coord[3*ndim],
381 weight );
382 }
383 // normal
384 if ( nnol_side==1 ) {
385 normal_dir[0] = -1.;
386 }
387 else if ( nnol_side==2 ) {
388 array_subtract( &tar_side_coord[1*ndim],
389 &tar_side_coord[0*ndim], vec1, ndim );
390 normal_dir[0] = -vec1[1];
391 normal_dir[1] = +vec1[0];
392 }
393 else if ( nnol_side==3 ) {
394 array_subtract( &tar_side_coord[1*ndim],
395 &tar_side_coord[0*ndim], vec1, ndim );
396 array_subtract( &tar_side_coord[2*ndim],
397 &tar_side_coord[0*ndim], vec2, ndim );
398 array_outproduct_3D( vec1, vec2, normal_dir );
399 }
400 else {
401 assert( nnol_side==4 );
402 array_subtract( &tar_side_coord[1*ndim],
403 &tar_side_coord[0*ndim], vec1, ndim );
404 array_subtract( &tar_side_coord[2*ndim],
405 &tar_side_coord[0*ndim], vec2, ndim );
406 array_outproduct_3D( vec1, vec2, normal_dir );
407 }
408 // check if normal is outward
409 array_subtract( average_tar_side_coord,
410 average_tar_coord, vec1, ndim );
411 tmp = array_inproduct( vec1, normal_dir, ndim );
412 if ( tmp<0. )
413 array_multiply( normal_dir, normal_dir, -1., ndim );
414 array_normalize( normal_dir, ndim );
415 // penetration
416 array_set( penetration_coord, 0., ndim );
417 for ( inol=0; inol<nnol_side; inol++ ) {
418 for ( idim=0; idim<ndim; idim++ ) {
419 penetration_coord[idim] += weight[inol] *
420 tar_side_coord[inol*ndim+idim];
421 }
422 }
423 array_subtract( cont_coord, penetration_coord, vec1, ndim );
424 penetration = array_inproduct( vec1, normal_dir, ndim );
425 if ( swit ) {
426 pri( "normal", normal_dir, ndim );
427 pri( "penetration_coord", penetration_coord, ndim );
428 pri( "penetration vector", vec1, ndim );
429 pri( "penetration normal", normal_dir, ndim );
430 pri( "penetration", penetration );
431 }
432 }
433 }
434 }
435 }
436 else {
437 assert( use_contact_geometry );
438 db( CONTACT_GEOMETRY, itar, contact_geometry, ddum,
439 ldum, VERSION_NORMAL, GET );
440 geometry( icont, ddum, contact_geometry, inside_target, rdum,
441 normal_dir, penetration, ddum, PLUS_DISPLACEMENT,
442 PROJECT_ON_EDGE, VERSION_NEW );
443 if ( inside_target ) {
444 if ( db_active_index( CONTACT_GEOMETRY_SWITCH,
445 itar, VERSION_NORMAL ) ) {
446 db( CONTACT_GEOMETRY_SWITCH, itar,
447 &contact_geometry_switch, ddum, ldum, VERSION_NORMAL, GET );
448 if ( contact_geometry_switch==-YES ) {
449 array_multiply( normal_dir, normal_dir, -1., ndim );
450 penetration = -penetration;
451 }
452 }
453 if ( penetration<0. )
454 inside_target = 1;
455 else
456 inside_target = 0;
457 if ( swit ) {
458 pri( "normal", normal_dir, ndim );
459 pri( "penetration", penetration );
460 pri( "inside_target", inside_target );
461 }
462 }
463 }
464 if ( inside_target ) {
465 if ( use_element_target ) {
466 fac = 2.;
467 if ( groundflow_pressure ) {
468 new_node_dof = db_dbl( NODE_DOF, icont, VERSION_NEW );
469 pres = -new_node_dof[pres_indx];
470 for ( inol=0; inol<nnol_side; inol++ ) {
471 inod = nodes[side_nodes[inol]];
472 new_node_dof = db_dbl( NODE_DOF, inod, VERSION_NEW );
473 pres += weight[inol] * new_node_dof[pres_indx];
474 }
475 pressure_force = pressure_penalty * pres;
476 }
477 if ( condif_temperature ) {
478 new_node_dof = db_dbl( NODE_DOF, icont, VERSION_NEW );
479 temp = -new_node_dof[temp_indx];
480 for ( inol=0; inol<nnol_side; inol++ ) {
481 inod = nodes[side_nodes[inol]];
482 new_node_dof = db_dbl( NODE_DOF, inod, VERSION_NEW );
483 temp += weight[inol] * new_node_dof[temp_indx];
484 }
485 temperature_force = temperature_penalty * temp;
486 }
487 }
488 else {
489 assert( use_contact_geometry );
490 fac = 1.;
491 array_set( tar_side_vel, 0., ndim );
492 }
493 // slip direction
494 array_subtract( tar_side_vel, cont_vel, slip_dir, ndim );
495 tmp = array_inproduct( normal_dir, slip_dir, ndim );
496 array_multiply( normal_dir, vec1, tmp, ndim );
497 array_subtract( slip_dir, vec1, slip_dir, ndim );
498 slip_size = array_size( slip_dir, ndim );
499 array_normalize( slip_dir, ndim );
500 // penetration force
501 normal_force = velocity_penalty * scalar_dabs(penetration);
502 // slip force
503 for ( idim=0; idim<ndim; idim++ ) {
504 ipuknwn = vel_indx/nder + idim;
505 node_force[idim] = node_rhside[ipuknwn] * slip_dir[idim];
506 }
507 if ( contact_stick==-NO ||
508 array_size(node_force,ndim)>=0.5*contact_friction*normal_force ) {
509 status = SLIP;
510 friction_force = contact_friction * normal_force;
511 friction_energy = contact_heat_generation * friction_force * slip_size;
512 }
513 else {
514 status = STICK;
515 friction_force = velocity_penalty * slip_size * dtime;
516 friction_energy = 0.;
517 }
518 if ( swit ) {
519 pri( "normal_force", normal_force );
520 pri( "slip_dir", slip_dir, ndim );
521 pri( "node_force", node_force, ndim );
522 if ( status==STICK ) pri( "status is STICK" );
523 else pri( "status is SLIP" );
524 pri( "friction_force", friction_force );
525 pri( "friction_energy", friction_energy );
526 pri( "pressure_force", pressure_force );
527 pri( "temperature_force", temperature_force );
528 }
529 // contacter contributions
530 node_lhside = db_dbl( NODE_LHSIDE, icont, VERSION_NORMAL );
531 node_rhside = db_dbl( NODE_RHSIDE, icont, VERSION_NORMAL );
532 if ( groundflow_pressure ) {
533 ipuknwn = pres_indx/nder;
534 node_lhside[ipuknwn] += pressure_penalty * fac;
535 node_rhside[ipuknwn] += pressure_force;
536 }
537 if ( condif_temperature ) {
538 ipuknwn = temp_indx/nder;
539 node_lhside[ipuknwn] += temperature_penalty * fac;
540 node_rhside[ipuknwn] += temperature_force;
541 if ( use_contact_geometry )
542 node_rhside[ipuknwn] += friction_energy;
543 else
544 node_rhside[ipuknwn] += 0.5 * friction_energy;
545 }
546 for ( idim=0; idim<ndim; idim++ ) {
547 ipuknwn = vel_indx/nder + idim;
548 // normal force on contacter
549 node_lhside[ipuknwn] += velocity_penalty * contact_relaxation *
550 scalar_dabs(normal_dir[idim]) * dtime * fac;
551 node_rhside[ipuknwn] += normal_force * normal_dir[idim];
552 // friction force on contacter
553 if ( status==STICK ) {
554 node_lhside[ipuknwn] += velocity_penalty *
555 scalar_dabs(slip_dir[idim]) * dtime * fac;
556 node_rhside[ipuknwn] += friction_force * slip_dir[idim];
557 }
558 else {
559 assert( status==SLIP );
560 node_rhside[ipuknwn] += friction_force * slip_dir[idim];
561 }
562 }
563 if ( swit ) {
564 pri( "cont_inod", icont );
565 pri( "node_lhside", node_lhside, npuknwn );
566 pri( "node_rhside", node_rhside, npuknwn );
567 }
568 // target contributions
569 if ( use_element_target ) {
570 parallel_sys_lock();
571 for ( inol=0; inol<nnol_side; inol++ ) {
572 inod = nodes[side_nodes[inol]];
573 node_lhside = db_dbl( NODE_LHSIDE, inod, VERSION_NORMAL );
574 node_rhside = db_dbl( NODE_RHSIDE, inod, VERSION_NORMAL );
575 if ( groundflow_pressure ) {
576 ipuknwn = pres_indx/nder;
577 node_lhside[ipuknwn] += pressure_penalty * weight[inol] * fac;
578 node_rhside[ipuknwn] -= pressure_force * weight[inol];
579 }
580 if ( condif_temperature ) {
581 ipuknwn = temp_indx/nder;
582 node_lhside[ipuknwn] += temperature_penalty * weight[inol] * fac;
583 node_rhside[ipuknwn] -= ( temperature_force +
584 0.5 * friction_energy ) * weight[inol];
585 }
586 for ( idim=0; idim<ndim; idim++ ) {
587 ipuknwn = vel_indx/nder + idim;
588 // normal force on target
589 node_lhside[ipuknwn] += contact_relaxation *
590 velocity_penalty * weight[inol] *
591 scalar_dabs(normal_dir[idim]) * dtime * fac;
592 node_rhside[ipuknwn] -=
593 normal_force * weight[inol] * normal_dir[idim];
594 // friction force on target
595 if ( status==STICK ) {
596 node_lhside[ipuknwn] +=
597 velocity_penalty * weight[inol] *
598 scalar_dabs(slip_dir[idim]) * dtime * fac;
599 node_rhside[ipuknwn] -=
600 friction_force * weight[inol] * slip_dir[idim];
601 }
602 else {
603 assert( status==SLIP );
604 node_rhside[ipuknwn] -=
605 friction_force * weight[inol] * slip_dir[idim];
606 }
607 }
608 if ( swit ) {
609 pri( "tar_inod", inod );
610 pri( "node_lhside", node_lhside, npuknwn );
611 pri( "node_rhside", node_rhside, npuknwn );
612 }
613 }
614 parallel_sys_unlock();
615 }
616 }
617 }
618 }
619 }
620 }
621 delete[] next_of_loop;
622 }
623
624 delete[] el;
625 delete[] nodes;
626 delete[] side_nodes;
627 delete[] tar_coord;
628 delete[] node_force;
629 delete[] normal_dir;
630 delete[] vec1;
631 delete[] vec2;
632 delete[] average_tar_coord;
633 delete[] average_tar_side_coord;
634 delete[] cont_coord;
635 delete[] tar_side_coord;
636 delete[] cont_vel;
637 delete[] average_tar_vel;
638 delete[] tar_side_vel;
639 delete[] slip_dir;
640 delete[] penetration_coord;
641 delete[] weight;
642
643 if ( swit ) pri( "Out routine CONTACT" );
644 }
645
646 }
647