v0.14.0
ftensor_and_adolc.cpp
Go to the documentation of this file.
1 
2 #ifndef WITH_ADOL_C
3 #error "MoFEM need to be compiled with ADOL-C"
4 #endif
5 
6 #include <ostream>
7 
8 #include <adolc/adolc.h>
9 
10 #define FTENSOR_DEBUG
11 #include <FTensor.hpp>
12 
13 #define BOOST_UBLAS_SHALLOW_ARRAY_ADAPTOR
14 
15 #include <boost/numeric/ublas/storage.hpp>
16 #include <boost/numeric/ublas/matrix.hpp>
17 #include <boost/numeric/ublas/matrix_proxy.hpp>
18 #include <boost/numeric/ublas/io.hpp>
19 
20 int main(int argc, char *argv[]) {
23 
24  boost::numeric::ublas::vector<double> input(3);
25  input.clear();
26  FTensor::Tensor2_symmetric<double *, 2> t2(&input[0], &input[1], &input[2]);
27  t2(0, 0) = 1;
28  t2(0, 1) = 0.5;
29  t2(1, 1) = 1;
30 
32  one(0, 0) = one(1, 1) = 1;
33  one(1, 0) = one(0, 1) = 1;
34 
35  for (int ii = 0; ii != 2; ii++) {
36  for (int jj = 0; jj != 2; jj++) {
37  std::cout << t2(ii, jj) << " ";
38  }
39  std::cout << std::endl;
40  }
41  std::cout << std::endl;
42 
43  double r;
44 
45  int tag = 1;
46  int keep = 1;
47 
48  trace_on(tag, keep);
50  // FTensor::Tensor2_symmetric<adouble,2> a_t2_tmp;
51  a_t2(I, J) <<= t2(I, J);
52  // a_t2_tmp = a_t2;
53  adouble a_r = a_t2(I, J) * a_t2(I, J);
54  // a_t2(0,0)*a_t2(0,0)+a_t2(1,1)*a_t2(1,1)+a_t2(0,1)*a_t2(0,1);
55  a_r >>= r;
56  trace_off();
57 
58  for (int ii = 0; ii != 2; ii++) {
59  for (int jj = 0; jj != 2; jj++) {
60  std::cout << a_t2(ii, jj) << " ";
61  }
62  std::cout << std::endl;
63  }
64  std::cout << std::endl;
65 
66  std::cout << a_r << " " << r << std::endl << std::endl;
67 
68  boost::numeric::ublas::vector<double> output(3);
69  FTensor::Tensor2_symmetric<double *, 2> jac_t2(&output[0], &output[1],
70  &output[2]);
71 
72  /*int ret_jac =*/::gradient(tag, 3, &input[0], &output[0]);
73 
74  std::cout << input << std::endl;
75  std::cout << output << std::endl;
76 
77  for (int ii = 0; ii != 2; ii++) {
78  for (int jj = 0; jj != 2; jj++) {
79  std::cout << jac_t2(ii, jj) << " ";
80  if (jac_t2(ii, jj) != 2) {
81  std::cerr << "Wrong value should be 2\n" << std::endl;
82  exit(-1);
83  }
84  }
85  std::cout << std::endl;
86  }
87  std::cout << std::endl;
88 
89  boost::numeric::ublas::matrix<double> Hessian(3, 3);
90  Hessian.clear();
91  double *H[3];
92  for (int nn = 0; nn != 3; nn++) {
93  H[nn] = &Hessian(nn, 0);
94  }
95  /*int ret_val_hessian =*/::hessian(tag, 3, &input[0], H);
96  std::cout << Hessian << std::endl;
97 
99  &Hessian(0, 0), &Hessian(0, 1), &Hessian(0, 2), &Hessian(1, 0),
100  &Hessian(1, 1), &Hessian(1, 2), &Hessian(2, 0), &Hessian(2, 1),
101  &Hessian(2, 2));
102  for (int ii = 0; ii != 2; ii++) {
103  for (int jj = 0; jj != 2; jj++) {
104  for (int II = 0; II != 2; II++) {
105  for (int JJ = 0; JJ != 2; JJ++) {
106  std::cout << "(" << ii << "," << jj << "," << II << "," << JJ << ") "
107  << t4(ii, jj, II, JJ) << std::endl;
108  }
109  }
110  }
111  }
112  std::cout << std::endl;
113 
114  if (t4(0, 0, 0, 0) != 2 || t4(1, 1, 1, 1) != 2 || t4(1, 0, 1, 0) != 4 ||
115  t4(0, 1, 0, 1) != 4 || t4(0, 1, 1, 0) != 4 || t4(1, 0, 0, 1) != 4) {
116  std::cerr << "Wrong value\n" << std::endl;
117  exit(-1);
118  }
119 
120  return 0;
121 }
main
int main(int argc, char *argv[])
Definition: ftensor_and_adolc.cpp:20
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
sdf.r
int r
Definition: sdf.py:8
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
FTensor::Index
Definition: Index.hpp:23
adouble
FTensor::Ddg
Definition: Ddg_value.hpp:7
FTensor.hpp
Tensors class implemented by Walter Landry.
H
double H
Hardening.
Definition: plastic.cpp:124