v0.14.0
Loading...
Searching...
No Matches
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
20int main(int argc, char *argv[]) {
21 FTensor::Index<'I', 2> I;
22 FTensor::Index<'J', 2> J;
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}
static Index< 'J', 3 > J
Tensors class implemented by Walter Landry.
int main()
constexpr IntegrationType I
double H
Hardening.
Definition plastic.cpp:175