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 data.hpp 45 \brief Generates and manages data for the Poisson example, including 46 all mesh and discretization data, matrices, etc. 47 */ 48 49 #ifndef ROL_PDEOPT_PDE_FEM_H 50 #define ROL_PDEOPT_PDE_FEM_H 51 52 #include "Teuchos_GlobalMPISession.hpp" 53 #include "Teuchos_TimeMonitor.hpp" 54 55 #include "Tpetra_MultiVector.hpp" 56 #include "Tpetra_Vector.hpp" 57 #include "Tpetra_CrsGraph.hpp" 58 #include "Tpetra_CrsMatrix.hpp" 59 #include "Tpetra_Version.hpp" 60 #include "Tpetra_RowMatrixTransposer.hpp" 61 #include "MatrixMarket_Tpetra.hpp" 62 63 #include "Intrepid_HGRAD_QUAD_C2_FEM.hpp" 64 #include "Intrepid_HGRAD_QUAD_C1_FEM.hpp" 65 #include "Intrepid_DefaultCubatureFactory.hpp" 66 #include "Intrepid_FunctionSpaceTools.hpp" 67 #include "Intrepid_CellTools.hpp" 68 69 #include "Amesos2.hpp" 70 #include "./dofmanager.hpp" 71 72 // Global Timers. 73 #ifdef ROL_TIMERS 74 ROL::Ptr<Teuchos::Time> FactorTime_example_PDEOPT_TOOLS_PDEFEM_GLOB = 75 Teuchos::TimeMonitor::getNewCounter("ROL: Factorization Time in PDEFEM"); 76 ROL::Ptr<Teuchos::Time> LUSubstitutionTime_example_PDEOPT_TOOLS_PDEFEM_GLOB = 77 Teuchos::TimeMonitor::getNewCounter("ROL: LU Substitution Time in PDEFEM"); 78 ROL::Ptr<Teuchos::Time> SolverUpdateTime_example_PDEOPT_TOOLS_PDEFEM_GLOB = 79 Teuchos::TimeMonitor::getNewCounter("ROL: Solver Update Time in PDEFEM"); 80 ROL::Ptr<Teuchos::Time> LocalAssemblyTime_example_PDEOPT_TOOLS_PDEFEM_GLOB = 81 Teuchos::TimeMonitor::getNewCounter("ROL: Local Assembly Time in PDEFEM"); 82 ROL::Ptr<Teuchos::Time> ConstraintDerivativeTime_example_PDEOPT_TOOLS_PDEFEM_GLOB = 83 Teuchos::TimeMonitor::getNewCounter("ROL: Constraint Derivative Application Time in PDEFEM"); 84 #endif 85 86 template<class Real> 87 class PDE_FEM { 88 private: 89 using GO = typename Tpetra::Map<>::global_ordinal_type; 90 91 protected: 92 93 ROL::Ptr<MeshManager<Real> > meshMgr_; 94 ROL::Ptr<DofManager<Real> > dofMgr_; 95 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > basisPtrs_; 96 ROL::Ptr<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > > coord_iface_; 97 Intrepid::FieldContainer<GO> cellDofs_; 98 99 Teuchos::RCP<Teuchos::ParameterList> parlist_; 100 ROL::Ptr<const Teuchos::Comm<int> > commPtr_; 101 ROL::Ptr<std::ostream> outStream_; 102 103 int numLocalDofs_; 104 Intrepid::FieldContainer<GO> ctn_; 105 Intrepid::FieldContainer<GO> cte_; 106 107 int myRank_; 108 int numProcs_; 109 110 Real alpha_; 111 int basisOrder_; 112 113 ROL::Ptr<const Tpetra::Map<> > myCellMap_; 114 ROL::Ptr<const Tpetra::Map<> > myOverlapMap_; 115 ROL::Ptr<const Tpetra::Map<> > myUniqueMap_; 116 ROL::Ptr<Tpetra::CrsGraph<> > matGraph_; 117 ROL::Ptr<Tpetra::CrsMatrix<> > matA_; 118 ROL::Ptr<Tpetra::CrsMatrix<> > matA_dirichlet_; 119 ROL::Ptr<Tpetra::CrsMatrix<> > matA_dirichlet_trans_; 120 ROL::Ptr<Tpetra::CrsMatrix<> > matM_; 121 ROL::Ptr<Tpetra::CrsMatrix<> > matM_dirichlet_; 122 ROL::Ptr<Tpetra::CrsMatrix<> > matM_dirichlet_trans_; 123 ROL::Ptr<Tpetra::MultiVector<> > vecUd_; 124 ROL::Ptr<Tpetra::MultiVector<> > vecF_; 125 ROL::Ptr<Tpetra::MultiVector<> > vecF_overlap_; 126 ROL::Ptr<Tpetra::MultiVector<> > vecF_dirichlet_; 127 128 Teuchos::Array<Real> myCellMeasure_; 129 Teuchos::Array<GO> myCellIds_; 130 // Elements on Boundary 131 std::vector<Teuchos::Array<GO> > myBoundaryCellIds_; 132 // DBC 133 Teuchos::Array<GO> myDirichletDofs_; 134 // BC Sides 135 std::vector<GO> my_dbc_; 136 std::vector<GO> my_nbc_; 137 138 //Point load on Bundary! 139 std::vector<GO> myPointLoadDofs_; 140 std::vector<Real> myPointLoadVals_; 141 142 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > solverA_; 143 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > solverA_trans_; 144 145 shards::CellTopology cellType_; 146 147 int sideDim_; 148 int spaceDim_; 149 int numNodesPerCell_; 150 int numCubPoints_; 151 int lfs_; 152 153 GO totalNumCells_; 154 GO totalNumDofs_; 155 GO numCells_; 156 157 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPoints_; 158 ROL::Ptr<Intrepid::FieldContainer<Real> > cubWeights_; 159 ROL::Ptr<Intrepid::FieldContainer<Real> > cellNodes_; 160 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJac_; 161 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacInv_; 162 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacDet_; 163 ROL::Ptr<Intrepid::FieldContainer<Real> > cellWeightedMeasure_; 164 ROL::Ptr<Intrepid::FieldContainer<Real> > valReference_; 165 ROL::Ptr<Intrepid::FieldContainer<Real> > gradReference_; 166 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysical_; 167 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysical_; 168 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysicalWeighted_; 169 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysicalWeighted_; 170 ROL::Ptr<Intrepid::FieldContainer<Real> > gradgradMats_; 171 ROL::Ptr<Intrepid::FieldContainer<Real> > valvalMats_; 172 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPointsPhysical_; 173 ROL::Ptr<Intrepid::FieldContainer<Real> > dataF_; 174 ROL::Ptr<Intrepid::FieldContainer<Real> > datavalVecF_; 175 ROL::Ptr<Intrepid::FieldContainer<Real> > dofPoints_; 176 ROL::Ptr<Intrepid::FieldContainer<Real> > dofPointsPhysical_; 177 ROL::Ptr<Intrepid::FieldContainer<Real> > dataUd_; 178 179 bool verbose_; 180 181 protected: 182 public: 183 184 // constructor 185 //PDE_FEM() {} PDE_FEM()186 PDE_FEM() {} 187 188 // destructor ~PDE_FEM()189 virtual ~PDE_FEM() {} 190 Initialize(const ROL::Ptr<const Teuchos::Comm<int>> & comm,const Teuchos::RCP<Teuchos::ParameterList> & parlist,const ROL::Ptr<std::ostream> & outStream)191 virtual void Initialize(const ROL::Ptr<const Teuchos::Comm<int> > &comm, 192 const Teuchos::RCP<Teuchos::ParameterList> &parlist, 193 const ROL::Ptr<std::ostream> &outStream) { 194 commPtr_ = comm; 195 parlist_ = parlist; 196 outStream_ = outStream; 197 myRank_ = comm->getRank(); 198 numProcs_ = comm->getSize(); 199 200 verbose_ = parlist->sublist("PDE FEM").get("Verbose Output",false); 201 if ( verbose_ ) { 202 *outStream_ << "Total number of processors: " << numProcs_ << std::endl; 203 } 204 205 basisOrder_ = parlist->sublist("PDE FEM").get<int>("Order of FE Discretization"); 206 } 207 SetDiscretization(const ROL::Ptr<MeshManager<Real>> & meshMgr,const std::vector<ROL::Ptr<Intrepid::Basis<Real,Intrepid::FieldContainer<Real>>>> & basisPtrs)208 void SetDiscretization(const ROL::Ptr<MeshManager<Real> > &meshMgr, 209 const std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > &basisPtrs) { 210 meshMgr_ = meshMgr; 211 totalNumCells_ = meshMgr_->getNumCells(); 212 213 basisPtrs_ = basisPtrs; 214 215 // Retrieve some basic cell information. 216 cellType_ = (basisPtrs_[0])->getBaseCellTopology(); // get the cell type from any basis 217 spaceDim_ = cellType_.getDimension(); // retrieve spatial dimension 218 numNodesPerCell_ = cellType_.getNodeCount(); // retrieve number of nodes per cell 219 220 sideDim_ = spaceDim_ - 1; 221 222 coord_iface_ = ROL::dynamicPtrCast<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > >(basisPtrs_[0]); 223 dofMgr_ = ROL::makePtr<DofManager<Real>>(meshMgr_, basisPtrs_); 224 } 225 SetParallelStructure(void)226 virtual void SetParallelStructure(void) { 227 int cellSplit = parlist_->sublist("Geometry").get<int>("Partition type"); 228 /****************************************************/ 229 /*** Build parallel communication infrastructure. ***/ 230 /****************************************************/ 231 // Partition the cells in the mesh. We use a basic quasi-equinumerous partitioning, 232 // where the remainder, if any, is assigned to the last processor. 233 Teuchos::Array<GO> myGlobIds_; 234 Teuchos::Array<GO> cellOffsets_(numProcs_, 0); 235 GO cellsPerProc = totalNumCells_ / numProcs_; 236 numCells_ = cellsPerProc; 237 switch(cellSplit) { 238 case 0: 239 if (myRank_ == 0) { // remainder in the first 240 numCells_ += totalNumCells_ % numProcs_; 241 } 242 for (int i=1; i<numProcs_; ++i) { 243 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i==1))*(totalNumCells_ % numProcs_); 244 } 245 break; 246 case 1: 247 if (myRank_ == numProcs_-1) { // remainder in the last 248 numCells_ += totalNumCells_ % numProcs_; 249 } 250 for (int i=1; i<numProcs_; ++i) { 251 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc; 252 } 253 break; 254 case 2: 255 if (myRank_ < (totalNumCells_%numProcs_)) { // spread remainder, starting from the first 256 numCells_++; 257 } 258 for (int i=1; i<numProcs_; ++i) { 259 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i-1<(totalNumCells_%numProcs_))); 260 } 261 break; 262 } 263 264 cellDofs_ = *(dofMgr_->getCellDofs()); 265 numLocalDofs_ = cellDofs_.dimension(1); 266 if ( verbose_ ) { 267 *outStream_ << "Cell offsets across processors: " << cellOffsets_ 268 << std::endl; 269 } 270 for (int i=0; i<numCells_; ++i) { 271 myCellIds_.push_back(cellOffsets_[myRank_]+i); 272 for (int j=0; j<numLocalDofs_; ++j) { 273 myGlobIds_.push_back( cellDofs_(cellOffsets_[myRank_]+i,j) ); 274 } 275 } 276 std::sort(myGlobIds_.begin(), myGlobIds_.end()); 277 myGlobIds_.erase( std::unique(myGlobIds_.begin(), myGlobIds_.end()), myGlobIds_.end() ); 278 279 // Build maps. 280 myOverlapMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), myGlobIds_, 0, commPtr_); 281 //std::cout << std::endl << myOverlapMap_->getNodeElementList()<<std::endl; 282 /** One can also use the non-member function: 283 myOverlapMap_ = Tpetra::createNonContigMap<int,int>(myGlobIds_, commPtr_); 284 to build the overlap map. 285 **/ 286 myUniqueMap_ = Tpetra::createOneToOne(myOverlapMap_); 287 //std::cout << std::endl << myUniqueMap_->getNodeElementList() << std::endl; 288 myCellMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), 289 this->myCellIds_, 0, this->commPtr_); 290 } 291 SetUpLocalIntrepidArrays(void)292 virtual void SetUpLocalIntrepidArrays(void) { 293 /****************************************************/ 294 /*** Set up local discretization data and arrays. ***/ 295 /****************************************************/ 296 // Cubature data. 297 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 298 int cubDegree = 4; // set cubature degree, e.g., 2 299 ROL::Ptr<Intrepid::Cubature<Real> > cellCub = cubFactory.create(cellType_, cubDegree); // create default cubature 300 numCubPoints_ = cellCub->getNumPoints(); // retrieve number of cubature points 301 302 lfs_ = dofMgr_->getLocalFieldSize(0); 303 // Discretization data. 304 cubPoints_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_, spaceDim_); 305 cubWeights_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_); 306 cubPointsPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_); 307 dofPoints_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs_, spaceDim_); 308 dofPointsPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs_, spaceDim_); 309 cellNodes_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numNodesPerCell_, spaceDim_); 310 cellJac_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 311 cellJacInv_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 312 cellJacDet_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 313 cellWeightedMeasure_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 314 valReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs_, numCubPoints_); 315 gradReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs_, numCubPoints_, spaceDim_); 316 valPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs_, numCubPoints_); 317 gradPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs_, numCubPoints_, spaceDim_); 318 valPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs_, numCubPoints_); 319 gradPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs_, numCubPoints_, spaceDim_); 320 321 // Geometric definition of the cells in the mesh, based on the cell-to-node map and the domain partition. 322 Intrepid::FieldContainer<Real> &nodes = *meshMgr_->getNodes(); 323 ctn_ = *meshMgr_->getCellToNodeMap(); 324 for (int i=0; i<numCells_; ++i) { 325 for (int j=0; j<numNodesPerCell_; ++j) { 326 for (int k=0; k<spaceDim_; ++k) { 327 (*cellNodes_)(i, j, k) = nodes(ctn_(myCellIds_[i],j), k); 328 } 329 } 330 } 331 332 //Compute Intrepid Arrays 333 cellCub->getCubature(*cubPoints_, *cubWeights_); // retrieve cubature points and weights 334 (*basisPtrs_[0]).getValues(*gradReference_, *cubPoints_, Intrepid::OPERATOR_GRAD); // evaluate grad operator at cubature points 335 (*basisPtrs_[0]).getValues(*valReference_, *cubPoints_, Intrepid::OPERATOR_VALUE); // evaluate value operator at cubature points 336 337 Intrepid::CellTools<Real>::setJacobian(*cellJac_, *cubPoints_, *cellNodes_, cellType_); // compute cell Jacobians 338 Intrepid::CellTools<Real>::setJacobianInv(*cellJacInv_, *cellJac_); // compute inverses of cell Jacobians 339 Intrepid::CellTools<Real>::setJacobianDet(*cellJacDet_, *cellJac_); // compute determinants of cell Jacobians 340 341 Intrepid::FunctionSpaceTools::computeCellMeasure<Real>(*cellWeightedMeasure_, // compute weighted cell measure 342 *cellJacDet_, 343 *cubWeights_); 344 345 Intrepid::FunctionSpaceTools::HGRADtransformGRAD<Real>(*gradPhysical_, // transform reference gradients into physical space 346 *cellJacInv_, 347 *gradReference_); 348 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*gradPhysicalWeighted_, // multiply with weighted measure 349 *cellWeightedMeasure_, 350 *gradPhysical_); 351 352 Intrepid::FunctionSpaceTools::HGRADtransformVALUE<Real>(*valPhysical_, // transform reference values into physical space 353 *valReference_); 354 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*valPhysicalWeighted_, // multiply with weighted measure 355 *cellWeightedMeasure_, 356 *valPhysical_); 357 Intrepid::CellTools<Real>::mapToPhysicalFrame(*cubPointsPhysical_, // map reference cubature points to physical space 358 *cubPoints_, 359 *cellNodes_, 360 cellType_); 361 362 ComputeCellTotalMeasures(); 363 } 364 ComputeCellTotalMeasures(void)365 virtual void ComputeCellTotalMeasures(void) 366 { 367 for (int i=0; i<numCells_; ++i){ 368 Real temp = 0.0; 369 for(int j=0; j<numCubPoints_; ++j){ 370 temp += (*cellWeightedMeasure_)(i, j); 371 } 372 myCellMeasure_.push_back(temp); 373 } 374 std::cout<<"First cell total measure: "<<myCellMeasure_[0]<<std::endl; 375 } 376 377 ComputeLocalSystemMats(void)378 virtual void ComputeLocalSystemMats(void) { } 379 ComputeLocalForceVec(void)380 virtual void ComputeLocalForceVec(void) { } 381 SetUpTrueDataOnNodes(void)382 virtual void SetUpTrueDataOnNodes(void) { } 383 AssembleSystemGraph(void)384 virtual void AssembleSystemGraph(void) { 385 /****************************************/ 386 /*** Assemble global graph structure. ***/ 387 /****************************************/ 388 389 Teuchos::ArrayRCP<const GO> cellDofsArrayRCP = cellDofs_.getData(); 390 Teuchos::Array<size_t> graphEntriesPerRow(myUniqueMap_->getNodeNumElements()); 391 for (GO i=0; i<numCells_; ++i) { 392 for (int j=0; j<numLocalDofs_; ++j) { 393 GO gid = cellDofs_(myCellIds_[i],j); 394 if(myUniqueMap_->isNodeGlobalElement(gid)) 395 graphEntriesPerRow[myUniqueMap_->getLocalElement(gid)] += numLocalDofs_; 396 } 397 } 398 matGraph_ = ROL::makePtr<Tpetra::CrsGraph<>>(myUniqueMap_, graphEntriesPerRow()); 399 for (GO i=0; i<numCells_; ++i) { 400 for (int j=0; j<numLocalDofs_; ++j) { 401 matGraph_->insertGlobalIndices(cellDofs_(myCellIds_[i],j), cellDofsArrayRCP(myCellIds_[i]*numLocalDofs_, numLocalDofs_)); 402 } 403 } 404 matGraph_->fillComplete(); 405 } 406 AssembleSystemMats(void)407 virtual void AssembleSystemMats(void) { 408 /****************************************/ 409 /*** Assemble global data structures. ***/ 410 /****************************************/ 411 // Assemble matrices. 412 // Stiffness matrix A. 413 matA_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 414 matA_dirichlet_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 415 int numLocalMatEntries = numLocalDofs_ * numLocalDofs_; 416 Teuchos::ArrayRCP<const GO> cellDofsArrayRCP = cellDofs_.getData(); 417 Teuchos::ArrayRCP<const Real> gradgradArrayRCP = gradgradMats_->getData(); 418 for (int i=0; i<numCells_; ++i) { 419 for (int j=0; j<numLocalDofs_; ++j) { 420 matA_->sumIntoGlobalValues(cellDofs_(myCellIds_[i],j), 421 cellDofsArrayRCP(myCellIds_[i] * numLocalDofs_, numLocalDofs_), 422 gradgradArrayRCP(i*numLocalMatEntries+j*numLocalDofs_, numLocalDofs_)); 423 matA_dirichlet_->sumIntoGlobalValues(cellDofs_(myCellIds_[i],j), 424 cellDofsArrayRCP(myCellIds_[i] * numLocalDofs_, numLocalDofs_), 425 gradgradArrayRCP(i*numLocalMatEntries+j*numLocalDofs_, numLocalDofs_)); 426 } 427 } 428 matA_->fillComplete(); 429 matA_dirichlet_->fillComplete(); 430 // Mass matrix M. 431 /* 432 matM_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 433 matM_dirichlet_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 434 Teuchos::ArrayRCP<const Real> valvalArrayRCP = valvalMats_->getData(); 435 for (int i=0; i<numCells_; ++i) { 436 for (int j=0; j<numLocalDofs_; ++j) { 437 matM_->sumIntoGlobalValues(cellDofs_(myCellIds_[i],j), 438 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs_, numLocalDofs_), 439 valvalArrayRCP(i*numLocalMatEntries+j*numLocalDofs_, numLocalDofs_)); 440 matM_dirichlet_->sumIntoGlobalValues(cellDofs_(myCellIds_[i],j), 441 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs_, numLocalDofs_), 442 valvalArrayRCP(i*numLocalMatEntries+j*numLocalDofs_, numLocalDofs_)); 443 } 444 } 445 matM_->fillComplete(); 446 matM_dirichlet_->fillComplete(); 447 */ 448 } 449 AssembleRHSVector(void)450 virtual void AssembleRHSVector(void) { 451 // Assemble vectors. 452 // vecF_ requires assembly using vecF_overlap_ and redistribution 453 vecF_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 454 vecF_dirichlet_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 455 vecF_overlap_ = ROL::makePtr<Tpetra::MultiVector<>>(myOverlapMap_, 1, true); 456 // assembly on the overlap map 457 for (int i=0; i<numCells_; ++i) { 458 for (int j=0; j<numLocalDofs_; ++j) { 459 vecF_overlap_->sumIntoGlobalValue(cellDofs_(myCellIds_[i],j), 460 0, 461 (*datavalVecF_)[i*numLocalDofs_+j]); 462 } 463 } 464 465 //Assemble the point loads! 466 for (unsigned i=0; i<myPointLoadDofs_.size(); ++i) { 467 vecF_overlap_->sumIntoGlobalValue(myPointLoadDofs_[i], 468 0, 469 myPointLoadVals_[i]); 470 } 471 472 // change map 473 Tpetra::Export<> exporter(vecF_overlap_->getMap(), vecF_->getMap()); // redistribution 474 vecF_->doExport(*vecF_overlap_, exporter, Tpetra::ADD); // from the overlap map to the unique map 475 vecF_dirichlet_->doExport(*vecF_overlap_, exporter, Tpetra::ADD); // from the overlap map to the unique map 476 } 477 478 AssembleDataVector(void)479 virtual void AssembleDataVector(void) { 480 // vecUd_ does not require assembly 481 vecUd_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getDomainMap(), 1, true); 482 for (int i=0; i<numCells_; ++i) { 483 for (int j=0; j<numLocalDofs_; ++j) { 484 if (vecUd_->getMap()->isNodeGlobalElement(cellDofs_(myCellIds_[i],j))) { 485 vecUd_->replaceGlobalValue(cellDofs_(myCellIds_[i],j), 486 0, 487 (*dataUd_)[i*numLocalDofs_+j]); 488 } 489 } 490 } 491 } 492 493 // find the local index of a global cell find_local_index(GO globalCell)494 virtual int find_local_index(GO globalCell) 495 { 496 return myCellMap_->getLocalElement(globalCell); 497 /* 498 for(int i=0; i<numCells_; i++) 499 { 500 if(myCellIds_[i] == globalCell) 501 return i; 502 } 503 return -1; 504 */ 505 } 506 507 // check to see if a globaldof is on dirichlet boundary check_myGlobalDof_on_boundary(GO globalDof)508 virtual bool check_myGlobalDof_on_boundary(GO globalDof) { 509 if (std::binary_search(myDirichletDofs_.begin(), myDirichletDofs_.end(), globalDof)) { 510 return true; 511 } 512 return false; 513 } 514 515 //create myBoundaryCellIds_ and myDirichletDofs_ SetUpMyDBCInfo(bool ifUseCoordsToSpecifyBC,std::vector<GO> dbc_side)516 virtual void SetUpMyDBCInfo(bool ifUseCoordsToSpecifyBC, std::vector<GO> dbc_side) 517 { 518 if(ifUseCoordsToSpecifyBC){ 519 my_dbc_.resize(4); 520 my_dbc_ = {0, 1, 2, 3}; 521 }else{ 522 my_dbc_ = dbc_side; 523 } 524 //Print to check my BC info 525 if ( verbose_ ) { 526 *outStream_ << "My dbc numbers: "; 527 for(unsigned i=0; i<my_dbc_.size(); ++i) { 528 *outStream_ << my_dbc_[i]; 529 } 530 *outStream_ << std::endl; 531 } 532 // 533 ROL::Ptr<std::vector<std::vector<Intrepid::FieldContainer<GO> > > > dirichletSideSets = meshMgr_->getSideSets(); 534 std::vector<std::vector<Intrepid::FieldContainer<GO> > > &dss = *dirichletSideSets; 535 Teuchos::Array<GO> mySortedCellIds_(myCellIds_); 536 std::sort(mySortedCellIds_.begin(), mySortedCellIds_.end()); 537 mySortedCellIds_.erase( std::unique(mySortedCellIds_.begin(), mySortedCellIds_.end()), mySortedCellIds_.end() ); 538 539 myBoundaryCellIds_.resize(dss[0].size()); 540 for (int i=0; i<static_cast<int>(dss[0].size()); ++i) { 541 for (int j=0; j<dss[0][i].dimension(0); ++j) { 542 if (std::binary_search(mySortedCellIds_.begin(), mySortedCellIds_.end(), dss[0][i](j))) { 543 myBoundaryCellIds_[i].push_back(dss[0][i](j)); 544 } 545 } 546 } 547 548 cte_ = *(meshMgr_->getCellToEdgeMap()); 549 Intrepid::FieldContainer<GO> &nodeDofs = *(dofMgr_->getNodeDofs()); 550 //Intrepid::FieldContainer<int> &edgeDofs = *(dofMgr_->getEdgeDofs()); 551 std::vector<std::vector<int> > dofTags = (basisPtrs_[0])->getAllDofTags(); 552 int numDofsPerNode = 0; 553 int numDofsPerEdge = 0; 554 for (int j=0; j<(basisPtrs_[0])->getCardinality(); ++j) { 555 if (dofTags[j][0] == 0) { 556 numDofsPerNode = dofTags[j][3]; 557 } 558 if (dofTags[j][0] == 1) { 559 numDofsPerEdge = dofTags[j][3]; 560 } 561 } 562 // vector field 563 int nfields = basisPtrs_.size(); 564 numDofsPerNode = numDofsPerNode * nfields; 565 numDofsPerEdge = numDofsPerEdge * nfields; 566 567 Intrepid::FieldContainer<Real> &nodes = *meshMgr_->getNodes(); 568 int n_dbc = my_dbc_.size(); 569 for (int i=0; i<static_cast<int>(myBoundaryCellIds_.size()); ++i) { 570 bool isdbc = false; 571 for(int i_dbc = 0; i_dbc < n_dbc; i_dbc++) { 572 if(i == my_dbc_[i_dbc]) { 573 isdbc = true; 574 break; 575 } 576 } 577 if(!isdbc) 578 continue; 579 580 for (int j=0; j<myBoundaryCellIds_[i].size(); ++j) { 581 /* 582 int ifDBCCell = check_DBC_By_Coords(myBoundaryCellIds_[i][j], i); 583 if(ifDBCCell < 1) 584 continue; 585 */ 586 std::vector<Real> x(spaceDim_); 587 const CellTopologyData * ctd = cellType_.getCellTopologyData(); 588 Teuchos::ArrayView<unsigned> locNodes(const_cast<unsigned *>(ctd->subcell[spaceDim_-1][i].node), cellType_.getVertexCount(spaceDim_-1, i)); 589 for (int l=0; l<static_cast<int>(cellType_.getVertexCount(spaceDim_-1, i)); ++l) { 590 x[0]=nodes(ctn_(myBoundaryCellIds_[i][j], locNodes[l]), 0); 591 x[1]=nodes(ctn_(myBoundaryCellIds_[i][j], locNodes[l]), 1); 592 // use coordinates to check if a DOF is DBC DOF 593 int ifDBCNode = check_DBC_Coords_Range( x ); 594 if(ifDBCNode < 0){ 595 continue; 596 } 597 else if(ifDBCNode == 0){ 598 if ( verbose_ ) { 599 *outStream_ << "DBC node: " 600 << ctn_(myBoundaryCellIds_[i][j], locNodes[l]) 601 << ", fixing X direction" << std::endl; 602 } 603 myDirichletDofs_.push_back(nodeDofs(ctn_(myBoundaryCellIds_[i][j], locNodes[l]), 0)); 604 } 605 else if(ifDBCNode == 1 && numDofsPerNode >= 2){ 606 if ( verbose_ ) { 607 *outStream_ << "DBC node: " 608 << ctn_(myBoundaryCellIds_[i][j], locNodes[l]) 609 << ", fixing Y direction" << std::endl; 610 } 611 myDirichletDofs_.push_back(nodeDofs(ctn_(myBoundaryCellIds_[i][j], locNodes[l]), 1)); 612 } 613 else { 614 if ( verbose_ ) { 615 *outStream_ << "DBC node: " 616 << ctn_(myBoundaryCellIds_[i][j], locNodes[l]) 617 << ", fixing ALL direction" << std::endl; 618 } 619 for (int k=0; k<numDofsPerNode; ++k) { 620 myDirichletDofs_.push_back(nodeDofs(ctn_(myBoundaryCellIds_[i][j], locNodes[l]), k)); 621 } 622 } 623 } 624 // edge dofs are NOT in use 625 /* 626 for (int k=0; k<numDofsPerEdge; ++k) { 627 myDirichletDofs_.push_back(edgeDofs(cte_(myBoundaryCellIds_[i][j], i), k)); 628 } 629 */ 630 } 631 } 632 std::sort(myDirichletDofs_.begin(), myDirichletDofs_.end()); 633 myDirichletDofs_.erase( std::unique(myDirichletDofs_.begin(), myDirichletDofs_.end()), myDirichletDofs_.end() ); 634 } 635 check_DBC_Coords_Range(const std::vector<Real> & x) const636 virtual int check_DBC_Coords_Range( const std::vector<Real> &x ) const { 637 // return value : 638 // -1: not a DBC node 639 // 0: x direction fixed 640 // 1: y direction fixed 641 // 5: both x, y direction are fixed 642 return 5; 643 } 644 // 645 // 646 // process_loading_information(const Teuchos::RCP<Teuchos::ParameterList> & parlist)647 virtual void process_loading_information(const Teuchos::RCP<Teuchos::ParameterList> &parlist) {} 648 // 649 //note that the point load is only allowed to applied on the boundary, not inside the body! 2D only check_and_Apply_PointLoad_By_Coords(int globalCellNum,int pointload_bc)650 virtual void check_and_Apply_PointLoad_By_Coords(int globalCellNum, int pointload_bc) { 651 Intrepid::FieldContainer<Real> &nodes = *meshMgr_->getNodes(); 652 const CellTopologyData * ctd = cellType_.getCellTopologyData(); 653 Teuchos::ArrayView<unsigned> locNodes(const_cast<unsigned *>(ctd->subcell[spaceDim_-1][pointload_bc].node), cellType_.getVertexCount(spaceDim_-1, pointload_bc)); 654 std::vector<Real> x1(spaceDim_); 655 std::vector<Real> x2(spaceDim_); 656 std::vector<int > localNodeNum(2); 657 x1[0]=nodes(ctn_(globalCellNum, locNodes[0]), 0); 658 x1[1]=nodes(ctn_(globalCellNum, locNodes[0]), 1); 659 x2[0]=nodes(ctn_(globalCellNum, locNodes[1]), 0); 660 x2[1]=nodes(ctn_(globalCellNum, locNodes[1]), 1); 661 ApplyPointLoad(pointload_bc, globalCellNum, localNodeNum, x1, x2); 662 } 663 ApplyPointLoad(const int pointload_bc,const GO globalCellNum,const std::vector<int> & localNodeNum,const std::vector<Real> & coord1,const std::vector<Real> & coord2)664 virtual void ApplyPointLoad(const int pointload_bc, 665 const GO globalCellNum, 666 const std::vector<int> &localNodeNum, 667 const std::vector<Real> &coord1, 668 const std::vector<Real> &coord2) { 669 //Intrepid::FieldContainer<int> &nodeDofs = *(dofMgr_->getNodeDofs()); 670 //bool isLoadPosContainedInCurrentSegment = false; 671 //int whichNodeIsCloserToPos = -1; 672 return; 673 } 674 // 675 // 676 // EnforceDBC(void)677 virtual void EnforceDBC(void) { 678 // Apply Dirichlet conditions. 679 // zero out row and column, make matrix symmetric 680 //ROL::Ptr<Tpetra::Details::DefaultTypes::node_type> node = matA_->getNode(); 681 //matA_dirichlet_ = matA_->clone(node); 682 //matM_dirichlet_ = matM_->clone(node); 683 //vecF_dirichlet_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 684 //Tpetra::deep_copy(*vecF_dirichlet_, *vecF_); 685 686 matA_dirichlet_->resumeFill(); 687 //matM_dirichlet_->resumeFill(); 688 689 GO gDof; 690 for(int i=0; i<numCells_; i++) { 691 for(int j=0; j<numLocalDofs_; j++) { 692 gDof = cellDofs_(myCellIds_[i], j); 693 if (myUniqueMap_->isNodeGlobalElement(gDof) && check_myGlobalDof_on_boundary(gDof)) { 694 size_t numRowEntries = matA_dirichlet_->getNumEntriesInGlobalRow(gDof); 695 Teuchos::Array<GO> indices(numRowEntries, 0); 696 Teuchos::Array<Real> values(numRowEntries, 0); 697 Teuchos::Array<Real> canonicalValues(numRowEntries, 0); 698 Teuchos::Array<Real> zeroValues(numRowEntries, 0); 699 matA_dirichlet_->getGlobalRowCopy(gDof, indices, values, numRowEntries); 700 //matM_dirichlet_->getGlobalRowCopy(gDof, indices, values, numRowEntries); 701 for (int k=0; k<indices.size(); k++) { 702 if (indices[k] == gDof) { 703 canonicalValues[k] = 1.0; 704 } 705 } 706 matA_dirichlet_->replaceGlobalValues(gDof, indices, canonicalValues); 707 //matM_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 708 vecF_dirichlet_->replaceGlobalValue (gDof, 0, 0); 709 } 710 711 if (!check_myGlobalDof_on_boundary(gDof)) { 712 size_t numDBCDofs = myDirichletDofs_.size(); 713 Teuchos::Array<GO> indices(numDBCDofs, 0); 714 Teuchos::Array<Real> zeroValues(numDBCDofs, 0); 715 for (int k=0; k<indices.size(); k++) { 716 indices[k] = myDirichletDofs_[k]; 717 } 718 matA_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 719 //matM_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 720 } 721 } 722 } 723 matA_dirichlet_->fillComplete(); 724 //matM_dirichlet_->fillComplete(); 725 726 // GenerateTransposedMats(); 727 // ConstructSolvers(); 728 // ConstructAdjointSolvers(); 729 } 730 731 // virtual void MatrixRemoveDBC(void) { 732 // // Apply Dirichlet conditions. 733 // // zero out row and column, make matrix symmetric 734 // //ROL::Ptr<Tpetra::Details::DefaultTypes::node_type> node = matA_->getNode(); 735 // //matA_dirichlet_ = matA_->clone(node); 736 // //matM_dirichlet_ = matM_->clone(node); 737 // 738 // matA_dirichlet_->resumeFill(); 739 // matM_dirichlet_->resumeFill(); 740 // 741 // GO gDof; 742 // for(int i=0; i<numCells_; i++) { 743 // for(int j=0; j<numLocalDofs_; j++) { 744 // gDof = cellDofs_(myCellIds_[i], j); 745 // if (myUniqueMap_->isNodeGlobalElement(gDof) && check_myGlobalDof_on_boundary(gDof)) { 746 // size_t numRowEntries = matA_dirichlet_->getNumEntriesInGlobalRow(gDof); 747 // Teuchos::Array<GO> indices(numRowEntries, 0); 748 // Teuchos::Array<Real> values(numRowEntries, 0); 749 // Teuchos::Array<Real> canonicalValues(numRowEntries, 0); 750 // Teuchos::Array<Real> zeroValues(numRowEntries, 0); 751 // matA_dirichlet_->getGlobalRowCopy(gDof, indices, values, numRowEntries); 752 // matM_dirichlet_->getGlobalRowCopy(gDof, indices, values, numRowEntries); 753 // for (int k=0; k<indices.size(); k++) { 754 // if (indices[k] == gDof) { 755 // canonicalValues[k] = 1.0; 756 // } 757 // } 758 // matA_dirichlet_->replaceGlobalValues(gDof, indices, canonicalValues); 759 // matM_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 760 // } 761 // 762 // if (!check_myGlobalDof_on_boundary(gDof)) { 763 // size_t numDBCDofs = myDirichletDofs_.size(); 764 // Teuchos::Array<GO> indices(numDBCDofs, 0); 765 // Teuchos::Array<Real> zeroValues(numDBCDofs, 0); 766 // for (int k=0; k<indices.size(); k++) { 767 // indices[k] = myDirichletDofs_[k]; 768 // } 769 // matA_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 770 // matM_dirichlet_->replaceGlobalValues(gDof, indices, zeroValues); 771 // } 772 // } 773 // } 774 // matA_dirichlet_->fillComplete(); 775 // matM_dirichlet_->fillComplete(); 776 // } 777 MatrixRemoveDBC(void)778 virtual void MatrixRemoveDBC(void) { 779 // Apply Dirichlet conditions. 780 // zero out row and column, make matrix symmetric 781 matA_dirichlet_->resumeFill(); 782 //matM_dirichlet_->resumeFill(); 783 784 for (GO i=0; i<myDirichletDofs_.size(); ++i) { 785 if (myUniqueMap_->isNodeGlobalElement(myDirichletDofs_[i])) { 786 size_t numRowEntries = matA_dirichlet_->getNumEntriesInGlobalRow(myDirichletDofs_[i]); 787 Teuchos::Array<GO> indices(numRowEntries, 0); 788 Teuchos::Array<Real> values(numRowEntries, 0); 789 Teuchos::Array<Real> canonicalValues(numRowEntries, 0); 790 Teuchos::Array<Real> zeroValues(numRowEntries, 0); 791 matA_dirichlet_->getGlobalRowCopy(myDirichletDofs_[i], indices, values, numRowEntries); 792 //matM_dirichlet_->getGlobalRowCopy(myDirichletDofs_[i], indices, values, numRowEntries); 793 for (int j=0; j < static_cast<int>(indices.size()); ++j) { 794 if (myDirichletDofs_[i] == indices[j]) { 795 canonicalValues[j] = 1.0; 796 } 797 } 798 matA_dirichlet_->replaceGlobalValues(myDirichletDofs_[i], indices, canonicalValues); 799 //matM_dirichlet_->replaceGlobalValues(myDirichletDofs_[i], indices, zeroValues); 800 for (int j=0; j < static_cast<int>(indices.size()); ++j) { 801 Teuchos::Array<GO> ind(1, myDirichletDofs_[i]); 802 Teuchos::Array<Real> valA(1, canonicalValues[j]); 803 Teuchos::Array<Real> valM(1, zeroValues[j]); 804 matA_dirichlet_->replaceGlobalValues(indices[j], ind, valA); 805 //matM_dirichlet_->replaceGlobalValues(indices[j], ind, valM); 806 } 807 } 808 } 809 matA_dirichlet_->fillComplete(); 810 //matM_dirichlet_->fillComplete(); 811 812 // GenerateTransposedMats(); 813 // ConstructSolvers(); 814 // ConstructAdjointSolvers(); 815 } 816 VectorRemoveDBC(void)817 virtual void VectorRemoveDBC(void) { 818 // Apply Dirichlet conditions. 819 //vecF_dirichlet_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 820 //Tpetra::deep_copy(*vecF_dirichlet_, *vecF_); 821 822 GO gDof(0); 823 for(GO i=0; i<numCells_; i++) { 824 for(int j=0; j<numLocalDofs_; j++) { 825 gDof = cellDofs_(myCellIds_[i], j); 826 if (myUniqueMap_->isNodeGlobalElement(gDof) && check_myGlobalDof_on_boundary(gDof)) { 827 vecF_dirichlet_->replaceGlobalValue (gDof, 0, 0); 828 } 829 } 830 } 831 } 832 GenerateTransposedMats()833 virtual void GenerateTransposedMats() { 834 // Create matrix transposes. 835 Tpetra::RowMatrixTransposer<> transposerA(matA_dirichlet_); 836 //Tpetra::RowMatrixTransposer<> transposerM(matM_dirichlet_); 837 matA_dirichlet_trans_ = transposerA.createTranspose(); 838 //matM_dirichlet_trans_ = transposerM.createTranspose(); 839 } 840 841 ConstructSolvers(void)842 virtual void ConstructSolvers(void) { 843 // Construct solver using Amesos2 factory. 844 #ifdef ROL_TIMERS 845 Teuchos::TimeMonitor LocalTimer(*FactorTime_example_PDEOPT_TOOLS_PDEFEM_GLOB); 846 #endif 847 try { 848 solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_); 849 //solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("PardisoMKL", matA_dirichlet_); 850 //solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("SuperLU", matA_dirichlet_); 851 } 852 catch (std::invalid_argument& e) { 853 std::cout << e.what() << std::endl; 854 } 855 // Perform factorization. 856 solverA_->numericFactorization(); 857 } 858 ConstructAdjointSolvers()859 virtual void ConstructAdjointSolvers() { 860 // Construct solver using Amesos2 factory. 861 try{ 862 solverA_trans_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_trans_); 863 //solverA_trans_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("PardisoMKL", matA_dirichlet_trans_); 864 //solverA_trans_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("SuperLU", matA_dirichlet_trans_); 865 } catch (std::invalid_argument& e) { 866 std::cout << e.what() << std::endl; 867 } 868 solverA_trans_->numericFactorization(); 869 } 870 getMatA(const bool & transpose=false) const871 ROL::Ptr<Tpetra::CrsMatrix<> > getMatA(const bool &transpose = false) const { 872 if (transpose) { 873 //return matA_dirichlet_trans_; 874 return matA_dirichlet_; 875 } 876 else { 877 return matA_dirichlet_; 878 } 879 } 880 881 getMatB(const bool & transpose=false) const882 ROL::Ptr<Tpetra::CrsMatrix<> > getMatB(const bool &transpose = false) const { 883 if (transpose) { 884 //return matM_dirichlet_trans_; 885 return matM_dirichlet_; 886 } 887 else { 888 return matM_dirichlet_; 889 } 890 } 891 892 getMatM() const893 ROL::Ptr<Tpetra::CrsMatrix<> > getMatM() const { 894 return matM_; 895 } 896 897 getMatR() const898 ROL::Ptr<Tpetra::CrsMatrix<> > getMatR() const { 899 return matM_; 900 } 901 902 getVecUd() const903 ROL::Ptr<Tpetra::MultiVector<> > getVecUd() const { 904 return vecUd_; 905 } 906 907 getVecF() const908 ROL::Ptr<Tpetra::MultiVector<> > getVecF() const { 909 return vecF_dirichlet_; 910 } 911 912 getSolver(const bool & transpose=false) const913 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > getSolver(const bool &transpose = false) const { 914 if (transpose) { 915 //return solverA_trans_; 916 return solverA_; 917 } 918 else { 919 return solverA_; 920 } 921 } 922 funcTarget(const Real & x1,const Real & x2) const923 virtual Real funcTarget(const Real &x1, const Real &x2) const { return 0.0; } 924 925 printMeshData(std::ostream & outStream) const926 void printMeshData(std::ostream &outStream) const { 927 928 ROL::Ptr<Intrepid::FieldContainer<Real> > nodesPtr = meshMgr_->getNodes(); 929 ROL::Ptr<Intrepid::FieldContainer<GO> > cellToNodeMapPtr = meshMgr_->getCellToNodeMap(); 930 Intrepid::FieldContainer<Real> &nodes = *nodesPtr; 931 Intrepid::FieldContainer<GO> &cellToNodeMap = *cellToNodeMapPtr; 932 outStream << "Number of nodes = " << meshMgr_->getNumNodes() << std::endl; 933 outStream << "Number of cells = " << meshMgr_->getNumCells() << std::endl; 934 outStream << "Number of edges = " << meshMgr_->getNumEdges() << std::endl; 935 // Print mesh to file. 936 if ((myRank_ == 0)) { 937 std::ofstream meshfile; 938 meshfile.open("cell_to_node_quad.txt"); 939 for (GO i=0; i<cellToNodeMap.dimension(0); ++i) { 940 for (GO j=0; j<cellToNodeMap.dimension(1); ++j) { 941 meshfile << cellToNodeMap(i,j) << " "; 942 } 943 meshfile << std::endl; 944 } 945 meshfile.close(); 946 947 meshfile.open("cell_to_node_tri.txt"); 948 for (GO i=0; i<cellToNodeMap.dimension(0); ++i) { 949 for (GO j=0; j<3; ++j) { 950 meshfile << cellToNodeMap(i,j) << " "; 951 } 952 meshfile << std::endl; 953 for (GO j=2; j<5; ++j) { 954 meshfile << cellToNodeMap(i,j%4) << " "; 955 } 956 meshfile << std::endl; 957 } 958 meshfile.close(); 959 960 meshfile.open("nodes.txt"); 961 meshfile.precision(16); 962 for (GO i=0; i<nodes.dimension(0); ++i) { 963 for (GO j=0; j<nodes.dimension(1); ++j) { 964 meshfile << std::scientific << nodes(i,j) << " "; 965 } 966 meshfile << std::endl; 967 } 968 meshfile.close(); 969 /* This somewhat clunky output is for gnuplot. 970 meshfile.open("mesh.txt"); 971 for (int i=0; i<cellToNodeMap.dimension(0); ++i) { 972 meshfile << nodes(cellToNodeMap(i,0), 0) << " " << nodes(cellToNodeMap(i,0), 1) << std::endl; 973 meshfile << nodes(cellToNodeMap(i,1), 0) << " " << nodes(cellToNodeMap(i,1), 1) << std::endl; 974 meshfile << nodes(cellToNodeMap(i,2), 0) << " " << nodes(cellToNodeMap(i,2), 1) << std::endl; 975 meshfile << nodes(cellToNodeMap(i,3), 0) << " " << nodes(cellToNodeMap(i,3), 1) << std::endl; 976 meshfile << nodes(cellToNodeMap(i,0), 0) << " " << nodes(cellToNodeMap(i,0), 1) << std::endl; 977 meshfile << nodes(cellToNodeMap(i,1), 0) << " " << nodes(cellToNodeMap(i,1), 1) << std::endl; 978 meshfile << nodes(cellToNodeMap(i,2), 0) << " " << nodes(cellToNodeMap(i,2), 1) << std::endl; 979 } 980 meshfile.close(); 981 */ 982 } //myRank 0 print 983 } // prinf function end 984 985 outputTpetraData() const986 void outputTpetraData() const { 987 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<> > matWriter; 988 matWriter.writeSparseFile("stiffness_mat", matA_); 989 matWriter.writeSparseFile("dirichlet_mat", matA_dirichlet_); 990 //matWriter.writeSparseFile("mass_mat", matM_); 991 matWriter.writeDenseFile("Ud_vec", vecUd_); 992 } 993 994 outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<>> & vec,const std::string & filename) const995 void outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<> > &vec, 996 const std::string &filename) const { 997 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<> > vecWriter; 998 vecWriter.writeDenseFile(filename, vec); 999 } 1000 1001 // ACCESSOR FUNCTIONS GetMeshManager(void)1002 ROL::Ptr<MeshManager<Real> >& GetMeshManager(void) { 1003 TEUCHOS_TEST_FOR_EXCEPTION(meshMgr_==ROL::nullPtr, std::logic_error, 1004 ">>> (PDE_FEM::GetMeshManager): MeshManager not initialized!"); 1005 return meshMgr_; 1006 } 1007 GetBasisPtr(const int ind)1008 ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > >& GetBasisPtr(const int ind) { 1009 TEUCHOS_TEST_FOR_EXCEPTION(ind > spaceDim_-1 || ind < 0, std::logic_error, 1010 ">>> (PDE_FEM::GetBasisPtr): ind out of bounds!"); 1011 TEUCHOS_TEST_FOR_EXCEPTION(basisPtrs_.size()==0, std::logic_error, 1012 ">>> (PDE_FEM::GetBasisPtr): BasisPntrs not initialized!"); 1013 return basisPtrs_[ind]; 1014 } 1015 GetBasisOrder(void) const1016 int GetBasisOrder(void) const { 1017 return basisOrder_; 1018 } 1019 GetSpaceDim(void) const1020 int GetSpaceDim(void) const { 1021 return spaceDim_; 1022 } 1023 GetNumCells(void) const1024 GO GetNumCells(void) const { 1025 return numCells_; 1026 } 1027 GetNumLocalDofs(void) const1028 int GetNumLocalDofs(void) const { 1029 return numLocalDofs_; 1030 } 1031 GetNumCubPoints(void) const1032 int GetNumCubPoints(void) const { 1033 return numCubPoints_; 1034 } 1035 GetLocalFieldSize(void) const1036 int GetLocalFieldSize(void) const { 1037 return lfs_; 1038 } 1039 1040 }; // class PDE_FEM 1041 1042 #endif 1043