diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cecf040a7b..3c8f2dac67 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -41,6 +41,7 @@ set (ALBANY_INCLUDE_DIRS ${Albany_SOURCE_DIR}/src/problems ${Albany_SOURCE_DIR}/src/responses ${Albany_SOURCE_DIR}/src/disc/stk + ${Albany_SOURCE_DIR}/src/disc/extruded ${Albany_SOURCE_DIR}/src/disc ${Albany_SOURCE_DIR}/src/utility ) @@ -234,9 +235,9 @@ list (APPEND SOURCES disc/Albany_MeshSpecs.cpp disc/Albany_DOFManager.cpp disc/Albany_DiscretizationUtils.cpp - disc/Albany_ExtrudedMesh.cpp - disc/Albany_ExtrudedConnManager.cpp - disc/Albany_ExtrudedDiscretization.cpp + disc/extruded/Albany_ExtrudedConnManager.cpp + disc/extruded/Albany_ExtrudedDiscretization.cpp + disc/extruded/Albany_ExtrudedMesh.cpp ) list (APPEND HEADERS disc/Albany_DiscretizationUtils.hpp @@ -247,9 +248,9 @@ list (APPEND HEADERS disc/Albany_MeshSpecs.hpp disc/Albany_ConnManager.hpp disc/Albany_DOFManager.hpp - disc/Albany_ExtrudedMesh.hpp - disc/Albany_ExtrudedConnManager.hpp - disc/Albany_ExtrudedDiscretization.hpp + disc/extruded/Albany_ExtrudedConnManager.hpp + disc/extruded/Albany_ExtrudedDiscretization.hpp + disc/extruded/Albany_ExtrudedMesh.hpp ) #stk diff --git a/src/disc/Albany_AbstractDiscretization.cpp b/src/disc/Albany_AbstractDiscretization.cpp index 2e87ec9568..fc938100ea 100644 --- a/src/disc/Albany_AbstractDiscretization.cpp +++ b/src/disc/Albany_AbstractDiscretization.cpp @@ -50,4 +50,294 @@ writeSolutionMV (const Thyra_MultiVector& soln, writeMeshDatabaseToFile(time, force_write_solution); } +auto AbstractDiscretization:: +get_dof_mgr (const std::string& part_name, + const FE_Type fe_type, + const int order, + const int dof_dim) + -> dof_mgr_ptr_t& +{ + // NOTE: we assume order<10, and dof_dim<10, which is virtually never going to change + int type_order_dim = 100*static_cast(fe_type) + 10*order + dof_dim; + + std::string key = part_name + "_" + std::to_string(type_order_dim); + return m_key_to_dof_mgr[key]; +} + +void AbstractDiscretization::buildSideSetsViews () +{ + GlobalSideSetList globalSideSetViews; + std::map> allLocalDOFViews; + + // (Kokkos Refactor) Convert sideSets to sideSetViews + + // 1) Compute view extents (num_local_worksets, max_sideset_length, max_sides) and local workset counter (current_local_index) + std::map num_local_worksets; + std::map max_sideset_length; + std::map max_sides; + std::map current_local_index; + for (const auto& ss : m_sideSets) { + for (const auto& [ss_key, ss_val] : ss) { + // Initialize values if this is the first time seeing a sideset key + if (num_local_worksets.find(ss_key) == num_local_worksets.end()) + num_local_worksets[ss_key] = 0; + if (max_sideset_length.find(ss_key) == max_sideset_length.end()) + max_sideset_length[ss_key] = 0; + if (max_sides.find(ss_key) == max_sides.end()) + max_sides[ss_key] = 0; + if (current_local_index.find(ss_key) == current_local_index.end()) + current_local_index[ss_key] = 0; + + // Update extents for given workset/sideset + num_local_worksets[ss_key]++; + max_sideset_length[ss_key] = std::max(max_sideset_length[ss_key], (int) ss_val.size()); + for (size_t j = 0; j < ss_val.size(); ++j) + max_sides[ss_key] = std::max(max_sides[ss_key], (int) ss_val[j].side_pos); + } + } + + // 2) Construct GlobalSideSetList (map of GlobalSideSetInfo) + for (const auto& ss_it : num_local_worksets) { + std::string ss_key = ss_it.first; + + max_sides[ss_key]++; // max sides is the largest local ID + 1 and needs to be incremented once for each key here + + globalSideSetViews[ss_key].num_local_worksets = num_local_worksets[ss_key]; + globalSideSetViews[ss_key].max_sideset_length = max_sideset_length[ss_key]; + globalSideSetViews[ss_key].side_GID = Kokkos::DualView("side_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].elem_GID = Kokkos::DualView("elem_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].ws_elem_idx = Kokkos::DualView("ws_elem_idx", num_local_worksets[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].elem_ebIndex = Kokkos::DualView("elem_ebIndex", num_local_worksets[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].side_pos = Kokkos::DualView("side_pos", num_local_worksets[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].max_sides = max_sides[ss_key]; + globalSideSetViews[ss_key].numCellsOnSide = Kokkos::DualView("numCellsOnSide", num_local_worksets[ss_key], max_sides[ss_key]); + globalSideSetViews[ss_key].cellsOnSide = Kokkos::DualView("cellsOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); + globalSideSetViews[ss_key].sideSetIdxOnSide = Kokkos::DualView("sideSetIdxOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); + } + + // 3) Populate global views + for (const auto& ss : m_sideSets) { + for (const auto& [ss_key, ss_val] : ss) { + int current_index = current_local_index[ss_key]; + int numSides = max_sides[ss_key]; + + int max_cells_on_side = 0; + std::vector numCellsOnSide(numSides); + std::vector> cellsOnSide(numSides); + std::vector> sideSetIdxOnSide(numSides); + for (size_t j = 0; j < ss_val.size(); ++j) { + int cell = ss_val[j].ws_elem_idx; + int side = ss_val[j].side_pos; + + cellsOnSide[side].push_back(cell); + sideSetIdxOnSide[side].push_back(j); + } + for (int side = 0; side < numSides; ++side) { + numCellsOnSide[side] = cellsOnSide[side].size(); + max_cells_on_side = std::max(max_cells_on_side, numCellsOnSide[side]); + } + + for (int side = 0; side < numSides; ++side) { + globalSideSetViews[ss_key].numCellsOnSide.h_view(current_index, side) = numCellsOnSide[side]; + for (int j = 0; j < numCellsOnSide[side]; ++j) { + globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = cellsOnSide[side][j]; + globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = sideSetIdxOnSide[side][j]; + } + for (int j = numCellsOnSide[side]; j < max_sideset_length[ss_key]; ++j) { + globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = -1; + globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = -1; + } + } + + for (size_t j = 0; j < ss_val.size(); ++j) { + globalSideSetViews[ss_key].side_GID.h_view(current_index, j) = ss_val[j].side_GID; + globalSideSetViews[ss_key].elem_GID.h_view(current_index, j) = ss_val[j].elem_GID; + globalSideSetViews[ss_key].ws_elem_idx.h_view(current_index, j) = ss_val[j].ws_elem_idx; + globalSideSetViews[ss_key].elem_ebIndex.h_view(current_index, j) = ss_val[j].elem_ebIndex; + globalSideSetViews[ss_key].side_pos.h_view(current_index, j) = ss_val[j].side_pos; + } + + globalSideSetViews[ss_key].side_GID.modify_host(); + globalSideSetViews[ss_key].elem_GID.modify_host(); + globalSideSetViews[ss_key].ws_elem_idx.modify_host(); + globalSideSetViews[ss_key].elem_ebIndex.modify_host(); + globalSideSetViews[ss_key].side_pos.modify_host(); + globalSideSetViews[ss_key].numCellsOnSide.modify_host(); + globalSideSetViews[ss_key].cellsOnSide.modify_host(); + globalSideSetViews[ss_key].sideSetIdxOnSide.modify_host(); + + globalSideSetViews[ss_key].side_GID.sync_device(); + globalSideSetViews[ss_key].elem_GID.sync_device(); + globalSideSetViews[ss_key].ws_elem_idx.sync_device(); + globalSideSetViews[ss_key].elem_ebIndex.sync_device(); + globalSideSetViews[ss_key].side_pos.sync_device(); + globalSideSetViews[ss_key].numCellsOnSide.sync_device(); + globalSideSetViews[ss_key].cellsOnSide.sync_device(); + globalSideSetViews[ss_key].sideSetIdxOnSide.sync_device(); + + current_local_index[ss_key]++; + } + } + + // 4) Reset current_local_index + std::map::iterator counter_it = current_local_index.begin(); + while (counter_it != current_local_index.end()) { + std::string counter_key = counter_it->first; + current_local_index[counter_key] = 0; + counter_it++; + } + + // 5) Populate map of LocalSideSetInfos + for (size_t i = 0; i < m_sideSets.size(); ++i) { + LocalSideSetInfoList& lssList = m_sideSetViews[i]; + + for (const auto& [ss_key,ss_val] : m_sideSets[i]) { + int current_index = current_local_index[ss_key]; + std::pair range(0, ss_val.size()); + + lssList[ss_key].size = ss_val.size(); + lssList[ss_key].side_GID = Kokkos::subview(globalSideSetViews[ss_key].side_GID, current_index, range ); + lssList[ss_key].elem_GID = Kokkos::subview(globalSideSetViews[ss_key].elem_GID, current_index, range ); + lssList[ss_key].ws_elem_idx = Kokkos::subview(globalSideSetViews[ss_key].ws_elem_idx, current_index, range ); + lssList[ss_key].elem_ebIndex = Kokkos::subview(globalSideSetViews[ss_key].elem_ebIndex, current_index, range ); + lssList[ss_key].side_pos = Kokkos::subview(globalSideSetViews[ss_key].side_pos, current_index, range ); + lssList[ss_key].numSides = globalSideSetViews[ss_key].max_sides; + lssList[ss_key].numCellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].numCellsOnSide, current_index, Kokkos::ALL() ); + lssList[ss_key].cellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].cellsOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); + lssList[ss_key].sideSetIdxOnSide = Kokkos::subview(globalSideSetViews[ss_key].sideSetIdxOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); + + current_local_index[ss_key]++; + } + } + + // 6) Determine size of global DOFView structure and allocate + std::map total_sideset_idx; + std::map sideset_idx_offset; + unsigned int maxSideNodes = 0; + const auto& cell_layers_data = getMeshStruct()->local_cell_layers_data; + if (!cell_layers_data.is_null()) { + const Teuchos::RCP cell_topo = Teuchos::rcp(new CellTopologyData(getMeshStruct()->meshSpecs[0]->ctd)); + const int numLayers = cell_layers_data->numLayers; + const int numComps = getDOFManager()->getNumFields(); + + // Determine maximum number of side nodes + for (unsigned int elem_side = 0; elem_side < cell_topo->side_count; ++elem_side) { + const CellTopologyData_Subcell& side = cell_topo->side[elem_side]; + const unsigned int numSideNodes = side.topology->node_count; + maxSideNodes = std::max(maxSideNodes, numSideNodes); + } + + // Determine total number of sideset indices per each sideset name + for (auto& ssList : m_sideSets) { + for (auto& ss_it : ssList) { + std::string ss_key = ss_it.first; + std::vector ss_val = ss_it.second; + + if (sideset_idx_offset.find(ss_key) == sideset_idx_offset.end()) + sideset_idx_offset[ss_key] = 0; + if (total_sideset_idx.find(ss_key) == total_sideset_idx.end()) + total_sideset_idx[ss_key] = 0; + + total_sideset_idx[ss_key] += ss_val.size(); + } + } + + // Allocate total localDOFView for each sideset name + for (auto& ss_it : num_local_worksets) { + std::string ss_key = ss_it.first; + allLocalDOFViews[ss_key] = Kokkos::DualView(ss_key + " localDOFView", total_sideset_idx[ss_key], maxSideNodes, numLayers+1, numComps); + } + } + + // Not all mesh structs that come through here are extruded mesh structs. + // If the mesh isn't extruded, we won't need to do any of the following work. + if (not cell_layers_data.is_null()) { + // Get topo data + auto ctd = getMeshStruct()->meshSpecs[0]->ctd; + + // Ensure we have ONE cell per layer. + const auto topo_hexa = shards::getCellTopologyData>(); + const auto topo_wedge = shards::getCellTopologyData>(); + TEUCHOS_TEST_FOR_EXCEPTION ( + ctd.name!=topo_hexa->name && + ctd.name!=topo_wedge->name, std::runtime_error, + "Extruded meshes only allowed if there is one element per layer (hexa or wedges).\n" + " - current topology name: " << ctd.name << "\n"); + + const auto& sol_dof_mgr = getDOFManager(); + const auto& elem_dof_lids = sol_dof_mgr->elem_dof_lids().host(); + + // Build a LayeredMeshNumbering for cells, so we can get the LIDs of elems over the column + const auto numLayers = cell_layers_data->numLayers; + const int top = cell_layers_data->top_side_pos; + const int bot = cell_layers_data->bot_side_pos; + + // 7) Populate localDOFViews for GatherVerticallyContractedSolution + for (int ws=0; ws>& wsldofViews = m_wsLocalDOFViews[ws]; + + const auto& elem_lids = getElementLIDs_host(ws); + + // Loop over the sides that form the boundary condition + // const Teuchos::ArrayRCP >& wsElNodeID_i = wsElNodeID[i]; + for (const auto& [ss_key,ss_val] : m_sideSets[ws]) { + Kokkos::DualView& globalDOFView = allLocalDOFViews[ss_key]; + + for (unsigned int sideSet_idx = 0; sideSet_idx < ss_val.size(); ++sideSet_idx) { + const auto& side = ss_val[sideSet_idx]; + + // Get the data that corresponds to the side + const int ws_elem_idx = side.ws_elem_idx; + const int side_pos = side.side_pos; + + // Check if this sideset is the top or bot of the mesh. If not, the data structure + // for coupling vertical dofs is not needed. + if (side_pos!=top && side_pos!=bot) + break; + + const int elem_LID = elem_lids(ws_elem_idx); + const int basal_elem_LID = cell_layers_data->getColumnId(elem_LID); + + for (int eq=0; eqgetGIDFieldOffsetsSide(eq,top,side_pos); + const auto& sol_bot_offsets = sol_dof_mgr->getGIDFieldOffsetsSide(eq,bot,side_pos); + const int numSideNodes = sol_top_offsets.size(); + + for (int j=0; jgetId(basal_elem_LID,il); + globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il, eq) = + elem_dof_lids(layer_elem_LID,sol_bot_offsets[j]); + } + + // Add top side in last layer + const int il = numLayers-1; + const LO layer_elem_LID = cell_layers_data->getId(basal_elem_LID,il); + globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il+1, eq) = + elem_dof_lids(layer_elem_LID,sol_top_offsets[j]); + } + } + } + + globalDOFView.modify_host(); + globalDOFView.sync_device(); + + // Set workset-local sub-view + std::pair range(sideset_idx_offset[ss_key], sideset_idx_offset[ss_key]+ss_val.size()); + wsldofViews[ss_key] = Kokkos::subview(globalDOFView, range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + + sideset_idx_offset[ss_key] += ss_val.size(); + } + } + } else { + // We still need this view to be present (even if of size 0), so create them + std::map> dummy; + for (int ws=0; ws>& - getLocalDOFViews(const int workset) const = 0; + const std::map>& + getLocalDOFViews(const int workset) const + { + return m_wsLocalDOFViews.at(workset); + } //! Get Dof Manager of field field_name Teuchos::RCP @@ -219,12 +221,14 @@ class AbstractDiscretization const strmap_t>& getSideSetDiscretizations() const { return sideSetDiscretizations; } //! Get the map side_id->side_set_elem_id - virtual const std::map>& - getSideToSideSetCellMap() const = 0; + const strmap_t>& getSideToSideSetCellMap() const { + return m_side_to_ss_cell; + } //! Get the map side_node_id->side_set_cell_node_id - virtual const std::map>>& - getSideNodeNumerationMap() const = 0; + const strmap_t>>& getSideNodeNumerationMap() const { + return m_side_nodes_to_ss_cell_nodes; + } //! Get MeshStruct virtual Teuchos::RCP @@ -262,8 +266,7 @@ class AbstractDiscretization getNumDim() const = 0; //! Get number of total DOFs per node - virtual int - getNumEq() const = 0; + int getNumEq() const { return m_neq; } //! Get Numbering for layered mesh (mesh structured in one direction) virtual Teuchos::RCP> @@ -358,19 +361,51 @@ class AbstractDiscretization virtual void adapt (const Teuchos::RCP& adaptData) = 0; protected: + dof_mgr_ptr_t& + get_dof_mgr (const std::string& part_name, + const FE_Type fe_type, + const int order, + const int dof_dim); + + // From std::vector build corresponding kokkos structures + void buildSideSetsViews (); + strmap_t> sideSetDiscretizations; //! Jacobian matrix operator factory Teuchos::RCP m_jac_factory; - // One dof mgr per dof per part // Notice that the dof mgr on a side is not the restriction // of the volume dof mgr to that side, since local ids are different. + // Note: the double map works as map[field_name][part_name] = dof_mgr strmap_t> m_dof_managers; - // Dof manager for a scalar node field + // Dof manager for a scalar node field (part_name->dof_mgr) strmap_t m_node_dof_managers; + // Store a all dof mgrs based on a key that encodes all params used to create it. + // This helps to build only one copy of dof mgrs with same specs + std::map m_key_to_dof_mgr; + + // For each ss, map the side_GID into vec, where vec[i]=k if the i-th + // node of the side corresponds to the k-th node of the cell in the side mesh + strmap_t>> m_side_nodes_to_ss_cell_nodes; + + // For each ss, map the side_GID to the cell_GID in the side mesh + strmap_t> m_side_to_ss_cell; + + //! side sets stored as std::map(string ID, SideArray classes) per workset + //! (std::vector across worksets) + std::vector m_sideSets; + std::map m_sideSetViews; + + //! Number of equations (and unknowns) per node + // TODO: this should soon be removed, in favor of more granular description of each dof/unknown + int m_neq; + + //! GatherVerticallyContractedSolution connectivity + std::map>> m_wsLocalDOFViews; + // Workset information WorksetArray m_workset_sizes; // size of each ws WorksetArray m_wsEBNames; // name of elem block that ws belongs diff --git a/src/disc/Albany_AbstractMeshStruct.hpp b/src/disc/Albany_AbstractMeshStruct.hpp index c56f1e1daf..c80e6339c3 100644 --- a/src/disc/Albany_AbstractMeshStruct.hpp +++ b/src/disc/Albany_AbstractMeshStruct.hpp @@ -37,11 +37,14 @@ struct AbstractMeshStruct { //! Internal mesh specs type needed virtual std::string meshLibName() const = 0; - virtual void setFieldData (const Teuchos::RCP& comm) = 0; + virtual void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) = 0; virtual void setBulkData(const Teuchos::RCP& comm) = 0; bool isBulkDataSet () const { return m_bulk_data_set; } + bool isFieldDataSet () const { return m_field_data_set; } virtual LO get_num_local_elements () const = 0; virtual LO get_num_local_nodes () const = 0; @@ -58,6 +61,7 @@ struct AbstractMeshStruct { std::map> sideSetMeshStructs; protected: + bool m_field_data_set = false; bool m_bulk_data_set = false; }; diff --git a/src/disc/Albany_DiscretizationFactory.cpp b/src/disc/Albany_DiscretizationFactory.cpp index 7ade403d39..99483f0d74 100644 --- a/src/disc/Albany_DiscretizationFactory.cpp +++ b/src/disc/Albany_DiscretizationFactory.cpp @@ -100,7 +100,7 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis //FixME very hacky! needed for printing 2d mesh Teuchos::RCP meshStruct2D; meshStruct2D = Teuchos::rcp(new AsciiSTKMesh2D(disc_params, comm, numParams)); - meshStruct2D->setFieldData(comm); + meshStruct2D->setFieldData(comm,Teuchos::null,{}); meshStruct2D->setBulkData(comm); Ioss::Init::Initializer io; Teuchos::RCP mesh_data = Teuchos::rcp(new stk::io::StkMeshIoBroker(MPI_COMM_WORLD)); @@ -188,12 +188,10 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis Teuchos::RCP params_ss; int sideDim = ms->numDim - 1; - for (int i(0); isideSetMeshStructs[ss_name]; - // If this is the basalside of an extruded mesh, we already created the mesh object + // Extruded meshes can create some side meshes, so check if mesh is already there if (ss_mesh.is_null()) { params_ss = Teuchos::rcp(new Teuchos::ParameterList(ssd_list.sublist(ss_name))); @@ -285,14 +283,7 @@ DiscretizationFactory::setMeshStructFieldData( const std::map >& side_set_sis) { TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: setMeshStructFieldData"); - meshStruct->setFieldData(comm); - meshStruct->get_field_accessor()->addStateStructs(sis); - for (auto& [ss_name, ss_mesh] : meshStruct->sideSetMeshStructs) { - ss_mesh->setFieldData(comm); - } - for (auto& [ss_name, ss_sis] : side_set_sis) { - meshStruct->sideSetMeshStructs.at(ss_name)->get_field_accessor()->addStateStructs(ss_sis); - } + meshStruct->setFieldData(comm,sis,side_set_sis); } void DiscretizationFactory:: diff --git a/src/disc/Albany_ExtrudedConnManager.cpp b/src/disc/extruded/Albany_ExtrudedConnManager.cpp similarity index 100% rename from src/disc/Albany_ExtrudedConnManager.cpp rename to src/disc/extruded/Albany_ExtrudedConnManager.cpp diff --git a/src/disc/Albany_ExtrudedConnManager.hpp b/src/disc/extruded/Albany_ExtrudedConnManager.hpp similarity index 100% rename from src/disc/Albany_ExtrudedConnManager.hpp rename to src/disc/extruded/Albany_ExtrudedConnManager.hpp diff --git a/src/disc/Albany_ExtrudedDiscretization.cpp b/src/disc/extruded/Albany_ExtrudedDiscretization.cpp similarity index 70% rename from src/disc/Albany_ExtrudedDiscretization.cpp rename to src/disc/extruded/Albany_ExtrudedDiscretization.cpp index 020a4aa061..d4b05a97e0 100644 --- a/src/disc/Albany_ExtrudedDiscretization.cpp +++ b/src/disc/extruded/Albany_ExtrudedDiscretization.cpp @@ -38,12 +38,13 @@ ExtrudedDiscretization (const Teuchos::RCP& discPara const std::map>& sideSetEquations) : m_comm(comm) , m_basal_disc (basal_disc) - , m_neq (neq) , m_sideSetEquations(sideSetEquations) , m_rigid_body_modes(rigidBodyModes) , m_extruded_mesh(extruded_mesh) , m_disc_params (discParams) { + m_neq = neq; + sideSetDiscretizations["basalside"] = basal_disc; } @@ -345,9 +346,17 @@ void ExtrudedDiscretization::createDOFManagers() // dof part name is "", we get the part stored in the stk mesh struct // for the element block, where we REQUIRE that there is only ONE element block. - strmap_t> name_to_partAndDim; - name_to_partAndDim[solution_dof_name()] = std::make_pair("",m_neq); - name_to_partAndDim[nodes_dof_name()] = std::make_pair("",1); + Teuchos::RCP dof_mgr; + + // Solution dof mgr + dof_mgr = create_dof_mgr("",solution_dof_name(),FE_Type::HGRAD,1,m_neq); + m_dof_managers[solution_dof_name()][""] = dof_mgr; + + // Nodes dof mgr + dof_mgr = create_dof_mgr("",nodes_dof_name(),FE_Type::HGRAD,1,1); + m_dof_managers[nodes_dof_name()][""] = dof_mgr; + m_node_dof_managers[""] = dof_mgr; + for (const auto& sis : m_extruded_mesh->get_field_accessor()->getNodalParameterSIS()) { const auto& dims = sis->dim; int dof_dim = -1; @@ -360,25 +369,11 @@ void ExtrudedDiscretization::createDOFManagers() "Error! Unsupported layout for nodal parameter '" + sis->name + ".\n"); } - name_to_partAndDim[sis->name] = std::make_pair(sis->meshPart,dof_dim); - } + dof_mgr = create_dof_mgr(sis->meshPart,sis->name,FE_Type::HGRAD,1,dof_dim); + m_dof_managers[sis->name][sis->meshPart] = dof_mgr; - for (const auto& it : name_to_partAndDim) { - const auto& field_name = it.first; - const auto& part_name = it.second.first; - const auto& dof_dim = it.second.second; - - // NOTE: for now we hard code P1. In the future, we must be able to - // store this info somewhere and retrieve it here. - auto dof_mgr = create_dof_mgr(part_name,field_name,FE_Type::HGRAD,1,dof_dim); - m_dof_managers[field_name][part_name] = dof_mgr; - m_node_dof_managers[part_name] = Teuchos::null; - } - - // For each part, also make a Node dof manager - for (auto& it : m_node_dof_managers) { - const auto& part_name = it.first; - it.second = create_dof_mgr(part_name, nodes_dof_name(), FE_Type::HGRAD,1,1); + dof_mgr = create_dof_mgr(sis->meshPart,sis->name,FE_Type::HGRAD,1,1); + m_node_dof_managers[sis->meshPart] = dof_mgr; } } @@ -775,278 +770,7 @@ ExtrudedDiscretization::computeSideSets() } } - // ============================================================= - // (Kokkos Refactor) Convert sideSets to sideSetViews - - // 1) Compute view extents (num_local_worksets, max_sideset_length, max_sides) and local workset counter (current_local_index) - std::map num_local_worksets; - std::map max_sideset_length; - std::map max_sides; - std::map current_local_index; - for (size_t i = 0; i < m_sideSets.size(); ++i) { - for (const auto& ss_it : m_sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - // Initialize values if this is the first time seeing a sideset key - if (num_local_worksets.find(ss_key) == num_local_worksets.end()) - num_local_worksets[ss_key] = 0; - if (max_sideset_length.find(ss_key) == max_sideset_length.end()) - max_sideset_length[ss_key] = 0; - if (max_sides.find(ss_key) == max_sides.end()) - max_sides[ss_key] = 0; - if (current_local_index.find(ss_key) == current_local_index.end()) - current_local_index[ss_key] = 0; - - // Update extents for given workset/sideset - num_local_worksets[ss_key]++; - max_sideset_length[ss_key] = std::max(max_sideset_length[ss_key], (int) ss_val.size()); - for (size_t j = 0; j < ss_val.size(); ++j) - max_sides[ss_key] = std::max(max_sides[ss_key], (int) ss_val[j].side_pos); - } - } - - // 2) Construct GlobalSideSetList (map of GlobalSideSetInfo) - for (const auto& ss_it : num_local_worksets) { - std::string ss_key = ss_it.first; - - max_sides[ss_key]++; // max sides is the largest local ID + 1 and needs to be incremented once for each key here - - globalSideSetViews[ss_key].num_local_worksets = num_local_worksets[ss_key]; - globalSideSetViews[ss_key].max_sideset_length = max_sideset_length[ss_key]; - globalSideSetViews[ss_key].side_GID = Kokkos::DualView("side_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].elem_GID = Kokkos::DualView("elem_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].ws_elem_idx = Kokkos::DualView("ws_elem_idx", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].elem_ebIndex = Kokkos::DualView("elem_ebIndex", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].side_pos = Kokkos::DualView("side_pos", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].max_sides = max_sides[ss_key]; - globalSideSetViews[ss_key].numCellsOnSide = Kokkos::DualView("numCellsOnSide", num_local_worksets[ss_key], max_sides[ss_key]); - globalSideSetViews[ss_key].cellsOnSide = Kokkos::DualView("cellsOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].sideSetIdxOnSide = Kokkos::DualView("sideSetIdxOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); - } - - // 3) Populate global views - for (size_t i = 0; i < m_sideSets.size(); ++i) { - for (const auto& ss_it : m_sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - int current_index = current_local_index[ss_key]; - int numSides = max_sides[ss_key]; - - int max_cells_on_side = 0; - std::vector numCellsOnSide(numSides); - std::vector> cellsOnSide(numSides); - std::vector> sideSetIdxOnSide(numSides); - for (size_t j = 0; j < ss_val.size(); ++j) { - int cell = ss_val[j].ws_elem_idx; - int side = ss_val[j].side_pos; - - cellsOnSide[side].push_back(cell); - sideSetIdxOnSide[side].push_back(j); - } - for (int side = 0; side < numSides; ++side) { - numCellsOnSide[side] = cellsOnSide[side].size(); - max_cells_on_side = std::max(max_cells_on_side, numCellsOnSide[side]); - } - - for (int side = 0; side < numSides; ++side) { - globalSideSetViews[ss_key].numCellsOnSide.h_view(current_index, side) = numCellsOnSide[side]; - for (int j = 0; j < numCellsOnSide[side]; ++j) { - globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = cellsOnSide[side][j]; - globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = sideSetIdxOnSide[side][j]; - } - for (int j = numCellsOnSide[side]; j < max_sideset_length[ss_key]; ++j) { - globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = -1; - globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = -1; - } - } - - for (size_t j = 0; j < ss_val.size(); ++j) { - globalSideSetViews[ss_key].side_GID.h_view(current_index, j) = ss_val[j].side_GID; - globalSideSetViews[ss_key].elem_GID.h_view(current_index, j) = ss_val[j].elem_GID; - globalSideSetViews[ss_key].ws_elem_idx.h_view(current_index, j) = ss_val[j].ws_elem_idx; - globalSideSetViews[ss_key].elem_ebIndex.h_view(current_index, j) = ss_val[j].elem_ebIndex; - globalSideSetViews[ss_key].side_pos.h_view(current_index, j) = ss_val[j].side_pos; - } - - globalSideSetViews[ss_key].side_GID.modify_host(); - globalSideSetViews[ss_key].elem_GID.modify_host(); - globalSideSetViews[ss_key].ws_elem_idx.modify_host(); - globalSideSetViews[ss_key].elem_ebIndex.modify_host(); - globalSideSetViews[ss_key].side_pos.modify_host(); - globalSideSetViews[ss_key].numCellsOnSide.modify_host(); - globalSideSetViews[ss_key].cellsOnSide.modify_host(); - globalSideSetViews[ss_key].sideSetIdxOnSide.modify_host(); - - globalSideSetViews[ss_key].side_GID.sync_device(); - globalSideSetViews[ss_key].elem_GID.sync_device(); - globalSideSetViews[ss_key].ws_elem_idx.sync_device(); - globalSideSetViews[ss_key].elem_ebIndex.sync_device(); - globalSideSetViews[ss_key].side_pos.sync_device(); - globalSideSetViews[ss_key].numCellsOnSide.sync_device(); - globalSideSetViews[ss_key].cellsOnSide.sync_device(); - globalSideSetViews[ss_key].sideSetIdxOnSide.sync_device(); - - current_local_index[ss_key]++; - } - } - - // 4) Reset current_local_index - std::map::iterator counter_it = current_local_index.begin(); - while (counter_it != current_local_index.end()) { - std::string counter_key = counter_it->first; - current_local_index[counter_key] = 0; - counter_it++; - } - - // 5) Populate map of LocalSideSetInfos - for (size_t i = 0; i < m_sideSets.size(); ++i) { - LocalSideSetInfoList& lssList = sideSetViews[i]; - - for (const auto& ss_it : m_sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - int current_index = current_local_index[ss_key]; - std::pair range(0, ss_val.size()); - - lssList[ss_key].size = ss_val.size(); - lssList[ss_key].side_GID = Kokkos::subview(globalSideSetViews[ss_key].side_GID, current_index, range ); - lssList[ss_key].elem_GID = Kokkos::subview(globalSideSetViews[ss_key].elem_GID, current_index, range ); - lssList[ss_key].ws_elem_idx = Kokkos::subview(globalSideSetViews[ss_key].ws_elem_idx, current_index, range ); - lssList[ss_key].elem_ebIndex = Kokkos::subview(globalSideSetViews[ss_key].elem_ebIndex, current_index, range ); - lssList[ss_key].side_pos = Kokkos::subview(globalSideSetViews[ss_key].side_pos, current_index, range ); - lssList[ss_key].numSides = globalSideSetViews[ss_key].max_sides; - lssList[ss_key].numCellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].numCellsOnSide, current_index, Kokkos::ALL() ); - lssList[ss_key].cellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].cellsOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); - lssList[ss_key].sideSetIdxOnSide = Kokkos::subview(globalSideSetViews[ss_key].sideSetIdxOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); - - current_local_index[ss_key]++; - } - } - - // 6) Determine size of global DOFView structure and allocate - std::map total_sideset_idx; - std::map sideset_idx_offset; - unsigned int maxSideNodes = 0; - const auto& cell_layers_data = m_extruded_mesh->cell_layers_lid(); - if (!cell_layers_data.is_null()) { - const Teuchos::RCP cell_topo = Teuchos::rcp(new CellTopologyData(m_extruded_mesh->meshSpecs[0]->ctd)); - const int numLayers = cell_layers_data->numLayers; - const int numComps = getDOFManager()->getNumFields(); - - // Determine maximum number of side nodes - for (unsigned int elem_side = 0; elem_side < cell_topo->side_count; ++elem_side) { - const CellTopologyData_Subcell& side = cell_topo->side[elem_side]; - const unsigned int numSideNodes = side.topology->node_count; - maxSideNodes = std::max(maxSideNodes, numSideNodes); - } - - // Determine total number of sideset indices per each sideset name - for (auto& ssList : m_sideSets) { - for (auto& ss_it : ssList) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - if (sideset_idx_offset.find(ss_key) == sideset_idx_offset.end()) - sideset_idx_offset[ss_key] = 0; - if (total_sideset_idx.find(ss_key) == total_sideset_idx.end()) - total_sideset_idx[ss_key] = 0; - - total_sideset_idx[ss_key] += ss_val.size(); - } - } - - // Allocate total localDOFView for each sideset name - for (auto& ss_it : num_local_worksets) { - std::string ss_key = ss_it.first; - allLocalDOFViews[ss_key] = Kokkos::DualView(ss_key + " localDOFView", total_sideset_idx[ss_key], maxSideNodes, numLayers+1, numComps); - } - } - - // Get topo data - auto ctd = m_extruded_mesh->meshSpecs[0]->ctd; - - // Ensure we have ONE cell per layer. - const auto topo_hexa = shards::getCellTopologyData>(); - const auto topo_wedge = shards::getCellTopologyData>(); - TEUCHOS_TEST_FOR_EXCEPTION ( - ctd.name!=topo_hexa->name && - ctd.name!=topo_wedge->name, std::runtime_error, - "Extruded meshes only allowed if there is one element per layer (hexa or wedges).\n" - " - current topology name: " << ctd.name << "\n"); - - const auto& sol_dof_mgr = getDOFManager(); - const auto& elem_dof_lids = sol_dof_mgr->elem_dof_lids().host(); - - // Build a LayeredMeshNumbering for cells, so we can get the LIDs of elems over the column - const auto numLayers = cell_layers_data->numLayers; - const int top = cell_layers_data->top_side_pos; - const int bot = cell_layers_data->bot_side_pos; - - // 7) Populate localDOFViews for GatherVerticallyContractedSolution - for (int ws=0; ws>& wsldofViews = wsLocalDOFViews[ws]; - - const auto& elem_lids = getElementLIDs_host(ws); - - // Loop over the sides that form the boundary condition - // const Teuchos::ArrayRCP >& wsElNodeID_i = wsElNodeID[i]; - for (auto& ss_it : m_sideSets[ws]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - Kokkos::DualView& globalDOFView = allLocalDOFViews[ss_key]; - - for (unsigned int sideSet_idx = 0; sideSet_idx < ss_val.size(); ++sideSet_idx) { - const auto& side = ss_val[sideSet_idx]; - - // Get the data that corresponds to the side - const int ws_elem_idx = side.ws_elem_idx; - const int side_pos = side.side_pos; - - // Check if this sideset is the top or bot of the mesh. If not, the data structure - // for coupling vertical dofs is not needed. - if (side_pos!=top && side_pos!=bot) - break; - - const int elem_LID = elem_lids(ws_elem_idx); - const int basal_elem_LID = cell_layers_data->getColumnId(elem_LID); - - for (int eq=0; eqgetGIDFieldOffsetsSide(eq,top,side_pos); - const auto& sol_bot_offsets = sol_dof_mgr->getGIDFieldOffsetsSide(eq,bot,side_pos); - const int numSideNodes = sol_top_offsets.size(); - - for (int j=0; jgetId(basal_elem_LID,il); - globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il, eq) = - elem_dof_lids(layer_elem_LID,sol_bot_offsets[j]); - } - - // Add top side in last layer - const int il = numLayers-1; - const LO layer_elem_LID = cell_layers_data->getId(basal_elem_LID,il); - globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il+1, eq) = - elem_dof_lids(layer_elem_LID,sol_top_offsets[j]); - } - } - } - - globalDOFView.modify_host(); - globalDOFView.sync_device(); - - // Set workset-local sub-view - std::pair range(sideset_idx_offset[ss_key], sideset_idx_offset[ss_key]+ss_val.size()); - wsldofViews[ss_key] = Kokkos::subview(globalDOFView, range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - - sideset_idx_offset[ss_key] += ss_val.size(); - } - } + buildSideSetsViews(); } void @@ -1160,8 +884,8 @@ buildCellSideNodeNumerationMaps() // ONLY for basalside and upperside, since that's where we are likely to load data from mesh for (int ws=0; wsgetColumnId(s.elem_GID); @@ -1200,17 +924,23 @@ Teuchos::RCP ExtrudedDiscretization:: create_dof_mgr (const std::string& part_name, const std::string& field_name, - const FE_Type /* fe_type */, - const int /* order */, - const int dof_dim) const + const FE_Type fe_type, + const int order, + const int dof_dim) { + auto& dof_mgr = get_dof_mgr(part_name,fe_type,order,dof_dim); + if (Teuchos::nonnull(dof_mgr)) { + // Not the first time we build a DOFManager for a field with these specs + return dof_mgr; + } + const auto& ebn = m_extruded_mesh->meshSpecs()[0]->ebName;; std::vector elem_blocks = {ebn}; // Create conn and dof managers auto conn_mgr_h = m_basal_disc->getDOFManager(field_name)->getAlbanyConnManager(); auto conn_mgr = Teuchos::rcp(new ExtrudedConnManager(conn_mgr_h,m_extruded_mesh)); - auto dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,m_comm,part_name)); + dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,m_comm,part_name)); const auto& topo = conn_mgr->get_topology(); Teuchos::RCP fp; @@ -1223,9 +953,9 @@ create_dof_mgr (const std::string& part_name, fp = Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis)); } // NOTE: we add $dof_dim copies of the field pattern to the dof mgr, - // and call the fields ${field_name}_n, n=0,..,$dof_dim-1 + // and call the fields cmp_n, n=0,..,$dof_dim-1 for (int i=0; iaddField(field_name + "_" + std::to_string(i),fp); + dof_mgr->addField("cmp " + std::to_string(i),fp); } dof_mgr->build(); diff --git a/src/disc/Albany_ExtrudedDiscretization.hpp b/src/disc/extruded/Albany_ExtrudedDiscretization.hpp similarity index 82% rename from src/disc/Albany_ExtrudedDiscretization.hpp rename to src/disc/extruded/Albany_ExtrudedDiscretization.hpp index 919bc1c992..f2b6225655 100644 --- a/src/disc/Albany_ExtrudedDiscretization.hpp +++ b/src/disc/extruded/Albany_ExtrudedDiscretization.hpp @@ -41,25 +41,6 @@ class ExtrudedDiscretization : public AbstractDiscretization const NodeSetGIDsList& getNodeSetGIDs() const override { return m_nodeSetGIDs; } const NodeSetCoordList& getNodeSetCoords() const override { return m_nodeSetCoords; } - //! Get Side set lists (typedef in Albany_AbstractDiscretization.hpp) - const SideSetList& getSideSets(const int workset) const override - { - return m_sideSets[workset]; - } - - //! Get Side set lists (typedef in Albany_AbstractDiscretization.hpp) - const LocalSideSetInfoList& getSideSetViews(const int workset) const override - { - return sideSetViews.at(workset); - } - - //! Get local DOF views for GatherVerticallyContractedSolution - const std::map>& - getLocalDOFViews(const int workset) const override - { - return wsLocalDOFViews.at(workset); - } - //! Get connectivity map from elementGID to workset WsLIDList& getElemGIDws() override { return m_elemGIDws; } const WsLIDList& getElemGIDws() const override { return m_elemGIDws; } @@ -72,17 +53,6 @@ class ExtrudedDiscretization : public AbstractDiscretization Teuchos::RCP getMeshStruct() const override { return m_extruded_mesh; } - const std::map>& getSideToSideSetCellMap() const override - { - throw NotYetImplemented("ExtrudedDiscretization::getSideToSideSetCellMap"); - } - - const std::map>>& - getSideNodeNumerationMap() const override - { - throw NotYetImplemented("ExtrudedDiscretization::getSideNodeNumerationMap"); - } - //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const override { return false; } @@ -99,9 +69,6 @@ class ExtrudedDiscretization : public AbstractDiscretization //! Get number of spatial dimensions int getNumDim() const override { return m_extruded_mesh->meshSpecs[0]->numDim; } - //! Get number of total DOFs per node - int getNumEq() const override { return m_neq; } - // --- Get/set solution/residual/field vectors to/from mesh --- // Teuchos::RCP getSolutionField (const bool overlapped = false) const override; @@ -147,7 +114,7 @@ class ExtrudedDiscretization : public AbstractDiscretization void setFieldData() override; - protected: +protected: void getSolutionField(Thyra_Vector& result, bool overlapped) const; @@ -193,7 +160,7 @@ class ExtrudedDiscretization : public AbstractDiscretization const std::string& field_name, const FE_Type fe_type, const int order, - const int dof_dim) const; + const int dof_dim); // ==================== Members =================== // @@ -202,9 +169,6 @@ class ExtrudedDiscretization : public AbstractDiscretization Teuchos::RCP m_basal_disc; - //! Number of equations (and unknowns) per node - const int m_neq; - //! Equations that are defined only on some side sets of the mesh std::map> m_sideSetEquations; @@ -213,15 +177,6 @@ class ExtrudedDiscretization : public AbstractDiscretization NodeSetGIDsList m_nodeSetGIDs; NodeSetCoordList m_nodeSetCoords; - //! side sets stored as std::map(string ID, SideArray classes) per workset - std::vector m_sideSets; - GlobalSideSetList globalSideSetViews; - std::map sideSetViews; - - //! GatherVerticallyContractedSolution connectivity - std::map> allLocalDOFViews; - std::map>> wsLocalDOFViews; - // Coordinates spliced together, as [x0,y0,z0,x1,y1,z1,...] Teuchos::ArrayRCP m_nodes_coordinates; @@ -238,8 +193,6 @@ class ExtrudedDiscretization : public AbstractDiscretization Teuchos::RCP m_disc_params; // Sideset discretizations - strmap_t> sideToSideSetCellMap; - strmap_t>> sideNodeNumerationMap; strmap_t> projectors; strmap_t> ov_projectors; }; diff --git a/src/disc/Albany_ExtrudedMesh.cpp b/src/disc/extruded/Albany_ExtrudedMesh.cpp similarity index 87% rename from src/disc/Albany_ExtrudedMesh.cpp rename to src/disc/extruded/Albany_ExtrudedMesh.cpp index 444ba6cfa6..13ceee706c 100644 --- a/src/disc/Albany_ExtrudedMesh.cpp +++ b/src/disc/extruded/Albany_ExtrudedMesh.cpp @@ -93,24 +93,35 @@ get_basal_part_name (const std::string& extruded_part_name) const } void ExtrudedMesh:: -setFieldData (const Teuchos::RCP& comm) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) { - // Ensure field data is set on basal mesh - m_basal_mesh->setFieldData(comm); - - // Register surface height and mesh thickness in the 2d mesh + // Ensure surface_height and thickness are valid states in the basal mesh std::string thickness_name = m_params->get("Thickness Field Name","thickness"); std::string surface_height_name = m_params->get("Surface Height Field Name","surface_height"); - StateInfoStruct mesh_sis; + auto& basal_sis = side_set_sis["basal_side"]; + if (basal_sis.is_null()) { + basal_sis = Teuchos::rcp(new StateInfoStruct()); + } auto NDTEN = StateStruct::MeshFieldEntity::NodalDataToElemNode; StateStruct::FieldDims dims = { static_cast(m_basal_mesh->meshSpecs[0]->worksetSize), m_basal_mesh->meshSpecs[0]->ctd.node_count }; - mesh_sis.emplace_back(Teuchos::rcp(new StateStruct(surface_height_name,NDTEN,dims,""))); - mesh_sis.emplace_back(Teuchos::rcp(new StateStruct(thickness_name,NDTEN,dims,""))); - m_basal_mesh->get_field_accessor()->addStateStructs(mesh_sis); + + if (basal_sis->find(surface_height_name).is_null()) { + basal_sis->emplace_back(Teuchos::rcp(new StateStruct(surface_height_name,NDTEN,dims,""))); + } + if (basal_sis->find(thickness_name).is_null()) { + basal_sis->emplace_back(Teuchos::rcp(new StateStruct(thickness_name,NDTEN,dims,""))); + } + + // Ensure field data is set on basal mesh + m_basal_mesh->setFieldData(comm,basal_sis,{}); + + // TODO: correctly register volume states } void ExtrudedMesh:: diff --git a/src/disc/Albany_ExtrudedMesh.hpp b/src/disc/extruded/Albany_ExtrudedMesh.hpp similarity index 95% rename from src/disc/Albany_ExtrudedMesh.hpp rename to src/disc/extruded/Albany_ExtrudedMesh.hpp index c1c47de307..04aa2f78dc 100644 --- a/src/disc/Albany_ExtrudedMesh.hpp +++ b/src/disc/extruded/Albany_ExtrudedMesh.hpp @@ -63,7 +63,9 @@ class ExtrudedMesh : public AbstractMeshStruct { return m_basal_mesh->get_field_accessor(); } - void setFieldData (const Teuchos::RCP& comm) override; + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) override; void setBulkData(const Teuchos::RCP& comm) override; diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index a47079660b..80b5990187 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -67,9 +67,9 @@ OmegahDiscretization (const Teuchos::RCP& discParams, , m_mesh_struct(mesh) , m_comm (comm) , m_side_set_equations(sideSetEquations) - , m_neq (neq) { m_num_time_deriv = m_disc_params->get("Number Of Time Derivatives",0); + m_neq = neq; // TODO: get solution names from param list m_sol_names.resize(m_num_time_deriv+1,solution_dof_name()); @@ -86,9 +86,12 @@ updateMesh () { printf ("TODO: change name to the method?\n"); + // Make sure we don't reuse old dof mgrs (if adapting) + m_key_to_dof_mgr.clear(); + // Create DOF managers - auto sol_dof_mgr = create_dof_mgr(solution_dof_name(),"",FE_Type::HGRAD,1,m_neq); - auto node_dof_mgr = create_dof_mgr(nodes_dof_name(),"",FE_Type::HGRAD,1,1); + auto sol_dof_mgr = create_dof_mgr("",FE_Type::HGRAD,1,m_neq); + auto node_dof_mgr = create_dof_mgr("",FE_Type::HGRAD,1,1); m_dof_managers[solution_dof_name()][""] = sol_dof_mgr; m_dof_managers[nodes_dof_name()][""] = node_dof_mgr; @@ -122,6 +125,8 @@ updateMesh () m_wsPhysIndex[i] = ms->ebNameToIndex[m_wsEBNames[i]]; } + m_mesh_struct->get_field_accessor()->createStateArrays(m_workset_sizes); + m_ws_elem_coords.resize(num_ws); auto coords_h = m_mesh_struct->coords_host(); auto node_gids = hostRead(mesh.globals(0)); @@ -157,10 +162,10 @@ updateMesh () elms_in_prior_worksets += m_workset_sizes[ws]; } - m_side_sets.resize(num_ws); + m_sideSets.resize(num_ws); for (int ws=0; wsname + "\n" " - input dims: (" + util::join(st->dim,",") + ")\n"); } - auto dof_mgr = create_dof_mgr (st->name,st->meshPart,FE_Type::HGRAD,1,numComps); + auto dof_mgr = create_dof_mgr (st->meshPart,FE_Type::HGRAD,1,numComps); m_dof_managers[st->name][st->meshPart] = dof_mgr; if (m_node_dof_managers.find(st->meshPart)==m_node_dof_managers.end()) { - auto node_dof_mgr = create_dof_mgr (nodes_dof_name(),st->meshPart,FE_Type::HGRAD,1,1); + auto node_dof_mgr = create_dof_mgr (st->meshPart,FE_Type::HGRAD,1,1); m_node_dof_managers[st->meshPart] = node_dof_mgr; } } @@ -328,17 +333,22 @@ setField (const Thyra_Vector& field_vector, Teuchos::RCP OmegahDiscretization:: -create_dof_mgr (const std::string& field_name, - const std::string& part_name, +create_dof_mgr (const std::string& part_name, const FE_Type fe_type, const int order, - const int dof_dim) const + const int dof_dim) { + auto& dof_mgr = get_dof_mgr(part_name,fe_type,order,dof_dim); + if (Teuchos::nonnull(dof_mgr)) { + // Not the first time we build a DOFManager for a field with these specs + return dof_mgr; + } + const auto& mesh_specs = m_mesh_struct->meshSpecs[0]; // Create conn and dof managers auto conn_mgr = Teuchos::rcp(new OmegahConnManager(m_mesh_struct)); - auto dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,m_comm,part_name)); + dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,m_comm,part_name)); shards::CellTopology topo (&mesh_specs->ctd); Teuchos::RCP fp; @@ -353,7 +363,7 @@ create_dof_mgr (const std::string& field_name, // NOTE: we add $dof_dim copies of the field pattern to the dof mgr, // and call the fields ${field_name}_n, n=0,..,$dof_dim-1 for (int i=0; iaddField(field_name + "_" + std::to_string(i),fp); + dof_mgr->addField("cmp_" + std::to_string(i),fp); } dof_mgr->build(); diff --git a/src/disc/omegah/Albany_OmegahDiscretization.hpp b/src/disc/omegah/Albany_OmegahDiscretization.hpp index c3cfb117e6..059f94f6fa 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.hpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.hpp @@ -46,24 +46,6 @@ class OmegahDiscretization : public AbstractDiscretization return m_node_set_coords; } - //! Get Side set lists - const SideSetList& - getSideSets(const int ws) const override { - return m_side_sets[ws]; - } - - //! Get Side set view lists - const LocalSideSetInfoList& - getSideSetViews(const int ws) const override { - return m_side_set_views.at(ws); - } - - //! Get local DOF views for GatherVerticallyContractedSolution - const strmap_t>& - getLocalDOFViews(const int workset) const override { - return m_ws_local_dof_views.at(workset); - } - //! Get coordinates (overlap map). const Teuchos::ArrayRCP& getCoordinates() const override { return m_nodes_coordinates; } @@ -73,18 +55,6 @@ class OmegahDiscretization : public AbstractDiscretization TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::printCoords"); } - //! Get the map side_id->side_set_elem_id - const std::map>& - getSideToSideSetCellMap() const override { - TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::getSideToSideSetCellMap"); - } - - //! Get the map side_node_id->side_set_cell_node_id - const std::map>>& - getSideNodeNumerationMap() const override { - TEUCHOS_TEST_FOR_EXCEPTION(true,NotYetImplemented,"OmegahDiscretization::getSideNodeNumerationMap"); - } - //! Get MeshStruct Teuchos::RCP getMeshStruct() const override { @@ -119,12 +89,6 @@ class OmegahDiscretization : public AbstractDiscretization return m_mesh_struct->meshSpecs[0]->numDim; } - //! Get number of total DOFs per node - int - getNumEq() const override { - return m_neq; - } - // --- Get/set solution/residual/field vectors to/from mesh --- // Teuchos::RCP getSolutionField(bool /* overlapped */ = false) const override { @@ -182,10 +146,9 @@ class OmegahDiscretization : public AbstractDiscretization Teuchos::RCP create_dof_mgr (const std::string& part_name, - const std::string& field_name, const FE_Type fe_type, const int order, - const int dof_dim) const; + const int dof_dim); void computeNodeSets (); void computeGraphs (); @@ -204,10 +167,6 @@ class OmegahDiscretization : public AbstractDiscretization // TODO: move stuff below in base class? Teuchos::RCP m_comm; - std::map m_side_set_views; - std::vector m_side_sets; - std::map>> m_ws_local_dof_views; - NodeSetList m_node_sets; NodeSetCoordList m_node_set_coords; @@ -217,10 +176,6 @@ class OmegahDiscretization : public AbstractDiscretization //! Equations that are defined only on some side sets of the mesh std::map> m_side_set_equations; - // Number of equations (and unknowns) per node - // TODO: this should soon be removed, in favor of more granular description of each dof/unknown - const int m_neq; - // TODO: I don't think this belongs in the disc class. We do store it for STK meshes, // mainly b/c we need to know how many components the soln vector has. // But I think we can get rid of it. In principle, we should handle time derivatives diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index 88782c68b2..36ec65df00 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -36,9 +36,17 @@ OmegahGenericMesh (const Teuchos::RCP& params) } void OmegahGenericMesh:: -setFieldData (const Teuchos::RCP& comm) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) { loadRequiredInputFields(comm); + get_field_accessor()->addStateStructs(sis); + + // We don't yet support side discretizations for omegah meshes + (void) side_set_sis; + + m_field_data_set = true; } void OmegahGenericMesh:: diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.hpp b/src/disc/omegah/Albany_OmegahGenericMesh.hpp index 529306547e..32073d54dc 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.hpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.hpp @@ -20,7 +20,9 @@ class OmegahGenericMesh : public AbstractMeshStruct { // ------------- Override from base class ------------- // std::string meshLibName () const override { return "Omega_h"; } - void setFieldData (const Teuchos::RCP& comm) override; + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) override; void setBulkData (const Teuchos::RCP& comm) override; LO get_num_local_nodes () const override; diff --git a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp index af6f4623db..865ef9d0ce 100644 --- a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp @@ -194,6 +194,15 @@ ExtrudedSTKMeshStruct(const Teuchos::RCP& params, this->initializeSideSetMeshSpecs(comm); } +void ExtrudedSTKMeshStruct:: +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) +{ + basalMeshStruct->setFieldData(comm,side_set_sis["basalside"],{}); + GenericSTKMeshStruct::setFieldData(comm,sis,side_set_sis); +} + void ExtrudedSTKMeshStruct:: setBulkData (const Teuchos::RCP& comm) { diff --git a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp index 0d5466093a..d0ddd8e0cf 100644 --- a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp @@ -26,6 +26,9 @@ class ExtrudedSTKMeshStruct : public GenericSTKMeshStruct ~ExtrudedSTKMeshStruct() = default; + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis); void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp index 8908b00593..d08542e671 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp @@ -75,7 +75,9 @@ GenericSTKMeshStruct::GenericSTKMeshStruct( } void GenericSTKMeshStruct:: -setFieldData (const Teuchos::RCP& comm) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) { TEUCHOS_TEST_FOR_EXCEPTION(!metaData->is_initialized(), std::logic_error, "[GenericSTKMeshStruct::SetupFieldData] metaData->initialize(numDim) not yet called" << std::endl); @@ -156,6 +158,14 @@ setFieldData (const Teuchos::RCP& comm) writeCoordsToMMFile = params->get("Write Coordinates to MatrixMarket", false); transferSolutionToCoords = params->get("Transfer Solution to Coordinates", false); + + fieldContainer->addStateStructs(sis); + for (auto& [ss_name, ss_mesh] : sideSetMeshStructs) { + if (not ss_mesh->isFieldDataSet()) { + ss_mesh->setFieldData(comm,side_set_sis[ss_name],{}); + } + } + m_field_data_set = true; } void GenericSTKMeshStruct::setAllPartsIO() diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp index 22ebac143f..4152a2853f 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp @@ -39,7 +39,9 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct int getNumParams() const {return num_params; } - void setFieldData (const Teuchos::RCP& comm) override; + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) override; void printParts(stk::mesh::MetaData *metaData); diff --git a/src/disc/stk/Albany_IossSTKMeshStruct.cpp b/src/disc/stk/Albany_IossSTKMeshStruct.cpp index 0cdc7919e4..071e1d4ff2 100644 --- a/src/disc/stk/Albany_IossSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_IossSTKMeshStruct.cpp @@ -276,9 +276,11 @@ IossSTKMeshStruct::~IossSTKMeshStruct() } void IossSTKMeshStruct:: -setFieldData (const Teuchos::RCP& comm) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) { - GenericSTKMeshStruct::setFieldData(comm); + GenericSTKMeshStruct::setFieldData(comm,sis,side_set_sis); if(mesh_data->is_bulk_data_null()) mesh_data->set_bulk_data(*bulkData); diff --git a/src/disc/stk/Albany_IossSTKMeshStruct.hpp b/src/disc/stk/Albany_IossSTKMeshStruct.hpp index a4a68aa6d6..62311e3018 100644 --- a/src/disc/stk/Albany_IossSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_IossSTKMeshStruct.hpp @@ -28,7 +28,9 @@ namespace Albany { ~IossSTKMeshStruct(); - void setFieldData (const Teuchos::RCP& comm); + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis); void setBulkData (const Teuchos::RCP& comm); diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index fab2247918..e87aa6cdc9 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -67,16 +67,17 @@ STKDiscretization::STKDiscretization( metaData(stkMeshStruct_->metaData), bulkData(stkMeshStruct_->bulkData), comm(comm_), - neq(neq_), sideSetEquations(sideSetEquations_), rigidBodyModes(rigidBodyModes_), stkMeshStruct(stkMeshStruct_), discParams(discParams_) { + m_neq = neq_; + if (stkMeshStruct->sideSetMeshStructs.size() > 0) { for (auto it : stkMeshStruct->sideSetMeshStructs) { auto stk_mesh = Teuchos::rcp_dynamic_cast(it.second,true); - auto side_disc = Teuchos::rcp(new STKDiscretization(discParams, neq, stk_mesh, comm)); + auto side_disc = Teuchos::rcp(new STKDiscretization(discParams, m_neq, stk_mesh, comm)); sideSetDiscretizations.insert(std::make_pair(it.first, side_disc)); sideSetDiscretizationsSTK.insert(std::make_pair(it.first, side_disc)); } @@ -907,9 +908,20 @@ void STKDiscretization::computeVectorSpaces() TEUCHOS_TEST_FOR_EXCEPTION (stkMeshStruct->ebNames_.size() < 1,std::logic_error, "Error! Albany problems must have at least 1 element block. Your mesh has 0 element blocks.\n"); - strmap_t> name_to_partAndDim; - name_to_partAndDim[solution_dof_name()] = std::make_pair("",neq); - name_to_partAndDim[nodes_dof_name()] = std::make_pair("",1); + // Make sure we don't reuse old dof mgrs (if adapting) + m_key_to_dof_mgr.clear(); + + Teuchos::RCP dof_mgr; + + // Solution dof mgr + dof_mgr = create_dof_mgr("",FE_Type::HGRAD,1,m_neq); + m_dof_managers[solution_dof_name()][""] = dof_mgr; + + // Nodes dof mgr + dof_mgr = create_dof_mgr("",FE_Type::HGRAD,1,1); + m_dof_managers[nodes_dof_name()][""] = dof_mgr; + m_node_dof_managers[""] = dof_mgr; + for (const auto& sis : stkMeshStruct->getFieldContainer()->getNodalParameterSIS()) { const auto& dims = sis->dim; int dof_dim = -1; @@ -922,25 +934,13 @@ void STKDiscretization::computeVectorSpaces() "Error! Unsupported layout for nodal parameter '" + sis->name + ".\n"); } - name_to_partAndDim[sis->name] = std::make_pair(sis->meshPart,dof_dim); - } - - for (const auto& it : name_to_partAndDim) { - const auto& field_name = it.first; - const auto& part_name = it.second.first; - const auto& dof_dim = it.second.second; - // NOTE: for now we hard code P1. In the future, we must be able to // store this info somewhere and retrieve it here. - auto dof_mgr = create_dof_mgr(part_name,field_name,FE_Type::HGRAD,1,dof_dim); - m_dof_managers[field_name][part_name] = dof_mgr; - m_node_dof_managers[part_name] = Teuchos::null; - } + dof_mgr = create_dof_mgr(sis->meshPart,FE_Type::HGRAD,1,dof_dim); + m_dof_managers[sis->name][sis->meshPart] = dof_mgr; - // For each part, also make a Node dof manager - for (auto& it : m_node_dof_managers) { - const auto& part_name = it.first; - it.second = create_dof_mgr(part_name, nodes_dof_name(), FE_Type::HGRAD,1,1); + dof_mgr = create_dof_mgr(sis->meshPart,FE_Type::HGRAD,1,1); + m_node_dof_managers[sis->meshPart] = dof_mgr; } const int meshDim = stkMeshStruct->numDim; @@ -959,7 +959,7 @@ STKDiscretization::computeGraphs() // as well as what eqn are on each sideset std::vector volumeEqns; std::map> ss_to_eqns; - for (int k=0; k < neq; ++k) { + for (int k=0; k < m_neq; ++k) { if (sideSetEquations.find(k) == sideSetEquations.end()) { volumeEqns.push_back(k); } @@ -971,7 +971,7 @@ STKDiscretization::computeGraphs() const int num_elems = sol_dof_mgr->cell_indexer()->getNumLocalElements(); // Handle the simple case, and return immediately - if (numVolumeEqns==neq) { + if (numVolumeEqns==m_neq) { // This is the easy case: couple everything with everything for (int icell=0; icellgetElementGIDs(icell); @@ -1069,7 +1069,7 @@ STKDiscretization::computeGraphs() for (const auto& ss_name : it.second) { for (int ws=0; wsnumLayers; const LO basal_elem_LID = cell_layers_data_lid->getColumnId(elem_LID); - for (int eq=0; eqgetGIDFieldOffsets(eq); const int num_col_gids = eq_offsets.size(); cols.resize(num_col_gids); @@ -1113,7 +1113,7 @@ STKDiscretization::computeGraphs() // contained in the other, but that starts to get too involved. Given that it's not // a common scenario (need 2+ ss eqn defined on 2 different sidesets), and that we // might have to redo this when we assemble by blocks, we just don't bother. - for (int col_eq=0; col_eqgetGIDFieldOffsetsSide(col_eq,side_pos); const int num_col_gids = col_eq_offsets.size(); cols.resize(num_col_gids); @@ -1288,23 +1288,23 @@ STKDiscretization::computeSideSets() { TEUCHOS_FUNC_TIME_MONITOR("STKDiscretization: computeSideSets"); // Clean up existing sideset structure if remeshing - for (size_t i = 0; i < sideSets.size(); ++i) { - sideSets[i].clear(); // empty the ith map + for (auto& ss : m_sideSets) { + ss.clear(); // empty the ith map } int numBuckets = m_wsEBNames.size(); - sideSets.resize(numBuckets); // Need a sideset list per workset + m_sideSets.resize(numBuckets); // Need a sideset list per workset - for (const auto& ss : stkMeshStruct->ssPartVec) { + for (const auto& [ss_name,ss_part] : stkMeshStruct->ssPartVec) { // Make sure the sideset exist even if no sides are owned - for (int i=0; ilocally_owned_part()); std::vector sides; @@ -1313,7 +1313,7 @@ STKDiscretization::computeSideSets() bulkData->buckets(metaData->side_rank()), sides); - *out << "STKDisc: sideset " << ss.first << " has size " << sides.size() + *out << "STKDisc: sideset " << ss_name << " has size " << sides.size() << " on Proc 0." << std::endl; // loop over the sides to see what they are, then fill in the data holder @@ -1328,7 +1328,7 @@ STKDiscretization::computeSideSets() bulkData->num_elements(sidee) != 1, std::logic_error, "STKDisc: cannot figure out side set topology for side set " - << ss.first << std::endl); + << ss_name << std::endl); stk::mesh::Entity elem = bulkData->begin_elements(sidee)[0]; @@ -1356,293 +1356,12 @@ STKDiscretization::computeSideSets() stkMeshStruct->meshSpecs[0]->ebNameToIndex[m_wsEBNames[workset]]; // Get or create the vector of side structs for this side set on this workset - auto& ss_vec = sideSets[workset][ss.first]; + auto& ss_vec = m_sideSets[workset][ss_name]; ss_vec.push_back(sStruct); } } - // ============================================================= - // (Kokkos Refactor) Convert sideSets to sideSetViews - - // 1) Compute view extents (num_local_worksets, max_sideset_length, max_sides) and local workset counter (current_local_index) - std::map num_local_worksets; - std::map max_sideset_length; - std::map max_sides; - std::map current_local_index; - for (size_t i = 0; i < sideSets.size(); ++i) { - for (const auto& ss_it : sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - // Initialize values if this is the first time seeing a sideset key - if (num_local_worksets.find(ss_key) == num_local_worksets.end()) - num_local_worksets[ss_key] = 0; - if (max_sideset_length.find(ss_key) == max_sideset_length.end()) - max_sideset_length[ss_key] = 0; - if (max_sides.find(ss_key) == max_sides.end()) - max_sides[ss_key] = 0; - if (current_local_index.find(ss_key) == current_local_index.end()) - current_local_index[ss_key] = 0; - - // Update extents for given workset/sideset - num_local_worksets[ss_key]++; - max_sideset_length[ss_key] = std::max(max_sideset_length[ss_key], (int) ss_val.size()); - for (size_t j = 0; j < ss_val.size(); ++j) - max_sides[ss_key] = std::max(max_sides[ss_key], (int) ss_val[j].side_pos); - } - } - - // 2) Construct GlobalSideSetList (map of GlobalSideSetInfo) - for (const auto& ss_it : num_local_worksets) { - std::string ss_key = ss_it.first; - - max_sides[ss_key]++; // max sides is the largest local ID + 1 and needs to be incremented once for each key here - - globalSideSetViews[ss_key].num_local_worksets = num_local_worksets[ss_key]; - globalSideSetViews[ss_key].max_sideset_length = max_sideset_length[ss_key]; - globalSideSetViews[ss_key].side_GID = Kokkos::DualView("side_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].elem_GID = Kokkos::DualView("elem_GID", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].ws_elem_idx = Kokkos::DualView("ws_elem_idx", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].elem_ebIndex = Kokkos::DualView("elem_ebIndex", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].side_pos = Kokkos::DualView("side_pos", num_local_worksets[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].max_sides = max_sides[ss_key]; - globalSideSetViews[ss_key].numCellsOnSide = Kokkos::DualView("numCellsOnSide", num_local_worksets[ss_key], max_sides[ss_key]); - globalSideSetViews[ss_key].cellsOnSide = Kokkos::DualView("cellsOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); - globalSideSetViews[ss_key].sideSetIdxOnSide = Kokkos::DualView("sideSetIdxOnSide", num_local_worksets[ss_key], max_sides[ss_key], max_sideset_length[ss_key]); - } - - // 3) Populate global views - for (size_t i = 0; i < sideSets.size(); ++i) { - for (const auto& ss_it : sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - int current_index = current_local_index[ss_key]; - int numSides = max_sides[ss_key]; - - int max_cells_on_side = 0; - std::vector numCellsOnSide(numSides); - std::vector> cellsOnSide(numSides); - std::vector> sideSetIdxOnSide(numSides); - for (size_t j = 0; j < ss_val.size(); ++j) { - int cell = ss_val[j].ws_elem_idx; - int side = ss_val[j].side_pos; - - cellsOnSide[side].push_back(cell); - sideSetIdxOnSide[side].push_back(j); - } - for (int side = 0; side < numSides; ++side) { - numCellsOnSide[side] = cellsOnSide[side].size(); - max_cells_on_side = std::max(max_cells_on_side, numCellsOnSide[side]); - } - - for (int side = 0; side < numSides; ++side) { - globalSideSetViews[ss_key].numCellsOnSide.h_view(current_index, side) = numCellsOnSide[side]; - for (int j = 0; j < numCellsOnSide[side]; ++j) { - globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = cellsOnSide[side][j]; - globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = sideSetIdxOnSide[side][j]; - } - for (int j = numCellsOnSide[side]; j < max_sideset_length[ss_key]; ++j) { - globalSideSetViews[ss_key].cellsOnSide.h_view(current_index, side, j) = -1; - globalSideSetViews[ss_key].sideSetIdxOnSide.h_view(current_index, side, j) = -1; - } - } - - for (size_t j = 0; j < ss_val.size(); ++j) { - globalSideSetViews[ss_key].side_GID.h_view(current_index, j) = ss_val[j].side_GID; - globalSideSetViews[ss_key].elem_GID.h_view(current_index, j) = ss_val[j].elem_GID; - globalSideSetViews[ss_key].ws_elem_idx.h_view(current_index, j) = ss_val[j].ws_elem_idx; - globalSideSetViews[ss_key].elem_ebIndex.h_view(current_index, j) = ss_val[j].elem_ebIndex; - globalSideSetViews[ss_key].side_pos.h_view(current_index, j) = ss_val[j].side_pos; - } - - globalSideSetViews[ss_key].side_GID.modify_host(); - globalSideSetViews[ss_key].elem_GID.modify_host(); - globalSideSetViews[ss_key].ws_elem_idx.modify_host(); - globalSideSetViews[ss_key].elem_ebIndex.modify_host(); - globalSideSetViews[ss_key].side_pos.modify_host(); - globalSideSetViews[ss_key].numCellsOnSide.modify_host(); - globalSideSetViews[ss_key].cellsOnSide.modify_host(); - globalSideSetViews[ss_key].sideSetIdxOnSide.modify_host(); - - globalSideSetViews[ss_key].side_GID.sync_device(); - globalSideSetViews[ss_key].elem_GID.sync_device(); - globalSideSetViews[ss_key].ws_elem_idx.sync_device(); - globalSideSetViews[ss_key].elem_ebIndex.sync_device(); - globalSideSetViews[ss_key].side_pos.sync_device(); - globalSideSetViews[ss_key].numCellsOnSide.sync_device(); - globalSideSetViews[ss_key].cellsOnSide.sync_device(); - globalSideSetViews[ss_key].sideSetIdxOnSide.sync_device(); - - current_local_index[ss_key]++; - } - } - - // 4) Reset current_local_index - std::map::iterator counter_it = current_local_index.begin(); - while (counter_it != current_local_index.end()) { - std::string counter_key = counter_it->first; - current_local_index[counter_key] = 0; - counter_it++; - } - - // 5) Populate map of LocalSideSetInfos - for (size_t i = 0; i < sideSets.size(); ++i) { - LocalSideSetInfoList& lssList = sideSetViews[i]; - - for (const auto& ss_it : sideSets[i]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - int current_index = current_local_index[ss_key]; - std::pair range(0, ss_val.size()); - - lssList[ss_key].size = ss_val.size(); - lssList[ss_key].side_GID = Kokkos::subview(globalSideSetViews[ss_key].side_GID, current_index, range ); - lssList[ss_key].elem_GID = Kokkos::subview(globalSideSetViews[ss_key].elem_GID, current_index, range ); - lssList[ss_key].ws_elem_idx = Kokkos::subview(globalSideSetViews[ss_key].ws_elem_idx, current_index, range ); - lssList[ss_key].elem_ebIndex = Kokkos::subview(globalSideSetViews[ss_key].elem_ebIndex, current_index, range ); - lssList[ss_key].side_pos = Kokkos::subview(globalSideSetViews[ss_key].side_pos, current_index, range ); - lssList[ss_key].numSides = globalSideSetViews[ss_key].max_sides; - lssList[ss_key].numCellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].numCellsOnSide, current_index, Kokkos::ALL() ); - lssList[ss_key].cellsOnSide = Kokkos::subview(globalSideSetViews[ss_key].cellsOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); - lssList[ss_key].sideSetIdxOnSide = Kokkos::subview(globalSideSetViews[ss_key].sideSetIdxOnSide, current_index, Kokkos::ALL(), Kokkos::ALL() ); - - current_local_index[ss_key]++; - } - } - - // 6) Determine size of global DOFView structure and allocate - std::map total_sideset_idx; - std::map sideset_idx_offset; - unsigned int maxSideNodes = 0; - const auto& cell_layers_data = stkMeshStruct->local_cell_layers_data; - if (!cell_layers_data.is_null()) { - const Teuchos::RCP cell_topo = Teuchos::rcp(new CellTopologyData(stkMeshStruct->meshSpecs[0]->ctd)); - const int numLayers = cell_layers_data->numLayers; - const int numComps = getDOFManager()->getNumFields(); - - // Determine maximum number of side nodes - for (unsigned int elem_side = 0; elem_side < cell_topo->side_count; ++elem_side) { - const CellTopologyData_Subcell& side = cell_topo->side[elem_side]; - const unsigned int numSideNodes = side.topology->node_count; - maxSideNodes = std::max(maxSideNodes, numSideNodes); - } - - // Determine total number of sideset indices per each sideset name - for (auto& ssList : sideSets) { - for (auto& ss_it : ssList) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - if (sideset_idx_offset.find(ss_key) == sideset_idx_offset.end()) - sideset_idx_offset[ss_key] = 0; - if (total_sideset_idx.find(ss_key) == total_sideset_idx.end()) - total_sideset_idx[ss_key] = 0; - - total_sideset_idx[ss_key] += ss_val.size(); - } - } - - // Allocate total localDOFView for each sideset name - for (auto& ss_it : num_local_worksets) { - std::string ss_key = ss_it.first; - allLocalDOFViews[ss_key] = Kokkos::DualView(ss_key + " localDOFView", total_sideset_idx[ss_key], maxSideNodes, numLayers+1, numComps); - } - } - - // Not all mesh structs that come through here are extruded mesh structs. - // If the mesh isn't extruded, we won't need to do any of the following work. - if (not cell_layers_data.is_null()) { - // Get topo data - auto ctd = stkMeshStruct->meshSpecs[0]->ctd; - - // Ensure we have ONE cell per layer. - const auto topo_hexa = shards::getCellTopologyData>(); - const auto topo_wedge = shards::getCellTopologyData>(); - TEUCHOS_TEST_FOR_EXCEPTION ( - ctd.name!=topo_hexa->name && - ctd.name!=topo_wedge->name, std::runtime_error, - "Extruded meshes only allowed if there is one element per layer (hexa or wedges).\n" - " - current topology name: " << ctd.name << "\n"); - - const auto& sol_dof_mgr = getDOFManager(); - const auto& elem_dof_lids = sol_dof_mgr->elem_dof_lids().host(); - - // Build a LayeredMeshNumbering for cells, so we can get the LIDs of elems over the column - const auto numLayers = cell_layers_data->numLayers; - const int top = cell_layers_data->top_side_pos; - const int bot = cell_layers_data->bot_side_pos; - - // 7) Populate localDOFViews for GatherVerticallyContractedSolution - for (int ws=0; ws>& wsldofViews = wsLocalDOFViews[ws]; - - const auto& elem_lids = getElementLIDs_host(ws); - - // Loop over the sides that form the boundary condition - // const Teuchos::ArrayRCP >& wsElNodeID_i = wsElNodeID[i]; - for (auto& ss_it : sideSets[ws]) { - std::string ss_key = ss_it.first; - std::vector ss_val = ss_it.second; - - Kokkos::DualView& globalDOFView = allLocalDOFViews[ss_key]; - - for (unsigned int sideSet_idx = 0; sideSet_idx < ss_val.size(); ++sideSet_idx) { - const auto& side = ss_val[sideSet_idx]; - - // Get the data that corresponds to the side - const int ws_elem_idx = side.ws_elem_idx; - const int side_pos = side.side_pos; - - // Check if this sideset is the top or bot of the mesh. If not, the data structure - // for coupling vertical dofs is not needed. - if (side_pos!=top && side_pos!=bot) - break; - - const int elem_LID = elem_lids(ws_elem_idx); - const int basal_elem_LID = cell_layers_data->getColumnId(elem_LID); - - for (int eq=0; eqgetGIDFieldOffsetsSide(eq,top,side_pos); - const auto& sol_bot_offsets = sol_dof_mgr->getGIDFieldOffsetsSide(eq,bot,side_pos); - const int numSideNodes = sol_top_offsets.size(); - - for (int j=0; jgetId(basal_elem_LID,il); - globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il, eq) = - elem_dof_lids(layer_elem_LID,sol_bot_offsets[j]); - } - - // Add top side in last layer - const int il = numLayers-1; - const LO layer_elem_LID = cell_layers_data->getId(basal_elem_LID,il); - globalDOFView.h_view(sideSet_idx + sideset_idx_offset[ss_key], j, il+1, eq) = - elem_dof_lids(layer_elem_LID,sol_top_offsets[j]); - } - } - } - - globalDOFView.modify_host(); - globalDOFView.sync_device(); - - // Set workset-local sub-view - std::pair range(sideset_idx_offset[ss_key], sideset_idx_offset[ss_key]+ss_val.size()); - wsldofViews[ss_key] = Kokkos::subview(globalDOFView, range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - - sideset_idx_offset[ss_key] += ss_val.size(); - } - } - } else { - // We still need this view to be present (even if of size 0), so create them - std::map> dummy; - for (int ws=0; wscell_indexer(); std::vector> offsets (sol_dof_mgr->getNumFields()); - for (int eq=0; eqgetGIDFieldOffsets(eq); } @@ -1935,8 +1654,8 @@ STKDiscretization::buildSideSetProjectors() // Recall, if node_map(i,j)=k, then on side_gid=i, the j-th side node corresponds // to the k-th node in the side set cell numeration. - const auto& node_numeration_map = sideNodeNumerationMap.at(it.first); - const auto& side_cell_gid_map = sideToSideSetCellMap.at(it.first); + const auto& node_numeration_map = m_side_nodes_to_ss_cell_nodes.at(it.first); + const auto& side_cell_gid_map = m_side_to_ss_cell.at(it.first); for (auto side : sides) { TEUCHOS_TEST_FOR_EXCEPTION (bulkData->num_elements(side)!=1, std::logic_error, "Error! Found a side with not exactly one element attached.\n" @@ -1956,7 +1675,7 @@ STKDiscretization::buildSideSetProjectors() const auto ss_elem_dof_gids = ss_dofMgr->getElementGIDs(ss_elem_lid); const auto& permutation = node_numeration_map.at(side_gid); - for (int eq=0; eqgetGIDFieldOffsetsSubcell(eq,sideDim,pos); const auto& ss_offsets = ss_dofMgr->getGIDFieldOffsets(eq); for (int i=0; iupdateMesh(); + for (const auto& [ss_name, ss_disc] : sideSetDiscretizationsSTK) { + ss_disc->updateMesh(); stkMeshStruct->buildCellSideNodeNumerationMap( - it.first, - sideToSideSetCellMap[it.first], - sideNodeNumerationMap[it.first]); + ss_name, + m_side_to_ss_cell[ss_name], + m_side_nodes_to_ss_cell_nodes[ss_name]); } if (sideSetDiscretizations.size()>0) { @@ -2073,10 +1792,10 @@ void STKDiscretization::setFieldData() if (Teuchos::nonnull(mSTKFieldContainer)) { solutionFieldContainer = Teuchos::rcp(new MultiSTKFieldContainer( - params, stkMeshStruct->metaData, stkMeshStruct->bulkData, neq, numDim, solution_vector, num_params)); + params, stkMeshStruct->metaData, stkMeshStruct->bulkData, m_neq, numDim, solution_vector, num_params)); } else if (Teuchos::nonnull(oSTKFieldContainer)) { solutionFieldContainer = Teuchos::rcp(new OrdinarySTKFieldContainer( - params, stkMeshStruct->metaData, stkMeshStruct->bulkData, neq, numDim, num_params)); + params, stkMeshStruct->metaData, stkMeshStruct->bulkData, m_neq, numDim, num_params)); } else { ALBANY_ABORT ("Error! Failed to cast the AbstractSTKFieldContainer to a concrete type.\n"); } @@ -2085,11 +1804,17 @@ void STKDiscretization::setFieldData() Teuchos::RCP STKDiscretization:: create_dof_mgr (const std::string& part_name, - const std::string& field_name, const FE_Type fe_type, const int order, - const int dof_dim) const + const int dof_dim) { + auto& dof_mgr = get_dof_mgr(part_name,fe_type,order,dof_dim); + if (Teuchos::nonnull(dof_mgr)) { + // Not the first time we build a DOFManager for a field with these specs + return dof_mgr; + } + // Teuchos::RCP dof_mgr; + // Figure out which element blocks this part belongs to std::vector elem_blocks; const auto& ebn = stkMeshStruct->ebNames_; @@ -2127,7 +1852,7 @@ create_dof_mgr (const std::string& part_name, // Create conn and dof managers auto conn_mgr = Teuchos::rcp(new STKConnManager(metaData,bulkData,elem_blocks)); - auto dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,comm,part_name)); + dof_mgr = Teuchos::rcp(new DOFManager(conn_mgr,comm,part_name)); const auto& topo = stkMeshStruct->elementBlockCT_.at(elem_blocks[0]); Teuchos::RCP fp; @@ -2140,9 +1865,9 @@ create_dof_mgr (const std::string& part_name, fp = Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis)); } // NOTE: we add $dof_dim copies of the field pattern to the dof mgr, - // and call the fields ${field_name}_n, n=0,..,$dof_dim-1 + // and call the fields cmp_N, N=0,..,$dof_dim-1 for (int i=0; iaddField(field_name + "_" + std::to_string(i),fp); + dof_mgr->addField("cmp_" + std::to_string(i),fp); } dof_mgr->build(); @@ -2222,7 +1947,7 @@ adapt (const Teuchos::RCP& adaptData) discParams->set("1D Elements",factor*ne_x); auto sis = Teuchos::rcp(new StateInfoStruct(getMeshStruct()->get_field_accessor()->getAllSIS())); stkMeshStruct = Teuchos::rcp(new TmplSTKMeshStruct<1>(discParams,comm,num_params)); - stkMeshStruct->setFieldData(comm); + stkMeshStruct->setFieldData(comm,sis,{}); stkMeshStruct->getFieldContainer()->addStateStructs(sis); this->setFieldData(); stkMeshStruct->setBulkData(comm); diff --git a/src/disc/stk/Albany_STKDiscretization.hpp b/src/disc/stk/Albany_STKDiscretization.hpp index aa1df8fbf7..d2c4cd2858 100644 --- a/src/disc/stk/Albany_STKDiscretization.hpp +++ b/src/disc/stk/Albany_STKDiscretization.hpp @@ -70,27 +70,6 @@ class STKDiscretization : public AbstractDiscretization return nodeSetCoords; } - //! Get Side set lists (typedef in Albany_AbstractDiscretization.hpp) - const SideSetList& - getSideSets(const int workset) const - { - return sideSets[workset]; - } - - //! Get Side set lists (typedef in Albany_AbstractDiscretization.hpp) - const LocalSideSetInfoList& - getSideSetViews(const int workset) const - { - return sideSetViews.at(workset); - } - - //! Get local DOF views for GatherVerticallyContractedSolution - const std::map>& - getLocalDOFViews(const int workset) const - { - return wsLocalDOFViews.at(workset); - } - //! Get connectivity map from elementGID to workset WsLIDList& getElemGIDws() @@ -123,18 +102,6 @@ class STKDiscretization : public AbstractDiscretization return stkMeshStruct; } - const std::map>& - getSideToSideSetCellMap() const - { - return sideToSideSetCellMap; - } - - const std::map>>& - getSideNodeNumerationMap() const - { - return sideNodeNumerationMap; - } - //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const @@ -165,13 +132,6 @@ class STKDiscretization : public AbstractDiscretization return stkMeshStruct->numDim; } - //! Get number of total DOFs per node - int - getNumEq() const - { - return neq; - } - const stk::mesh::MetaData& getSTKMetaData() const { @@ -308,10 +268,9 @@ class STKDiscretization : public AbstractDiscretization // If node_as_elements=true, build the ConnMgr as if nodes are the "cells". Teuchos::RCP create_dof_mgr (const std::string& part_name, - const std::string& field_name, const FE_Type fe_type, const int order, - const int dof_dim) const; + const int dof_dim); // ==================== Members =================== // @@ -324,9 +283,6 @@ class STKDiscretization : public AbstractDiscretization //! Teuchos communicator Teuchos::RCP comm; - //! Number of equations (and unknowns) per node - const int neq; - //! Equations that are defined only on some side sets of the mesh std::map> sideSetEquations; @@ -338,16 +294,6 @@ class STKDiscretization : public AbstractDiscretization NodeSetGIDsList nodeSetGIDs; NodeSetCoordList nodeSetCoords; - //! side sets stored as std::map(string ID, SideArray classes) per workset - //! (std::vector across worksets) - std::vector sideSets; - GlobalSideSetList globalSideSetViews; - std::map sideSetViews; - - //! GatherVerticallyContractedSolution connectivity - std::map> allLocalDOFViews; - std::map>> wsLocalDOFViews; - mutable Teuchos::ArrayRCP coordinates; Teuchos::RCP coordMV; @@ -370,8 +316,6 @@ class STKDiscretization : public AbstractDiscretization // Sideset discretizations std::map> sideSetDiscretizationsSTK; - std::map> sideToSideSetCellMap; - std::map>> sideNodeNumerationMap; std::map> projectors; std::map> ov_projectors; diff --git a/src/disc/stk/Albany_TmplSTKMeshStruct.hpp b/src/disc/stk/Albany_TmplSTKMeshStruct.hpp index 67b2b31612..85ae1bf6ed 100644 --- a/src/disc/stk/Albany_TmplSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_TmplSTKMeshStruct.hpp @@ -92,7 +92,9 @@ class TmplSTKMeshStruct : public GenericSTKMeshStruct { ~TmplSTKMeshStruct() = default; //! Sets mesh generation parameters - void setFieldData (const Teuchos::RCP& comm); + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis); void setBulkData (const Teuchos::RCP& comm); diff --git a/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp b/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp index e12aa61cfd..48dac6198e 100644 --- a/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp +++ b/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp @@ -298,7 +298,9 @@ TmplSTKMeshStruct::TmplSTKMeshStruct( template void TmplSTKMeshStruct:: -setFieldData (const Teuchos::RCP& comm) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis, + std::map > side_set_sis) { // Create global mesh: Dim-D structured, rectangular std::vector h_dim[traits_type::size]; @@ -318,7 +320,7 @@ setFieldData (const Teuchos::RCP& comm) x[idx][i] = x[idx][i - 1] + h_dim[idx][i - 1]; // place the coordinates of the element nodes } - GenericSTKMeshStruct::setFieldData(comm); + GenericSTKMeshStruct::setFieldData(comm,sis,side_set_sis); } template diff --git a/tests/unit/disc/generic/DummyMesh.hpp b/tests/unit/disc/generic/DummyMesh.hpp index 5efb57da8f..86901c70c5 100644 --- a/tests/unit/disc/generic/DummyMesh.hpp +++ b/tests/unit/disc/generic/DummyMesh.hpp @@ -180,7 +180,9 @@ struct DummyMesh : public AbstractMeshStruct { //! Internal mesh specs type needed std::string meshLibName() const override { return "dummy"; } - void setFieldData (const Teuchos::RCP& /* comm */) override {} + void setFieldData (const Teuchos::RCP& /* comm */, + const Teuchos::RCP& /* sis */, + std::map > /* side_set_sis */) override {} void setBulkData(const Teuchos::RCP& /* comm */) override {}