v0.14.0
triangle_felippa_rule.c
Go to the documentation of this file.
1 # include <stdlib.h>
2 # include <stdio.h>
3 # include <time.h>
4 # include <std::string.h>
5 # include <math.h>
6 
7 # include "triangle_felippa_rule.h"
8 
9 /******************************************************************************/
10 
11 void comp_next ( int n, int k, int a[], int *more, int *h, int *t )
12 
13 /******************************************************************************/
14 /*
15  Purpose:
16 
17  COMP_NEXT computes the compositions of the integer N into K parts.
18 
19  Discussion:
20 
21  A composition of the integer N into K parts is an ordered sequence
22  of K nonnegative integers which sum to N. The compositions (1,2,1)
23  and (1,1,2) are considered to be distinct.
24 
25  The routine computes one composition on each call until there are no more.
26  For instance, one composition of 6 into 3 parts is
27  3+2+1, another would be 6+0+0.
28 
29  On the first call to this routine, set MORE = FALSE. The routine
30  will compute the first element in the sequence of compositions, and
31  return it, as well as setting MORE = TRUE. If more compositions
32  are desired, call again, and again. Each time, the routine will
33  return with a new composition.
34 
35  However, when the LAST composition in the sequence is computed
36  and returned, the routine will reset MORE to FALSE, signaling that
37  the end of the sequence has been reached.
38 
39  This routine originally used a STATICE statement to maintain the
40  variables H and T. I have decided (based on an wasting an
41  entire morning trying to track down a problem) that it is safer
42  to pass these variables as arguments, even though the user should
43  never alter them. This allows this routine to safely shuffle
44  between several ongoing calculations.
45 
46 
47  There are 28 compositions of 6 into three parts. This routine will
48  produce those compositions in the following order:
49 
50  I A
51  - ---------
52  1 6 0 0
53  2 5 1 0
54  3 4 2 0
55  4 3 3 0
56  5 2 4 0
57  6 1 5 0
58  7 0 6 0
59  8 5 0 1
60  9 4 1 1
61  10 3 2 1
62  11 2 3 1
63  12 1 4 1
64  13 0 5 1
65  14 4 0 2
66  15 3 1 2
67  16 2 2 2
68  17 1 3 2
69  18 0 4 2
70  19 3 0 3
71  20 2 1 3
72  21 1 2 3
73  22 0 3 3
74  23 2 0 4
75  24 1 1 4
76  25 0 2 4
77  26 1 0 5
78  27 0 1 5
79  28 0 0 6
80 
81  Licensing:
82 
83  This code is distributed under the GNU LGPL license.
84 
85  Modified:
86 
87  02 July 2008
88 
89  Author:
90 
91  Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
92  C version by John Burkardt.
93 
94  Reference:
95 
96  Albert Nijenhuis, Herbert Wilf,
97  Combinatorial Algorithms for Computers and Calculators,
98  Second Edition,
99  Academic Press, 1978,
100  ISBN: 0-12-519260-6,
101  LC: QA164.N54.
102 
103  Parameters:
104 
105  Input, int N, the integer whose compositions are desired.
106 
107  Input, int K, the number of parts in the composition.
108 
109  Input/output, int A[K], the parts of the composition.
110 
111  Input/output, int *MORE.
112  Set MORE = FALSE on first call. It will be reset to TRUE on return
113  with a new composition. Each new call returns another composition until
114  MORE is set to FALSE when the last composition has been computed
115  and returned.
116 
117  Input/output, int *H, *T, two internal parameters needed for the
118  computation. The user should allocate space for these in the calling
119  program, include them in the calling sequence, but never alter them!
120 */
121 {
122  int i;
123 
124  if ( !( *more ) )
125  {
126  *t = n;
127  *h = 0;
128  a[0] = n;
129  for ( i = 1; i < k; i++ )
130  {
131  a[i] = 0;
132  }
133  }
134  else
135  {
136  if ( 1 < *t )
137  {
138  *h = 0;
139  }
140  *h = *h + 1;
141  *t = a[*h-1];
142  a[*h-1] = 0;
143  a[0] = *t - 1;
144  a[*h] = a[*h] + 1;
145  }
146 
147  *more = ( a[k-1] != n );
148 
149  return;
150 }
151 /******************************************************************************/
152 
153 double *monomial_value ( int m, int n, int e[], double x[] )
154 
155 /******************************************************************************/
156 /*
157  Purpose:
158 
159  MONOMIAL_VALUE evaluates a monomial.
160 
161  Discussion:
162 
163  This routine evaluates a monomial of the form
164 
165  product ( 1 <= i <= m ) x(i)^e(i)
166 
167  The combination 0.0^0 is encountered is treated as 1.0.
168 
169  Licensing:
170 
171  This code is distributed under the GNU LGPL license.
172 
173  Modified:
174 
175  17 August 2014
176 
177  Author:
178 
179  John Burkardt
180 
181  Parameters:
182 
183  Input, int M, the spatial dimension.
184 
185  Input, int N, the number of evaluation points.
186 
187  Input, int E[M], the exponents.
188 
189  Input, double X[M*N], the point coordinates.
190 
191  Output, double MONOMIAL_VALUE[N], the monomial values.
192 */
193 {
194  int i;
195  int j;
196  double *v;
197 
198  v = ( double * ) malloc ( n * sizeof ( double ) );
199 
200  for ( j = 0; j < n; j++ )
201  {
202  v[j] = 1.0;
203  }
204 
205  for ( i = 0; i < m; i++ )
206  {
207  if ( 0 != e[i] )
208  {
209  for ( j = 0; j < n; j++ )
210  {
211  v[j] = v[j] * pow ( x[i+j*m], e[i] );
212  }
213  }
214  }
215 
216  return v;
217 }
218 /******************************************************************************/
219 
220 void r8vec_copy ( int n, double a1[], double a2[] )
221 
222 /******************************************************************************/
223 /*
224  Purpose:
225 
226  R8VEC_COPY copies an R8VEC.
227 
228  Discussion:
229 
230  An R8VEC is a vector of R8's.
231 
232  Licensing:
233 
234  This code is distributed under the GNU LGPL license.
235 
236  Modified:
237 
238  03 July 2005
239 
240  Author:
241 
242  John Burkardt
243 
244  Parameters:
245 
246  Input, int N, the number of entries in the vectors.
247 
248  Input, double A1[N], the vector to be copied.
249 
250  Input, double A2[N], the copy of A1.
251 */
252 {
253  int i;
254 
255  for ( i = 0; i < n; i++ )
256  {
257  a2[i] = a1[i];
258  }
259  return;
260 }
261 /******************************************************************************/
262 
263 double r8vec_dot_product ( int n, double a1[], double a2[] )
264 
265 /******************************************************************************/
266 /*
267  Purpose:
268 
269  R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's.
270 
271  Licensing:
272 
273  This code is distributed under the GNU LGPL license.
274 
275  Modified:
276 
277  26 July 2007
278 
279  Author:
280 
281  John Burkardt
282 
283  Parameters:
284 
285  Input, int N, the number of entries in the vectors.
286 
287  Input, double A1[N], A2[N], the two vectors to be considered.
288 
289  Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors.
290 */
291 {
292  int i;
293  double value;
294 
295  value = 0.0;
296  for ( i = 0; i < n; i++ )
297  {
298  value = value + a1[i] * a2[i];
299  }
300  return value;
301 }
302 /******************************************************************************/
303 
304 void subcomp_next ( int n, int k, int a[], int *more, int *h, int *t )
305 
306 /******************************************************************************/
307 /*
308  Purpose:
309 
310  SUBCOMP_NEXT computes the next subcomposition of N into K parts.
311 
312  Discussion:
313 
314  A composition of the integer N into K parts is an ordered sequence
315  of K nonnegative integers which sum to a value of N.
316 
317  A subcomposition of the integer N into K parts is a composition
318  of M into K parts, where 0 <= M <= N.
319 
320  Licensing:
321 
322  This code is distributed under the GNU LGPL license.
323 
324  Modified:
325 
326  27 March 2009
327 
328  Author:
329 
330  Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
331  C version by John Burkardt.
332 
333  Reference:
334 
335  Albert Nijenhuis, Herbert Wilf,
336  Combinatorial Algorithms for Computers and Calculators,
337  Second Edition,
338  Academic Press, 1978,
339  ISBN: 0-12-519260-6,
340  LC: QA164.N54.
341 
342  Parameters:
343 
344  Input, int N, the integer whose subcompositions are desired.
345 
346  Input, int K, the number of parts in the subcomposition.
347 
348  Input/output, int A[K], the parts of the subcomposition.
349 
350  Input/output, int *MORE, set by the user to start the computation,
351  and by the routine to terminate it.
352 
353  Input/output, int *H, *T, two internal parameters needed for the
354  computation. The user should allocate space for these in the calling
355  program, include them in the calling sequence, but never alter them!
356 */
357 {
358  int i;
359  static int more2 = 0;
360  static int n2 = 0;
361 /*
362  The first computation.
363 */
364  if ( !( *more ) )
365  {
366  n2 = 0;
367 
368  for ( i = 0; i < k; i++ )
369  {
370  a[i] = 0;
371  }
372  more2 = 0;
373  *h = 0;
374  *t = 0;
375 
376  *more = 1;
377  }
378 /*
379  Do the next element at the current value of N.
380 */
381  else if ( more2 )
382  {
383  comp_next ( n2, k, a, &more2, h, t );
384  }
385  else
386  {
387  more2 = 0;
388  n2 = n2 + 1;
389 
390  comp_next ( n2, k, a, &more2, h, t );
391  }
392 /*
393  Termination occurs if MORE2 = FALSE and N2 = N.
394 */
395  if ( !more2 && n2 == n )
396  {
397  *more = 0;
398  }
399 
400  return;
401 }
402 /******************************************************************************/
403 
404 //void timestamp ( ) //NOTE: it is defined in gm_rule.c
405 
406 /******************************************************************************/
407 /*
408  Purpose:
409 
410  TIMESTAMP prints the current YMDHMS date as a time stamp.
411 
412  Example:
413 
414  31 May 2001 09:45:54 AM
415 
416  Licensing:
417 
418  This code is distributed under the GNU LGPL license.
419 
420  Modified:
421 
422  24 September 2003
423 
424  Author:
425 
426  John Burkardt
427 
428  Parameters:
429 
430  None
431 */
432 /*{
433 # define TIME_SIZE 40
434 
435  static char time_buffer[TIME_SIZE];
436  const struct tm *tm;
437  size_t len;
438  time_t now;
439 
440  now = time ( NULL );
441  tm = localtime ( &now );
442 
443  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
444 
445  fprintf ( stdout, "%s\n", time_buffer );
446 
447  return;
448 # undef TIME_SIZE
449 }*/
450 
451 /******************************************************************************/
452 
453 double triangle_unit_monomial ( int expon[2] )
454 
455 /******************************************************************************/
456 /*
457  Purpose:
458 
459  TRIANGLE_UNIT_MONOMIAL integrates a monomial over the unit triangle.
460 
461  Discussion:
462 
463  This routine integrates a monomial of the form
464 
465  product ( 1 <= dim <= 2 ) x(dim)^expon(dim)
466 
467  where the exponents are nonnegative integers. Note that
468  if the combination 0^0 is encountered, it should be treated
469  as 1.
470 
471  Integral ( over unit triangle ) x^m y^n dx dy = m! * n! / ( m + n + 2 )!
472 
473  The integration region is:
474 
475  0 <= X
476  0 <= Y
477  X + Y <= 1.
478 
479  Licensing:
480 
481  This code is distributed under the GNU LGPL license.
482 
483  Modified:
484 
485  18 April 2009
486 
487  Author:
488 
489  John Burkardt
490 
491  Parameters:
492 
493  Input, int EXPON[2], the exponents.
494 
495  Output, double TRIANGLE_UNIT_MONOMIAL, the integral of the monomial.
496 */
497 {
498  int i;
499  int k;
500  double value;
501 /*
502  The first computation ends with VALUE = 1.0;
503 */
504  value = 1.0;
505 
506 /* k = 0;
507 
508  The first loop simply computes 1 so we short circuit it!
509 
510  for ( i = 1; i <= expon[0]; i++ )
511  {
512  k = k + 1;
513  value = value * ( double ) ( i ) / ( double ) ( k );
514  }*/
515 
516  k = expon[0];
517 
518  for ( i = 1; i <= expon[1]; i++ )
519  {
520  k = k + 1;
521  value = value * ( double ) ( i ) / ( double ) ( k );
522  }
523 
524  k = k + 1;
525  value = value / ( double ) ( k );
526 
527  k = k + 1;
528  value = value / ( double ) ( k );
529 
530  return value;
531 }
532 /******************************************************************************/
533 
534 void triangle_unit_o01 ( double w[], double xy[] )
535 
536 /******************************************************************************/
537 /*
538  Purpose:
539 
540  TRIANGLE_UNIT_O01 returns a 1 point quadrature rule for the unit triangle.
541 
542  Discussion:
543 
544  This rule is precise for monomials through degree 1.
545 
546  The integration region is:
547 
548  0 <= X
549  0 <= Y
550  X + Y <= 1.
551 
552  Licensing:
553 
554  This code is distributed under the GNU LGPL license.
555 
556  Modified:
557 
558  16 April 2009
559 
560  Author:
561 
562  John Burkardt
563 
564  Reference:
565 
566  Carlos Felippa,
567  A compendium of FEM integration formulas for symbolic work,
568  Engineering Computation,
569  Volume 21, Number 8, 2004, pages 867-890.
570 
571  Parameters:
572 
573  Output, double W[1], the weights.
574 
575  Output, double XY[2*1], the abscissas.
576 */
577 {
578  int order = 1;
579 
580  double w_save[1] = {
581  1.0 };
582  double xy_save[2*1] = {
583  0.33333333333333333333, 0.33333333333333333333 };
584 
585  r8vec_copy ( order, w_save, w );
586  r8vec_copy ( 2*order, xy_save, xy );
587 
588  return;
589 }
590 /******************************************************************************/
591 
592 void triangle_unit_o03 ( double w[], double xy[] )
593 
594 /******************************************************************************/
595 /*
596  Purpose:
597 
598  TRIANGLE_UNIT_O03 returns a 3 point quadrature rule for the unit triangle.
599 
600  Discussion:
601 
602  This rule is precise for monomials through degree 2.
603 
604  The integration region is:
605 
606  0 <= X
607  0 <= Y
608  X + Y <= 1.
609 
610  Licensing:
611 
612  This code is distributed under the GNU LGPL license.
613 
614  Modified:
615 
616  16 April 2009
617 
618  Author:
619 
620  John Burkardt
621 
622  Reference:
623 
624  Carlos Felippa,
625  A compendium of FEM integration formulas for symbolic work,
626  Engineering Computation,
627  Volume 21, Number 8, 2004, pages 867-890.
628 
629  Parameters:
630 
631  Output, double W[3], the weights.
632 
633  Output, double XY[2*3], the abscissas.
634 */
635 {
636  int order = 3;
637 
638  double w_save[3] = {
639  0.33333333333333333333,
640  0.33333333333333333333,
641  0.33333333333333333333 };
642  double xy_save[2*3] = {
643  0.66666666666666666667, 0.16666666666666666667,
644  0.16666666666666666667, 0.66666666666666666667,
645  0.16666666666666666667, 0.16666666666666666667 };
646 
647  r8vec_copy ( order, w_save, w );
648  r8vec_copy ( 2*order, xy_save, xy );
649 
650  return;
651 }
652 /******************************************************************************/
653 
654 void triangle_unit_o03b ( double w[], double xy[] )
655 
656 /******************************************************************************/
657 /*
658  Purpose:
659 
660  TRIANGLE_UNIT_O03B returns a 3 point quadrature rule for the unit triangle.
661 
662  Discussion:
663 
664  This rule is precise for monomials through degree 2.
665 
666  The integration region is:
667 
668  0 <= X
669  0 <= Y
670  X + Y <= 1.
671 
672  Licensing:
673 
674  This code is distributed under the GNU LGPL license.
675 
676  Modified:
677 
678  16 April 2009
679 
680  Author:
681 
682  John Burkardt
683 
684  Reference:
685 
686  Carlos Felippa,
687  A compendium of FEM integration formulas for symbolic work,
688  Engineering Computation,
689  Volume 21, Number 8, 2004, pages 867-890.
690 
691  Parameters:
692 
693  Output, double W[3], the weights.
694 
695  Output, double XY[2*3], the abscissas.
696 */
697 {
698  int order = 3;
699 
700  double w_save[3] = {
701  0.33333333333333333333,
702  0.33333333333333333333,
703  0.33333333333333333333 };
704  double xy_save[2*3] = {
705  0.0, 0.5,
706  0.5, 0.0,
707  0.5, 0.5 };
708 
709  r8vec_copy ( order, w_save, w );
710  r8vec_copy ( 2*order, xy_save, xy );
711 
712  return;
713 }
714 /******************************************************************************/
715 
716 void triangle_unit_o06 ( double w[], double xy[] )
717 
718 /******************************************************************************/
719 /*
720  Purpose:
721 
722  TRIANGLE_UNIT_O06 returns a 6 point quadrature rule for the unit triangle.
723 
724  Discussion:
725 
726  This rule is precise for monomials through degree 4.
727 
728  The integration region is:
729 
730  0 <= X
731  0 <= Y
732  X + Y <= 1.
733 
734  Licensing:
735 
736  This code is distributed under the GNU LGPL license.
737 
738  Modified:
739 
740  16 April 2009
741 
742  Author:
743 
744  John Burkardt
745 
746  Reference:
747 
748  Carlos Felippa,
749  A compendium of FEM integration formulas for symbolic work,
750  Engineering Computation,
751  Volume 21, Number 8, 2004, pages 867-890.
752 
753  Parameters:
754 
755  Output, double W[6], the weights.
756 
757  Output, double XY[2*6], the abscissas.
758 */
759 {
760  int order = 6;
761 
762  double w_save[6] = {
763  0.22338158967801146570,
764  0.22338158967801146570,
765  0.22338158967801146570,
766  0.10995174365532186764,
767  0.10995174365532186764,
768  0.10995174365532186764 };
769  double xy_save[2*6] = {
770  0.10810301816807022736, 0.44594849091596488632,
771  0.44594849091596488632, 0.10810301816807022736,
772  0.44594849091596488632, 0.44594849091596488632,
773  0.81684757298045851308, 0.091576213509770743460,
774  0.091576213509770743460, 0.81684757298045851308,
775  0.091576213509770743460, 0.091576213509770743460 };
776 
777  r8vec_copy ( order, w_save, w );
778  r8vec_copy ( 2*order, xy_save, xy );
779 
780  return;
781 }
782 /******************************************************************************/
783 
784 void triangle_unit_o06b ( double w[], double xy[] )
785 
786 /******************************************************************************/
787 /*
788  Purpose:
789 
790  TRIANGLE_UNIT_O06B returns a 6 point quadrature rule for the unit triangle.
791 
792  Discussion:
793 
794  This rule is precise for monomials through degree 3.
795 
796  The integration region is:
797 
798  0 <= X
799  0 <= Y
800  X + Y <= 1.
801 
802  Licensing:
803 
804  This code is distributed under the GNU LGPL license.
805 
806  Modified:
807 
808  16 April 2009
809 
810  Author:
811 
812  John Burkardt
813 
814  Reference:
815 
816  Carlos Felippa,
817  A compendium of FEM integration formulas for symbolic work,
818  Engineering Computation,
819  Volume 21, Number 8, 2004, pages 867-890.
820 
821  Parameters:
822 
823  Output, double W[6], the weights.
824 
825  Output, double XY[2*6], the abscissas.
826 */
827 {
828  int order = 6;
829 
830  double w_save[6] = {
831  0.30000000000000000000,
832  0.30000000000000000000,
833  0.30000000000000000000,
834  0.033333333333333333333,
835  0.033333333333333333333,
836  0.033333333333333333333 };
837  double xy_save[2*6] = {
838  0.66666666666666666667, 0.16666666666666666667,
839  0.16666666666666666667, 0.66666666666666666667,
840  0.16666666666666666667, 0.16666666666666666667,
841  0.0, 0.5,
842  0.5, 0.0,
843  0.5, 0.5 };
844 
845  r8vec_copy ( order, w_save, w );
846  r8vec_copy ( 2*order, xy_save, xy );
847 
848  return;
849 }
850 /******************************************************************************/
851 
852 void triangle_unit_o07 ( double w[], double xy[] )
853 
854 /******************************************************************************/
855 /*
856  Purpose:
857 
858  TRIANGLE_UNIT_O07 returns a 7 point quadrature rule for the unit triangle.
859 
860  Discussion:
861 
862  This rule is precise for monomials through degree 5.
863 
864  The integration region is:
865 
866  0 <= X
867  0 <= Y
868  X + Y <= 1.
869 
870  Licensing:
871 
872  This code is distributed under the GNU LGPL license.
873 
874  Modified:
875 
876  16 April 2009
877 
878  Author:
879 
880  John Burkardt
881 
882  Reference:
883 
884  Carlos Felippa,
885  A compendium of FEM integration formulas for symbolic work,
886  Engineering Computation,
887  Volume 21, Number 8, 2004, pages 867-890.
888 
889  Parameters:
890 
891  Output, double W[7], the weights.
892 
893  Output, double XY[2*7], the abscissas.
894 */
895 {
896  int order = 7;
897 
898  double w_save[7] = {
899  0.12593918054482715260,
900  0.12593918054482715260,
901  0.12593918054482715260,
902  0.13239415278850618074,
903  0.13239415278850618074,
904  0.13239415278850618074,
905  0.22500000000000000000 };
906  double xy_save[2*7] = {
907  0.79742698535308732240, 0.10128650732345633880,
908  0.10128650732345633880, 0.79742698535308732240,
909  0.10128650732345633880, 0.10128650732345633880,
910  0.059715871789769820459, 0.47014206410511508977,
911  0.47014206410511508977, 0.059715871789769820459,
912  0.47014206410511508977, 0.47014206410511508977,
913  0.33333333333333333333, 0.33333333333333333333 };
914 
915  r8vec_copy ( order, w_save, w );
916  r8vec_copy ( 2*order, xy_save, xy );
917 
918  return;
919 }
920 /******************************************************************************/
921 
922 void triangle_unit_o12 ( double w[], double xy[] )
923 
924 /******************************************************************************/
925 /*
926  Purpose:
927 
928  TRIANGLE_UNIT_O12 returns a 12 point quadrature rule for the unit triangle.
929 
930  Discussion:
931 
932  This rule is precise for monomials through degree 6.
933 
934  The integration region is:
935 
936  0 <= X
937  0 <= Y
938  X + Y <= 1.
939 
940  Licensing:
941 
942  This code is distributed under the GNU LGPL license.
943 
944  Modified:
945 
946  19 April 2009
947 
948  Author:
949 
950  John Burkardt
951 
952  Reference:
953 
954  Carlos Felippa,
955  A compendium of FEM integration formulas for symbolic work,
956  Engineering Computation,
957  Volume 21, Number 8, 2004, pages 867-890.
958 
959  Parameters:
960 
961  Output, double W[12], the weights.
962 
963  Output, double XY[2*12], the abscissas.
964 */
965 {
966  int order = 12;
967 
968  double w_save[12] = {
969  0.050844906370206816921,
970  0.050844906370206816921,
971  0.050844906370206816921,
972  0.11678627572637936603,
973  0.11678627572637936603,
974  0.11678627572637936603,
975  0.082851075618373575194,
976  0.082851075618373575194,
977  0.082851075618373575194,
978  0.082851075618373575194,
979  0.082851075618373575194,
980  0.082851075618373575194 };
981  double xy_save[2*12] = {
982  0.87382197101699554332, 0.063089014491502228340,
983  0.063089014491502228340, 0.87382197101699554332,
984  0.063089014491502228340, 0.063089014491502228340,
985  0.50142650965817915742, 0.24928674517091042129,
986  0.24928674517091042129, 0.50142650965817915742,
987  0.24928674517091042129, 0.24928674517091042129,
988  0.053145049844816947353, 0.31035245103378440542,
989  0.31035245103378440542, 0.053145049844816947353,
990  0.053145049844816947353, 0.63650249912139864723,
991  0.31035245103378440542, 0.63650249912139864723,
992  0.63650249912139864723, 0.053145049844816947353,
993  0.63650249912139864723, 0.31035245103378440542 };
994 
995  r8vec_copy ( order, w_save, w );
996  r8vec_copy ( 2*order, xy_save, xy );
997 
998  return;
999 }
1000 /******************************************************************************/
1001 
1003 
1004 /******************************************************************************/
1005 /*
1006  Purpose:
1007 
1008  TRIANGLE_UNIT_VOLUME returns the "volume" of the unit triangle in 2D.
1009 
1010  Discussion:
1011 
1012  The "volume" of a triangle is usually called its area.
1013 
1014  The integration region is:
1015 
1016  0 <= X,
1017  0 <= Y,
1018  X + Y <= 1.
1019 
1020  Licensing:
1021 
1022  This code is distributed under the GNU LGPL license.
1023 
1024  Modified:
1025 
1026  13 March 2008
1027 
1028  Author:
1029 
1030  John Burkardt
1031 
1032  Parameters:
1033 
1034  Output, double TRIANGLE_UNIT_VOLUME, the volume.
1035 */
1036 {
1037  double volume;
1038 
1039  volume = 1.0 / 2.0;
1040 
1041  return volume;
1042 }
triangle_unit_o01
void triangle_unit_o01(double w[], double xy[])
Definition: triangle_felippa_rule.c:534
a2
constexpr double a2
Definition: hcurl_check_approx_in_2d.cpp:39
triangle_unit_o06b
void triangle_unit_o06b(double w[], double xy[])
Definition: triangle_felippa_rule.c:784
comp_next
void comp_next(int n, int k, int a[], int *more, int *h, int *t)
Definition: triangle_felippa_rule.c:11
triangle_unit_o03
void triangle_unit_o03(double w[], double xy[])
Definition: triangle_felippa_rule.c:592
r8vec_dot_product
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: triangle_felippa_rule.c:263
monomial_value
double * monomial_value(int m, int n, int e[], double x[])
Definition: triangle_felippa_rule.c:153
order
constexpr int order
Definition: dg_projection.cpp:18
triangle_unit_o06
void triangle_unit_o06(double w[], double xy[])
Definition: triangle_felippa_rule.c:716
triangle_unit_o12
void triangle_unit_o12(double w[], double xy[])
Definition: triangle_felippa_rule.c:922
triangle_unit_o07
void triangle_unit_o07(double w[], double xy[])
Definition: triangle_felippa_rule.c:852
a1
constexpr double a1
Definition: hcurl_check_approx_in_2d.cpp:38
a
constexpr double a
Definition: approx_sphere.cpp:30
double
h
double h
Definition: photon_diffusion.cpp:60
triangle_felippa_rule.h
subcomp_next
void subcomp_next(int n, int k, int a[], int *more, int *h, int *t)
Definition: triangle_felippa_rule.c:304
triangle_unit_o03b
void triangle_unit_o03b(double w[], double xy[])
Definition: triangle_felippa_rule.c:654
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
r8vec_copy
void r8vec_copy(int n, double a1[], double a2[])
Definition: triangle_felippa_rule.c:220
convert.n
n
Definition: convert.py:82
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
triangle_unit_monomial
double triangle_unit_monomial(int expon[2])
Definition: triangle_felippa_rule.c:453
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
triangle_unit_volume
double triangle_unit_volume()
Definition: triangle_felippa_rule.c:1002
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20