diff --git a/packages/panzer/adapters-stk/src/Panzer_STK_LocalMeshUtilities.cpp b/packages/panzer/adapters-stk/src/Panzer_STK_LocalMeshUtilities.cpp index 6be6923ac1b3..84ad1b20fd40 100644 --- a/packages/panzer/adapters-stk/src/Panzer_STK_LocalMeshUtilities.cpp +++ b/packages/panzer/adapters-stk/src/Panzer_STK_LocalMeshUtilities.cpp @@ -256,12 +256,15 @@ buildNodeToCellMatrix(const Teuchos::RCP > & comm, RCP cell_to_node; { PANZER_FUNC_TIME_MONITOR_DIFF("Build matrix",BuildMatrix); - // The matrix is indexed by (global cell, global node) = local node - cell_to_node = rcp(new crs_type(cell_map,0)); // fill in the cell to node matrix const unsigned int num_local_cells = owned_cells_to_nodes.extent(0); const unsigned int num_nodes_per_cell = owned_cells_to_nodes.extent(1); + + // The matrix is indexed by (global cell, global node) = local node + cell_to_node = rcp(new crs_type(cell_map,num_nodes_per_cell, + Tpetra::StaticProfile)); + std::vector local_node_indexes(num_nodes_per_cell); std::vector global_node_indexes(num_nodes_per_cell); for(unsigned int i=0;i::invalid(),indices,0,comm_)); } - RCP ghostedGraph = rcp(new GraphType(ghostedTargetMap,ghostedSourceMap,0)); // Now insert the non-zero pattern per row + // count number of entries per row; required by CrsGraph constructor + std::vector nEntriesPerRow(ghostedTargetMap->getNodeNumElements(),0); std::vector elementBlockIds; targetGlobalIndexer_->getElementBlockIds(elementBlockIds); std::vector::const_iterator blockItr; + for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { + std::string blockId = *blockItr; + const std::vector & elements = targetGlobalIndexer_->getElementBlock(blockId); + + std::vector row_gids; + std::vector col_gids; + + for(std::size_t elmt=0;elmtgetElementGIDs(elements[elmt],row_gids); + sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids); + for(std::size_t row=0;rowgetLocalElement(row_gids[row]); + nEntriesPerRow[lid] += col_gids.size(); + } + } + } + + Teuchos::ArrayView nEntriesPerRowView(nEntriesPerRow); + RCP ghostedGraph = rcp(new GraphType(ghostedTargetMap,ghostedSourceMap,nEntriesPerRowView,Tpetra::StaticProfile)); + for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { std::string blockId = *blockItr; const std::vector & elements = targetGlobalIndexer_->getElementBlock(blockId); diff --git a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp index ad42eeb46c83..94c9cd309d5b 100644 --- a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp +++ b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp @@ -1001,8 +1001,6 @@ buildTpetraGhostedGraph(int i,int j) const RCP map_i = getGhostedMap(i); RCP map_j = getGhostedMap(j); - RCP graph = rcp(new CrsGraphType(map_i,map_j,0)); - std::vector elementBlockIds; Teuchos::RCP rowProvider, colProvider; @@ -1013,8 +1011,37 @@ buildTpetraGhostedGraph(int i,int j) const gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the // same element blocks - // graph information about the mesh + // Count number of entries in each row of graph; needed for graph constructor + std::vector nEntriesPerRow(map_i->getNodeNumElements(), 0); std::vector::const_iterator blockItr; + for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { + std::string blockId = *blockItr; + // grab elements for this block + const std::vector & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the + // same elements in each element block + + // get information about number of indicies + std::vector row_gids; + std::vector col_gids; + + // loop over the elemnts + for(std::size_t elmt=0;elmtgetElementGIDs(elements[elmt],row_gids); + colProvider->getElementGIDs(elements[elmt],col_gids); + for(std::size_t row=0;rowgetLocalElement(row_gids[row]); + nEntriesPerRow[lid] += col_gids.size(); + } + } + } + Teuchos::ArrayView nEntriesPerRowView(nEntriesPerRow); + RCP graph = rcp(new CrsGraphType(map_i,map_j, nEntriesPerRowView, + Tpetra::StaticProfile)); + + + + // graph information about the mesh for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { std::string blockId = *blockItr; diff --git a/packages/panzer/disc-fe/src/lof/Panzer_TpetraLinearObjFactory_impl.hpp b/packages/panzer/disc-fe/src/lof/Panzer_TpetraLinearObjFactory_impl.hpp index 6e8b459c1a44..3bd0682ee710 100644 --- a/packages/panzer/disc-fe/src/lof/Panzer_TpetraLinearObjFactory_impl.hpp +++ b/packages/panzer/disc-fe/src/lof/Panzer_TpetraLinearObjFactory_impl.hpp @@ -664,7 +664,6 @@ buildGhostedGraph() const // build the map and allocate the space for the graph Teuchos::RCP rMap = getGhostedMap(); Teuchos::RCP cMap = getGhostedColMap(); - Teuchos::RCP graph = Teuchos::rcp(new CrsGraphType(rMap,cMap,0)); std::vector elementBlockIds; gidProvider_->getElementBlockIds(elementBlockIds); @@ -675,6 +674,9 @@ buildGhostedGraph() const const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors(); // graph information about the mesh + // Count number of entries per graph row; needed for graph constructor + std::vector nEntriesPerRow(rMap->getNodeNumElements(), 0); + std::vector::const_iterator blockItr; for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { std::string blockId = *blockItr; @@ -686,6 +688,44 @@ buildGhostedGraph() const std::vector gids; std::vector col_gids; + // loop over the elemnts + for(std::size_t i=0;igetElementGIDs(elements[i],gids); + + colGidProvider->getElementGIDs(elements[i],col_gids); + if (han) { + const std::vector& aes = conn_mgr->getAssociatedNeighbors(elements[i]); + for (typename std::vector::const_iterator eit = aes.begin(); + eit != aes.end(); ++eit) { + std::vector other_col_gids; + colGidProvider->getElementGIDs(*eit, other_col_gids); + col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end()); + } + } + + for(std::size_t j=0;jgetLocalElement(gids[j]); + nEntriesPerRow[lid] += col_gids.size(); + } + } + } + + Teuchos::ArrayView nEntriesPerRowView(nEntriesPerRow); + Teuchos::RCP graph = Teuchos::rcp(new CrsGraphType(rMap,cMap, + nEntriesPerRowView, + Tpetra::StaticProfile)); + + // Now insert entries into the graph + for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) { + std::string blockId = *blockItr; + + // grab elements for this block + const std::vector & elements = gidProvider_->getElementBlock(blockId); + + // get information about number of indicies + std::vector gids; + std::vector col_gids; + // loop over the elemnts for(std::size_t i=0;igetElementGIDs(elements[i],gids);