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