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 /*! \file traction.hpp 45 \brief Implements traction loads for the structural topology 46 optimization problem. 47 */ 48 49 #ifndef ROL_PDEOPT_ELASTICITY_TRACTION_HPP 50 #define ROL_PDEOPT_ELASTICITY_TRACTION_HPP 51 52 #include "ROL_Ptr.hpp" 53 #include "Teuchos_ParameterList.hpp" 54 #include "Intrepid_FieldContainer.hpp" 55 56 #include "../../../TOOLS/fe.hpp" 57 58 #include <vector> 59 60 template<class Real> 61 class Traction { 62 private: 63 std::vector<int> sidesets_, sideids_; 64 std::vector<Real> loadMagnitude_; 65 std::vector<Real> loadPolar_, loadAzimuth_; 66 std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > bdryCellNodes_; 67 std::vector<std::vector<std::vector<int> > > bdryCellLocIds_; 68 int dim_, offset_; 69 70 protected: DegreesToRadians(const Real deg) const71 Real DegreesToRadians(const Real deg) const { 72 return deg * static_cast<Real>(M_PI) / static_cast<Real>(180); 73 } 74 75 public: ~Traction()76 virtual ~Traction() {} 77 Traction(void)78 Traction(void) {}; 79 Traction(Teuchos::ParameterList & parlist,std::string ex="Default")80 Traction(Teuchos::ParameterList & parlist, std::string ex = "Default") { 81 if (ex == "2D Wheel") { 82 dim_ = 2; 83 offset_ = 0; 84 sidesets_.push_back(2); 85 sideids_.push_back(0); 86 loadMagnitude_.push_back(static_cast<Real>(1)); 87 loadPolar_.push_back(static_cast<Real>(270)); 88 loadAzimuth_.push_back(static_cast<Real>(0)); 89 } 90 else if (ex == "2D Cantilever with 1 Load" || 91 ex == "2D Truss") { 92 dim_ = 2; 93 offset_ = 2; 94 sidesets_.clear(); 95 sideids_.clear(); 96 loadMagnitude_.clear(); 97 loadPolar_.clear(); 98 loadAzimuth_.clear(); 99 } 100 else if (ex == "2D Cantilever with 3 Loads") { 101 dim_ = 2; 102 offset_ = 0; 103 sidesets_.push_back(1); 104 sidesets_.push_back(3); 105 sidesets_.push_back(5); 106 sideids_.push_back(0); 107 sideids_.push_back(0); 108 sideids_.push_back(0); 109 loadMagnitude_.push_back(static_cast<Real>(1)); 110 loadMagnitude_.push_back(static_cast<Real>(1)); 111 loadMagnitude_.push_back(static_cast<Real>(1)); 112 loadPolar_.push_back(static_cast<Real>(270)); 113 loadPolar_.push_back(static_cast<Real>(270)); 114 loadPolar_.push_back(static_cast<Real>(270)); 115 loadAzimuth_.push_back(static_cast<Real>(0)); 116 loadAzimuth_.push_back(static_cast<Real>(0)); 117 loadAzimuth_.push_back(static_cast<Real>(0)); 118 } 119 else if (ex == "2D Beams") { 120 dim_ = 2; 121 offset_ = 4; 122 sidesets_.clear(); 123 sideids_.clear(); 124 loadMagnitude_.clear(); 125 loadPolar_.clear(); 126 loadAzimuth_.clear(); 127 } 128 else if (ex == "2D Carrier Plate") { 129 dim_ = 2; 130 offset_ = 0; 131 sidesets_.push_back(2); 132 sidesets_.push_back(4); 133 sideids_.push_back(2); 134 sideids_.push_back(2); 135 loadMagnitude_.push_back(static_cast<Real>(1.0)); 136 loadMagnitude_.push_back(static_cast<Real>(0.5)); 137 loadPolar_.push_back(static_cast<Real>(270)); 138 loadPolar_.push_back(static_cast<Real>(270)); 139 loadAzimuth_.push_back(static_cast<Real>(0)); 140 loadAzimuth_.push_back(static_cast<Real>(0)); 141 } 142 else if (ex == "3D Cantilever") { 143 dim_ = 3; 144 offset_ = 3; 145 sidesets_.clear(); 146 sideids_.clear(); 147 loadMagnitude_.clear(); 148 loadPolar_.clear(); 149 loadAzimuth_.clear(); 150 } 151 else { 152 if (parlist.isSublist("Traction")) { 153 // Grab problem dimension 154 dim_ = parlist.get("Problem Dimension",2); 155 156 offset_ = 0; 157 if (parlist.isSublist("Load")) { 158 Teuchos::Array<Real> loadMag 159 = Teuchos::getArrayFromStringParameter<double>(parlist.sublist("Load"), "Magnitude"); 160 offset_ = dim_*static_cast<int>(loadMag.size()); 161 } 162 163 // Grab Sidesets 164 Teuchos::Array<int> sidesets 165 = Teuchos::getArrayFromStringParameter<int>(parlist.sublist("Traction"), "Sidesets"); 166 sidesets_ = sidesets.toVector(); 167 168 // Grab Side IDs 169 Teuchos::Array<int> sideids 170 = Teuchos::getArrayFromStringParameter<int>(parlist.sublist("Traction"), "Side IDs"); 171 sideids_ = sideids.toVector(); 172 173 // Grab Magnitudes 174 Teuchos::Array<Real> magnitude 175 = Teuchos::getArrayFromStringParameter<double>(parlist.sublist("Traction"), "Magnitude"); 176 loadMagnitude_ = magnitude.toVector(); 177 178 // Grab Polar Angle 179 Teuchos::Array<Real> polar 180 = Teuchos::getArrayFromStringParameter<double>(parlist.sublist("Traction"), "Polar Angle"); 181 loadPolar_ = polar.toVector(); 182 183 if (dim_ == 3) { 184 // Grab Azimuth Angle 185 Teuchos::Array<Real> azimuth 186 = Teuchos::getArrayFromStringParameter<double>(parlist.sublist("Traction"), "Azimuth Angle"); 187 loadAzimuth_ = azimuth.toVector(); 188 } 189 } 190 } 191 192 int polarSize = loadPolar_.size(); 193 for (int i=0; i<polarSize; ++i) { 194 loadPolar_[i] = DegreesToRadians(loadPolar_[i]); 195 } 196 197 int azimuthSize = loadAzimuth_.size(); 198 for (int i=0; i<azimuthSize; ++i) { 199 loadAzimuth_[i] = DegreesToRadians(loadAzimuth_[i]); 200 } 201 } 202 isNull(void) const203 bool isNull(void) const { 204 return (sidesets_.size()==0); 205 } 206 setCellNodes(const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & bdryCellNodes,const std::vector<std::vector<std::vector<int>>> & bdryCellLocIds)207 void setCellNodes(const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > &bdryCellNodes, 208 const std::vector<std::vector<std::vector<int> > > &bdryCellLocIds) { 209 bdryCellNodes_ = bdryCellNodes; 210 bdryCellLocIds_ = bdryCellLocIds; 211 } 212 evaluate(const int dir,const int sideset,const std::vector<Real> & param) const213 virtual Real evaluate(const int dir, 214 const int sideset, 215 const std::vector<Real> ¶m) const { 216 // const int numSideSets = sidesets_.size(); 217 const int paramSize = param.size(); 218 Real val(0); 219 Real loadMagNoise(0), loadAngNoise0(0); 220 if (paramSize > 0 + dim_*sideset) { 221 loadMagNoise = param[offset_ + 0 + dim_*sideset]; 222 } 223 if (paramSize > 1 + dim_*sideset) { 224 loadAngNoise0 = DegreesToRadians(param[offset_ + 1 + dim_*sideset]); 225 } 226 const Real loadMagnitude = loadMagnitude_[sideset] + loadMagNoise; 227 const Real loadAngle0 = loadPolar_[sideset] + loadAngNoise0; 228 229 if (dim_==2) { 230 if (dir==0) { 231 val = loadMagnitude*std::cos(loadAngle0); 232 } 233 if (dir==1) { 234 val = loadMagnitude*std::sin(loadAngle0); 235 } 236 } 237 if (dim_==3) { 238 Real loadAngNoise1(0); 239 if (paramSize > 2 + dim_*sideset) { 240 loadAngNoise1 = DegreesToRadians(param[offset_ + 2 + dim_*sideset]); 241 } 242 const Real loadAngle1 = loadAzimuth_[sideset] + loadAngNoise1; 243 244 if (dir==0) { 245 val = loadMagnitude*std::sin(loadAngle0)*std::cos(loadAngle1); 246 } 247 if (dir==1) { 248 val = loadMagnitude*std::sin(loadAngle0)*std::sin(loadAngle1); 249 } 250 if (dir==2) { 251 val = loadMagnitude*std::cos(loadAngle0); 252 } 253 } 254 return val; 255 } 256 compute(std::vector<std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>>> & traction,const std::vector<std::vector<ROL::Ptr<FE<Real>>>> & fe,const std::vector<Real> & param,const Real scale=1) const257 void compute(std::vector<std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > > &traction, 258 const std::vector<std::vector<ROL::Ptr<FE<Real> > > > &fe, 259 const std::vector<Real> ¶m, 260 const Real scale = 1) const { 261 traction.clear(); 262 traction.resize(bdryCellLocIds_.size()); 263 const int numSideSets = sidesets_.size(); 264 if (numSideSets > 0) { 265 for (int i = 0; i < numSideSets; ++i) { 266 traction[sidesets_[i]].resize(bdryCellLocIds_[sidesets_[i]].size()); 267 const int numCellsSide = bdryCellLocIds_[sidesets_[i]][sideids_[i]].size(); 268 if (numCellsSide > 0) { 269 traction[sidesets_[i]][sideids_[i]].resize(dim_); 270 const int numCubPerSide = fe[0][0]->gradN()->dimension(2); 271 for (int k = 0; k < dim_; ++k) { 272 traction[sidesets_[i]][sideids_[i]][k] 273 = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCellsSide, numCubPerSide); 274 for (int c = 0; c < numCellsSide; ++c) { 275 for (int p = 0; p < numCubPerSide; ++p) { 276 (*traction[sidesets_[i]][sideids_[i]][k])(c,p) = scale*evaluate(k,i,param); 277 } 278 } 279 } 280 } 281 } 282 } 283 } 284 apply(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & R,const std::vector<std::vector<ROL::Ptr<FE<Real>>>> & fe,const std::vector<Real> & param,const Real scale=1) const285 void apply(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > &R, 286 const std::vector<std::vector<ROL::Ptr<FE<Real> > > > &fe, 287 const std::vector<Real> ¶m, 288 const Real scale = 1) const { 289 const int numSideSets = sidesets_.size(); 290 const int nf = R[0]->dimension(1); 291 if (numSideSets > 0) { 292 for (int i = 0; i < numSideSets; ++i) { 293 const int numCellsSide = bdryCellLocIds_[sidesets_[i]][sideids_[i]].size(); 294 if (numCellsSide > 0) { 295 const int numCubPerSide = fe[0][0]->gradN()->dimension(2); 296 for (int k = 0; k < dim_; ++k) { 297 Intrepid::FieldContainer<Real> traction(numCellsSide, numCubPerSide); 298 for (int c = 0; c < numCellsSide; ++c) { 299 for (int p = 0; p < numCubPerSide; ++p) { 300 traction(c,p) = scale*evaluate(k,i,param); 301 } 302 } 303 Intrepid::FieldContainer<Real> tractionResidual(numCellsSide, nf); 304 Intrepid::FunctionSpaceTools::integrate<Real>(tractionResidual, 305 traction, 306 *(fe[sidesets_[i]][sideids_[i]]->NdetJ()), 307 Intrepid::COMP_CPP, false); 308 // Add Robin residual to volume residual 309 for (int l = 0; l < numCellsSide; ++l) { 310 int cidx = bdryCellLocIds_[sidesets_[i]][sideids_[i]][l]; 311 for (int m = 0; m < nf; ++m) { 312 (*R[k])(cidx,m) += tractionResidual(l,m); 313 } 314 } 315 } 316 } 317 } 318 } 319 } 320 }; 321 322 #endif 323