v0.9.1
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
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,
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
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
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
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
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,
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
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
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
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
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
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
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
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
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
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
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 }
void triangle_unit_o03b(double w[], double xy[])
void triangle_unit_o06(double w[], double xy[])
void triangle_unit_o06b(double w[], double xy[])
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
double * monomial_value(int m, int n, int e[], double x[])
void subcomp_next(int n, int k, int a[], int *more, int *h, int *t)
constexpr double a2
double r8vec_dot_product(int n, double a1[], double a2[])
void comp_next(int n, int k, int a[], int *more, int *h, int *t)
void triangle_unit_o07(double w[], double xy[])
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
double triangle_unit_monomial(int expon[2])
double triangle_unit_volume()
void triangle_unit_o12(double w[], double xy[])
constexpr int order
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
constexpr double a1
void triangle_unit_o01(double w[], double xy[])
void triangle_unit_o03(double w[], double xy[])
void r8vec_copy(int n, double a1[], double a2[])
double w(const double g, const double t)
Definition: ContactOps.hpp:160
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26