Context Navigation

network.c @ bf9faf6

Visit:

stereo-2025

Last change on this file since bf9faf6 was bf9faf6, checked in by Olly Betts <olly@…>, 8 months ago

Fix component counting bugs

The component count was wrong in some cases, and we calculate the number
of loops using this component count, so the loop count would be wrong by
the same amount in these cases.

Property mode set to 100644

File size: 22.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: git/src/network.c @ bf9faf6

Download in other formats: