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_H 46 #define ROL_PDEOPT_ELASTICITYSIMP_H 47 48 #include "Intrepid_HGRAD_QUAD_C2_FEM.hpp" 49 #include "Intrepid_HGRAD_QUAD_C1_FEM.hpp" 50 #include "../TOOLS/elasticity.hpp" 51 #include "../TOOLS/materials.hpp" 52 53 template<class Real> 54 class ElasticitySIMP : public Elasticity <Real> { 55 protected: 56 Real initDensity_; 57 Real minDensity_; 58 int powerP_; 59 60 std::vector<ROL::Ptr<Material_SIMP<Real> > >SIMPmaterial_; 61 ROL::Ptr<Tpetra::MultiVector<> > myDensity_; 62 ROL::Ptr<Tpetra::MultiVector<> > myCellArea_; 63 ROL::Ptr<Intrepid::FieldContainer<Real> > CBMat0_; 64 ROL::Ptr<Intrepid::FieldContainer<Real> > gradgradMats0_; 65 66 public: ElasticitySIMP()67 ElasticitySIMP() {} ~ElasticitySIMP()68 ~ElasticitySIMP() {} 69 ElasticitySIMP_Initialize(const ROL::Ptr<const Teuchos::Comm<int>> & comm,const Teuchos::RCP<Teuchos::ParameterList> & parlist,const ROL::Ptr<std::ostream> & outStream)70 virtual void ElasticitySIMP_Initialize(const ROL::Ptr<const Teuchos::Comm<int> > &comm, 71 const Teuchos::RCP<Teuchos::ParameterList> &parlist, 72 const ROL::Ptr<std::ostream> &outStream) { 73 this->Initialize(comm, parlist, outStream); 74 // new material parameters 75 Teuchos::ParameterList &list = this->parlist_->sublist("ElasticitySIMP"); 76 powerP_ = list.get<int>("SIMP Power"); 77 initDensity_ = list.get<Real>("Initial Density"); 78 minDensity_ = list.get<Real>("Minimum Density"); 79 } 80 81 SetSIMPParallelStructure()82 virtual void SetSIMPParallelStructure() { 83 PDE_FEM<Real>::SetParallelStructure(); 84 this->myCellMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), 85 this->myCellIds_, 0, this->commPtr_); 86 } 87 getDomainMapA() const88 ROL::Ptr<const Tpetra::Map<> > getDomainMapA() const { 89 return this->matA_->getDomainMap(); 90 } 91 getCellMap() const92 ROL::Ptr<const Tpetra::Map<> > getCellMap() const { 93 return this->myCellMap_; 94 } 95 getMaterialDensity() const96 ROL::Ptr<Tpetra::MultiVector<> > getMaterialDensity() const { 97 return myDensity_; 98 } 99 100 // return cell measures getCellAreas()101 ROL::Ptr<Tpetra::MultiVector<> > getCellAreas() { 102 myCellArea_ = ROL::makePtr<Tpetra::MultiVector<>>(this->myCellMap_, 1, true); 103 for (int i=0; i<this->numCells_; i++){ 104 myCellArea_ -> replaceGlobalValue(this->myCellIds_[i], 0, this->myCellMeasure_[i]); 105 } 106 return myCellArea_; 107 } 108 // ApplyBCToVec(const ROL::Ptr<Tpetra::MultiVector<>> & u)109 void ApplyBCToVec (const ROL::Ptr<Tpetra::MultiVector<> > & u) { 110 Real zero(0.0); 111 // u is myOverlapMap_ 112 for (int i=0; i<this->myDirichletDofs_.size(); ++i) { 113 if (this->myOverlapMap_->isNodeGlobalElement(this->myDirichletDofs_[i])) 114 u->replaceGlobalValue(this->myDirichletDofs_[i], 0, zero); 115 } 116 } 117 resetMaterialDensity(const Real val)118 void resetMaterialDensity (const Real val) { 119 myDensity_ = ROL::makePtr<Tpetra::MultiVector<>>(this->myCellMap_, 1, true); 120 myDensity_->putScalar(val); 121 renewMaterialVector(); 122 } 123 updateMaterialDensity(const ROL::Ptr<const Tpetra::MultiVector<>> & newDensity)124 void updateMaterialDensity (const ROL::Ptr<const Tpetra::MultiVector<> > & newDensity) { 125 myDensity_ = ROL::makePtr<Tpetra::MultiVector<>>(this->myCellMap_, 1, true); 126 Tpetra::deep_copy(*myDensity_, *newDensity); 127 renewMaterialVector(); 128 } 129 renewMaterialVector()130 void renewMaterialVector () { 131 Teuchos::ArrayRCP<const Real> densData = myDensity_->get1dView(); 132 for(int i=0; i<this->numCells_; i++) { 133 Real dens = densData[myDensity_->getMap()->getLocalElement(this->myCellIds_[i])]; 134 SIMPmaterial_[i]->setDensity(dens); 135 } 136 } 137 CreateMaterial()138 virtual void CreateMaterial() { 139 for(int i=0; i<this->numCells_; i++) { 140 ROL::Ptr<Material_SIMP<Real> > CellMaterial = ROL::makePtr<Material_SIMP<Real>>(); 141 CellMaterial->InitializeSIMP(this->spaceDim_, this->planeStrain_, this->E_, 142 this->poissonR_, initDensity_, powerP_, minDensity_); 143 this->materialTensorDim_ = CellMaterial->GetMaterialTensorDim(); 144 SIMPmaterial_.push_back(CellMaterial); 145 } 146 resetMaterialDensity (initDensity_); 147 } 148 149 ComputeLocalSystemMats(bool ifInitial)150 virtual void ComputeLocalSystemMats(bool ifInitial) { 151 int full_lfs = this->lfs_ * this->spaceDim_; 152 if(!ifInitial) { 153 #ifdef ROL_TIMERS 154 Teuchos::TimeMonitor LocalTimer(*LocalAssemblyTime_example_PDEOPT_TOOLS_PDEFEM_GLOB); 155 #endif 156 renewMaterialVector(); 157 //this->gradgradMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, full_lfs); 158 Construct_CBmats(ifInitial); 159 Intrepid::FunctionSpaceTools::integrate<Real>(*this->gradgradMats_, // compute local grad.grad (stiffness) matrices 160 *this->CBMat_, 161 *this->BMatWeighted_, 162 Intrepid::COMP_CPP); 163 return; 164 } 165 CreateMaterial(); 166 this->gradgradMats0_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, full_lfs); 167 this->gradgradMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, full_lfs); 168 this->valvalMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, full_lfs); 169 170 this->BMat_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->materialTensorDim_); 171 this->BMatWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->materialTensorDim_); 172 CBMat0_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->materialTensorDim_); 173 this->CBMat_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->materialTensorDim_); 174 this->NMat_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->spaceDim_); 175 this->NMatWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numCells_, full_lfs, this->numCubPoints_, this->spaceDim_); 176 this->Construct_N_B_mats(); 177 Construct_CBmats(ifInitial); 178 Intrepid::FunctionSpaceTools::integrate<Real>(*this->gradgradMats0_, // compute local grad.grad (stiffness) matrices 179 *CBMat0_, 180 *this->BMatWeighted_, 181 Intrepid::COMP_CPP); 182 Intrepid::FunctionSpaceTools::integrate<Real>(*this->gradgradMats_, // compute local grad.grad (stiffness) matrices 183 *this->CBMat_, 184 *this->BMatWeighted_, 185 Intrepid::COMP_CPP); 186 Intrepid::FunctionSpaceTools::integrate<Real>(*this->valvalMats_, // compute local val.val (mass) matrices 187 *this->NMat_, 188 *this->NMatWeighted_, 189 Intrepid::COMP_CPP); 190 } 191 192 // test matrices test_mats()193 virtual void test_mats() { 194 ROL::Ptr<Intrepid::FieldContainer<Real> > test_Jaco_Mat; 195 ROL::Ptr<Intrepid::FieldContainer<Real> > test_Grad_Mat; 196 ROL::Ptr<Intrepid::FieldContainer<Real> > test_N_Mat; 197 ROL::Ptr<Intrepid::FieldContainer<Real> > test_B_Mat; 198 ROL::Ptr<Intrepid::FieldContainer<Real> > test_K0_Mat; 199 ROL::Ptr<Intrepid::FieldContainer<Real> > test_K_Mat; 200 ROL::Ptr<Intrepid::FieldContainer<Real> > test_M_Mat; 201 ROL::Ptr<Intrepid::FieldContainer<Real> > test_F_Vec; 202 203 test_Jaco_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->spaceDim_, this->spaceDim_); 204 test_Grad_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->lfs_, this->spaceDim_); 205 test_N_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, this->spaceDim_); 206 test_B_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, this->materialTensorDim_); 207 test_K0_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, this->numLocalDofs_); 208 test_K_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, this->numLocalDofs_); 209 test_M_Mat = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, this->numLocalDofs_); 210 test_F_Vec = ROL::makePtr<Intrepid::FieldContainer<Real>>(this->numLocalDofs_, 1); 211 212 for(int i=0; i<this->spaceDim_; i++) { 213 for(int j=0; j<this->spaceDim_; j++) { 214 (*test_Jaco_Mat)(i, j) = (*this->cellJac_)(0, 0, i, j); 215 } 216 } 217 for(int i=0; i<this->numLocalDofs_; i++) { 218 for(int j=0; j<this->spaceDim_; j++) { 219 if(i<this->lfs_) 220 (*test_Grad_Mat)(i, j) = (*this->gradReference_)(i, 0, j); 221 222 (*test_N_Mat)(i, j) = (*this->NMat_)(0, i, 0, j); 223 } 224 for(int j=0; j<this->materialTensorDim_; j++) { 225 (*test_B_Mat)(i, j) = (*this->BMat_)(0, i, 0, j); 226 } 227 for(int j=0; j<this->numLocalDofs_; j++) { 228 (*test_K0_Mat)(i, j) = (*this->gradgradMats0_)(0, i, j); 229 (*test_K_Mat)(i, j) = (*this->gradgradMats_)(0, i, j); 230 (*test_M_Mat)(i, j) = (*this->valvalMats_)(0, i, j); 231 } 232 (*test_F_Vec)(i, 0) = (*this->datavalVecF_)(0, i); 233 } 234 std::cout<<*(SIMPmaterial_[0]->GetMaterialTensor())<<std::endl; 235 std::cout<<*test_Jaco_Mat<<std::endl; 236 std::cout<<*test_Grad_Mat<<std::endl; 237 std::cout<<*test_N_Mat<<std::endl; 238 std::cout<<*test_B_Mat<<std::endl; 239 std::cout<<*test_M_Mat<<std::endl; 240 std::cout<<*test_F_Vec<<std::endl; 241 std::cout<<*test_K0_Mat<<std::endl; 242 std::cout<<*test_K_Mat<<std::endl; 243 } 244 245 Construct_CBmats(const bool ifInitial)246 void Construct_CBmats(const bool ifInitial) { 247 if(ifInitial) 248 CBMat0_->initialize(0.0); 249 250 this->CBMat_->initialize(0.0); 251 Real SIMPScale; 252 for (int i=0; i<this->numCells_; ++i) { 253 254 SIMPScale = SIMPmaterial_[i]->getSIMPScaleFactor(); 255 ROL::Ptr<Intrepid::FieldContainer<Real> > materialMat = SIMPmaterial_[i]->GetMaterialTensor(); 256 for (int j=0; j<this->numCubPoints_; ++j) { 257 for (int m=0; m<this->lfs_*this->spaceDim_; m++) { 258 for (int n=0; n<this->materialTensorDim_; n++) { 259 for (int k=0; k<this->materialTensorDim_; k++) { 260 if(ifInitial) 261 (*CBMat0_)(i, m, j, n) += (*this->BMat_)(i, m, j, k) * (*materialMat)(k, n); 262 (*this->CBMat_)(i, m, j, n) += SIMPScale * (*this->BMat_)(i, m, j, k) * (*materialMat)(k, n); 263 } 264 } 265 } 266 } 267 } 268 } 269 270 271 }; // class ElasticitySIMPData 272 273 #endif 274