ROL
ROL_CantileverBeam.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
51#ifndef ROL_CANTILEVERBEAM_HPP
52#define ROL_CANTILEVERBEAM_HPP
53
54#include "ROL_StdVector.hpp"
55#include "ROL_TestProblem.hpp"
56#include "ROL_Bounds.hpp"
57#include "ROL_StdObjective.hpp"
58#include "ROL_StdConstraint.hpp"
59
60namespace ROL {
61namespace ZOO {
62
63template<class Real>
65private:
66 int nseg_;
67 Real len_;
68 std::vector<Real> l_;
69
70public:
71 Objective_CantileverBeam(const int nseg = 5)
72 : nseg_(nseg), len_(500) {
73 l_.clear(); l_.resize(nseg, len_/static_cast<Real>(nseg_));
74 }
75
76 Real value( const std::vector<Real> &x, Real &tol ) {
77 Real val(0);
78 for (int i = 0; i < nseg_; ++i) {
79 val += l_[i]*x[i]*x[i+nseg_];
80 }
81 return val;
82 }
83
84 void gradient( std::vector<Real> &g, const std::vector<Real> &x, Real &tol ) {
85 for (int i = 0; i < nseg_; ++i) {
86 g[i] = l_[i]*x[i+nseg_];
87 g[i+nseg_] = l_[i]*x[i];
88 }
89 }
90
91 void hessVec( std::vector<Real> &hv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
92 for (int i = 0; i < nseg_; ++i) {
93 hv[i] = l_[i]*v[i+nseg_];
94 hv[i+nseg_] = l_[i]*v[i];
95 }
96 }
97}; // class Objective_CantileverBeam
98
99
100template<class Real>
102private:
103 int nseg_;
105 std::vector<Real> l_, suml_, M_, dyp_;
106
107public:
108 Constraint_CantileverBeam(const int nseg = 5)
109 : nseg_(nseg), len_(500), P_(50e3), E_(200e5), Sigma_max_(14e3), tip_max_(2.5), suml_(0) {
110 const Real half(0.5);
111 l_.clear();
112 l_.resize(nseg_, len_/static_cast<Real>(nseg_));
113 suml_.resize(nseg_);
114 M_.resize(nseg_);
115 dyp_.resize(nseg_);
116 for (int i = 0; i < nseg_; ++i) {
117 suml_[i] = static_cast<Real>(i+1)*l_[i];
118 M_[i] = P_*(len_ - suml_[i] + l_[i]);
119 dyp_[i] = (len_ - suml_[i] + half*l_[i])*static_cast<Real>(nseg_-i-1);
120 }
121 }
122
123 void value( std::vector<Real> &c, const std::vector<Real> &x, Real &tol ) {
124 const Real one(1), two(2), three(3), twelve(12), twenty(20);
125 std::vector<Real> y1(nseg_), yp(nseg_);
126 Real Inertia(0), sigma(0), sumy1(0), sumypl(0);
127 for (int i = 0; i < nseg_; ++i) {
128 // stress constraint
129 Inertia = x[i]*std::pow(x[i+nseg_],3)/twelve;
130 sigma = M_[i]*(x[i+nseg_]/two)/Inertia;
131 c[i] = sigma / Sigma_max_ - one;
132 // lateral slope constraint
133 c[i+nseg_] = x[i+nseg_] - twenty*x[i];
134 // tip deflection constraint
135 y1[i] = P_*std::pow(l_[i],2)/(two*E_*Inertia) * (len_ - suml_[i] + two/three*l_[i]);
136 yp[i] = P_*l_[i]/(E_*Inertia) * (len_ - suml_[i] + l_[i]/two);
137 sumy1 += y1[i];
138 }
139 // tip deflection constraint
140 for (int i = 1; i < nseg_; ++i) {
141 yp[i] += yp[i-1];
142 sumypl += yp[i-1]*l_[i];
143 }
144 Real yN = sumy1 + sumypl;
145 c[2*nseg_] = yN / tip_max_ - one;
146 }
147
148 void applyJacobian( std::vector<Real> &jv, const std::vector<Real> &v,
149 const std::vector<Real> &x, Real &tol ) {
150 const Real two(2), three(3), six(6), twelve(12), twenty(20);
151 // const Real one(1); // Unused
152 Real Inertia(0), dyN(0), sumW(0), sumH(0);
153 for (int i = 0; i < nseg_; ++i) {
154 // stress constraint
155 jv[i] = -six/Sigma_max_*M_[i]*v[i]/std::pow(x[i]*x[i+nseg_],2)
156 -twelve/Sigma_max_*M_[i]*v[i+nseg_]/(x[i]*std::pow(x[i+nseg_],3));
157 // lateral slope constraint
158 jv[i+nseg_] = v[i+nseg_] - twenty*v[i];
159 // tip deflection constraint
160 Inertia = x[i]*std::pow(x[i+nseg_],3)/twelve;
161 dyN = -P_/E_ * std::pow(l_[i],2)/Inertia*((len_ - suml_[i] + two/three*l_[i])/two + dyp_[i]);
162 sumW += v[i]*(dyN/x[i])/tip_max_;
163 sumH += three*v[i+nseg_]*(dyN/x[i+nseg_])/tip_max_;
164 }
165 // tip deflection constraint
166 jv[2*nseg_] = sumW+sumH;
167 }
168
169 void applyAdjointJacobian( std::vector<Real> &ajv, const std::vector<Real> &v,
170 const std::vector<Real> &x, Real &tol ) {
171 const Real two(2), three(3), six(6), twelve(12), twenty(20);
172// const Real one(1); // Unused
173 Real Inertia(0), dyN(0);
174 for (int i = 0; i < nseg_; ++i) {
175 // stress constraint
176 ajv[i] = -six/Sigma_max_*M_[i]*v[i]/std::pow(x[i]*x[i+nseg_],2);
177 ajv[i+nseg_] = -twelve/Sigma_max_*M_[i]*v[i]/(x[i]*std::pow(x[i+nseg_],3));
178 // lateral slope constraint
179 ajv[i] += -twenty*v[i+nseg_];
180 ajv[i+nseg_] += v[i+nseg_];
181 // tip deflection constraint
182 Inertia = x[i]*std::pow(x[i+nseg_],3)/twelve;
183 dyN = -P_/E_ * std::pow(l_[i],2)/Inertia*((len_ - suml_[i] + two/three*l_[i])/two + dyp_[i]);
184 ajv[i] += v[2*nseg_]*(dyN/x[i])/tip_max_;
185 ajv[i+nseg_] += three*v[2*nseg_]*(dyN/x[i+nseg_])/tip_max_;
186 }
187 }
188
189// void applyAdjointHessian( std::vector<Real> &ahuv, const std::vector<Real> &u,
190// const std::vector<Real> &v, const std::vector<Real> &x,
191// Real &tol ) {
192// }
193
194}; // class Constraint_CantileverBeam
195
196
197
198template<class Real>
199class getCantileverBeam : public TestProblem<Real> {
200private:
201 int nseg_;
202
203public:
204 getCantileverBeam(int nseg = 5) : nseg_(nseg) {}
205
206 Ptr<Objective<Real> > getObjective( void ) const {
207 return makePtr<Objective_CantileverBeam<Real>>(nseg_);
208 }
209
210 Ptr<Constraint<Real> > getInequalityConstraint( void ) const {
211 return makePtr<Constraint_CantileverBeam<Real>>(nseg_);
212 }
213
214 Ptr<BoundConstraint<Real>> getBoundConstraint( void ) const {
215 // Lower bound is zero
216 Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(2*nseg_);
217 for (int i = 0; i < nseg_; ++i) {
218 (*lp)[i] = static_cast<Real>(1);
219 (*lp)[i+nseg_] = static_cast<Real>(5);
220 }
221
222 // No upper bound
223 Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(2*nseg_,ROL_INF<Real>());
224
225 Ptr<Vector<Real>> l = makePtr<StdVector<Real>>(lp);
226 Ptr<Vector<Real>> u = makePtr<StdVector<Real>>(up);
227
228 return makePtr<Bounds<Real>>(l,u);
229 }
230
231 Ptr<Vector<Real>> getInitialGuess( void ) const {
232 Ptr<std::vector<Real> > x0p = makePtr<std::vector<Real>>(2*nseg_,0);
233 for (int i = 0; i < nseg_; ++i) {
234 (*x0p)[i] = static_cast<Real>(5);
235 (*x0p)[i+nseg_] = static_cast<Real>(40);
236 }
237
238 return makePtr<StdVector<Real>>(x0p);
239 }
240
241 Ptr<Vector<Real>> getSolution( const int i = 0 ) const {
242 //Ptr<std::vector<Real> > xp = makePtr<std::vector<Real>>(2);
243 //(*xp)[0] = 3.0;
244 //(*xp)[1] = std::sqrt(3.0);
245
246 Ptr<std::vector<Real>> xp = makePtr<std::vector<Real>>(2*nseg_);
247 for (int i = 0; i < nseg_; ++i) {
248 (*xp)[i] = 3.0;
249 (*xp)[i+nseg_] = std::sqrt(3.0);
250 }
251
252 return makePtr<StdVector<Real>>(xp);
253 }
254
255 Ptr<Vector<Real>> getInequalityMultiplier( void ) const {
256 Ptr<std::vector<Real> > lp = makePtr<std::vector<Real>>(2*nseg_+1,0.0);
257 return makePtr<StdVector<Real>>(lp);
258 }
259
260 Ptr<BoundConstraint<Real>> getSlackBoundConstraint(void) const {
261 // No lower bound
262 Ptr<std::vector<Real> > lp = makePtr<std::vector<Real>>(2*nseg_+1,ROL_NINF<Real>());
263
264 // Upper bound is zero
265 Ptr<std::vector<Real> > up = makePtr<std::vector<Real>>(2*nseg_+1,0);
266
267 Ptr<Vector<Real> > l = makePtr<StdVector<Real>>(lp);
268 Ptr<Vector<Real> > u = makePtr<StdVector<Real>>(up);
269
270 return makePtr<Bounds<Real>>(l,u);
271 }
272};
273
274
275} // namespace ZOO
276} // namespace ROL
277
278#endif // ROL_CANTILEVERBEAM_HPP
279
Contains definitions of test objective functions.
Defines the equality constraint operator interface for StdVectors.
Specializes the ROL::Objective interface for objective functions that operate on ROL::StdVector's.
void applyAdjointJacobian(std::vector< Real > &ajv, 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 applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void gradient(std::vector< Real > &g, const std::vector< Real > &x, Real &tol)
Real value(const std::vector< Real > &x, Real &tol)
void hessVec(std::vector< Real > &hv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Ptr< Vector< Real > > getSolution(const int i=0) const
Ptr< Objective< Real > > getObjective(void) const
Ptr< Constraint< Real > > getInequalityConstraint(void) const
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Ptr< Vector< Real > > getInequalityMultiplier(void) const
Ptr< BoundConstraint< Real > > getSlackBoundConstraint(void) const