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
23#include "ROL_StdObjective.hpp"
24#include "ROL_StdConstraint.hpp"
25#include "ROL_TestProblem.hpp"
26
27namespace ROL {
28namespace 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
Contains definitions of test objective functions.
Defines the equality constraint operator interface for StdVectors.
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 ,...
Specializes the ROL::Objective interface for objective functions that operate on ROL::StdVector's.
virtual void hessVec(std::vector< Real > &hv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, 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)
Real displacement(const Real w, const Real t, const int deriv=0, const int comp1=0, const int comp2=0) const
Real stress(const Real w, const Real t, const int deriv=0, const int comp1=0, const int comp2=0) const
void gradient(std::vector< Real > &g, const std::vector< Real > &x, Real &tol)
Real value(const std::vector< Real > &x, Real &tol)
Ptr< Objective< Real > > getObjective(void) const
Ptr< Vector< Real > > getSolution(const int i=0) const
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Ptr< Constraint< Real > > getInequalityConstraint(void) const
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< Vector< Real > > getInequalityMultiplier(void) const
Ptr< BoundConstraint< Real > > getSlackBoundConstraint(void) const