?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix coding style
[survex.git] / src / network.c
blobacfc179c77eca0cf5a27703aeb819c3ed3464889
1 /* network.c
2 * Survex network reduction - find patterns and apply network reductions
3 * Copyright (C) 1991-2002,2005 Olly Betts
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
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #if 0
21 #define DEBUG_INVALID 1
22 #define VALIDATE 1
23 #define DUMP_NETWORK 1
24 #endif
25
26 #ifdef HAVE_CONFIG_H
27 # include <config.h>
28 #endif
29
30 #include "validate.h"
31 #include "debug.h"
32 #include "cavern.h"
33 #include "message.h"
34 #include "netbits.h"
35 #include "network.h"
36 #include "out.h"
37
38 /* type field isn't vital - join3 is unused except for deltastar, so
39 * we can set its value to indicate which type this is:
40 * join3 == NULL for noose, join3 == join1 for ||, otherwise D* */
41 #ifdef EXPLICIT_STACKRED_TYPE
42 #define SET_NOOSE(SR) (SR)->type = 1
43 #define IS_NOOSE(SR) ((SR)->type == 1)
44 #define SET_PARALLEL(SR) (SR)->type = 0
45 #define IS_PARALLEL(SR) ((SR)->type == 0)
46 #define SET_DELTASTAR(SR) (SR)->type = 2
47 #define IS_DELTASTAR(SR) ((SR)->type == 2)
48 #else
49 #define IS_NOOSE(SR) ((SR)->join3 == NULL)
50 #define SET_NOOSE(SR) (SR)->join3 = NULL
51 #define IS_PARALLEL(SR) ((SR)->join3 == (SR)->join1)
52 #define SET_PARALLEL(SR) (SR)->join3 = (SR)->join1
53 #define IS_DELTASTAR(SR) (!IS_NOOSE(SR) && !IS_PARALLEL(SR))
54 #define SET_DELTASTAR(SR) NOP
55 #endif
56
57 typedef struct StackRed {
58 struct Link *join1, *join2, *join3;
59 #ifdef EXPLICIT_STACKRED_TYPE
60 int type; /* 1 => noose, 0 => parallel legs, 2 => delta-star */
61 #endif
62 struct StackRed *next;
63 } stackRed;
64
65 static stackRed *ptrRed; /* Ptr to TRaverse linked list for C*-*< , -*=*- */
66
67 /* can be altered by -z<letters> on command line */
68 unsigned long optimize = BITA('l') | BITA('p') | BITA('d');
69 /* Lollipops, Parallel legs, Iterate mx, Delta* */
70
71 extern void
72 remove_subnets(void)
73 {
74 node *stn, *stn2, *stn3, *stn4;
75 int dirn, dirn2, dirn3, dirn4;
76 stackRed *trav;
77 linkfor *newleg, *newleg2;
78 bool fMore = fTrue;
79
80 ptrRed = NULL;
81
82 out_current_action(msg(/*Simplifying network*/129));
83
84 while (fMore) {
85 fMore = fFalse;
86 if (optimize & BITA('l')) {
87 #if PRINT_NETBITS
88 printf("replacing lollipops\n");
89 #endif
90 /* _
91 * ( )
92 * * stn
93 * |
94 * * stn2
95 * stn /|
96 * 4 * * stn3 --> stn4 *-* stn3
97 * : : : :
98 */
99 /* NB can have non-fixed 0 nodes */
100 FOR_EACH_STN(stn, stnlist) {
101 if (!fixed(stn) && three_node(stn)) {
102 dirn = -1;
103 if (stn->leg[1]->l.to == stn) dirn++;
104 if (stn->leg[0]->l.to == stn) dirn += 2;
105 if (dirn < 0) continue;
106
107 stn2 = stn->leg[dirn]->l.to;
108 if (fixed(stn2)) continue;
109
110 SVX_ASSERT(three_node(stn2));
111
112 dirn2 = reverse_leg_dirn(stn->leg[dirn]);
113 dirn2 = (dirn2 + 1) % 3;
114 stn3 = stn2->leg[dirn2]->l.to;
115 if (stn2 == stn3) continue; /* dumb-bell - leave alone */
116
117 dirn3 = reverse_leg_dirn(stn2->leg[dirn2]);
118
119 trav = osnew(stackRed);
120 newleg2 = (linkfor*)osnew(linkrev);
121
122 newleg = copy_link(stn3->leg[dirn3]);
123
124 dirn2 = (dirn2 + 1) % 3;
125 stn4 = stn2->leg[dirn2]->l.to;
126 dirn4 = reverse_leg_dirn(stn2->leg[dirn2]);
127 #if 0
128 printf("Noose found with stn...stn4 = \n");
129 print_prefix(stn->name); putnl();
130 print_prefix(stn2->name); putnl();
131 print_prefix(stn3->name); putnl();
132 print_prefix(stn4->name); putnl();
133 #endif
134
135 addto_link(newleg, stn2->leg[dirn2]);
136
137 /* remove stn and stn2 */
138 remove_stn_from_list(&stnlist, stn);
139 remove_stn_from_list(&stnlist, stn2);
140
141 /* stack noose and replace with a leg between stn3 and stn4 */
142 trav->join1 = stn3->leg[dirn3];
143 newleg->l.to = stn4;
144 newleg->l.reverse = dirn4 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
145
146 trav->join2 = stn4->leg[dirn4];
147 newleg2->l.to = stn3;
148 newleg2->l.reverse = dirn3 | FLAG_REPLACEMENTLEG;
149
150 stn3->leg[dirn3] = newleg;
151 stn4->leg[dirn4] = newleg2;
152
153 trav->next = ptrRed;
154 SET_NOOSE(trav);
155 #if PRINT_NETBITS
156 printf("remove noose\n");
157 #endif
158 ptrRed = trav;
159 fMore = fTrue;
160 }
161 }
162 }
163
164 if (optimize & BITA('p')) {
165 #if PRINT_NETBITS
166 printf("replacing parallel legs\n");
167 #endif
168 FOR_EACH_STN(stn, stnlist) {
169 /*
170 * :
171 * * stn3
172 * | :
173 * * stn * stn3
174 * ( ) -> |
175 * * stn2 * stn4
176 * | :
177 * * stn4
178 * :
179 */
180 if (!fixed(stn) && three_node(stn)) {
181 stn2 = stn->leg[0]->l.to;
182 if (stn2 == stn->leg[1]->l.to) {
183 dirn = 2;
184 } else if (stn2 == stn->leg[2]->l.to) {
185 dirn = 1;
186 } else {
187 if (stn->leg[1]->l.to != stn->leg[2]->l.to) continue;
188 stn2 = stn->leg[1]->l.to;
189 dirn = 0;
190 }
191
192 /* stn == stn2 => noose */
193 if (fixed(stn2) || stn == stn2) continue;
194
195 SVX_ASSERT(three_node(stn2));
196
197 stn3 = stn->leg[dirn]->l.to;
198 /* 3 parallel legs (=> nothing else) so leave */
199 if (stn3 == stn2) continue;
200
201 dirn3 = reverse_leg_dirn(stn->leg[dirn]);
202 dirn2 = (0 + 1 + 2 - reverse_leg_dirn(stn->leg[(dirn + 1) % 3])
203 - reverse_leg_dirn(stn->leg[(dirn + 2) % 3]));
204
205 stn4 = stn2->leg[dirn2]->l.to;
206 dirn4 = reverse_leg_dirn(stn2->leg[dirn2]);
207
208 trav = osnew(stackRed);
209
210 newleg = copy_link(stn->leg[(dirn + 1) % 3]);
211 /* use newleg2 for scratch */
212 newleg2 = copy_link(stn->leg[(dirn + 2) % 3]);
213 {
214 #ifdef NO_COVARIANCES
215 vars sum;
216 var prod;
217 delta temp, temp2;
218 addss(&sum, &newleg->v, &newleg2->v);
219 SVX_ASSERT2(!fZeros(&sum), "loop of zero variance found");
220 mulss(&prod, &newleg->v, &newleg2->v);
221 mulsd(&temp, &newleg2->v, &newleg->d);
222 mulsd(&temp2, &newleg->v, &newleg2->d);
223 adddd(&temp, &temp, &temp2);
224 divds(&newleg->d, &temp, &sum);
225 sdivvs(&newleg->v, &prod, &sum);
226 #else
227 svar inv1, inv2, sum;
228 delta temp, temp2;
229 /* if leg one is an equate, we can just ignore leg two
230 * whatever it is */
231 if (invert_svar(&inv1, &newleg->v)) {
232 if (invert_svar(&inv2, &newleg2->v)) {
233 addss(&sum, &inv1, &inv2);
234 if (!invert_svar(&newleg->v, &sum)) {
235 BUG("matrix singular in parallel legs replacement");
236 }
237
238 mulsd(&temp, &inv1, &newleg->d);
239 mulsd(&temp2, &inv2, &newleg2->d);
240 adddd(&temp, &temp, &temp2);
241 mulsd(&newleg->d, &newleg->v, &temp);
242 } else {
243 /* leg two is an equate, so just ignore leg 1 */
244 linkfor *tmpleg;
245 tmpleg = newleg;
246 newleg = newleg2;
247 newleg2 = tmpleg;
248 }
249 }
250 #endif
251 }
252 osfree(newleg2);
253 newleg2 = (linkfor*)osnew(linkrev);
254
255 addto_link(newleg, stn2->leg[dirn2]);
256 addto_link(newleg, stn3->leg[dirn3]);
257
258 #if 0
259 printf("Parallel found with stn...stn4 = \n");
260 (dump_node)(stn); (dump_node)(stn2); (dump_node)(stn3); (dump_node)(stn4);
261 printf("dirns = %d %d %d %d\n", dirn, dirn2, dirn3, dirn4);
262 #endif
263 SVX_ASSERT2(stn3->leg[dirn3]->l.to == stn, "stn3 end of || doesn't recip");
264 SVX_ASSERT2(stn4->leg[dirn4]->l.to == stn2, "stn4 end of || doesn't recip");
265 SVX_ASSERT2(stn->leg[(dirn+1)%3]->l.to == stn2 && stn->leg[(dirn + 2) % 3]->l.to == stn2, "|| legs aren't");
266
267 /* remove stn and stn2 (already discarded triple parallel) */
268 /* so stn!=stn4 <=> stn2!=stn3 */
269 remove_stn_from_list(&stnlist, stn);
270 remove_stn_from_list(&stnlist, stn2);
271
272 /* stack parallel and replace with a leg between stn3 and stn4 */
273 trav->join1 = stn3->leg[dirn3];
274 newleg->l.to = stn4;
275 newleg->l.reverse = dirn4 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
276
277 trav->join2 = stn4->leg[dirn4];
278 newleg2->l.to = stn3;
279 newleg2->l.reverse = dirn3 | FLAG_REPLACEMENTLEG;
280
281 stn3->leg[dirn3] = newleg;
282 stn4->leg[dirn4] = newleg2;
283
284 trav->next = ptrRed;
285 SET_PARALLEL(trav);
286 #if PRINT_NETBITS
287 printf("remove parallel\n");
288 #endif
289 ptrRed = trav;
290 fMore = fTrue;
291 }
292 }
293 }
294
295 if (optimize & BITA('d')) {
296 node *stn5, *stn6;
297 int dirn5, dirn6, dirn0;
298 linkfor *legAB, *legBC, *legCA;
299 #if PRINT_NETBITS
300 printf("replacing deltas with stars\n");
301 #endif
302 FOR_EACH_STN(stn, stnlist) {
303 /* printf("*");*/
304 /*
305 * :
306 * * stn5 :
307 * | * stn5
308 * * stn2 |
309 * / \ -> O stnZ
310 * stn *---* stn3 / \
311 * / \ stn4 * * stn6
312 * stn4 * * stn6 : :
313 * : :
314 */
315 if (!fixed(stn) && three_node(stn)) {
316 for (dirn0 = 0; ; dirn0++) {
317 if (dirn0 >= 3) goto nodeltastar; /* continue outer loop */
318 dirn = dirn0;
319 stn2 = stn->leg[dirn]->l.to;
320 if (fixed(stn2) || stn2 == stn) continue;
321 dirn2 = reverse_leg_dirn(stn->leg[dirn]);
322 dirn2 = (dirn2 + 1) % 3;
323 stn3 = stn2->leg[dirn2]->l.to;
324 if (fixed(stn3) || stn3 == stn || stn3 == stn2)
325 goto nextdirn2;
326 dirn3 = reverse_leg_dirn(stn2->leg[dirn2]);
327 dirn3 = (dirn3 + 1) % 3;
328 if (stn3->leg[dirn3]->l.to == stn) {
329 legAB = copy_link(stn->leg[dirn]);
330 legBC = copy_link(stn2->leg[dirn2]);
331 legCA = copy_link(stn3->leg[dirn3]);
332 dirn = 0 + 1 + 2 - dirn - reverse_leg_dirn(stn3->leg[dirn3]);
333 dirn2 = (dirn2 + 1) % 3;
334 dirn3 = (dirn3 + 1) % 3;
335 } else if (stn3->leg[(dirn3 + 1) % 3]->l.to == stn) {
336 legAB = copy_link(stn->leg[dirn]);
337 legBC = copy_link(stn2->leg[dirn2]);
338 legCA = copy_link(stn3->leg[(dirn3 + 1) % 3]);
339 dirn = (0 + 1 + 2 - dirn
340 - reverse_leg_dirn(stn3->leg[(dirn3 + 1) % 3]));
341 dirn2 = (dirn2 + 1) % 3;
342 break;
343 } else {
344 nextdirn2:;
345 dirn2 = (dirn2 + 1) % 3;
346 stn3 = stn2->leg[dirn2]->l.to;
347 if (fixed(stn3) || stn3 == stn || stn3 == stn2) continue;
348 dirn3 = reverse_leg_dirn(stn2->leg[dirn2]);
349 dirn3 = (dirn3 + 1) % 3;
350 if (stn3->leg[dirn3]->l.to == stn) {
351 legAB = copy_link(stn->leg[dirn]);
352 legBC = copy_link(stn2->leg[dirn2]);
353 legCA = copy_link(stn3->leg[dirn3]);
354 dirn = (0 + 1 + 2 - dirn
355 - reverse_leg_dirn(stn3->leg[dirn3]));
356 dirn2 = (dirn2 + 2) % 3;
357 dirn3 = (dirn3 + 1) % 3;
358 break;
359 } else if (stn3->leg[(dirn3 + 1) % 3]->l.to == stn) {
360 legAB = copy_link(stn->leg[dirn]);
361 legBC = copy_link(stn2->leg[dirn2]);
362 legCA = copy_link(stn3->leg[(dirn3 + 1) % 3]);
363 dirn = (0 + 1 + 2 - dirn
364 - reverse_leg_dirn(stn3->leg[(dirn3 + 1) % 3]));
365 dirn2 = (dirn2 + 2) % 3;
366 break;
367 }
368 }
369 }
370
371 SVX_ASSERT(three_node(stn2));
372 SVX_ASSERT(three_node(stn3));
373
374 stn4 = stn->leg[dirn]->l.to;
375 stn5 = stn2->leg[dirn2]->l.to;
376 stn6 = stn3->leg[dirn3]->l.to;
377
378 if (stn4 == stn2 || stn4 == stn3 || stn5 == stn3) break;
379
380 dirn4 = reverse_leg_dirn(stn->leg[dirn]);
381 dirn5 = reverse_leg_dirn(stn2->leg[dirn2]);
382 dirn6 = reverse_leg_dirn(stn3->leg[dirn3]);
383 #if 0
384 printf("delta-star, stn ... stn6 are:\n");
385 (dump_node)(stn);
386 (dump_node)(stn2);
387 (dump_node)(stn3);
388 (dump_node)(stn4);
389 (dump_node)(stn5);
390 (dump_node)(stn6);
391 #endif
392 SVX_ASSERT(stn4->leg[dirn4]->l.to == stn);
393 SVX_ASSERT(stn5->leg[dirn5]->l.to == stn2);
394 SVX_ASSERT(stn6->leg[dirn6]->l.to == stn3);
395
396 trav = osnew(stackRed);
397 {
398 linkfor *legAZ, *legBZ, *legCZ;
399 node *stnZ;
400 prefix *nameZ;
401 svar invAB, invBC, invCA, tmp, sum, inv;
402 var vtmp;
403 svar sumAZBZ, sumBZCZ, sumCZAZ;
404 delta temp, temp2;
405
406 /* FIXME: ought to handle cases when some legs are
407 * equates, but handle as a special case maybe? */
408 if (!invert_svar(&invAB, &legAB->v)) break;
409 if (!invert_svar(&invBC, &legBC->v)) break;
410 if (!invert_svar(&invCA, &legCA->v)) break;
411
412 addss(&sum, &legBC->v, &legCA->v);
413 addss(&tmp, &sum, &legAB->v);
414 if (!invert_svar(&inv, &tmp)) {
415 /* impossible - loop of zero variance */
416 BUG("loop of zero variance found");
417 }
418
419 legAZ = osnew(linkfor);
420 legBZ = osnew(linkfor);
421 legCZ = osnew(linkfor);
422
423 /* AZBZ */
424 /* done above: addvv(&sum, &legBC->v, &legCA->v); */
425 mulss(&vtmp, &sum, &inv);
426 smulvs(&sumAZBZ, &vtmp, &legAB->v);
427
428 adddd(&temp, &legBC->d, &legCA->d);
429 divds(&temp2, &temp, &sum);
430 mulsd(&temp, &invAB, &legAB->d);
431 subdd(&temp, &temp2, &temp);
432 mulsd(&legBZ->d, &sumAZBZ, &temp);
433
434 /* leg vectors after transform are determined up to
435 * a constant addition, so arbitrarily fix AZ = 0 */
436 legAZ->d[2] = legAZ->d[1] = legAZ->d[0] = 0;
437
438 /* BZCZ */
439 addss(&sum, &legCA->v, &legAB->v);
440 mulss(&vtmp, &sum, &inv);
441 smulvs(&sumBZCZ, &vtmp, &legBC->v);
442
443 /* CZAZ */
444 addss(&sum, &legAB->v, &legBC->v);
445 mulss(&vtmp, &sum, &inv);
446 smulvs(&sumCZAZ, &vtmp, &legCA->v);
447
448 adddd(&temp, &legAB->d, &legBC->d);
449 divds(&temp2, &temp, &sum);
450 mulsd(&temp, &invCA, &legCA->d);
451 /* NB: swapped arguments to negate answer for legCZ->d */
452 subdd(&temp, &temp, &temp2);
453 mulsd(&legCZ->d, &sumCZAZ, &temp);
454
455 osfree(legAB);
456 osfree(legBC);
457 osfree(legCA);
458
459 /* Now add two, subtract third, and scale by 0.5 */
460 addss(&sum, &sumAZBZ, &sumCZAZ);
461 subss(&sum, &sum, &sumBZCZ);
462 mulsc(&legAZ->v, &sum, 0.5);
463
464 addss(&sum, &sumBZCZ, &sumAZBZ);
465 subss(&sum, &sum, &sumCZAZ);
466 mulsc(&legBZ->v, &sum, 0.5);
467
468 addss(&sum, &sumCZAZ, &sumBZCZ);
469 subss(&sum, &sum, &sumAZBZ);
470 mulsc(&legCZ->v, &sum, 0.5);
471
472 nameZ = osnew(prefix);
473 nameZ->pos = osnew(pos);
474 nameZ->ident = NULL;
475 nameZ->shape = 3;
476 stnZ = osnew(node);
477 stnZ->name = nameZ;
478 nameZ->stn = stnZ;
479 nameZ->up = NULL;
480 nameZ->min_export = nameZ->max_export = 0;
481 unfix(stnZ);
482 add_stn_to_list(&stnlist, stnZ);
483 legAZ->l.to = stnZ;
484 legAZ->l.reverse = 0 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
485 legBZ->l.to = stnZ;
486 legBZ->l.reverse = 1 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
487 legCZ->l.to = stnZ;
488 legCZ->l.reverse = 2 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
489 stnZ->leg[0] = (linkfor*)osnew(linkrev);
490 stnZ->leg[1] = (linkfor*)osnew(linkrev);
491 stnZ->leg[2] = (linkfor*)osnew(linkrev);
492 stnZ->leg[0]->l.to = stn4;
493 stnZ->leg[0]->l.reverse = dirn4;
494 stnZ->leg[1]->l.to = stn5;
495 stnZ->leg[1]->l.reverse = dirn5;
496 stnZ->leg[2]->l.to = stn6;
497 stnZ->leg[2]->l.reverse = dirn6;
498 addto_link(legAZ, stn4->leg[dirn4]);
499 addto_link(legBZ, stn5->leg[dirn5]);
500 addto_link(legCZ, stn6->leg[dirn6]);
501 /* stack stuff */
502 trav->join1 = stn4->leg[dirn4];
503 trav->join2 = stn5->leg[dirn5];
504 trav->join3 = stn6->leg[dirn6];
505 trav->next = ptrRed;
506 SET_DELTASTAR(trav);
507 #if PRINT_NETBITS
508 printf("remove delta*\n");
509 #endif
510 ptrRed = trav;
511 fMore = fTrue;
512
513 remove_stn_from_list(&stnlist, stn);
514 remove_stn_from_list(&stnlist, stn2);
515 remove_stn_from_list(&stnlist, stn3);
516 stn4->leg[dirn4] = legAZ;
517 stn5->leg[dirn5] = legBZ;
518 stn6->leg[dirn6] = legCZ;
519 }
520
521 }
522 nodeltastar:;
523 }
524 }
525
526 }
527 }
528
529 extern void
530 replace_subnets(void)
531 {
532 stackRed *ptrOld;
533 node *stn2, *stn3, *stn4;
534 int dirn2, dirn3, dirn4;
535
536 /* help to catch bad accesses */
537 stn2 = stn3 = stn4 = NULL;
538 dirn2 = dirn3 = dirn4 = 0;
539
540 out_current_action(msg(/*Calculating network*/130));
541
542 while (ptrRed != NULL) {
543 /* printf("replace_subnets() type %d\n", ptrRed->type);*/
544
545 #if PRINT_NETBITS
546 printf("replace_subnets\n");
547 if (IS_NOOSE(ptrRed)) printf("isnoose\n");
548 if (IS_PARALLEL(ptrRed)) printf("isparallel\n");
549 if (IS_DELTASTAR(ptrRed)) printf("isdelta*\n");
550 #endif
551
552 if (!IS_DELTASTAR(ptrRed)) {
553 linkfor *leg;
554 leg = ptrRed->join1; leg = reverse_leg(leg);
555 stn3 = leg->l.to; dirn3 = reverse_leg_dirn(leg);
556 leg = ptrRed->join2; leg = reverse_leg(leg);
557 stn4 = leg->l.to; dirn4 = reverse_leg_dirn(leg);
558
559 SVX_ASSERT(fixed(stn3));
560 SVX_ASSERT(fixed(stn4));
561 SVX_ASSERT(data_here(stn3->leg[dirn3]));
562 }
563
564 if (IS_NOOSE(ptrRed)) {
565 /* noose (hanging-loop) */
566 node *stn;
567 delta e;
568 linkfor *leg;
569 int zero;
570
571 SVX_ASSERT(fixed(stn3));
572 SVX_ASSERT(fixed(stn4));
573
574 leg = stn3->leg[dirn3];
575 stn2 = ptrRed->join1->l.to;
576 dirn2 = reverse_leg_dirn(ptrRed->join1);
577
578 zero = fZeros(&leg->v);
579 if (!zero) {
580 delta tmp;
581 subdd(&e, &POSD(stn4), &POSD(stn3));
582 subdd(&tmp, &e, &leg->d);
583 divds(&e, &tmp, &leg->v);
584 }
585 if (data_here(ptrRed->join1)) {
586 adddd(&POSD(stn2), &POSD(stn3), &ptrRed->join1->d);
587 if (!zero) {
588 delta tmp;
589 mulsd(&tmp, &ptrRed->join1->v, &e);
590 adddd(&POSD(stn2), &POSD(stn2), &tmp);
591 }
592 } else {
593 subdd(&POSD(stn2), &POSD(stn3), &stn2->leg[dirn2]->d);
594 if (!zero) {
595 delta tmp;
596 mulsd(&tmp, &stn2->leg[dirn2]->v, &e);
597 adddd(&POSD(stn2), &POSD(stn2), &tmp);
598 }
599 }
600 fix(stn2);
601 dirn2 = (dirn2 + 2) % 3; /* point back at stn again */
602 stn = stn2->leg[dirn2]->l.to;
603 #if 0
604 printf("Replacing noose with stn...stn4 = \n");
605 print_prefix(stn->name); putnl();
606 print_prefix(stn2->name); putnl();
607 print_prefix(stn3->name); putnl();
608 print_prefix(stn4->name); putnl();
609 #endif
610 if (data_here(stn2->leg[dirn2]))
611 adddd(&POSD(stn), &POSD(stn2), &stn2->leg[dirn2]->d);
612 else
613 subdd(&POSD(stn), &POSD(stn2), &reverse_leg(stn2->leg[dirn2])->d);
614
615 /* the "rope" of the noose is a new articulation */
616 stn2->leg[dirn2]->l.reverse |= FLAG_ARTICULATION;
617 reverse_leg(stn2->leg[dirn2])->l.reverse |= FLAG_ARTICULATION;
618
619 fix(stn);
620
621 add_stn_to_list(&stnlist, stn);
622 add_stn_to_list(&stnlist, stn2);
623
624 osfree(stn3->leg[dirn3]);
625 stn3->leg[dirn3] = ptrRed->join1;
626 osfree(stn4->leg[dirn4]);
627 stn4->leg[dirn4] = ptrRed->join2;
628 } else if (IS_PARALLEL(ptrRed)) {
629 /* parallel legs */
630 node *stn;
631 delta e, e2;
632 linkfor *leg;
633 int dirn;
634
635 stn = ptrRed->join1->l.to;
636 stn2 = ptrRed->join2->l.to;
637 SVX_ASSERT(fixed(stn3));
638 SVX_ASSERT(fixed(stn4));
639
640 dirn = reverse_leg_dirn(ptrRed->join1);
641 dirn2 = reverse_leg_dirn(ptrRed->join2);
642
643 leg = stn3->leg[dirn3];
644
645 if (leg->l.reverse & FLAG_ARTICULATION) {
646 ptrRed->join1->l.reverse |= FLAG_ARTICULATION;
647 stn->leg[dirn]->l.reverse |= FLAG_ARTICULATION;
648 ptrRed->join2->l.reverse |= FLAG_ARTICULATION;
649 stn2->leg[dirn2]->l.reverse |= FLAG_ARTICULATION;
650 }
651
652 if (fZeros(&leg->v))
653 e[0] = e[1] = e[2] = 0.0;
654 else {
655 delta tmp;
656 subdd(&e, &POSD(stn4), &POSD(stn3));
657 subdd(&tmp, &e, &leg->d);
658 divds(&e, &tmp, &leg->v);
659 }
660
661 if (data_here(ptrRed->join1)) {
662 leg = ptrRed->join1;
663 adddd(&POSD(stn), &POSD(stn3), &leg->d);
664 } else {
665 leg = stn->leg[dirn];
666 subdd(&POSD(stn), &POSD(stn3), &leg->d);
667 }
668 mulsd(&e2, &leg->v, &e);
669 adddd(&POSD(stn), &POSD(stn), &e2);
670
671 if (data_here(ptrRed->join2)) {
672 leg = ptrRed->join2;
673 adddd(&POSD(stn2), &POSD(stn4), &leg->d);
674 } else {
675 leg = stn2->leg[dirn2];
676 subdd(&POSD(stn2), &POSD(stn4), &leg->d);
677 }
678 mulsd(&e2, &leg->v, &e);
679 subdd(&POSD(stn2), &POSD(stn2), &e2);
680 fix(stn);
681 fix(stn2);
682 #if 0
683 printf("Replacing parallel with stn...stn4 = \n");
684 print_prefix(stn->name); putnl();
685 print_prefix(stn2->name); putnl();
686 print_prefix(stn3->name); putnl();
687 print_prefix(stn4->name); putnl();
688 #endif
689
690 add_stn_to_list(&stnlist, stn);
691 add_stn_to_list(&stnlist, stn2);
692
693 osfree(stn3->leg[dirn3]);
694 stn3->leg[dirn3] = ptrRed->join1;
695 osfree(stn4->leg[dirn4]);
696 stn4->leg[dirn4] = ptrRed->join2;
697 } else if (IS_DELTASTAR(ptrRed)) {
698 node *stnZ;
699 node *stn[3];
700 int dirn[3];
701 linkfor *legs[3];
702 int i;
703 linkfor *leg;
704
705 legs[0] = ptrRed->join1;
706 legs[1] = ptrRed->join2;
707 legs[2] = ptrRed->join3;
708
709 /* work out ends as we don't bother stacking them */
710 leg = reverse_leg(legs[0]);
711 stn[0] = leg->l.to;
712 dirn[0] = reverse_leg_dirn(leg);
713 stnZ = stn[0]->leg[dirn[0]]->l.to;
714 SVX_ASSERT(fixed(stnZ));
715 stn[1] = stnZ->leg[1]->l.to;
716 dirn[1] = reverse_leg_dirn(stnZ->leg[1]);
717 stn[2] = stnZ->leg[2]->l.to;
718 dirn[2] = reverse_leg_dirn(stnZ->leg[2]);
719 /*print_prefix(stnZ->name);printf(" %p\n",(void*)stnZ);*/
720
721 for (i = 0; i < 3; i++) {
722 SVX_ASSERT2(fixed(stn[i]), "stn not fixed for D*");
723
724 leg = stn[i]->leg[dirn[i]];
725
726 SVX_ASSERT2(data_here(leg), "data not on leg for D*");
727 SVX_ASSERT2(leg->l.to == stnZ, "bad sub-network for D*");
728
729 stn2 = legs[i]->l.to;
730
731 if (data_here(legs[i])) {
732 adddd(&POSD(stn2), &POSD(stn[i]), &legs[i]->d);
733 } else {
734 subdd(&POSD(stn2), &POSD(stn[i]), &reverse_leg(legs[i])->d);
735 }
736
737 if (!fZeros(&leg->v)) {
738 delta e, tmp;
739 subdd(&e, &POSD(stnZ), &POSD(stn[i]));
740 subdd(&e, &e, &leg->d);
741 divds(&tmp, &e, &leg->v);
742 if (data_here(legs[i])) {
743 mulsd(&e, &legs[i]->v, &tmp);
744 } else {
745 mulsd(&e, &reverse_leg(legs[i])->v, &tmp);
746 }
747 adddd(&POSD(stn2), &POSD(stn2), &e);
748 }
749 fix(stn2);
750 add_stn_to_list(&stnlist, stn2);
751 osfree(leg);
752 stn[i]->leg[dirn[i]] = legs[i];
753 /* transfer the articulation status of the radial legs */
754 if (stnZ->leg[i]->l.reverse & FLAG_ARTICULATION) {
755 legs[i]->l.reverse |= FLAG_ARTICULATION;
756 reverse_leg(legs[i])->l.reverse |= FLAG_ARTICULATION;
757 }
758 osfree(stnZ->leg[i]);
759 stnZ->leg[i] = NULL;
760 }
761 /*printf("---%f %f %f\n",POS(stnZ, 0), POS(stnZ, 1), POS(stnZ, 2));*/
762 remove_stn_from_list(&stnlist, stnZ);
763 osfree(stnZ->name);
764 osfree(stnZ);
765 } else {
766 BUG("ptrRed has unknown type");
767 }
768
769 ptrOld = ptrRed;
770 ptrRed = ptrRed->next;
771 osfree(ptrOld);
772 }
773 }
Survex is a software suite to aid mapping cave systems