ROL
ROL_Cantilever.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
15 #ifndef USE_HESSVEC
16 #define USE_HESSVEC 1
17 #endif
18 
19 #ifndef ROL_CANTILEVER_HPP
20 #define ROL_CANTILEVER_HPP
21 
22 #include "ROL_ScaledStdVector.hpp"
23 #include "ROL_StdObjective.hpp"
24 #include "ROL_StdConstraint.hpp"
25 #include "ROL_TestProblem.hpp"
26 
27 namespace ROL {
28 namespace ZOO {
29 
30  template<class Real>
31  class Objective_Cantilever : public StdObjective<Real> {
32  public:
34 
35  Real value( const std::vector<Real> &x, Real &tol ) {
36  return x[0]*x[1];
37  }
38 
39  void gradient( std::vector<Real> &g, const std::vector<Real> &x, Real &tol ) {
40  g[0] = x[1];
41  g[1] = x[0];
42  }
43 #if USE_HESSVEC
44  void hessVec( std::vector<Real> &hv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
45  hv[0] = v[1];
46  hv[1] = v[0];
47  }
48 #endif
49  };
50 
51  template<class Real>
52  class Constraint_Cantilever : public StdConstraint<Real> {
53  private:
54  Real stress(const Real w, const Real t, const int deriv = 0, const int comp1 = 0, const int comp2 = 0) const {
55  const Real scale(600), X(500), Y(1000);
56  Real val(0);
57  if (deriv == 0) {
58  val = scale*(Y/(w*t*t) + X/(w*w*t));
59  }
60  else if (deriv == 1) {
61  if (comp1 == 0) {
62  const Real two(2);
63  val = scale*(-Y/(w*w*t*t) - two*X/(w*w*w*t));
64  }
65  else if (comp1 == 1) {
66  const Real two(2);
67  val = scale*(-two*Y/(w*t*t*t) - X/(w*w*t*t));
68  }
69  }
70  else if (deriv == 2) {
71  if (comp1 == 0 && comp2 == 0) {
72  const Real two(2), six(6);
73  val = scale*(two*Y/(w*w*w*t*t) + six*X/(w*w*w*w*t));
74  }
75  else if (comp1 == 1 && comp2 == 1) {
76  const Real two(2), six(6);
77  val = scale*(six*Y/(w*t*t*t*t) + two*X/(w*w*t*t*t));
78  }
79  else if (comp1 == 0 && comp2 == 1) {
80  const Real two(2);
81  val = scale*two*(Y/(w*w*t*t*t) + X/(w*w*w*t*t));
82  }
83  else if (comp1 == 1 && comp2 == 0) {
84  const Real two(2);
85  val = scale*two*(Y/(w*w*t*t*t) + X/(w*w*w*t*t));
86  }
87  }
88  return val;
89  }
90 
91  Real displacement(const Real w, const Real t, const int deriv = 0, const int comp1 = 0, const int comp2 = 0) const {
92  const Real four(4), L(100), E(2.9e7), X(500), Y(1000);
93  const Real C = four*std::pow(L,3)/E;
94  Real arg1 = std::pow(Y/(t*t),2), arg2 = std::pow(X/(w*w),2);
95  Real mag = std::sqrt(arg1 + arg2);
96  Real val(0);
97  if (deriv == 0) {
98  val = C/(w*t)*mag;
99  }
100  else if (deriv == 1) {
101  if (comp1 == 0) {
102  const Real three(3);
103  val = -C * (three * std::pow(X*t*t,2) + std::pow(Y*w*w,2))
104  / (std::pow(w,6)*std::pow(t,5)*mag);
105  }
106  else if (comp1 == 1) {
107  const Real three(3);
108  val = -C * (std::pow(X*t*t,2) + three*std::pow(Y*w*w,2))
109  / (std::pow(w,5)*std::pow(t,6)*mag);
110  }
111  }
112  else if (deriv == 2) {
113  if (comp1 == 0 && comp2 == 0) {
114  const Real two(2), six(6), nine(9);
115  val = C * two * mag * (std::pow(Y*w*w,4) + nine*std::pow(Y*X*w*w*t*t,2) + six*std::pow(X*t*t,4))
116  / (std::pow(w,3)*t*std::pow(std::pow(Y*w*w,2)+std::pow(X*t*t,2),2));
117  }
118  else if (comp1 == 1 && comp2 == 1) {
119  const Real two(2), six(6), nine(9);
120  val = C * two * mag * (six*std::pow(Y*w*w,4) + nine*std::pow(Y*X*w*w*t*t,2) + std::pow(X*t*t,4))
121  / (std::pow(t,3)*w*std::pow(std::pow(Y*w*w,2)+std::pow(X*t*t,2),2));
122  }
123  else if (comp1 == 0 && comp2 == 1) {
124  const Real two(2), three(3);
125  val = C * (three*std::pow(X*t*t,4) + two*std::pow(X*Y*t*t*w*w,2) + three*std::pow(Y*w*w,4))
126  / (std::pow(t*w,6)*mag*(std::pow(X*t*t,2) + std::pow(Y*w*w,2)));
127  }
128  else if (comp1 == 1 && comp2 == 0) {
129  const Real two(2), three(3);
130  val = C * (three*std::pow(X*t*t,4) + two*std::pow(X*Y*t*t*w*w,2) + three*std::pow(Y*w*w,4))
131  / (std::pow(t*w,6)*mag*(std::pow(X*t*t,2) + std::pow(Y*w*w,2)));
132  }
133  }
134  return val;
135  }
136  public:
138 
139  void value( std::vector<Real> &c, const std::vector<Real> &x, Real &tol ) {
140  const Real R(40000), D(2.2535), one(1);
141  Real s = stress(x[0],x[1],0)/R;
142  Real d = displacement(x[0],x[1],0)/D;
143  c[0] = s - one;
144  c[1] = d - one;
145  }
146 
147  void applyJacobian( std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
148  const Real R(40000), D(2.2535);
149  Real s0 = stress(x[0],x[1],1,0)/R, s1 = stress(x[0],x[1],1,1)/R;
150  Real d0 = displacement(x[0],x[1],1,0)/D, d1 = displacement(x[0],x[1],1,1)/D;
151  jv[0] = s0*v[0] + s1*v[1];
152  jv[1] = d0*v[0] + d1*v[1];
153  }
154 
155  void applyAdjointJacobian( std::vector<Real> &ajv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
156  const Real R(40000), D(2.2535);
157  Real s0 = stress(x[0],x[1],1,0)/R, s1 = stress(x[0],x[1],1,1)/R;
158  Real d0 = displacement(x[0],x[1],1,0)/D, d1 = displacement(x[0],x[1],1,1)/D;
159  ajv[0] = s0*v[0] + d0*v[1];
160  ajv[1] = s1*v[0] + d1*v[1];
161  }
162 #if USE_HESSVEC
163  void applyAdjointHessian( std::vector<Real> &ahuv, const std::vector<Real> &u, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
164  const Real R(40000), D(2.2535);
165  Real s00 = stress(x[0],x[1],2,0,0)/R, s01 = stress(x[0],x[1],2,0,1)/R;
166  Real s10 = stress(x[0],x[1],2,1,0)/R, s11 = stress(x[0],x[1],2,1,1)/R;
167  Real d00 = displacement(x[0],x[1],2,0,0)/D, d01 = displacement(x[0],x[1],2,0,1)/D;
168  Real d10 = displacement(x[0],x[1],2,1,0)/D, d11 = displacement(x[0],x[1],2,1,1)/D;
169  ahuv[0] = (s00*u[0] + d00*u[1])*v[0] + (s01*u[0] + d01*u[1])*v[1];
170  ahuv[1] = (s10*u[0] + d10*u[1])*v[0] + (s11*u[0] + d11*u[1])*v[1];
171  }
172 #endif
173  };
174 
175  template<class Real>
176  class getCantilever : public TestProblem<Real> {
177  public:
179 
180  Ptr<Objective<Real>> getObjective(void) const {
181  return makePtr<Objective_Cantilever<Real>>();
182  }
183 
184  Ptr<Vector<Real>> getInitialGuess(void) const {
185  int n = 2;
186  Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
187  Ptr<std::vector<Real>> xp = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
188  (*xp)[0] = static_cast<Real>(2.0);
189  (*xp)[1] = static_cast<Real>(2.0);
190  return makePtr<PrimalScaledStdVector<Real>>(xp,scale);
191  }
192 
193  Ptr<Vector<Real>> getSolution(const int i = 0) const {
194  int n = 2;
195  Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
196  Ptr<std::vector<Real>> xp = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
197  (*xp)[0] = static_cast<Real>(2.3520341271);
198  (*xp)[1] = static_cast<Real>(3.3262784077);
199  return makePtr<PrimalScaledStdVector<Real>>(xp,scale);
200  }
201 
202  Ptr<BoundConstraint<Real>> getBoundConstraint(void) const {
203  int n = 2;
204  Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
205  Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
206  Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(n,static_cast<Real>(4.0));
207  Ptr<Vector<Real>> l = makePtr<PrimalScaledStdVector<Real>>(lp,scale);
208  Ptr<Vector<Real>> u = makePtr<PrimalScaledStdVector<Real>>(up,scale);
209  return makePtr<Bounds<Real>>(l,u);
210  }
211 
212  Ptr<Constraint<Real>> getInequalityConstraint(void) const {
213  return makePtr<Constraint_Cantilever<Real>>();
214  }
215 
216  Ptr<Vector<Real>> getInequalityMultiplier(void) const {
217  Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(2,static_cast<Real>(1.0));
218  Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(2,static_cast<Real>(0.0));
219  return makePtr<DualScaledStdVector<Real>>(lp,scale);
220  }
221 
222  Ptr<BoundConstraint<Real>> getSlackBoundConstraint(void) const {
223  Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(2,static_cast<Real>(1.0));
224  Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(2,static_cast<Real>(0.0));
225  Ptr<Vector<Real>> u = makePtr<DualScaledStdVector<Real>>(up,scale);
226  return makePtr<Bounds<Real>>(*u,false);
227  }
228  };
229 
230 }// End ZOO Namespace
231 }// End ROL Namespace
232 
233 #endif
Ptr< Objective< Real > > getObjective(void) const
void gradient(std::vector< Real > &g, const std::vector< Real > &x, Real &tol)
void applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Ptr< BoundConstraint< Real > > getSlackBoundConstraint(void) const
Real displacement(const Real w, const Real t, const int deriv=0, const int comp1=0, const int comp2=0) const
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Real stress(const Real w, const Real t, const int deriv=0, const int comp1=0, const int comp2=0) const
Defines the equality constraint operator interface for StdVectors.
virtual void hessVec(std::vector< Real > &hv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Specializes the ROL::Objective interface for objective functions that operate on ROL::StdVector&#39;s.
Contains definitions of test objective functions.
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
Ptr< Vector< Real > > getSolution(const int i=0) const
Ptr< Constraint< Real > > getInequalityConstraint(void) const
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Real value(const std::vector< Real > &x, Real &tol)
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol) override
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< Vector< Real > > getInequalityMultiplier(void) const