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 Stefan-Boltzmann example, including 46 all mesh and discretization data, matrices, etc. 47 */ 48 49 #ifndef ROL_PDEOPT_STEFANBOLTZMANN_DATA_H 50 #define ROL_PDEOPT_STEFANBOLTZMANN_DATA_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 71 #include "../TOOLS/dofmanager.hpp" 72 73 template<class Real> 74 class StefanBoltzmannData { 75 76 private: 77 using GO = typename Tpetra::Map<>::global_ordinal_type; 78 79 ROL::Ptr<MeshManager<Real> > meshMgr_; 80 ROL::Ptr<DofManager<Real> > dofMgr_; 81 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > basisPtrs_; 82 83 ROL::Ptr<const Teuchos::Comm<int> > commPtr_; 84 int myRank_; 85 int numProcs_; 86 87 Real alpha_; 88 int basisOrder_; 89 90 ROL::Ptr<const Tpetra::Map<> > myOverlapMap_; 91 ROL::Ptr<const Tpetra::Map<> > myUniqueMap_; 92 ROL::Ptr<Tpetra::CrsGraph<> > matGraph_; 93 ROL::Ptr<Tpetra::CrsMatrix<> > matA_; 94 ROL::Ptr<Tpetra::CrsMatrix<> > matA_dirichlet_; 95 ROL::Ptr<Tpetra::CrsMatrix<> > matA_dirichlet_trans_; 96 ROL::Ptr<Tpetra::CrsMatrix<> > matM_; 97 ROL::Ptr<Tpetra::CrsMatrix<> > matM_dirichlet_; 98 ROL::Ptr<Tpetra::CrsMatrix<> > matM_dirichlet_trans_; 99 ROL::Ptr<Tpetra::MultiVector<> > vecUd_; 100 ROL::Ptr<Tpetra::MultiVector<> > vecF_; 101 ROL::Ptr<Tpetra::MultiVector<> > vecF_overlap_; 102 ROL::Ptr<Tpetra::MultiVector<> > vecF_dirichlet_; 103 104 Teuchos::Array<GO> myCellIds_; 105 Teuchos::Array<GO> myDirichletDofs_; 106 107 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > solverA_; 108 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > solverA_trans_; 109 110 shards::CellTopology cellType_; 111 int spaceDim_; 112 int numNodesPerCell_; 113 int numCubPoints_; 114 115 GO totalNumCells_; 116 GO totalNumDofs_; 117 GO numCells_; 118 119 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPoints_; 120 ROL::Ptr<Intrepid::FieldContainer<Real> > cubWeights_; 121 ROL::Ptr<Intrepid::FieldContainer<Real> > cellNodes_; 122 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJac_; 123 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacInv_; 124 ROL::Ptr<Intrepid::FieldContainer<Real> > cellJacDet_; 125 ROL::Ptr<Intrepid::FieldContainer<Real> > cellWeightedMeasure_; 126 ROL::Ptr<Intrepid::FieldContainer<Real> > valReference_; 127 ROL::Ptr<Intrepid::FieldContainer<Real> > gradReference_; 128 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysical_; 129 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysical_; 130 ROL::Ptr<Intrepid::FieldContainer<Real> > kappaGradPhysical_; 131 ROL::Ptr<Intrepid::FieldContainer<Real> > valPhysicalWeighted_; 132 ROL::Ptr<Intrepid::FieldContainer<Real> > gradPhysicalWeighted_; 133 ROL::Ptr<Intrepid::FieldContainer<Real> > gradgradMats_; 134 ROL::Ptr<Intrepid::FieldContainer<Real> > valvalMats_; 135 ROL::Ptr<Intrepid::FieldContainer<Real> > cubPointsPhysical_; 136 ROL::Ptr<Intrepid::FieldContainer<Real> > kappa_; 137 ROL::Ptr<Intrepid::FieldContainer<Real> > dataF_; 138 ROL::Ptr<Intrepid::FieldContainer<Real> > datavalVecF_; 139 ROL::Ptr<Intrepid::FieldContainer<Real> > dofPoints_; 140 ROL::Ptr<Intrepid::FieldContainer<Real> > dofPointsPhysical_; 141 ROL::Ptr<Intrepid::FieldContainer<Real> > dataUd_; 142 143 int sdim_; 144 145 public: 146 StefanBoltzmannData(const ROL::Ptr<const Teuchos::Comm<int>> & comm,const Teuchos::RCP<Teuchos::ParameterList> & parlist,const ROL::Ptr<std::ostream> & outStream)147 StefanBoltzmannData(const ROL::Ptr<const Teuchos::Comm<int> > &comm, 148 const Teuchos::RCP<Teuchos::ParameterList> &parlist, 149 const ROL::Ptr<std::ostream> &outStream) { 150 sdim_ = parlist->sublist("Problem").get("Stochastic Dimension",6); 151 std::vector<Real> par(sdim_,1.0); 152 153 /************************************/ 154 /*** Retrieve communication data. ***/ 155 /************************************/ 156 commPtr_ = comm; 157 myRank_ = commPtr_->getRank(); 158 numProcs_ = commPtr_->getSize(); 159 *outStream << "Total number of processors: " << numProcs_ << std::endl; 160 /************************************/ 161 /************************************/ 162 163 /*************************************/ 164 /*** Retrieve parameter list data. ***/ 165 /*************************************/ 166 alpha_ = parlist->sublist("Problem").get("Penalty", 1e-2); 167 basisOrder_ = parlist->sublist("Problem").get("Order of FE discretization", 1); 168 int cellSplit = parlist->sublist("Geometry").get("Partition type", 1); 169 /*************************************/ 170 /*************************************/ 171 172 /****************************************************************************/ 173 /*** Initialize mesh / finite element fields / degree-of-freedom manager. ***/ 174 /****************************************************************************/ 175 176 // Mesh manager. 177 meshMgr_ = ROL::makePtr<MeshManager_Rectangle<Real>>(*parlist); 178 printMeshData(*outStream); 179 // Finite element fields. 180 ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > basisPtr; 181 if (basisOrder_ == 1) { 182 basisPtr = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C1_FEM<Real, Intrepid::FieldContainer<Real> >>(); 183 } 184 else if (basisOrder_ == 2) { 185 basisPtr = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C2_FEM<Real, Intrepid::FieldContainer<Real> >>(); 186 } 187 basisPtrs_.resize(1, ROL::nullPtr); 188 basisPtrs_[0] = basisPtr; 189 // DOF coordinate interface. 190 ROL::Ptr<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > > coord_iface = 191 ROL::dynamicPtrCast<Intrepid::DofCoordsInterface<Intrepid::FieldContainer<Real> > >(basisPtrs_[0]); 192 // Degree-of-freedom manager. 193 dofMgr_ = ROL::makePtr<DofManager<Real>>(meshMgr_, basisPtrs_); 194 // Retrieve total number of cells in the mesh. 195 totalNumCells_ = meshMgr_->getNumCells(); 196 // Retrieve total number of degrees of freedom in the mesh. 197 totalNumDofs_ = dofMgr_->getNumDofs(); 198 199 /****************************************************************************/ 200 /****************************************************************************/ 201 202 203 /****************************************************/ 204 /*** Build parallel communication infrastructure. ***/ 205 /****************************************************/ 206 207 // Partition the cells in the mesh. We use a basic quasi-equinumerous partitioning, 208 // where the remainder, if any, is assigned to the last processor. 209 Teuchos::Array<GO> myGlobIds_; 210 Teuchos::Array<GO> cellOffsets_(numProcs_, 0); 211 GO cellsPerProc = totalNumCells_ / numProcs_; 212 numCells_ = cellsPerProc; 213 switch(cellSplit) { 214 case 0: 215 if (myRank_ == 0) { // remainder in the first 216 numCells_ += totalNumCells_ % numProcs_; 217 } 218 for (int i=1; i<numProcs_; ++i) { 219 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i==1))*(totalNumCells_ % numProcs_); 220 } 221 break; 222 case 1: 223 if (myRank_ == numProcs_-1) { // remainder in the last 224 numCells_ += totalNumCells_ % numProcs_; 225 } 226 for (int i=1; i<numProcs_; ++i) { 227 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc; 228 } 229 break; 230 case 2: 231 if (myRank_ < (totalNumCells_%numProcs_)) { // spread remainder, starting from the first 232 numCells_++; 233 } 234 for (int i=1; i<numProcs_; ++i) { 235 cellOffsets_[i] = cellOffsets_[i-1] + cellsPerProc + (static_cast<int>(i-1<(totalNumCells_%numProcs_))); 236 } 237 break; 238 } 239 Intrepid::FieldContainer<GO> &cellDofs = *(dofMgr_->getCellDofs()); 240 int numLocalDofs = cellDofs.dimension(1); 241 *outStream << "Cell offsets across processors: " << cellOffsets_ << std::endl; 242 for (GO i=0; i<numCells_; ++i) { 243 myCellIds_.push_back(cellOffsets_[myRank_]+i); 244 for (int j=0; j<numLocalDofs; ++j) { 245 myGlobIds_.push_back( cellDofs(cellOffsets_[myRank_]+i,j) ); 246 } 247 } 248 std::sort(myGlobIds_.begin(), myGlobIds_.end()); 249 myGlobIds_.erase( std::unique(myGlobIds_.begin(), myGlobIds_.end()), myGlobIds_.end() ); 250 251 // Build maps. 252 myOverlapMap_ = ROL::makePtr<Tpetra::Map<>>(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(), 253 myGlobIds_, 0, comm); 254 //std::cout << std::endl << myOverlapMap_->getNodeElementList(); 255 /** One can also use the non-member function: 256 myOverlapMap_ = Tpetra::createNonContigMap<int,int>(myGlobIds_, comm); 257 to build the overlap map. 258 **/ 259 myUniqueMap_ = Tpetra::createOneToOne(myOverlapMap_); 260 //std::cout << std::endl << myUniqueMap_->getNodeElementList() << std::endl; 261 262 /****************************************************/ 263 /****************************************************/ 264 265 266 /****************************************************/ 267 /*** Set up local discretization data and arrays. ***/ 268 /****************************************************/ 269 270 // Retrieve some basic cell information. 271 cellType_ = (basisPtrs_[0])->getBaseCellTopology(); // get the cell type from any basis 272 spaceDim_ = cellType_.getDimension(); // retrieve spatial dimension 273 numNodesPerCell_ = cellType_.getNodeCount(); // retrieve number of nodes per cell 274 275 // Cubature data. 276 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 277 int cubDegree = 4; // set cubature degree, e.g., 2 278 ROL::Ptr<Intrepid::Cubature<Real> > cellCub = cubFactory.create(cellType_, cubDegree); // create default cubature 279 numCubPoints_ = cellCub->getNumPoints(); // retrieve number of cubature points 280 281 int lfs = dofMgr_->getLocalFieldSize(0); 282 283 // Discretization data. 284 cubPoints_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_, spaceDim_); 285 cubWeights_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCubPoints_); 286 cubPointsPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_); 287 dofPoints_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs, spaceDim_); 288 dofPointsPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, spaceDim_); 289 cellNodes_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numNodesPerCell_, spaceDim_); 290 cellJac_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 291 cellJacInv_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_, spaceDim_, spaceDim_); 292 cellJacDet_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 293 cellWeightedMeasure_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 294 valReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs, numCubPoints_); 295 gradReference_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(lfs, numCubPoints_, spaceDim_); 296 valPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_); 297 gradPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 298 kappaGradPhysical_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 299 valPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_); 300 gradPhysicalWeighted_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, numCubPoints_, spaceDim_); 301 gradgradMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, lfs); 302 valvalMats_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs, lfs); 303 kappa_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 304 dataF_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, numCubPoints_); 305 datavalVecF_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs); 306 dataUd_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(numCells_, lfs); 307 308 // Geometric definition of the cells in the mesh, based on the cell-to-node map and the domain partition. 309 Intrepid::FieldContainer<Real> &nodes = *meshMgr_->getNodes(); 310 Intrepid::FieldContainer<GO> &ctn = *meshMgr_->getCellToNodeMap(); 311 for (GO i=0; i<numCells_; ++i) { 312 for (int j=0; j<numNodesPerCell_; ++j) { 313 for (int k=0; k<spaceDim_; ++k) { 314 (*cellNodes_)(i, j, k) = nodes(ctn(myCellIds_[i],j), k); 315 } 316 } 317 } 318 319 /****************************************************/ 320 /****************************************************/ 321 322 323 /****************************************************************/ 324 /*** Assemble cellwise contributions to vectors and matrices. ***/ 325 /****************************************************************/ 326 327 cellCub->getCubature(*cubPoints_, *cubWeights_); // retrieve cubature points and weights 328 (*basisPtrs_[0]).getValues(*gradReference_, *cubPoints_, Intrepid::OPERATOR_GRAD); // evaluate grad operator at cubature points 329 (*basisPtrs_[0]).getValues(*valReference_, *cubPoints_, Intrepid::OPERATOR_VALUE); // evaluate value operator at cubature points 330 331 Intrepid::CellTools<Real>::setJacobian(*cellJac_, *cubPoints_, *cellNodes_, cellType_); // compute cell Jacobians 332 Intrepid::CellTools<Real>::setJacobianInv(*cellJacInv_, *cellJac_); // compute inverses of cell Jacobians 333 Intrepid::CellTools<Real>::setJacobianDet(*cellJacDet_, *cellJac_); // compute determinants of cell Jacobians 334 335 Intrepid::FunctionSpaceTools::computeCellMeasure<Real>(*cellWeightedMeasure_, // compute weighted cell measure 336 *cellJacDet_, 337 *cubWeights_); 338 339 Intrepid::CellTools<Real>::mapToPhysicalFrame(*cubPointsPhysical_, // map reference cubature points to physical space 340 *cubPoints_, 341 *cellNodes_, 342 cellType_); 343 344 Intrepid::FunctionSpaceTools::HGRADtransformGRAD<Real>(*gradPhysical_, // transform reference gradients into physical space 345 *cellJacInv_, 346 *gradReference_); 347 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*gradPhysicalWeighted_, // multiply with weighted measure 348 *cellWeightedMeasure_, 349 *gradPhysical_); 350 for (GO i=0; i<numCells_; ++i) { // evaluate conductivity kappa at cubature points 351 for (int j=0; j<numCubPoints_; ++j) { 352 (*kappa_)(i, j) = funcKappa((*cubPointsPhysical_)(i, j, 0), 353 (*cubPointsPhysical_)(i, j, 1), 354 par); 355 } 356 } 357 Intrepid::FunctionSpaceTools::tensorMultiplyDataField<Real>(*kappaGradPhysical_, // multiply with conductivity kappa 358 *kappa_, 359 *gradPhysical_); 360 Intrepid::FunctionSpaceTools::integrate<Real>(*gradgradMats_, // compute local grad.(kappa)grad (stiffness) matrices 361 *kappaGradPhysical_, 362 *gradPhysicalWeighted_, 363 Intrepid::COMP_CPP); 364 365 Intrepid::FunctionSpaceTools::HGRADtransformVALUE<Real>(*valPhysical_, // transform reference values into physical space 366 *valReference_); 367 Intrepid::FunctionSpaceTools::multiplyMeasure<Real>(*valPhysicalWeighted_, // multiply with weighted measure 368 *cellWeightedMeasure_, 369 *valPhysical_); 370 Intrepid::FunctionSpaceTools::integrate<Real>(*valvalMats_, // compute local val.val (mass) matrices 371 *valPhysical_, 372 *valPhysicalWeighted_, 373 Intrepid::COMP_CPP); 374 375 for (GO i=0; i<numCells_; ++i) { // evaluate RHS function at cubature points 376 for (int j=0; j<numCubPoints_; ++j) { 377 (*dataF_)(i, j) = funcRHS((*cubPointsPhysical_)(i, j, 0), 378 (*cubPointsPhysical_)(i, j, 1), 379 par); 380 } 381 } 382 Intrepid::FunctionSpaceTools::integrate<Real>(*datavalVecF_, // compute local data.val vectors for RHS F 383 *dataF_, 384 *valPhysicalWeighted_, 385 Intrepid::COMP_CPP); 386 387 coord_iface->getDofCoords(*dofPoints_); // get coordinates of DOFs in reference cell 388 Intrepid::CellTools<Real>::mapToPhysicalFrame(*dofPointsPhysical_, // map reference DOF locations to physical space 389 *dofPoints_, 390 *cellNodes_, 391 cellType_); 392 for (GO i=0; i<numCells_; ++i) { // evaluate target function at these points 393 for (int j=0; j<lfs; ++j) { 394 (*dataUd_)(i, j) = funcTarget((*dofPointsPhysical_)(i, j, 0), (*dofPointsPhysical_)(i, j, 1)); 395 } 396 } 397 398 399 /****************************************************************/ 400 /****************************************************************/ 401 402 /****************************************/ 403 /*** Assemble global data structures. ***/ 404 /****************************************/ 405 406 // Assemble graph. 407 Teuchos::ArrayRCP<const GO> cellDofsArrayRCP = cellDofs.getData(); 408 Teuchos::Array<size_t> graphEntriesPerRow(myUniqueMap_->getNodeNumElements()); 409 for (GO i=0; i<numCells_; ++i) { 410 for (int j=0; j<numLocalDofs; ++j) { 411 GO gid = cellDofs(myCellIds_[i],j); 412 if(myUniqueMap_->isNodeGlobalElement(gid)) 413 graphEntriesPerRow[myUniqueMap_->getLocalElement(gid)] += numLocalDofs; 414 } 415 } 416 matGraph_ = ROL::makePtr<Tpetra::CrsGraph<>>(myUniqueMap_, graphEntriesPerRow()); 417 for (GO i=0; i<numCells_; ++i) { 418 for (int j=0; j<numLocalDofs; ++j) { 419 matGraph_->insertGlobalIndices(cellDofs(myCellIds_[i],j), cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs)); 420 } 421 } 422 matGraph_->fillComplete(); 423 424 // Assemble matrices. 425 // Stiffness matrix A. 426 matA_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 427 int numLocalMatEntries = numLocalDofs*numLocalDofs; 428 Teuchos::ArrayRCP<const Real> gradgradArrayRCP = gradgradMats_->getData(); 429 for (GO i=0; i<numCells_; ++i) { 430 for (int j=0; j<numLocalDofs; ++j) { 431 matA_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 432 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs), 433 gradgradArrayRCP(i*numLocalMatEntries+j*numLocalDofs, numLocalDofs)); 434 } 435 } 436 matA_->fillComplete(); 437 438 // Mass matrix M. 439 matM_ = ROL::makePtr<Tpetra::CrsMatrix<>>(matGraph_); 440 Teuchos::ArrayRCP<const Real> valvalArrayRCP = valvalMats_->getData(); 441 for (GO i=0; i<numCells_; ++i) { 442 for (int j=0; j<numLocalDofs; ++j) { 443 matM_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 444 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs), 445 valvalArrayRCP(i*numLocalMatEntries+j*numLocalDofs, numLocalDofs)); 446 } 447 } 448 matM_->fillComplete(); 449 450 // Assemble vectors. 451 // vecF_ requires assembly using vecF_overlap_ and redistribution 452 vecF_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 453 vecF_overlap_ = ROL::makePtr<Tpetra::MultiVector<>>(myOverlapMap_, 1, true); 454 for (GO i=0; i<numCells_; ++i) { // assembly on the overlap map 455 for (int j=0; j<numLocalDofs; ++j) { 456 vecF_overlap_->sumIntoGlobalValue(cellDofs(myCellIds_[i],j), 457 0, 458 (*datavalVecF_)[i*numLocalDofs+j]); 459 } 460 } 461 Tpetra::Export<> exporter(vecF_overlap_->getMap(), vecF_->getMap()); // redistribution: 462 vecF_->doExport(*vecF_overlap_, exporter, Tpetra::ADD); // from the overlap map to the unique map 463 // vecUd_ does not require assembly 464 vecUd_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getDomainMap(), 1, true); 465 for (GO i=0; i<numCells_; ++i) { 466 for (int j=0; j<numLocalDofs; ++j) { 467 if (vecUd_->getMap()->isNodeGlobalElement(cellDofs(myCellIds_[i],j))) { 468 vecUd_->replaceGlobalValue(cellDofs(myCellIds_[i],j), 469 0, 470 (*dataUd_)[i*numLocalDofs+j]); 471 } 472 } 473 } 474 475 // Apply Dirichlet conditions. 476 // Stiffness matrix with Dirichlet conditions: 477 // AD = [ A11 A12 ] where A = [ A11 A12 ] 478 // [ 0 I ] [ A21 A22 ] 479 // Mass matrix with Dirichlet conditions: 480 // MD = [ M11 M12 ] where M = [ M11 M12 ] 481 // [ 0 0 ] [ M21 M22 ] 482 // Vector F with Dirichlet conditions G: 483 // FD = [ F1 ] where F = [ F1 ] 484 // [ G ] [ F2 ] 485 matA_dirichlet_ = ROL::makePtr<Tpetra::CrsMatrix<>>(*matA_, Teuchos::DataAccess::Copy); 486 matM_dirichlet_ = ROL::makePtr<Tpetra::CrsMatrix<>>(*matA_, Teuchos::DataAccess::Copy); 487 vecF_dirichlet_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 488 Tpetra::deep_copy(*vecF_dirichlet_, *vecF_); 489 ROL::Ptr<std::vector<std::vector<Intrepid::FieldContainer<GO> > > > dirichletSideSets = meshMgr_->getSideSets(); 490 std::vector<std::vector<Intrepid::FieldContainer<GO> > > &dss = *dirichletSideSets; 491 Teuchos::Array<GO> mySortedCellIds_(myCellIds_); 492 std::sort(mySortedCellIds_.begin(), mySortedCellIds_.end()); 493 mySortedCellIds_.erase( std::unique(mySortedCellIds_.begin(), mySortedCellIds_.end()), mySortedCellIds_.end() ); 494 std::vector<Teuchos::Array<GO> > myDirichletCellIds_(dss[0].size()); 495 for (int i=0; i<static_cast<int>(dss[0].size()); ++i) { 496 for (int j=0; j<dss[0][i].dimension(0); ++j) { 497 if (std::binary_search(mySortedCellIds_.begin(), mySortedCellIds_.end(), dss[0][i](j))) { 498 myDirichletCellIds_[i].push_back(dss[0][i](j)); 499 } 500 } 501 } 502 Intrepid::FieldContainer<GO> &cte = *(meshMgr_->getCellToEdgeMap()); 503 Intrepid::FieldContainer<GO> &nodeDofs = *(dofMgr_->getNodeDofs()); 504 Intrepid::FieldContainer<GO> &edgeDofs = *(dofMgr_->getEdgeDofs()); 505 std::vector<std::vector<int> > dofTags = (basisPtrs_[0])->getAllDofTags(); 506 int numDofsPerNode = 0; 507 int numDofsPerEdge = 0; 508 for (int j=0; j<(basisPtrs_[0])->getCardinality(); ++j) { 509 if (dofTags[j][0] == 0) { 510 numDofsPerNode = dofTags[j][3]; 511 } 512 if (dofTags[j][0] == 1) { 513 numDofsPerEdge = dofTags[j][3]; 514 } 515 } 516 for (int i=0; i<static_cast<int>(myDirichletCellIds_.size()); ++i) { 517 for (int j=0; j<myDirichletCellIds_[i].size(); ++j) { 518 for (int k=0; k<numDofsPerNode; ++k) { 519 const CellTopologyData * ctd = cellType_.getCellTopologyData(); 520 Teuchos::ArrayView<unsigned> locNodes(const_cast<unsigned *>(ctd->subcell[spaceDim_-1][i].node), cellType_.getVertexCount(spaceDim_-1, i)); 521 for (int l=0; l<static_cast<int>(cellType_.getVertexCount(spaceDim_-1, i)); ++l) { 522 myDirichletDofs_.push_back(nodeDofs(ctn(myDirichletCellIds_[i][j], locNodes[l]), k)); 523 } 524 } 525 for (int k=0; k<numDofsPerEdge; ++k) { 526 myDirichletDofs_.push_back(edgeDofs(cte(myDirichletCellIds_[i][j], i), k)); 527 } 528 } 529 } 530 std::sort(myDirichletDofs_.begin(), myDirichletDofs_.end()); 531 myDirichletDofs_.erase( std::unique(myDirichletDofs_.begin(), myDirichletDofs_.end()), myDirichletDofs_.end() ); 532 matA_dirichlet_->resumeFill(); 533 matM_dirichlet_->resumeFill(); 534 for (int i=0; i<myDirichletDofs_.size(); ++i) { 535 if (myUniqueMap_->isNodeGlobalElement(myDirichletDofs_[i])) { 536 size_t numRowEntries = matA_dirichlet_->getNumEntriesInGlobalRow(myDirichletDofs_[i]); 537 Teuchos::Array<GO> indices(numRowEntries, 0); 538 Teuchos::Array<Real> values(numRowEntries, 0); 539 Teuchos::Array<Real> canonicalValues(numRowEntries, 0); 540 Teuchos::Array<Real> zeroValues(numRowEntries, 0); 541 matA_dirichlet_->getGlobalRowCopy(myDirichletDofs_[i], indices, values, numRowEntries); 542 matM_dirichlet_->getGlobalRowCopy(myDirichletDofs_[i], indices, values, numRowEntries); 543 for (int j=0; j<indices.size(); ++j) { 544 if (myDirichletDofs_[i] == indices[j]) { 545 canonicalValues[j] = 1.0; 546 } 547 } 548 matA_dirichlet_->replaceGlobalValues(myDirichletDofs_[i], indices, canonicalValues); 549 matM_dirichlet_->replaceGlobalValues(myDirichletDofs_[i], indices, zeroValues); 550 vecF_dirichlet_->replaceGlobalValue(myDirichletDofs_[i], 0, 0); 551 } 552 } 553 matA_dirichlet_->fillComplete(); 554 matM_dirichlet_->fillComplete(); 555 556 // Create matrix transposes. 557 Tpetra::RowMatrixTransposer<> transposerA(matA_dirichlet_); 558 Tpetra::RowMatrixTransposer<> transposerM(matM_dirichlet_); 559 matA_dirichlet_trans_ = transposerA.createTranspose(); 560 matM_dirichlet_trans_ = transposerM.createTranspose(); 561 562 updateA(par); 563 564 /*********************************/ 565 /*** Construct solver objects. ***/ 566 /*********************************/ 567 568 // Construct solver using Amesos2 factory. 569 try{ 570 solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_); 571 } catch (std::invalid_argument& e) { 572 std::cout << e.what() << std::endl; 573 } 574 try{ 575 solverA_trans_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_trans_); 576 } catch (std::invalid_argument& e) { 577 std::cout << e.what() << std::endl; 578 } 579 solverA_->numericFactorization(); 580 solverA_trans_->numericFactorization(); 581 582 583 /****************************************/ 584 /****************************************/ 585 586 //outputTpetraData(); 587 588 } 589 590 getMatA(const bool & transpose=false) const591 ROL::Ptr<Tpetra::CrsMatrix<> > getMatA(const bool &transpose = false) const { 592 if (transpose) { 593 return matA_dirichlet_trans_; 594 } 595 else { 596 return matA_dirichlet_; 597 } 598 } 599 600 getMatB(const bool & transpose=false) const601 ROL::Ptr<Tpetra::CrsMatrix<> > getMatB(const bool &transpose = false) const { 602 if (transpose) { 603 return matM_dirichlet_trans_; 604 } 605 else { 606 return matM_dirichlet_; 607 } 608 } 609 610 getMatM() const611 ROL::Ptr<Tpetra::CrsMatrix<> > getMatM() const { 612 return matM_; 613 } 614 615 getMatR() const616 ROL::Ptr<Tpetra::CrsMatrix<> > getMatR() const { 617 return matM_; 618 } 619 620 getVecUd() const621 ROL::Ptr<Tpetra::MultiVector<> > getVecUd() const { 622 return vecUd_; 623 } 624 625 getVecF() const626 ROL::Ptr<Tpetra::MultiVector<> > getVecF() const { 627 return vecF_dirichlet_; 628 } 629 630 getSolver(const bool & transpose=false) const631 ROL::Ptr<Amesos2::Solver< Tpetra::CrsMatrix<>, Tpetra::MultiVector<> > > getSolver(const bool &transpose = false) const { 632 if (transpose) { 633 return solverA_trans_; 634 } 635 else { 636 return solverA_; 637 } 638 } 639 640 funcRHS(const Real & x1,const Real & x2,const std::vector<Real> & par) const641 Real funcRHS(const Real &x1, const Real &x2, const std::vector<Real> &par) const { 642 Real xi1 = 0.5*(par[sdim_-2]+1.0), xi2 = 0.5*(par[sdim_-1]+1.0); 643 return std::exp(-1.e-2*(std::pow(x1-xi1,2.0)+std::pow(x2-xi2,2.0))); 644 // return 2.0*M_PI*M_PI*std::sin(M_PI*x1)*std::sin(M_PI*x2) + (1.0/(alpha_*128.0*M_PI*M_PI))*std::sin(8.0*M_PI*x1)*std::sin(8.0*M_PI*x2); 645 } 646 647 funcTarget(const Real & x1,const Real & x2) const648 Real funcTarget(const Real &x1, const Real &x2) const { 649 return 1.0; 650 // return std::sin(M_PI*x1)*std::sin(M_PI*x2) - std::sin(8.0*M_PI*x1)*std::sin(8.0*M_PI*x2); 651 } 652 653 funcStateSolution(const Real & x1,const Real & x2) const654 Real funcStateSolution(const Real &x1, const Real &x2) const { 655 return std::sin(M_PI*x1)*std::sin(M_PI*x2); 656 } 657 658 funcKappa(const Real & x1,const Real & x2,const std::vector<Real> & par) const659 Real funcKappa(const Real &x1, const Real &x2, const std::vector<Real> &par) const { 660 // Random diffusion coefficient from I. Babuska, F. Nobile, R. Tempone 2010. 661 // Simplified model for random stratified media. 662 Real Lc = 1.0/16.0, sqrtpi = std::sqrt(M_PI), xi = std::sqrt(sqrtpi*Lc); 663 Real sqrt3 = std::sqrt(3.0), f = 0.0, phi = 0.0; 664 Real val = 1.0 + sqrt3*par[0]*std::sqrt(sqrtpi*Lc*0.5); 665 for (int i = 1; i < sdim_-2; i++) { 666 f = floor((Real)(i+1)/2.0); 667 phi = ((i+1)%2 ? std::sin(f*M_PI*x1) : std::cos(f*M_PI*x1)); 668 val += xi*std::exp(-std::pow(f*M_PI*Lc,2.0)/8.0)*phi*sqrt3*par[i]; 669 } 670 return 0.5 + std::exp(val); 671 } 672 updateF(const std::vector<Real> & par)673 void updateF(const std::vector<Real> &par) { 674 Intrepid::FieldContainer<GO> &cellDofs = *(dofMgr_->getCellDofs()); 675 int numLocalDofs = cellDofs.dimension(1); 676 for (GO i=0; i<numCells_; ++i) { // evaluate RHS function at cubature points 677 for (int j=0; j<numCubPoints_; ++j) { 678 (*dataF_)(i, j) = funcRHS((*cubPointsPhysical_)(i, j, 0), 679 (*cubPointsPhysical_)(i, j, 1), 680 par); 681 } 682 } 683 Intrepid::FunctionSpaceTools::integrate<Real>(*datavalVecF_, // compute local data.val vectors for RHS F 684 *dataF_, 685 *valPhysicalWeighted_, 686 Intrepid::COMP_CPP); 687 688 // vecF_ requires assembly using vecF_overlap_ and redistribution 689 vecF_ = ROL::makePtr<Tpetra::MultiVector<>>(matA_->getRangeMap(), 1, true); 690 vecF_overlap_ = ROL::makePtr<Tpetra::MultiVector<>>(myOverlapMap_, 1, true); 691 for (GO i=0; i<numCells_; ++i) { // assembly on the overlap map 692 for (int j=0; j<numLocalDofs; ++j) { 693 vecF_overlap_->sumIntoGlobalValue(cellDofs(myCellIds_[i],j), 694 0, 695 (*datavalVecF_)[i*numLocalDofs+j]); 696 } 697 } 698 Tpetra::Export<> exporter(vecF_overlap_->getMap(), vecF_->getMap()); // redistribution: 699 vecF_->doExport(*vecF_overlap_, exporter, Tpetra::ADD); // from the overlap map to the unique map 700 701 Tpetra::deep_copy(*vecF_dirichlet_, *vecF_); 702 for (int i=0; i<myDirichletDofs_.size(); ++i) { 703 if (myUniqueMap_->isNodeGlobalElement(myDirichletDofs_[i])) { 704 vecF_dirichlet_->replaceGlobalValue(myDirichletDofs_[i], 0, 0); 705 } 706 } 707 } 708 updateA(const std::vector<Real> & par)709 void updateA(const std::vector<Real> &par) { 710 711 Intrepid::FieldContainer<GO> &cellDofs = *(dofMgr_->getCellDofs()); 712 Teuchos::ArrayRCP<const GO> cellDofsArrayRCP = cellDofs.getData(); 713 int numLocalDofs = cellDofs.dimension(1); 714 715 for (GO i=0; i<numCells_; ++i) { // evaluate conductivity kappa at cubature points 716 for (int j=0; j<numCubPoints_; ++j) { 717 (*kappa_)(i, j) = funcKappa((*cubPointsPhysical_)(i, j, 0), 718 (*cubPointsPhysical_)(i, j, 1), 719 par); 720 } 721 } 722 Intrepid::FunctionSpaceTools::tensorMultiplyDataField<Real>(*kappaGradPhysical_, // multiply with conductivity kappa 723 *kappa_, 724 *gradPhysical_); 725 Intrepid::FunctionSpaceTools::integrate<Real>(*gradgradMats_, // compute local grad.(kappa)grad (stiffness) matrices 726 *kappaGradPhysical_, 727 *gradPhysicalWeighted_, 728 Intrepid::COMP_CPP); 729 730 int numLocalMatEntries = numLocalDofs*numLocalDofs; 731 Teuchos::ArrayRCP<const Real> gradgradArrayRCP = gradgradMats_->getData(); 732 matA_->resumeFill(); 733 matA_->setAllToScalar(0); 734 for (int i=0; i<numCells_; ++i) { 735 for (int j=0; j<numLocalDofs; ++j) { 736 matA_->sumIntoGlobalValues(cellDofs(myCellIds_[i],j), 737 cellDofsArrayRCP(myCellIds_[i]*numLocalDofs, numLocalDofs), 738 gradgradArrayRCP(i*numLocalMatEntries+j*numLocalDofs, numLocalDofs)); 739 } 740 } 741 matA_->fillComplete(); 742 743 matA_dirichlet_ = ROL::makePtr<Tpetra::CrsMatrix<>>(*matA_, Teuchos::DataAccess::Copy); 744 matA_dirichlet_->resumeFill(); 745 for (int i=0; i<myDirichletDofs_.size(); ++i) { 746 if (myUniqueMap_->isNodeGlobalElement(myDirichletDofs_[i])) { 747 size_t numRowEntries = matA_dirichlet_->getNumEntriesInGlobalRow(myDirichletDofs_[i]); 748 Teuchos::Array<GO> indices(numRowEntries, 0); 749 Teuchos::Array<Real> values(numRowEntries, 0); 750 Teuchos::Array<Real> canonicalValues(numRowEntries, 0); 751 matA_dirichlet_->getGlobalRowCopy(myDirichletDofs_[i], indices, values, numRowEntries); 752 for (int j=0; j<indices.size(); ++j) { 753 if (myDirichletDofs_[i] == indices[j]) { 754 canonicalValues[j] = 1.0; 755 } 756 } 757 matA_dirichlet_->replaceGlobalValues(myDirichletDofs_[i], indices, canonicalValues); 758 } 759 } 760 matA_dirichlet_->fillComplete(); 761 762 Tpetra::RowMatrixTransposer<> transposerA(matA_dirichlet_); 763 matA_dirichlet_trans_ = transposerA.createTranspose(); 764 765 // Construct solver using Amesos2 factory. 766 try{ 767 solverA_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_); 768 } catch (std::invalid_argument& e) { 769 std::cout << e.what() << std::endl; 770 } 771 try{ 772 solverA_trans_ = Amesos2::create< Tpetra::CrsMatrix<>,Tpetra::MultiVector<> >("KLU2", matA_dirichlet_trans_); 773 } catch (std::invalid_argument& e) { 774 std::cout << e.what() << std::endl; 775 } 776 solverA_->numericFactorization(); 777 solverA_trans_->numericFactorization(); 778 779 } 780 computeStateError(const ROL::Ptr<const Tpetra::MultiVector<>> & soln) const781 Real computeStateError(const ROL::Ptr<const Tpetra::MultiVector<> > &soln) const { 782 783 ROL::Ptr<Tpetra::MultiVector<> > soln_overlap = 784 ROL::makePtr<Tpetra::MultiVector<>>(vecF_overlap_->getMap(), 1, true); 785 Tpetra::Import<> importer(vecUd_->getMap(), soln_overlap->getMap()); // redistribution: 786 soln_overlap->doImport(*soln, importer, Tpetra::REPLACE); // from the unique map to the overlap map 787 788 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 789 int cubDeg = 6; // set cubature degree, e.g., 6 790 ROL::Ptr<Intrepid::Cubature<Real> > cellCub = cubFactory.create(cellType_, cubDeg); // create cubature for error computation 791 int numCubPts = cellCub->getNumPoints(); // retrieve number of cubature points 792 int lfs = dofMgr_->getLocalFieldSize(0); 793 Intrepid::FieldContainer<Real> cubPts(numCubPts, spaceDim_); 794 Intrepid::FieldContainer<Real> cubWts(numCubPts); 795 Intrepid::FieldContainer<Real> cubPtsPhys(numCells_, numCubPts, spaceDim_); 796 Intrepid::FieldContainer<Real> jac(numCells_, numCubPts, spaceDim_, spaceDim_); 797 Intrepid::FieldContainer<Real> jacDet(numCells_, numCubPts); 798 Intrepid::FieldContainer<Real> valRef(lfs, numCubPts); 799 Intrepid::FieldContainer<Real> valPhys(numCells_, lfs, numCubPts); 800 Intrepid::FieldContainer<Real> wtMeas(numCells_, numCubPts); 801 Intrepid::FieldContainer<Real> inCoeffs(numCells_, lfs); 802 Intrepid::FieldContainer<Real> funcVals(numCells_, numCubPts); 803 Intrepid::FieldContainer<Real> funcValsWt(numCells_, numCubPts); 804 Intrepid::FieldContainer<Real> normSquaredError(numCells_); 805 806 cellCub->getCubature(cubPts, cubWts); // retrieve cubature points and weights 807 (*basisPtrs_[0]).getValues(valRef, cubPts, Intrepid::OPERATOR_VALUE); // evaluate value operator at cubature points 808 809 Intrepid::CellTools<Real>::setJacobian(jac, cubPts, *cellNodes_, cellType_); // compute cell Jacobians 810 Intrepid::CellTools<Real>::setJacobianDet(jacDet, jac); // compute determinants of cell Jacobians 811 812 Intrepid::FunctionSpaceTools::HGRADtransformVALUE<Real>(valPhys, // transform reference values into physical space 813 valRef); 814 815 Intrepid::FunctionSpaceTools::computeCellMeasure<Real>(wtMeas, // compute weighted cell measure 816 jacDet, 817 cubWts); 818 819 Intrepid::CellTools<Real>::mapToPhysicalFrame(cubPtsPhys, // map reference cubature points to physical space 820 cubPts, 821 *cellNodes_, 822 cellType_); 823 824 Intrepid::FieldContainer<GO> &cellDofs = *(dofMgr_->getCellDofs()); 825 Teuchos::ArrayRCP<const Real> soln_data = soln_overlap->get1dView(); // populate inCoeffs 826 for (GO i=0; i<numCells_; ++i) { 827 for (int j=0; j<lfs; ++j) { 828 inCoeffs(i, j) = soln_data[soln_overlap->getMap()->getLocalElement(cellDofs(myCellIds_[i],j))]; 829 } 830 } 831 832 Intrepid::FunctionSpaceTools::evaluate<Real>(funcVals, inCoeffs, valPhys); 833 834 for (GO i=0; i<numCells_; ++i) { // compute error 835 for (int j=0; j<numCubPts; ++j) { 836 funcVals(i, j) -= funcStateSolution(cubPtsPhys(i, j, 0), cubPtsPhys(i, j, 1)); 837 } 838 } 839 840 Intrepid::FunctionSpaceTools::scalarMultiplyDataData<Real>(funcValsWt, // multiply with weighted measure 841 wtMeas, 842 funcVals); 843 844 Intrepid::FunctionSpaceTools::integrate<Real>(normSquaredError, // compute norm squared of local error 845 funcVals, 846 funcValsWt, 847 Intrepid::COMP_CPP); 848 849 Real localErrorSum(0); 850 Real globalErrorSum(0); 851 for (GO i=0; i<numCells_; ++i) { 852 localErrorSum += normSquaredError(i); 853 } 854 ROL::Ptr<const Teuchos::Comm<int> > comm = soln_overlap->getMap()->getComm(); 855 Teuchos::reduceAll<int, Real>(*comm, Teuchos::REDUCE_SUM, 1, &localErrorSum, &globalErrorSum); 856 857 return globalErrorSum; 858 } 859 860 printMeshData(std::ostream & outStream) const861 void printMeshData(std::ostream &outStream) const { 862 ROL::Ptr<Intrepid::FieldContainer<Real> > nodesPtr = meshMgr_->getNodes(); 863 ROL::Ptr<Intrepid::FieldContainer<GO> > cellToNodeMapPtr = meshMgr_->getCellToNodeMap(); 864 Intrepid::FieldContainer<Real> &nodes = *nodesPtr; 865 Intrepid::FieldContainer<GO> &cellToNodeMap = *cellToNodeMapPtr; 866 outStream << "Number of nodes = " << meshMgr_->getNumNodes() << std::endl; 867 outStream << "Number of cells = " << meshMgr_->getNumCells() << std::endl; 868 outStream << "Number of edges = " << meshMgr_->getNumEdges() << std::endl; 869 // Print mesh to file. 870 if ((myRank_ == 0)) { 871 std::ofstream meshfile; 872 meshfile.open("cell_to_node_quad.txt"); 873 for (int i=0; i<cellToNodeMap.dimension(0); ++i) { 874 for (int j=0; j<cellToNodeMap.dimension(1); ++j) { 875 meshfile << cellToNodeMap(i,j) << " "; 876 } 877 meshfile << std::endl; 878 } 879 meshfile.close(); 880 meshfile.open("cell_to_node_tri.txt"); 881 for (int i=0; i<cellToNodeMap.dimension(0); ++i) { 882 for (int j=0; j<3; ++j) { 883 meshfile << cellToNodeMap(i,j) << " "; 884 } 885 meshfile << std::endl; 886 for (int j=2; j<5; ++j) { 887 meshfile << cellToNodeMap(i,j%4) << " "; 888 } 889 meshfile << std::endl; 890 } 891 meshfile.close(); 892 meshfile.open("nodes.txt"); 893 meshfile.precision(16); 894 for (int i=0; i<nodes.dimension(0); ++i) { 895 for (int j=0; j<nodes.dimension(1); ++j) { 896 meshfile << std::scientific << nodes(i,j) << " "; 897 } 898 meshfile << std::endl; 899 } 900 meshfile.close(); 901 /* This somewhat clunky output is for gnuplot. 902 meshfile.open("mesh.txt"); 903 for (int i=0; i<cellToNodeMap.dimension(0); ++i) { 904 meshfile << nodes(cellToNodeMap(i,0), 0) << " " << nodes(cellToNodeMap(i,0), 1) << std::endl; 905 meshfile << nodes(cellToNodeMap(i,1), 0) << " " << nodes(cellToNodeMap(i,1), 1) << std::endl; 906 meshfile << nodes(cellToNodeMap(i,2), 0) << " " << nodes(cellToNodeMap(i,2), 1) << std::endl; 907 meshfile << nodes(cellToNodeMap(i,3), 0) << " " << nodes(cellToNodeMap(i,3), 1) << std::endl; 908 meshfile << nodes(cellToNodeMap(i,0), 0) << " " << nodes(cellToNodeMap(i,0), 1) << std::endl; 909 meshfile << nodes(cellToNodeMap(i,1), 0) << " " << nodes(cellToNodeMap(i,1), 1) << std::endl; 910 meshfile << nodes(cellToNodeMap(i,2), 0) << " " << nodes(cellToNodeMap(i,2), 1) << std::endl; 911 } 912 meshfile.close(); 913 */ 914 } 915 } 916 917 outputTpetraData() const918 void outputTpetraData() const { 919 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<> > matWriter; 920 matWriter.writeSparseFile("stiffness_mat", matA_); 921 matWriter.writeSparseFile("dirichlet_mat", matA_dirichlet_); 922 matWriter.writeSparseFile("mass_mat", matM_); 923 matWriter.writeDenseFile("Ud_vec", vecUd_); 924 } 925 926 outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<>> & vec,const std::string & filename) const927 void outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<> > &vec, 928 const std::string &filename) const { 929 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<> > vecWriter; 930 vecWriter.writeDenseFile(filename, vec); 931 } 932 933 }; // class StefanBoltzmannData 934 935 #endif 936