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 44 45 #ifndef ROL_PDEOPT_ELASTICITYSIMP_CONSTRAINT_H 46 #define ROL_PDEOPT_ELASTICITYSIMP_CONSTRAINT_H 47 48 #include "ROL_Constraint_SimOpt.hpp" 49 #include "ROL_TpetraMultiVector.hpp" 50 #include "Amesos2.hpp" 51 #include "filter.hpp" 52 #include "data.hpp" 53 54 template<class Real> 55 class EqualityConstraint_PDEOPT_ElasticitySIMP : public ROL::Constraint_SimOpt<Real> { 56 private: 57 58 const ROL::Ptr<ElasticitySIMPOperators<Real> > data_; 59 const ROL::Ptr<DensityFilter<Real> > filter_; 60 61 public: 62 EqualityConstraint_PDEOPT_ElasticitySIMP(const ROL::Ptr<ElasticitySIMPOperators<Real>> & data,const ROL::Ptr<DensityFilter<Real>> & filter,const Teuchos::RCP<Teuchos::ParameterList> & parlist)63 EqualityConstraint_PDEOPT_ElasticitySIMP(const ROL::Ptr<ElasticitySIMPOperators<Real> > &data, 64 const ROL::Ptr<DensityFilter<Real> > &filter, 65 const Teuchos::RCP<Teuchos::ParameterList> &parlist) 66 : data_(data), filter_(filter) {} 67 68 using ROL::Constraint_SimOpt<Real>::value; 69 value(ROL::Vector<Real> & c,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)70 void value(ROL::Vector<Real> &c, 71 const ROL::Vector<Real> &u, 72 const ROL::Vector<Real> &z, 73 Real &tol) { 74 ROL::Ptr<Tpetra::MultiVector<> > cp 75 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(c)).getVector(); 76 ROL::Ptr<const Tpetra::MultiVector<> > up 77 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(u)).getVector(); 78 ROL::Ptr<const Tpetra::MultiVector<> > zp 79 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 80 81 data_->ApplyJacobian1ToVec(cp, up); 82 Real one(1); 83 cp->update(-one, *(data_->getVecF()), one); 84 } 85 solve(ROL::Vector<Real> & c,ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)86 void solve(ROL::Vector<Real> &c, 87 ROL::Vector<Real> &u, 88 const ROL::Vector<Real> &z, 89 Real &tol) { 90 ROL::Ptr<Tpetra::MultiVector<> > up 91 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(u)).getVector(); 92 ROL::Ptr<const Tpetra::MultiVector<> > zp 93 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 94 // Solve PDE 95 data_->ApplyInverseJacobian1ToVec(up, data_->getVecF(), false); 96 // Compute residual 97 EqualityConstraint_PDEOPT_ElasticitySIMP<Real>::value(c,u,z,tol); 98 } 99 applyJacobian_1(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)100 void applyJacobian_1(ROL::Vector<Real> &jv, 101 const ROL::Vector<Real> &v, 102 const ROL::Vector<Real> &u, 103 const ROL::Vector<Real> &z, 104 Real &tol) { 105 ROL::Ptr<Tpetra::MultiVector<> > jvp 106 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(jv)).getVector(); 107 ROL::Ptr<const Tpetra::MultiVector<> > vp 108 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 109 ROL::Ptr<const Tpetra::MultiVector<> > zp 110 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 111 112 data_->ApplyJacobian1ToVec(jvp, vp); 113 } 114 115 applyJacobian_2(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)116 void applyJacobian_2(ROL::Vector<Real> &jv, 117 const ROL::Vector<Real> &v, 118 const ROL::Vector<Real> &u, 119 const ROL::Vector<Real> &z, 120 Real &tol) { 121 ROL::Ptr<Tpetra::MultiVector<> > jvp 122 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(jv)).getVector(); 123 ROL::Ptr<const Tpetra::MultiVector<> > vp 124 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 125 ROL::Ptr<const Tpetra::MultiVector<> > up 126 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(u)).getVector(); 127 ROL::Ptr<const Tpetra::MultiVector<> > zp 128 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 129 130 ROL::Ptr<Tpetra::MultiVector<> > Fvp 131 = ROL::makePtr<Tpetra::MultiVector<>>(vp->getMap(), 1); 132 filter_->apply(Fvp, vp); 133 data_->ApplyJacobian2ToVec (jvp, up, Fvp); 134 } 135 136 applyAdjointJacobian_1(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)137 void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, 138 const ROL::Vector<Real> &v, 139 const ROL::Vector<Real> &u, 140 const ROL::Vector<Real> &z, 141 Real &tol) { 142 ROL::Ptr<Tpetra::MultiVector<> > ajvp 143 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(ajv)).getVector(); 144 ROL::Ptr<const Tpetra::MultiVector<> > vp 145 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 146 ROL::Ptr<const Tpetra::MultiVector<> > zp 147 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 148 149 data_->ApplyJacobian1ToVec(ajvp, vp); 150 } 151 152 applyAdjointJacobian_2(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)153 void applyAdjointJacobian_2(ROL::Vector<Real> &ajv, 154 const ROL::Vector<Real> &v, 155 const ROL::Vector<Real> &u, 156 const ROL::Vector<Real> &z, 157 Real &tol) { 158 ROL::Ptr<Tpetra::MultiVector<> > ajvp 159 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(ajv)).getVector(); 160 ROL::Ptr<const Tpetra::MultiVector<> > vp 161 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 162 ROL::Ptr<const Tpetra::MultiVector<> > up 163 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(u)).getVector(); 164 ROL::Ptr<const Tpetra::MultiVector<> > zp 165 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 166 167 ROL::Ptr<Tpetra::MultiVector<> > tmp 168 = ROL::makePtr<Tpetra::MultiVector<>>(ajvp->getMap(), 1); 169 data_->ApplyAdjointJacobian2ToVec (tmp, up, vp); 170 filter_->apply(ajvp, tmp); 171 } 172 173 applyAdjointHessian_11(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)174 void applyAdjointHessian_11(ROL::Vector<Real> &ahwv, 175 const ROL::Vector<Real> &w, 176 const ROL::Vector<Real> &v, 177 const ROL::Vector<Real> &u, 178 const ROL::Vector<Real> &z, 179 Real &tol) { 180 ahwv.zero(); 181 } 182 183 applyAdjointHessian_12(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)184 void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, 185 const ROL::Vector<Real> &w, 186 const ROL::Vector<Real> &v, 187 const ROL::Vector<Real> &u, 188 const ROL::Vector<Real> &z, 189 Real &tol) { 190 applyAdjointJacobian_2(ahwv, w, v, z, tol); 191 } 192 193 applyAdjointHessian_21(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)194 void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, 195 const ROL::Vector<Real> &w, 196 const ROL::Vector<Real> &v, 197 const ROL::Vector<Real> &u, 198 const ROL::Vector<Real> &z, 199 Real &tol) { 200 applyJacobian_2(ahwv, v, w, z, tol); 201 } 202 203 applyAdjointHessian_22(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)204 void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, 205 const ROL::Vector<Real> &w, 206 const ROL::Vector<Real> &v, 207 const ROL::Vector<Real> &u, 208 const ROL::Vector<Real> &z, 209 Real &tol) { 210 ROL::Ptr<Tpetra::MultiVector<> > ahwvp 211 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(ahwv)).getVector(); 212 ROL::Ptr<const Tpetra::MultiVector<> > wp 213 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(w)).getVector(); 214 ROL::Ptr<const Tpetra::MultiVector<> > vp 215 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 216 ROL::Ptr<const Tpetra::MultiVector<> > up 217 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(u)).getVector(); 218 ROL::Ptr<const Tpetra::MultiVector<> > zp 219 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 220 221 ROL::Ptr<Tpetra::MultiVector<> > Fvp 222 = ROL::makePtr<Tpetra::MultiVector<>>(vp->getMap(), 1); 223 filter_->apply(Fvp, vp); 224 ROL::Ptr<Tpetra::MultiVector<> > tmp 225 = ROL::makePtr<Tpetra::MultiVector<>>(ahwvp->getMap(), 1); 226 data_->ApplyAdjointHessian22ToVec (tmp, up, Fvp, wp); 227 filter_->apply(ahwvp, tmp); 228 } 229 230 applyInverseJacobian_1(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)231 void applyInverseJacobian_1(ROL::Vector<Real> &ijv, 232 const ROL::Vector<Real> &v, 233 const ROL::Vector<Real> &u, 234 const ROL::Vector<Real> &z, 235 Real &tol) { 236 ROL::Ptr<Tpetra::MultiVector<> > ijvp 237 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(ijv)).getVector(); 238 ROL::Ptr<const Tpetra::MultiVector<> > vp 239 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 240 ROL::Ptr<const Tpetra::MultiVector<> > zp 241 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 242 243 data_->ApplyInverseJacobian1ToVec (ijvp, vp, false); 244 } 245 246 applyInverseAdjointJacobian_1(ROL::Vector<Real> & iajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)247 void applyInverseAdjointJacobian_1(ROL::Vector<Real> &iajv, 248 const ROL::Vector<Real> &v, 249 const ROL::Vector<Real> &u, 250 const ROL::Vector<Real> &z, 251 Real &tol) { 252 ROL::Ptr<Tpetra::MultiVector<> > iajvp 253 = (dynamic_cast<ROL::TpetraMultiVector<Real>&>(iajv)).getVector(); 254 ROL::Ptr<const Tpetra::MultiVector<> > vp 255 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(v)).getVector(); 256 ROL::Ptr<const Tpetra::MultiVector<> > zp 257 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 258 259 data_->ApplyInverseJacobian1ToVec (iajvp, vp, true); 260 } 261 update_2(const ROL::Vector<Real> & z,bool flag=true,int iter=-1)262 void update_2(const ROL::Vector<Real> &z, bool flag = true, int iter = -1) { 263 ROL::Ptr<const Tpetra::MultiVector<> > zp 264 = (dynamic_cast<const ROL::TpetraMultiVector<Real>&>(z)).getVector(); 265 ROL::Ptr<Tpetra::MultiVector<> > Fzp 266 = ROL::makePtr<Tpetra::MultiVector<>>(zp->getMap(), 1); 267 filter_->apply(Fzp, zp); 268 data_->updateMaterialDensity(Fzp); 269 data_->constructSolverWithNewMaterial(); 270 } 271 272 }; 273 274 #endif 275