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