v0.14.0
Searching...
No Matches
gm_rule.c
Go to the documentation of this file.
1# include <stdlib.h>
2# include <stdio.h>
3# include <math.h>
4# include <time.h>
5
6# include "gm_rule.h"
7
8/******************************************************************************/
9
10void comp_next ( int n, int k, int a[], int *more, int *h, int *t )
11
12/******************************************************************************/
13/*
14 Purpose:
15
16 COMP_NEXT computes the compositions of the integer N into K parts.
17
18 Discussion:
19
20 A composition of the integer N into K parts is an ordered sequence
21 of K nonnegative integers which sum to N. The compositions (1,2,1)
22 and (1,1,2) are considered to be distinct.
23
24 The routine computes one composition on each call until there are no more.
25 For instance, one composition of 6 into 3 parts is
26 3+2+1, another would be 6+0+0.
27
28 On the first call to this routine, set MORE = FALSE. The routine
29 will compute the first element in the sequence of compositions, and
30 return it, as well as setting MORE = TRUE. If more compositions
31 are desired, call again, and again. Each time, the routine will
32 return with a new composition.
33
34 However, when the LAST composition in the sequence is computed
35 and returned, the routine will reset MORE to FALSE, signaling that
36 the end of the sequence has been reached.
37
38 This routine originally used a STATICE statement to maintain the
39 variables H and T. I have decided (based on an wasting an
40 entire morning trying to track down a problem) that it is safer
41 to pass these variables as arguments, even though the user should
42 never alter them. This allows this routine to safely shuffle
43 between several ongoing calculations.
44
45
46 There are 28 compositions of 6 into three parts. This routine will
47 produce those compositions in the following order:
48
49 I A
50 - ---------
51 1 6 0 0
52 2 5 1 0
53 3 4 2 0
54 4 3 3 0
55 5 2 4 0
56 6 1 5 0
57 7 0 6 0
58 8 5 0 1
59 9 4 1 1
60 10 3 2 1
61 11 2 3 1
62 12 1 4 1
63 13 0 5 1
64 14 4 0 2
65 15 3 1 2
66 16 2 2 2
67 17 1 3 2
68 18 0 4 2
69 19 3 0 3
70 20 2 1 3
71 21 1 2 3
72 22 0 3 3
73 23 2 0 4
74 24 1 1 4
75 25 0 2 4
76 26 1 0 5
77 27 0 1 5
78 28 0 0 6
79
80 Licensing:
81
83
84 Modified:
85
86 02 July 2008
87
88 Author:
89
90 Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
91 C version by John Burkardt.
92
93 Reference:
94
95 Albert Nijenhuis, Herbert Wilf,
96 Combinatorial Algorithms for Computers and Calculators,
97 Second Edition,
99 ISBN: 0-12-519260-6,
100 LC: QA164.N54.
101
102 Parameters:
103
104 Input, int N, the integer whose compositions are desired.
105
106 Input, int K, the number of parts in the composition.
107
108 Input/output, int A[K], the parts of the composition.
109
110 Input/output, int *MORE.
111 Set MORE = FALSE on first call. It will be reset to TRUE on return
112 with a new composition. Each new call returns another composition until
113 MORE is set to FALSE when the last composition has been computed
114 and returned.
115
116 Input/output, int *H, *T, two internal parameters needed for the
117 computation. The user should allocate space for these in the calling
118 program, include them in the calling sequence, but never alter them!
119*/
120{
121 int i;
122
123 if ( !( *more ) )
124 {
125 *t = n;
126 *h = 0;
127 a[0] = n;
128 for ( i = 1; i < k; i++ )
129 {
130 a[i] = 0;
131 }
132 }
133 else
134 {
135 if ( 1 < *t )
136 {
137 *h = 0;
138 }
139 *h = *h + 1;
140 *t = a[*h-1];
141 a[*h-1] = 0;
142 a[0] = *t - 1;
143 a[*h] = a[*h] + 1;
144 }
145
146 *more = ( a[k-1] != n );
147
148 return;
149}
150/******************************************************************************/
151
152void gm_rule_set ( int rule, int dim_num, int point_num, double w[],
153 double x[] )
154
155/******************************************************************************/
156/*
157 Purpose:
158
159 GM_RULE_SET sets a Grundmann-Moeller rule.
160
161 Discussion:
162
163 This is a revised version of the calculation which seeks to compute
164 the value of the weight in a cautious way that avoids intermediate
165 overflow. Thanks to John Peterson for pointing out the problem on
166 26 June 2008.
167
168 This rule returns weights and abscissas of a Grundmann-Moeller
169 quadrature rule for the DIM_NUM-dimensional unit simplex.
170
171 The dimension POINT_NUM can be determined by calling GM_RULE_SIZE.
172
173 Licensing:
174
176
177 Modified:
178
179 30 August 2012
180
181 Author:
182
183 John Burkardt
184
185 Reference:
186
187 Axel Grundmann, Michael Moeller,
188 Invariant Integration Formulas for the N-Simplex
189 by Combinatorial Methods,
190 SIAM Journal on Numerical Analysis,
191 Volume 15, Number 2, April 1978, pages 282-290.
192
193 Parameters:
194
195 Input, int RULE, the index of the rule.
196 0 <= RULE.
197
198 Input, int DIM_NUM, the spatial dimension.
199 1 <= DIM_NUM.
200
201 Input, int POINT_NUM, the number of points in the rule.
202
203 Output, double W[POINT_NUM], the weights.
204
205 Output, double X[DIM_NUM*POINT_NUM], the abscissas.
206*/
207{
208 int *beta;
209 int beta_sum;
210 int d;
211 int dim;
212 int h;
213 int i;
214 int j;
215 int j_hi;
216 int k;
217 int more;
218 int n;
219 double one_pm;
220 int s;
221 int t;
222 double weight;
223
224 s = rule;
225 d = 2 * s + 1;
226 k = 0;
227 n = dim_num;
228 one_pm = 1.0;
229
230 beta = ( int * ) malloc ( ( dim_num + 1 ) * sizeof ( int ) );
231
232 for ( i = 0; i <= s; i++ )
233 {
234 weight = ( double ) one_pm;
235
236 j_hi = i4_max ( n, i4_max ( d, d + n - i ) );
237
238 for ( j = 1; j <= j_hi; j++ )
239 {
240 if ( j <= n )
241 {
242 weight = weight * ( double ) ( j );
243 }
244 if ( j <= d )
245 {
246 weight = weight * ( double ) ( d + n - 2 * i );
247 }
248 if ( j <= 2 * s )
249 {
250 weight = weight / 2.0;
251 }
252 if ( j <= i )
253 {
254 weight = weight / ( double ) ( j );
255 }
256 if ( j <= d + n - i )
257 {
258 weight = weight / ( double ) ( j );
259 }
260 }
261
262 one_pm = - one_pm;
263
264 beta_sum = s - i;
265 more = 0;
266 h = 0;
267 t = 0;
268
269 for ( ; ; )
270 {
271 comp_next ( beta_sum, dim_num + 1, beta, &more, &h, &t );
272
273 w[k] = weight;
274 for ( dim = 0; dim < dim_num; dim++ )
275 {
276 x[dim+k*dim_num] = ( double ) ( 2 * beta[dim+1] + 1 )
277 / ( double ) ( d + n - 2 * i );
278 }
279 k = k + 1;
280
281 if ( !more )
282 {
283 break;
284 }
285 }
286 }
287
288 free ( beta );
289
290 return;
291}
292/******************************************************************************/
293
294int gm_rule_size ( int rule, int dim_num )
295
296/******************************************************************************/
297/*
298 Purpose:
299
300 GM_RULE_SIZE determines the size of a Grundmann-Moeller rule.
301
302 Discussion:
303
304 This rule returns the value of POINT_NUM, the number of points associated
305 with a GM rule of given index.
306
307 After calling this rule, the user can use the value of POINT_NUM to
308 allocate space for the weight vector as W(POINT_NUM) and the abscissa
309 vector as X(DIM_NUM,POINT_NUM), and then call GM_RULE_SET.
310
311 Licensing:
312
314
315 Modified:
316
317 30 August 2012
318
319 Author:
320
321 John Burkardt
322
323 Reference:
324
325 Axel Grundmann, Michael Moeller,
326 Invariant Integration Formulas for the N-Simplex
327 by Combinatorial Methods,
328 SIAM Journal on Numerical Analysis,
329 Volume 15, Number 2, April 1978, pages 282-290.
330
331 Parameters:
332
333 Input, int RULE, the index of the rule.
334 0 <= RULE.
335
336 Input, int DIM_NUM, the spatial dimension.
337 1 <= DIM_NUM.
338
339 Output, int GM_RULE_SIZE, the number of points in the rule.
340*/
341{
342 int arg1;
343 int point_num;
344
345 arg1 = dim_num + rule + 1;
346
347 point_num = i4_choose ( arg1, rule );
348
349 return point_num;
350}
351/******************************************************************************/
352
353int i4_choose ( int n, int k )
354
355/******************************************************************************/
356/*
357 Purpose:
358
359 I4_CHOOSE computes the binomial coefficient C(N,K).
360
361 Discussion:
362
363 The value is calculated in such a way as to avoid overflow and
364 roundoff. The calculation is done in integer arithmetic.
365
366 The formula used is:
367
368 C(N,K) = N! / ( K! * (N-K)! )
369
370 Licensing:
371
373
374 Modified:
375
376 03 May 2008
377
378 Author:
379
380 John Burkardt
381
382 Reference:
383
384 ML Wolfson, HV Wright,
385 Algorithm 160:
386 Combinatorial of M Things Taken N at a Time,
387 Communications of the ACM,
388 Volume 6, Number 4, April 1963, page 161.
389
390 Parameters:
391
392 Input, int N, K, are the values of N and K.
393
394 Output, int I4_CHOOSE, the number of combinations of N
395 things taken K at a time.
396*/
397{
398 int i;
399 int mn;
400 int mx;
401 int value;
402
403 mn = i4_min ( k, n - k );
404
405 if ( mn < 0 )
406 {
407 value = 0;
408 }
409 else if ( mn == 0 )
410 {
411 value = 1;
412 }
413 else
414 {
415 mx = i4_max ( k, n - k );
416 value = mx + 1;
417
418 for ( i = 2; i <= mn; i++ )
419 {
420 value = ( value * ( mx + i ) ) / i;
421 }
422 }
423
424 return value;
425}
426/******************************************************************************/
427
428int i4_huge ( void )
429
430/******************************************************************************/
431/*
432 Purpose:
433
434 I4_HUGE returns a "huge" I4.
435
436 Licensing:
437
439
440 Modified:
441
442 29 August 2006
443
444 Author:
445
446 John Burkardt
447
448 Parameters:
449
450 Output, int I4_HUGE, a "huge" integer.
451*/
452{
453 static int value = 2147483647;
454
455 return value;
456}
457/******************************************************************************/
458
459int i4_max ( int i1, int i2 )
460
461/******************************************************************************/
462/*
463 Purpose:
464
465 I4_MAX returns the maximum of two I4's.
466
467 Licensing:
468
470
471 Modified:
472
473 29 August 2006
474
475 Author:
476
477 John Burkardt
478
479 Parameters:
480
481 Input, int I1, I2, are two integers to be compared.
482
483 Output, int I4_MAX, the larger of I1 and I2.
484*/
485{
486 int value;
487
488 if ( i2 < i1 )
489 {
490 value = i1;
491 }
492 else
493 {
494 value = i2;
495 }
496 return value;
497}
498/******************************************************************************/
499
500int i4_min ( int i1, int i2 )
501
502/******************************************************************************/
503/*
504 Purpose:
505
506 I4_MIN returns the smaller of two I4's.
507
508 Licensing:
509
511
512 Modified:
513
514 29 August 2006
515
516 Author:
517
518 John Burkardt
519
520 Parameters:
521
522 Input, int I1, I2, two integers to be compared.
523
524 Output, int I4_MIN, the smaller of I1 and I2.
525*/
526{
527 int value;
528
529 if ( i1 < i2 )
530 {
531 value = i1;
532 }
533 else
534 {
535 value = i2;
536 }
537 return value;
538}
539/******************************************************************************/
540
541int i4_power ( int i, int j )
542
543/******************************************************************************/
544/*
545 Purpose:
546
547 I4_POWER returns the value of I^J.
548
549 Licensing:
550
552
553 Modified:
554
555 23 October 2007
556
557 Author:
558
559 John Burkardt
560
561 Parameters:
562
563 Input, int I, J, the base and the power. J should be nonnegative.
564
565 Output, int I4_POWER, the value of I^J.
566*/
567{
568 int k;
569 int value;
570
571 if ( j < 0 )
572 {
573 if ( i == 1 )
574 {
575 value = 1;
576 }
577 else if ( i == 0 )
578 {
579 fprintf ( stderr, "\n" );
580 fprintf ( stderr, "I4_POWER - Fatal error!\n" );
581 fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" );
582 exit ( 1 );
583 }
584 else
585 {
586 value = 0;
587 }
588 }
589 else if ( j == 0 )
590 {
591 if ( i == 0 )
592 {
593 fprintf ( stderr, "\n" );
594 fprintf ( stderr, "I4_POWER - Fatal error!\n" );
595 fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" );
596 exit ( 1 );
597 }
598 else
599 {
600 value = 1;
601 }
602 }
603 else if ( j == 1 )
604 {
605 value = i;
606 }
607 else
608 {
609 value = 1;
610 for ( k = 1; k <= j; k++ )
611 {
612 value = value * i;
613 }
614 }
615 return value;
616}
617/******************************************************************************/
618
619double *monomial_value ( int dim_num, int point_num, double x[], int expon[] )
620
621/******************************************************************************/
622/*
623 Purpose:
624
625 MONOMIAL_VALUE evaluates a monomial.
626
627 Discussion:
628
629 This routine evaluates a monomial of the form
630
631 product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
632
633 where the exponents are nonnegative integers. Note that
634 if the combination 0^0 is encountered, it should be treated
635 as 1.
636
637 Licensing:
638
640
641 Modified:
642
643 30 August 2012
644
645 Author:
646
647 John Burkardt
648
649 Parameters:
650
651 Input, int DIM_NUM, the spatial dimension.
652
653 Input, int POINT_NUM, the number of points at which the
654 monomial is to be evaluated.
655
656 Input, double X[DIM_NUM*POINT_NUM], the point coordinates.
657
658 Input, int EXPON[DIM_NUM], the exponents.
659
660 Output, double MONOMIAL_VALUE[POINT_NUM], the value of the monomial.
661*/
662{
663 int dim;
664 int point;
665 double *value;
666
667 value = ( double * ) malloc ( point_num * sizeof ( double ) );
668
669 for ( point = 0; point < point_num; point++ )
670 {
671 value[point] = 1.0;
672 }
673
674 for ( dim = 0; dim < dim_num; dim++ )
675 {
676 if ( 0 != expon[dim] )
677 {
678 for ( point = 0; point < point_num; point++ )
679 {
680 value[point] = value[point] * pow ( x[dim+point*dim_num], expon[dim] );
681 }
682 }
683 }
684
685 return value;
686}
687/******************************************************************************/
688
689double r8_abs ( double x )
690
691/******************************************************************************/
692/*
693 Purpose:
694
695 R8_ABS returns the absolute value of an R8.
696
697 Licensing:
698
700
701 Modified:
702
703 07 May 2006
704
705 Author:
706
707 John Burkardt
708
709 Parameters:
710
711 Input, double X, the quantity whose absolute value is desired.
712
713 Output, double R8_ABS, the absolute value of X.
714*/
715{
716 double value;
717
718 if ( 0.0 <= x )
719 {
720 value = + x;
721 }
722 else
723 {
724 value = - x;
725 }
726 return value;
727}
728/******************************************************************************/
729
730double r8_factorial ( int n )
731
732/******************************************************************************/
733/*
734 Purpose:
735
736 R8_FACTORIAL computes the factorial of N.
737
738 Discussion:
739
740 factorial ( N ) = product ( 1 <= I <= N ) I
741
742 Licensing:
743
745
746 Modified:
747
748 26 June 2008
749
750 Author:
751
752 John Burkardt
753
754 Parameters:
755
756 Input, int N, the argument of the factorial function.
757 If N is less than 1, the function value is returned as 1.
758
759 Output, double R8_FACTORIAL, the factorial of N.
760*/
761{
762 int i;
763 double value;
764
765 value = 1.0;
766
767 for ( i = 1; i <= n; i++ )
768 {
769 value = value * ( double ) ( i );
770 }
771
772 return value;
773}
774/******************************************************************************/
775
776double r8vec_dot_product ( int n, double a1[], double a2[] )
777
778/******************************************************************************/
779/*
780 Purpose:
781
782 R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's.
783
784 Licensing:
785
787
788 Modified:
789
790 26 July 2007
791
792 Author:
793
794 John Burkardt
795
796 Parameters:
797
798 Input, int N, the number of entries in the vectors.
799
800 Input, double A1[N], A2[N], the two vectors to be considered.
801
802 Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors.
803*/
804{
805 int i;
806 double value;
807
808 value = 0.0;
809 for ( i = 0; i < n; i++ )
810 {
811 value = value + a1[i] * a2[i];
812 }
813 return value;
814}
815/******************************************************************************/
816
817double *r8vec_uniform_01_new ( int n, int *seed )
818
819/******************************************************************************/
820/*
821 Purpose:
822
823 R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC.
824
825 Discussion:
826
827 This routine implements the recursion
828
829 seed = 16807 * seed mod ( 2^31 - 1 )
830 unif = seed / ( 2^31 - 1 )
831
832 The integer arithmetic never requires more than 32 bits,
833 including a sign bit.
834
835 Licensing:
836
838
839 Modified:
840
841 19 August 2004
842
843 Author:
844
845 John Burkardt
846
847 Reference:
848
849 Paul Bratley, Bennett Fox, Linus Schrage,
850 A Guide to Simulation,
851 Second Edition,
852 Springer, 1987,
853 ISBN: 0387964673,
854 LC: QA76.9.C65.B73.
855
856 Bennett Fox,
857 Algorithm 647:
858 Implementation and Relative Efficiency of Quasirandom
859 Sequence Generators,
860 ACM Transactions on Mathematical Software,
861 Volume 12, Number 4, December 1986, pages 362-376.
862
863 Pierre L'Ecuyer,
864 Random Number Generation,
865 in Handbook of Simulation,
866 edited by Jerry Banks,
867 Wiley, 1998,
868 ISBN: 0471134031,
869 LC: T57.62.H37.
870
871 Peter Lewis, Allen Goodman, James Miller,
872 A Pseudo-Random Number Generator for the System/360,
873 IBM Systems Journal,
874 Volume 8, Number 2, 1969, pages 136-143.
875
876 Parameters:
877
878 Input, int N, the number of entries in the vector.
879
880 Input/output, int *SEED, a seed for the random number generator.
881
882 Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values.
883*/
884{
885 int i;
886 int i4_huge = 2147483647;
887 int k;
888 double *r;
889
890 if ( *seed == 0 )
891 {
892 fprintf ( stderr, "\n" );
893 fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" );
894 fprintf ( stderr, " Input value of SEED = 0.\n" );
895 exit ( 1 );
896 }
897
898 r = ( double * ) malloc ( n * sizeof ( double ) );
899
900 for ( i = 0; i < n; i++ )
901 {
902 k = *seed / 127773;
903
904 *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
905
906 if ( *seed < 0 )
907 {
908 *seed = *seed + i4_huge;
909 }
910
911 r[i] = ( double ) ( *seed ) * 4.656612875E-10;
912 }
913
914 return r;
915}
916/******************************************************************************/
917
918double simplex_unit_monomial_int ( int dim_num, int expon[] )
919
920/******************************************************************************/
921/*
922 Purpose:
923
924 SIMPLEX_UNIT_MONOMIAL_INT integrates a monomial over a simplex.
925
926 Discussion:
927
928 This routine evaluates a monomial of the form
929
930 product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
931
932 where the exponents are nonnegative integers. Note that
933 if the combination 0^0 is encountered, it should be treated
934 as 1.
935
936 Licensing:
937
939
940 Modified:
941
942 30 August 2012
943
944 Author:
945
946 John Burkardt
947
948 Parameters:
949
950 Input, int DIM_NUM, the spatial dimension.
951
952 Input, int EXPON[DIM_NUM], the exponents.
953
954 Output, double SIMPLEX_UNIT_MONOMIAL_INT, the value of the integral
955 of the monomial.
956*/
957{
958 int dim;
959 int i;
960 int k;
961 double value;
962
963 value = 1.0;
964
965 k = 0;
966
967 for ( dim = 0; dim < dim_num; dim++ )
968 {
969 for ( i = 1; i <= expon[dim]; i++ )
970 {
971 k = k + 1;
972 value = value * ( double ) ( i ) / ( double ) ( k );
973 }
974 }
975
976 for ( dim = 0; dim < dim_num; dim++ )
977 {
978 k = k + 1;
979 value = value / ( double ) ( k );
980 }
981
982 return value;
983}
984/******************************************************************************/
985
986double simplex_unit_monomial_quadrature ( int dim_num, int expon[],
987 int point_num, double x[], double w[] )
988
989/******************************************************************************/
990/*
991 Purpose:
992
994
995 Licensing:
996
998
999 Modified:
1000
1001 30 August 2012
1002
1003 Author:
1004
1005 John Burkardt
1006
1007 Parameters:
1008
1009 Input, int DIM_NUM, the spatial dimension.
1010
1011 Input, int EXPON[DIM_NUM], the exponents.
1012
1013 Input, int POINT_NUM, the number of points in the rule.
1014
1015 Input, double X[DIM_NUM*POINT_NUM], the quadrature points.
1016
1017 Input, double W[POINT_NUM], the quadrature weights.
1018
1020*/
1021{
1022 double exact = 1.0;
1025 double scale;
1026 double *value;
1027 double volume;
1028/*
1029 Get the exact value of the integral of the unscaled monomial.
1030*/
1031 scale = simplex_unit_monomial_int ( dim_num, expon );
1032/*
1033 Evaluate the monomial at the quadrature points.
1034*/
1035 value = monomial_value ( dim_num, point_num, x, expon );
1036/*
1037 Compute the weighted sum and divide by the exact value.
1038*/
1039 volume = simplex_unit_volume ( dim_num );
1040 quad = volume * r8vec_dot_product ( point_num, w, value ) / scale;
1041
1042 free ( value );
1043/*
1044 Error:
1045*/
1047
1049}
1050/******************************************************************************/
1051
1052double *simplex_unit_sample ( int dim_num, int n, int *seed )
1053
1054/******************************************************************************/
1055/*
1056 Purpose:
1057
1058 SIMPLEX_UNIT_SAMPLE returns uniformly random points from a general simplex.
1059
1060 Discussion:
1061
1062 The interior of the unit DIM_NUM dimensional simplex is the set of
1063 points X(1:DIM_NUM) such that each X(I) is nonnegative, and
1064 sum(X(1:DIM_NUM)) <= 1.
1065
1066 This routine is valid for any spatial dimension DIM_NUM.
1067
1068 Licensing:
1069
1071
1072 Modified:
1073
1074 30 August 2012
1075
1076 Author:
1077
1078 John Burkardt
1079
1080 Reference:
1081
1082 Reuven Rubinstein,
1083 Monte Carlo Optimization, Simulation, and Sensitivity
1084 of Queueing Networks,
1085 Krieger, 1992,
1086 ISBN: 0894647644,
1087 LC: QA298.R79.
1088
1089 Parameters:
1090
1091 Input, int DIM_NUM, the dimension of the space.
1092
1093 Input, int N, the number of points.
1094
1095 Input/output, int *SEED, a seed for the random number generator.
1096
1097 Output, double UNIFORM_IN_SIMPLEX01_MAP[DIM_NUM*N], the points.
1098*/
1099{
1100 double *e;
1101 int i;
1102 int j;
1103 double total;
1104 double *x;
1105/*
1106 The construction begins by sampling DIM_NUM+1 points from the
1107 exponential distribution with parameter 1.
1108*/
1109 x = ( double * ) malloc ( dim_num * n * sizeof ( double ) );
1110
1111 for ( j = 0; j < n; j++ )
1112 {
1113 e = r8vec_uniform_01_new ( dim_num+1, seed );
1114
1115 for ( i = 0; i <= dim_num; i++ )
1116 {
1117 e[i] = -log ( e[i] );
1118 }
1119
1120 total = 0.0;
1121 for ( i = 0; i <= dim_num; i++ )
1122 {
1123 total = total + e[i];
1124 }
1125
1126 for ( i = 0; i < dim_num; i++ )
1127 {
1128 x[i+dim_num*j] = e[i] / total;
1129 }
1130 free ( e );
1131 }
1132
1133 return x;
1134}
1135/******************************************************************************/
1136
1137double *simplex_unit_to_general ( int dim_num, int point_num, double t[],
1138 double ref[] )
1139
1140/******************************************************************************/
1141/*
1142 Purpose:
1143
1144 SIMPLEX_UNIT_TO_GENERAL maps the unit simplex to a general simplex.
1145
1146 Discussion:
1147
1148 Given that the unit simplex has been mapped to a general simplex
1149 with vertices T, compute the images in T, under the same linear
1150 mapping, of points whose coordinates in the unit simplex are REF.
1151
1152 The vertices of the unit simplex are listed as suggested in the
1153 following:
1154
1155 (0,0,0,...,0)
1156 (1,0,0,...,0)
1157 (0,1,0,...,0)
1158 (0,0,1,...,0)
1159 (...........)
1160 (0,0,0,...,1)
1161
1162 Thanks to Andrei ("spiritualworlds") for pointing out a mistake in the
1163 previous implementation of this routine, 02 March 2008.
1164
1165 Licensing:
1166
1168
1169 Modified:
1170
1171 30 August 2012
1172
1173 Author:
1174
1175 John Burkardt
1176
1177 Parameters:
1178
1179 Input, int DIM_NUM, the spatial dimension.
1180
1181 Input, int POINT_NUM, the number of points to transform.
1182
1183 Input, double T[DIM_NUM*(DIM_NUM+1)], the vertices of the
1184 general simplex.
1185
1186 Input, double REF[DIM_NUM*POINT_NUM], points in the
1187 reference triangle.
1188
1189 Output, double SIMPLEX_UNIT_TO_GENERAL[DIM_NUM*POINT_NUM],
1190 corresponding points in the physical triangle.
1191*/
1192{
1193 int dim;
1194 double *phy;
1195 int point;
1196 int vertex;
1197
1198 phy = ( double * ) malloc ( dim_num * point_num * sizeof ( double ) );
1199//
1200// The image of each point is initially the image of the origin.
1201//
1202// Insofar as the pre-image differs from the origin in a given vertex
1203// direction, add that proportion of the difference between the images
1204// of the origin and the vertex.
1205//
1206 for ( point = 0; point < point_num; point++ )
1207 {
1208 for ( dim = 0; dim < dim_num; dim++ )
1209 {
1210 phy[dim+point*dim_num] = t[dim+0*dim_num];
1211
1212 for ( vertex = 1; vertex < dim_num + 1; vertex++ )
1213 {
1214 phy[dim+point*dim_num] = phy[dim+point*dim_num]
1215 + ( t[dim+vertex*dim_num] - t[dim+0*dim_num] ) * ref[vertex-1+point*dim_num];
1216 }
1217 }
1218 }
1219
1220 return phy;
1221}
1222/******************************************************************************/
1223
1224double simplex_unit_volume ( int dim_num )
1225
1226/******************************************************************************/
1227/*
1228 Purpose:
1229
1230 SIMPLEX_UNIT_VOLUME computes the volume of the unit simplex.
1231
1232 Discussion:
1233
1234 The formula is simple: volume = 1/DIM_NUM!.
1235
1236 Licensing:
1237
1239
1240 Modified:
1241
1242 04 September 2003
1243
1244 Author:
1245
1246 John Burkardt
1247
1248 Parameters:
1249
1250 Input, int DIM_NUM, the dimension of the space.
1251
1252 Output, double SIMPLEX_UNIT_VOLUME, the volume of the cone.
1253*/
1254{
1255 int i;
1256 double volume;
1257
1258 volume = 1.0;
1259 for ( i = 1; i <= dim_num; i++ )
1260 {
1261 volume = volume / ( ( double ) i );
1262 }
1263
1264 return volume;
1265}
1266/******************************************************************************/
1267
1268void timestamp ( void )
1269
1270/******************************************************************************/
1271/*
1272 Purpose:
1273
1274 TIMESTAMP prints the current YMDHMS date as a time stamp.
1275
1276 Example:
1277
1278 31 May 2001 09:45:54 AM
1279
1280 Licensing:
1281
1283
1284 Modified:
1285
1286 24 September 2003
1287
1288 Author:
1289
1290 John Burkardt
1291
1292 Parameters:
1293
1294 None
1295*/
1296{
1297# define TIME_SIZE 40
1298
1299 static char time_buffer[TIME_SIZE];
1300 const struct tm *tm;
1301 size_t len;
1302 time_t now;
1303
1304 (void)(len);
1305
1306 now = time ( NULL );
1307 tm = localtime ( &now );
1308
1309 len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
1310
1311 fprintf ( stdout, "%s\n", time_buffer );
1312
1313 return;
1314# undef TIME_SIZE
1315}
1316
constexpr double a
const int dim
FTensor::Index< 'n', SPACE_DIM > n
void timestamp(void)
Definition: gm_rule.c:1268
int i4_max(int i1, int i2)
Definition: gm_rule.c:459
double r8_factorial(int n)
Definition: gm_rule.c:730
double * simplex_unit_sample(int dim_num, int n, int *seed)
Definition: gm_rule.c:1052
double simplex_unit_monomial_int(int dim_num, int expon[])
Definition: gm_rule.c:918
double simplex_unit_volume(int dim_num)
Definition: gm_rule.c:1224
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
double r8_abs(double x)
Definition: gm_rule.c:689
double * r8vec_uniform_01_new(int n, int *seed)
Definition: gm_rule.c:817
int i4_power(int i, int j)
Definition: gm_rule.c:541
double simplex_unit_monomial_quadrature(int dim_num, int expon[], int point_num, double x[], double w[])
Definition: gm_rule.c:986
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: gm_rule.c:776
int i4_choose(int n, int k)
Definition: gm_rule.c:353
#define TIME_SIZE
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
double * monomial_value(int dim_num, int point_num, double x[], int expon[])
Definition: gm_rule.c:619
int i4_huge(void)
Definition: gm_rule.c:428
void comp_next(int n, int k, int a[], int *more, int *h, int *t)
Definition: gm_rule.c:10
double * simplex_unit_to_general(int dim_num, int point_num, double t[], double ref[])
Definition: gm_rule.c:1137
int i4_min(int i1, int i2)
Definition: gm_rule.c:500
constexpr double a2
constexpr double a1
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double h
double scale
Definition: plastic.cpp:170
constexpr double t
plate stiffness
Definition: plate.cpp:59