diff --git a/src/Albany_Application.cpp b/src/Albany_Application.cpp index 1cd5cc4952..b1b2ac11ef 100644 --- a/src/Albany_Application.cpp +++ b/src/Albany_Application.cpp @@ -3376,7 +3376,7 @@ Application::setScaleBCDofs( // nodeset specified via the sideset discretization. // /*const auto &sdn = - disc->getMeshStruct()->getMeshSpecs()[0]->sideSetMeshNames; + disc->getMeshStruct()->meshSpecs[0]->sideSetMeshNames; for (int isd = 0; isd < sdn.size(); ++isd) { const auto &sd = disc->getSideSetDiscretizations().at(sdn[isd]); for (auto iterator = sd->getNodeSets().begin(); diff --git a/src/SolutionManager.cpp b/src/SolutionManager.cpp index ab662c5b3a..14a9841662 100644 --- a/src/SolutionManager.cpp +++ b/src/SolutionManager.cpp @@ -84,7 +84,9 @@ SolutionManager::SolutionManager( // This call allocates the non-overlapped MV current_soln = Thyra::createMembers(owned_vs,num_time_deriv+1); - disc_->getSolutionMV(*current_soln); + if (disc_->hasRestartSolution()) { + disc_->getSolutionMV(*current_soln); + } // Create the CombineAndScatterManager for handling distributed memory linear // algebra communications diff --git a/src/corePDEs/problems/Albany_PopulateMesh.cpp b/src/corePDEs/problems/Albany_PopulateMesh.cpp index 2659e9ce4a..dbfc96c2ad 100644 --- a/src/corePDEs/problems/Albany_PopulateMesh.cpp +++ b/src/corePDEs/problems/Albany_PopulateMesh.cpp @@ -26,11 +26,6 @@ PopulateMesh::PopulateMesh (const Teuchos::RCP& params_, neq = 1; } -PopulateMesh::~PopulateMesh() -{ - // Nothing to be done here -} - void PopulateMesh::buildProblem (Teuchos::ArrayRCP> meshSpecs, StateManager& stateMgr) { diff --git a/src/corePDEs/problems/Albany_PopulateMesh.hpp b/src/corePDEs/problems/Albany_PopulateMesh.hpp index c2ffb2fb70..27e0c3855d 100644 --- a/src/corePDEs/problems/Albany_PopulateMesh.hpp +++ b/src/corePDEs/problems/Albany_PopulateMesh.hpp @@ -44,7 +44,7 @@ class PopulateMesh : public Albany::AbstractProblem const Teuchos::RCP& paramLib_); //! Destructor - ~PopulateMesh(); + ~PopulateMesh() = default; //! Return number of spatial dimensions virtual int spatialDimension() const { return 0; } @@ -66,16 +66,6 @@ class PopulateMesh : public Albany::AbstractProblem bool useSDBCs() const { return use_sdbcs_; } -private: - - //! Private to prohibit copying - PopulateMesh(const PopulateMesh&); - - //! Private to prohibit copying - PopulateMesh& operator=(const PopulateMesh&); - -public: - //! Main problem setup routine. Not directly called, but indirectly by following functions template Teuchos::RCP constructEvaluators (PHX::FieldManager& fm0, diff --git a/src/disc/Albany_AbstractMeshStruct.hpp b/src/disc/Albany_AbstractMeshStruct.hpp index 11291d9c67..f01fbfebff 100644 --- a/src/disc/Albany_AbstractMeshStruct.hpp +++ b/src/disc/Albany_AbstractMeshStruct.hpp @@ -36,34 +36,22 @@ struct AbstractMeshStruct { //! Internal mesh specs type needed virtual std::string meshType() const = 0; - virtual void setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}) = 0; - - virtual void setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}) = 0; - - void setFieldAndBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}) - { - setFieldData(commT, sis, worksetSize, side_set_sis); - setBulkData(commT, sis, worksetSize, side_set_sis); - } - - virtual Teuchos::ArrayRCP >& getMeshSpecs() = 0; - virtual const Teuchos::ArrayRCP >& getMeshSpecs() const = 0; + virtual void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) = 0; + + virtual void setBulkData(const Teuchos::RCP& comm) = 0; + + bool isBulkDataSet () const { return m_bulk_data_set; } Teuchos::RCP > global_cell_layers_data; Teuchos::RCP > local_cell_layers_data; Teuchos::ArrayRCP mesh_layers_ratio; + + Teuchos::ArrayRCP > meshSpecs; + std::map> sideSetMeshStructs; + protected: + + bool m_bulk_data_set = false; }; } // Namespace Albany diff --git a/src/disc/Albany_DiscretizationFactory.cpp b/src/disc/Albany_DiscretizationFactory.cpp index 65ea49eb98..8831f1a023 100644 --- a/src/disc/Albany_DiscretizationFactory.cpp +++ b/src/disc/Albany_DiscretizationFactory.cpp @@ -33,9 +33,9 @@ namespace Albany { DiscretizationFactory::DiscretizationFactory( const Teuchos::RCP& topLevelParams, - const Teuchos::RCP& commT_, + const Teuchos::RCP& comm_, const bool explicit_scheme_) : - commT(commT_), + comm(comm_), explicit_scheme(explicit_scheme_) { discParams = Teuchos::sublist(topLevelParams, "Discretization", true); @@ -51,8 +51,8 @@ DiscretizationFactory::DiscretizationFactory( Teuchos::ArrayRCP > DiscretizationFactory::createMeshSpecs() { // First, create the mesh struct - meshStruct = createMeshStruct(discParams, commT, num_params); - return meshStruct->getMeshSpecs(); + meshStruct = createMeshStruct(discParams, comm, num_params); + return meshStruct->meshSpecs; } Teuchos::RCP @@ -60,22 +60,23 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis Teuchos::RCP comm, const int numParams) { std::string& method = disc_params->get("Method", "STK1D"); + Teuchos::RCP mesh; if (method == "STK1D") { - return Teuchos::rcp(new TmplSTKMeshStruct<1>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new TmplSTKMeshStruct<1>(disc_params, comm, numParams)); } else if (method == "STK0D") { - return Teuchos::rcp(new TmplSTKMeshStruct<0>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new TmplSTKMeshStruct<0>(disc_params, comm, numParams)); } else if (method == "STK2D") { - return Teuchos::rcp(new TmplSTKMeshStruct<2>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new TmplSTKMeshStruct<2>(disc_params, comm, numParams)); } else if (method == "STK3D") { - return Teuchos::rcp(new TmplSTKMeshStruct<3>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new TmplSTKMeshStruct<3>(disc_params, comm, numParams)); } else if (method == "STK3D") { - return Teuchos::rcp(new TmplSTKMeshStruct<3>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new TmplSTKMeshStruct<3>(disc_params, comm, numParams)); } else if (method == "STK3DPoint") { - return Teuchos::rcp(new STK3DPointStruct(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new STK3DPointStruct(disc_params, comm, numParams)); } else if (method == "Ioss" || method == "Exodus" || method == "Pamgen") { #ifdef ALBANY_SEACAS - return Teuchos::rcp(new IossSTKMeshStruct(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new IossSTKMeshStruct(disc_params, comm, numParams)); #else TEUCHOS_TEST_FOR_EXCEPTION(method == "Ioss" || method == "Exodus" || method == "Pamgen", Teuchos::Exceptions::InvalidParameter, @@ -85,25 +86,25 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis } #ifdef ALBANY_OMEGAH else if (method == "Box1D") { - return Teuchos::rcp(new OmegahBoxMesh<1>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new OmegahBoxMesh<1>(disc_params, comm, numParams)); } else if (method == "Box2D") { - return Teuchos::rcp(new OmegahBoxMesh<2>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new OmegahBoxMesh<2>(disc_params, comm, numParams)); } else if (method == "Box3D") { - return Teuchos::rcp(new OmegahBoxMesh<3>(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new OmegahBoxMesh<3>(disc_params, comm, numParams)); } #endif else if (method == "Ascii") { - return Teuchos::rcp(new AsciiSTKMeshStruct(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new AsciiSTKMeshStruct(disc_params, comm, numParams)); } else if (method == "Ascii2D") { - return Teuchos::rcp(new AsciiSTKMesh2D(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new AsciiSTKMesh2D(disc_params, comm, numParams)); #ifdef ALBANY_SEACAS // Fails to compile without SEACAS } else if (method == "Hacky Ascii2D") { //FixME very hacky! needed for printing 2d mesh Teuchos::RCP meshStruct2D; meshStruct2D = Teuchos::rcp(new AsciiSTKMesh2D(disc_params, comm, numParams)); Teuchos::RCP sis = Teuchos::rcp(new StateInfoStruct); - meshStruct2D->setFieldAndBulkData(comm, - sis, meshStruct2D->getMeshSpecs()[0]->worksetSize); + meshStruct2D->setFieldData(comm, sis); + meshStruct2D->setBulkData(comm); Ioss::Init::Initializer io; Teuchos::RCP mesh_data = Teuchos::rcp(new stk::io::StkMeshIoBroker(MPI_COMM_WORLD)); mesh_data->set_bulk_data(*meshStruct2D->bulkData); @@ -112,7 +113,7 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis mesh_data->process_output_request(idx, 0.0); #endif // ALBANY_SEACAS } else if (method == "Gmsh") { - return Teuchos::rcp(new GmshSTKMeshStruct(disc_params, comm, numParams)); + mesh = Teuchos::rcp(new GmshSTKMeshStruct(disc_params, comm, numParams)); } else if (method == "Extruded") { Teuchos::RCP basalMesh; @@ -140,30 +141,84 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP dis } basalMesh = createMeshStruct(basal_params, comm, numParams); - return Teuchos::rcp(new ExtrudedSTKMeshStruct(disc_params, comm, basalMesh, numParams)); + mesh = Teuchos::rcp(new ExtrudedSTKMeshStruct(disc_params, comm, basalMesh, numParams)); } else if (method == "Cubit") { TEUCHOS_TEST_FOR_EXCEPTION(method == "Cubit", Teuchos::Exceptions::InvalidParameter, "Error: Discretization method " << method << " requested, but no longer supported as of 10/2017" << std::endl); - } else + } else { TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, std::endl << "Error! Unknown discretization method in DiscretizationFactory: " << method << "!" << std::endl << "Supplied parameter list is " << std::endl << *disc_params << "\nValid Methods are: STK1D, STK2D, STK3D, STK3DPoint, Ioss," << " Exodus, Ascii," << " Ascii2D, Extruded" << std::endl); + } - return Teuchos::null; -} + if (disc_params->isSublist ("Side Set Discretizations")) { + TEUCHOS_TEST_FOR_EXCEPTION (mesh->meshSpecs.size()!=1, std::logic_error, + "Error! So far, side set mesh is allowed only for meshes with 1 element block.\n"); + auto ms = mesh->meshSpecs[0]; -Teuchos::RCP -DiscretizationFactory::createDiscretization(unsigned int neq, - const Teuchos::RCP& sis, - const Teuchos::RCP& rigidBodyModes) -{ - return createDiscretization(neq, empty_side_set_equations, sis, empty_side_set_sis, rigidBodyModes); + const Teuchos::ParameterList& ssd_list = disc_params->sublist("Side Set Discretizations"); + const Teuchos::Array& sideSets = ssd_list.get >("Side Sets"); + + 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 + if (ss_mesh.is_null()) { + params_ss = Teuchos::rcp(new Teuchos::ParameterList(ssd_list.sublist(ss_name))); + + if (!params_ss->isParameter("Number Of Time Derivatives")) + params_ss->set("Number Of Time Derivatives",disc_params->get("Number Of Time Derivatives")); + + // Set sideset discretization workset size based on sideset mesh spec if a single workset is used + const auto &sideSetMeshSpecs = ms->sideSetMeshSpecs; + auto sideSetMeshSpecIter = sideSetMeshSpecs.find(ss_name); + TEUCHOS_TEST_FOR_EXCEPTION(sideSetMeshSpecIter == sideSetMeshSpecs.end(), std::runtime_error, + "Cannot find " << ss_name << " in sideSetMeshSpecs!\n"); + + std::string ss_method = params_ss->get("Method"); + if (ss_method=="SideSetSTK") { + ss_mesh = Teuchos::rcp(new SideSetSTKMeshStruct(*ms, params_ss, comm, numParams)); + + auto mesh_stk = Teuchos::rcp_dynamic_cast(mesh,true); + auto ss_mesh_stk = Teuchos::rcp_dynamic_cast(ss_mesh,true); + ss_mesh_stk->setParentMeshInfo(*mesh_stk, ss_name); + + // If requested, we ignore the side maps already stored in the imported side mesh (if any) + // This can be useful for side mesh of an extruded mesh, in the case it was constructed + // as side mesh of an extruded mesh with a different ordering and/or different number + // of layers. Notice that if that's the case, it probably is impossible to build a new + // set of maps, since there is no way to correctly map the side nodes to the cell nodes. + ss_mesh_stk->ignore_side_maps = params_ss->get("Ignore Side Maps", false); + } else { + // This can be the case if we restart from existing volume and side meshes + ss_mesh = createMeshStruct (params_ss,comm, numParams); + } + } + + auto ss_ms = ss_mesh->meshSpecs; + + // Checking that the side mesh has the correct dimension (in case they were loaded from file, + // and the user mistakenly gave the wrong file name) + TEUCHOS_TEST_FOR_EXCEPTION (sideDim!=ss_ms[0]->numDim, std::logic_error, + "Error! Mesh on side " << ss_name << " has the wrong dimension.\n"); + } + + auto stk_mesh = Teuchos::rcp_dynamic_cast(mesh); + if (stk_mesh) + stk_mesh->createSideMeshMaps(); + } + + return mesh; } Teuchos::RCP @@ -191,7 +246,7 @@ DiscretizationFactory::createDiscretization( setFieldData(it.second,{}); } } - setMeshStructBulkData(sis, side_set_sis); + setMeshStructBulkData(); completeDiscSetup(result); return result; @@ -200,7 +255,7 @@ DiscretizationFactory::createDiscretization( Teuchos::ArrayRCP > DiscretizationFactory::createMeshSpecs(Teuchos::RCP mesh) { meshStruct = mesh; - return meshStruct->getMeshSpecs(); + return meshStruct->meshSpecs; } void @@ -214,54 +269,48 @@ DiscretizationFactory::setMeshStructFieldData( const Teuchos::RCP& sis, const std::map >& side_set_sis) { - TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: setMeshStructFieldData"); - meshStruct->setFieldData(commT, sis, - meshStruct->getMeshSpecs()[0]->worksetSize, side_set_sis); -} - -void -DiscretizationFactory::setMeshStructBulkData( - const Teuchos::RCP& sis) { - setMeshStructBulkData(sis, empty_side_set_sis); + TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: setMeshStructFieldData"); + meshStruct->setFieldData(comm, sis); + for (auto& it : meshStruct->sideSetMeshStructs) { + auto this_ss_sis = side_set_sis.count(it.first)>0 ? side_set_sis.at(it.first) : Teuchos::null; + it.second->setFieldData(comm,this_ss_sis); + } } -void -DiscretizationFactory::setMeshStructBulkData( - const Teuchos::RCP& sis, - const std::map >& side_set_sis) +void DiscretizationFactory:: +setMeshStructBulkData() { - TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: setMeshStructBulkData"); - meshStruct->setBulkData(commT, sis, - meshStruct->getMeshSpecs()[0]->worksetSize, side_set_sis); -} - -Teuchos::RCP -DiscretizationFactory::createDiscretizationFromInternalMeshStruct( - const int neq, - const Teuchos::RCP& rigidBodyModes) { - return createDiscretizationFromInternalMeshStruct(neq, empty_side_set_equations, rigidBodyModes); + TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: setMeshStructBulkData"); + meshStruct->setBulkData(comm); + for (auto& it : meshStruct->sideSetMeshStructs) { + // For extruded meshes, the bulk data of the basal mesh + // should be set from inside the extruded mesh call, + // during the 'setBulkData' call above + if (not it.second->isBulkDataSet()) { + it.second->setBulkData(comm); + } + } } Teuchos::RCP DiscretizationFactory::createDiscretizationFromInternalMeshStruct( const int neq, const std::map >& sideSetEquations, - const Teuchos::RCP& rigidBodyModes) { - TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: createDiscretizationFromInternalMeshStruct"); - - if (!piroParams.is_null() && !rigidBodyModes.is_null()) - - rigidBodyModes->setPiroPL(piroParams); + const Teuchos::RCP& rigidBodyModes) +{ + TEUCHOS_FUNC_TIME_MONITOR("Albany_DiscrFactory: createDiscretizationFromInternalMeshStruct"); + if (!piroParams.is_null() && !rigidBodyModes.is_null()) + rigidBodyModes->setPiroPL(piroParams); Teuchos::RCP disc; if (meshStruct->meshType()=="STK") { auto ms = Teuchos::rcp_dynamic_cast(meshStruct); - disc = Teuchos::rcp(new STKDiscretization(discParams, neq, ms, commT, rigidBodyModes, sideSetEquations)); + disc = Teuchos::rcp(new STKDiscretization(discParams, neq, ms, comm, rigidBodyModes, sideSetEquations)); #ifdef ALBANY_OMEGAH } else if (meshStruct->meshType()=="Omega_h") { auto ms = Teuchos::rcp_dynamic_cast(meshStruct); - disc = Teuchos::rcp(new OmegahDiscretization(discParams, neq, ms, commT, rigidBodyModes, sideSetEquations)); + disc = Teuchos::rcp(new OmegahDiscretization(discParams, neq, ms, comm, rigidBodyModes, sideSetEquations)); #endif } return disc; @@ -269,7 +318,7 @@ DiscretizationFactory::createDiscretizationFromInternalMeshStruct( void DiscretizationFactory::setFieldData(Teuchos::RCP disc, - const Teuchos::RCP& sis) { + const Teuchos::RCP& sis) { if (meshStruct->meshType()=="STK") { auto stk_disc = Teuchos::rcp_dynamic_cast(disc); diff --git a/src/disc/Albany_DiscretizationFactory.hpp b/src/disc/Albany_DiscretizationFactory.hpp index 6da5a76583..302d3c588a 100644 --- a/src/disc/Albany_DiscretizationFactory.hpp +++ b/src/disc/Albany_DiscretizationFactory.hpp @@ -29,94 +29,70 @@ class DiscretizationFactory { //! Default constructor DiscretizationFactory( const Teuchos::RCP& topLevelParams, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, const bool explicit_scheme_ = false); //! Destructor ~DiscretizationFactory() {} - static Teuchos::RCP + Teuchos::RCP createMeshStruct (Teuchos::RCP disc_params, Teuchos::RCP comm, const int numParams); - Teuchos::RCP getMeshStruct() { - return meshStruct; - } + Teuchos::ArrayRCP > createMeshSpecs(); - Teuchos::ArrayRCP > createMeshSpecs(); + Teuchos::ArrayRCP > createMeshSpecs(Teuchos::RCP mesh); - Teuchos::ArrayRCP > createMeshSpecs(Teuchos::RCP mesh); - - Teuchos::RCP - createDiscretization(unsigned int num_equations, - const Teuchos::RCP& sis, - const Teuchos::RCP& rigidBodyModes = Teuchos::null); - Teuchos::RCP + Teuchos::RCP createDiscretization(unsigned int num_equations, const std::map >& sideSetEquations, - const Teuchos::RCP& sis, - const std::map >& side_set_sis, - const Teuchos::RCP& rigidBodyModes = Teuchos::null); + const Teuchos::RCP& sis, + const std::map >& side_set_sis, + const Teuchos::RCP& rigidBodyModes = Teuchos::null); void setMeshStructFieldData( - const Teuchos::RCP& sis); - - void - setMeshStructFieldData( - const Teuchos::RCP& sis, - const std::map >& side_set_sis); + const Teuchos::RCP& sis); void - setMeshStructBulkData( - const Teuchos::RCP& sis); + setMeshStructFieldData( + const Teuchos::RCP& sis, + const std::map >& side_set_sis); - void - setMeshStructBulkData( - const Teuchos::RCP& sis, - const std::map >& side_set_sis); + void setMeshStructBulkData(); - Teuchos::RCP createDiscretizationFromInternalMeshStruct( + Teuchos::RCP createDiscretizationFromInternalMeshStruct( const int neq, - const Teuchos::RCP& rigidBodyModes); + const Teuchos::RCP& rigidBodyModes); - Teuchos::RCP createDiscretizationFromInternalMeshStruct( + Teuchos::RCP createDiscretizationFromInternalMeshStruct( const int neq, const std::map >& sideSetEquations, - const Teuchos::RCP& rigidBodyModes); + const Teuchos::RCP& rigidBodyModes); /* This function overwrite previous discretization parameter list */ void setDiscretizationParameters(Teuchos::RCP disc_params); - private: - - //! Private to prohibit copying - DiscretizationFactory(const DiscretizationFactory&); - - //! Private to prohibit copying - DiscretizationFactory& operator=(const DiscretizationFactory&); + protected: const std::map > empty_side_set_equations; - const std::map > empty_side_set_sis; + const std::map > empty_side_set_sis; - void - setFieldData(Teuchos::RCP disc, - const Teuchos::RCP& sis); + void setFieldData (Teuchos::RCP disc, + const Teuchos::RCP& sis); void completeDiscSetup(Teuchos::RCP disc); - protected: - //! Parameter list specifying what element to create Teuchos::RCP discParams; //! Parameter list specifying solver parameters Teuchos::RCP piroParams; - Teuchos::RCP commT; + Teuchos::RCP comm; //The following are for Aeras hydrostatic problems int numLevels; @@ -125,7 +101,7 @@ class DiscretizationFactory { //Flag for explicit time-integration scheme, used in Aeras bool explicit_scheme; - Teuchos::RCP meshStruct; + Teuchos::RCP meshStruct; //Number of parameters int num_params = 0; diff --git a/src/disc/omegah/Albany_OmegahBoxMesh.hpp b/src/disc/omegah/Albany_OmegahBoxMesh.hpp index 9edaf55f10..4d06e20e07 100644 --- a/src/disc/omegah/Albany_OmegahBoxMesh.hpp +++ b/src/disc/omegah/Albany_OmegahBoxMesh.hpp @@ -13,10 +13,7 @@ class OmegahBoxMesh : public OmegahGenericMesh { OmegahBoxMesh (const Teuchos::RCP& params, const Teuchos::RCP& comm, const int numParams); - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) override; + void setBulkData (const Teuchos::RCP& comm) override; Omega_h::Read create_ns_tag (const std::string& name, diff --git a/src/disc/omegah/Albany_OmegahBoxMesh_Def.hpp b/src/disc/omegah/Albany_OmegahBoxMesh_Def.hpp index 5896819e0c..dc5fe26afc 100644 --- a/src/disc/omegah/Albany_OmegahBoxMesh_Def.hpp +++ b/src/disc/omegah/Albany_OmegahBoxMesh_Def.hpp @@ -119,8 +119,8 @@ OmegahBoxMesh (const Teuchos::RCP& params, case 3: ctd = shards::getCellTopologyData>(); break; } - this->m_mesh_specs.resize(1); - this->m_mesh_specs[0] = Teuchos::rcp(new MeshSpecsStruct(*ctd, Dim, + this->meshSpecs.resize(1); + this->meshSpecs[0] = Teuchos::rcp(new MeshSpecsStruct(*ctd, Dim, nsNames, ssNames, m_mesh.nelems(), ebName, ebNameToIndex)); } @@ -142,14 +142,10 @@ create_ns_tag (const std::string& name, } template -void OmegahBoxMesh::setBulkData( - const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void OmegahBoxMesh:: +setBulkData (const Teuchos::RCP& /* comm */) { - // We can finally extract the side set meshes and set the fields and bulk data in all of them - // this->setSideSetBulkData(comm, side_set_sis, worksetSize); + m_bulk_data_set = true; } template diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index c2a378e5ef..30e55266ef 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -49,7 +49,7 @@ updateMesh () // Compute workset information // NOTE: these arrays are all of size 1, for the foreseable future. // Still, make impl generic (where possible), in case things change. - const auto& ms = m_mesh_struct->getMeshSpecs()[0]; + const auto& ms = m_mesh_struct->meshSpecs[0]; const auto& mesh = m_mesh_struct->getOmegahMesh(); int nelems = mesh.nelems(); int max_ws_size = ms->worksetSize; @@ -115,7 +115,7 @@ updateMesh () void OmegahDiscretization:: computeNodeSets () { - const auto& nsNames = getMeshStruct()->getMeshSpecs()[0]->nsNames; + const auto& nsNames = getMeshStruct()->meshSpecs[0]->nsNames; using Omega_h::I32; using Omega_h::I8; @@ -271,7 +271,7 @@ create_dof_mgr (const std::string& field_name, const int order, const int dof_dim) const { - const auto& mesh_specs = m_mesh_struct->getMeshSpecs()[0]; + const auto& mesh_specs = m_mesh_struct->meshSpecs[0]; // Create conn and dof managers auto conn_mgr = Teuchos::rcp(new OmegahConnManager(m_mesh_struct)); diff --git a/src/disc/omegah/Albany_OmegahDiscretization.hpp b/src/disc/omegah/Albany_OmegahDiscretization.hpp index 2e72d8cc91..72fcd5b3e7 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.hpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.hpp @@ -131,7 +131,7 @@ class OmegahDiscretization : public AbstractDiscretization //! Get number of spatial dimensions int getNumDim() const override { - return m_mesh_struct->getMeshSpecs()[0]->numDim; + return m_mesh_struct->meshSpecs[0]->numDim; } //! Get number of total DOFs per node diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index 2b0573b23d..4149a3f45b 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -5,10 +5,8 @@ namespace Albany { void OmegahGenericMesh:: -setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) { m_field_accessor = Teuchos::rcp(new OmegahMeshFieldAccessor(m_mesh)); if (not sis.is_null()) { diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.hpp b/src/disc/omegah/Albany_OmegahGenericMesh.hpp index 4f811d5ed4..4dd1f5ab5d 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.hpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.hpp @@ -18,17 +18,8 @@ class OmegahGenericMesh : public AbstractMeshStruct { // ------------- Override from base class ------------- // std::string meshType () const override { return "Omega_h"; } - Teuchos::ArrayRCP>& - getMeshSpecs() override { return m_mesh_specs; } - - const Teuchos::ArrayRCP>& - getMeshSpecs() const override { return m_mesh_specs; } - - - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) override; + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) override; // ------------- Omegah specific methods -------------- // @@ -63,8 +54,6 @@ class OmegahGenericMesh : public AbstractMeshStruct { // Given a part name, returns its topology (in the form of an Omega_h enum std::map m_part_topo; - Teuchos::ArrayRCP > m_mesh_specs; - Teuchos::RCP m_field_accessor; DeviceView1d m_coords_d; diff --git a/src/disc/omegah/Albany_OmegahOshMesh.cpp b/src/disc/omegah/Albany_OmegahOshMesh.cpp index 76d8d7e14a..3ef3ca8a21 100644 --- a/src/disc/omegah/Albany_OmegahOshMesh.cpp +++ b/src/disc/omegah/Albany_OmegahOshMesh.cpp @@ -93,8 +93,8 @@ OmegahOshMesh (const Teuchos::RCP& params, this->declare_part(ebName,elem_topo); // Omega_h does not know what worksets are, so all elements are in one workset - this->m_mesh_specs.resize(1); - this->m_mesh_specs[0] = Teuchos::rcp(new MeshSpecsStruct(*ctd, m_mesh.dim(), + this->meshSpecs.resize(1); + this->meshSpecs[0] = Teuchos::rcp(new MeshSpecsStruct(*ctd, m_mesh.dim(), nsNames, ssNames, m_mesh.nelems(), ebName, ebNameToIndex)); } diff --git a/src/disc/omegah/Albany_OmegahOshMesh.hpp b/src/disc/omegah/Albany_OmegahOshMesh.hpp index fb57c7b698..4345e8f650 100644 --- a/src/disc/omegah/Albany_OmegahOshMesh.hpp +++ b/src/disc/omegah/Albany_OmegahOshMesh.hpp @@ -11,12 +11,10 @@ class OmegahOshMesh : public OmegahGenericMesh OmegahOshMesh (const Teuchos::RCP& params, const Teuchos::RCP& comm, const int numParams); - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) override + void setBulkData (const Teuchos::RCP& comm) { throw NotYetImplemented("OmegahOshMesh::setBulkData"); + m_bulk_data_set = true; } }; diff --git a/src/disc/omegah/OmegahConnManager.cpp b/src/disc/omegah/OmegahConnManager.cpp index 9df0ccb3b5..5958f7cfad 100644 --- a/src/disc/omegah/OmegahConnManager.cpp +++ b/src/disc/omegah/OmegahConnManager.cpp @@ -26,7 +26,7 @@ namespace Albany { Omega_h::Read getIsEntInPart(const OmegahGenericMesh& albanyMesh, const std::string& part_name) { const auto& mesh = albanyMesh.getOmegahMesh(); const int part_dim = albanyMesh.part_dim(part_name); - if(part_name == albanyMesh.getMeshSpecs()[0]->ebName) { + if(part_name == albanyMesh.meshSpecs[0]->ebName) { return Omega_h::Read(mesh.nents(part_dim), 1); } else { return mesh.get_array(part_dim, part_name); @@ -54,7 +54,7 @@ Omega_h::LOs numberEntsInPart(const OmegahGenericMesh& albanyMesh, const std::st LO getNumEntsInPart(const OmegahGenericMesh& albanyMesh, const std::string& part_name) { const auto& mesh = albanyMesh.getOmegahMesh(); const int part_dim = albanyMesh.part_dim(part_name); - if(part_name == albanyMesh.getMeshSpecs()[0]->ebName) { + if(part_name == albanyMesh.meshSpecs[0]->ebName) { return mesh.nents(part_dim); } else { const auto isInPart = getIsEntInPart(albanyMesh,part_name); @@ -65,7 +65,7 @@ LO getNumEntsInPart(const OmegahGenericMesh& albanyMesh, const std::string& part [[nodiscard]] std::vector getLocalElmIds(const OmegahGenericMesh& albanyMesh, const std::string& part_name) { const auto& mesh = albanyMesh.getOmegahMesh(); - if(part_name==albanyMesh.getMeshSpecs()[0]->ebName) { + if(part_name==albanyMesh.meshSpecs[0]->ebName) { std::vector localElmIds(mesh.nelems()); std::iota(localElmIds.begin(), localElmIds.end(), 0); return localElmIds; @@ -79,7 +79,7 @@ std::vector getLocalElmIds(const OmegahGenericMesh& albanyMesh, const std:: OmegahConnManager:: OmegahConnManager(const Teuchos::RCP& in_mesh) : OmegahConnManager(in_mesh, - in_mesh->getMeshSpecs()[0]->ebName) + in_mesh->meshSpecs[0]->ebName) { // Nothing to do here } diff --git a/src/disc/stk/Albany_AbstractSTKMeshStruct.hpp b/src/disc/stk/Albany_AbstractSTKMeshStruct.hpp index 89a6e16986..c519e276b3 100644 --- a/src/disc/stk/Albany_AbstractSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_AbstractSTKMeshStruct.hpp @@ -45,7 +45,7 @@ struct AbstractSTKMeshStruct : public AbstractMeshStruct virtual ~AbstractSTKMeshStruct() = default; public: - std::string meshType () const { return "STK"; } + std::string meshType () const override { return "STK"; } Teuchos::RCP metaData; Teuchos::RCP bulkData; @@ -133,9 +133,6 @@ struct AbstractSTKMeshStruct : public AbstractMeshStruct virtual double restartDataTime() const = 0; - virtual bool - useCompositeTet() = 0; - // Flag for transforming STK mesh; currently only needed for LandIce/Aeras // problems std::string transformType; @@ -166,9 +163,6 @@ struct AbstractSTKMeshStruct : public AbstractMeshStruct // Info for periodic BCs -- only for hand-coded STK meshes struct PeriodicBCStruct PBCStruct; - std::map> sideSetMeshStructs; - - bool fieldAndBulkDataSet; virtual void buildCellSideNodeNumerationMap( @@ -177,8 +171,8 @@ struct AbstractSTKMeshStruct : public AbstractMeshStruct std::map>& sideNodeMap) = 0; // Useful for loading side meshes from file - bool side_maps_present; - bool ignore_side_maps; + bool side_maps_present = false; + bool ignore_side_maps = false; protected: Teuchos::RCP fieldContainer; diff --git a/src/disc/stk/Albany_AsciiSTKMesh2D.cpp b/src/disc/stk/Albany_AsciiSTKMesh2D.cpp index d6dcb8f5f7..a099a25eef 100644 --- a/src/disc/stk/Albany_AsciiSTKMesh2D.cpp +++ b/src/disc/stk/Albany_AsciiSTKMesh2D.cpp @@ -32,11 +32,14 @@ #include "Albany_Utils.hpp" -Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams) : - GenericSTKMeshStruct (params, 2, numParams), - periodic (false) +namespace Albany { + +AsciiSTKMesh2D:: +AsciiSTKMesh2D (const Teuchos::RCP& params, + const Teuchos::RCP& comm, + const int numParams) + : GenericSTKMeshStruct (params, 2, numParams) + , periodic (false) { NumElemNodes = NumNodes = NumElems = NumBdEdges = 0; @@ -45,7 +48,7 @@ Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCPgetRank() == 0) + if (comm->getRank() == 0) { std::ifstream ifile; @@ -166,9 +169,9 @@ Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCPvalidateParameters(*getValidDiscretizationParameters(), 0); std::string ebn = "Element Block 0"; - if (commT->getSize() > 0){ // running in parallel, need to exchange the topology with PE 0 + if (comm->getSize() > 0){ // running in parallel, need to exchange the topology with PE 0 int enum_to_int = static_cast(etopology); - Teuchos::broadcast(*commT,0,1,&enum_to_int); + Teuchos::broadcast(*comm,0,1,&enum_to_int); etopology = static_cast(enum_to_int); } partVec.push_back(&metaData->declare_part_with_topology(ebn, etopology)); @@ -207,12 +210,12 @@ Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCP(*commT, 0, &numBdNodeTags); + Teuchos::broadcast(*comm, 0, &numBdNodeTags); std::vector bdNodeTagsArray(numBdNodeTags); std::set::iterator it=bdNodeTags.begin(); for (int k=0; it!=bdNodeTags.end(); ++it,++k) bdNodeTagsArray[k] = *it; - Teuchos::broadcast(*commT, 0, numBdNodeTags, bdNodeTagsArray.data()); + Teuchos::broadcast(*comm, 0, numBdNodeTags, bdNodeTagsArray.data()); // Adding boundary nodesets and sidesets separating different labels for (int k=0; k(*commT, 0, &numBdEdgeTags); + Teuchos::broadcast(*comm, 0, &numBdEdgeTags); std::vector bdEdgeTagsArray(numBdEdgeTags); it=bdEdgeTags.begin(); for (int k=0; it!=bdEdgeTags.end(); ++it,++k) bdEdgeTagsArray[k] = *it; - Teuchos::broadcast(*commT, 0, numBdEdgeTags, bdEdgeTagsArray.data()); + Teuchos::broadcast(*comm, 0, numBdEdgeTags, bdEdgeTagsArray.data()); // Adding boundary nodesets and sidesets separating different labels for (int k=0; kuniversal_part()); #endif - Teuchos::broadcast(*commT, 0, &NumElemNodes); + Teuchos::broadcast(*comm, 0, &NumElemNodes); if(NumElemNodes == 3) { stk::mesh::set_topology(*partVec[0], stk::topology::TRI_3_2D); } @@ -278,7 +281,7 @@ Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCPget("Workset Size", DEFAULT_WORKSET_SIZE); - Teuchos::broadcast(*commT, 0, &NumElems); + Teuchos::broadcast(*comm, 0, &NumElems); int worksetSize = this->computeWorksetSize(worksetSizeMax, NumElems); @@ -286,43 +289,26 @@ Albany::AsciiSTKMesh2D::AsciiSTKMesh2D (const Teuchos::RCPmeshSpecs[0] = Teuchos::rcp ( - new Albany::MeshSpecsStruct (ctd, numDim, nsNames, ssNames, + new MeshSpecsStruct (ctd, numDim, nsNames, ssNames, worksetSize, partVec[0]->name(), ebNameToIndex)); // Create a mesh specs object for EACH side set - this->initializeSideSetMeshSpecs(commT); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(commT); -} - -Albany::AsciiSTKMesh2D::~AsciiSTKMesh2D() {} - -void Albany::AsciiSTKMesh2D::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& /* side_set_sis */) -{ - this->SetupFieldData(commT, sis, worksetSize); + this->initializeSideSetMeshSpecs(comm); } -void Albany::AsciiSTKMesh2D::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void AsciiSTKMesh2D:: +setBulkData (const Teuchos::RCP& comm) { metaData->commit(); bulkData->modification_begin(); // Begin modifying the mesh Teuchos::RCP out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); - out->setProcRankAndSize(commT->getRank(), commT->getSize()); + out->setProcRankAndSize(comm->getRank(), comm->getSize()); out->setOutputToRootOnly(0); // Only proc 0 has loaded the file - if (commT->getRank()==0) + if (comm->getRank()==0) { stk::mesh::PartVector singlePartVec(1); unsigned int ebNo = 0; //element block #??? @@ -363,7 +349,7 @@ void Albany::AsciiSTKMesh2D::setBulkData( if(proc_rank_field){ int* p_rank = stk::mesh::field_data(*proc_rank_field, elem); if(p_rank) - p_rank[0] = commT->getRank(); + p_rank[0] = comm->getRank(); } } *out << "done!\n"; @@ -434,19 +420,19 @@ void Albany::AsciiSTKMesh2D::setBulkData( params->set("Rebalance Mesh", true); // Rebalance the mesh before starting the simulation if indicated - rebalanceInitialMeshT(commT); + rebalanceInitialMesh(comm); #endif // Loading required input fields from file - this->loadRequiredInputFields (commT); - - // Finally, perform the setup of the (possible) side set meshes (including extraction if of type SideSetSTKMeshStruct) - this->setSideSetFieldAndBulkData(commT, side_set_sis, worksetSize); + this->loadRequiredInputFields (comm); - fieldAndBulkDataSet = true; + m_bulk_data_set = true; } -Teuchos::RCP Albany::AsciiSTKMesh2D::getValidDiscretizationParameters() const { +Teuchos::RCP +AsciiSTKMesh2D:: +getValidDiscretizationParameters() const +{ Teuchos::RCP validPL = this->getValidGenericSTKParameters("Valid ASCII_DiscParams"); validPL->set("Ascii Input Mesh File Name", "greenland.msh", @@ -456,3 +442,5 @@ Teuchos::RCP Albany::AsciiSTKMesh2D::getValidDiscr return validPL; } + +} // namespace Albany diff --git a/src/disc/stk/Albany_AsciiSTKMesh2D.hpp b/src/disc/stk/Albany_AsciiSTKMesh2D.hpp index a970ea9ba4..f317d8702c 100644 --- a/src/disc/stk/Albany_AsciiSTKMesh2D.hpp +++ b/src/disc/stk/Albany_AsciiSTKMesh2D.hpp @@ -53,22 +53,13 @@ namespace Albany { public: - AsciiSTKMesh2D( - const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams); - - ~AsciiSTKMesh2D(); - - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); - - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + AsciiSTKMesh2D (const Teuchos::RCP& params, + const Teuchos::RCP& comm, + const int numParams); + + ~AsciiSTKMesh2D() = default; + + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return false; } diff --git a/src/disc/stk/Albany_AsciiSTKMeshStruct.cpp b/src/disc/stk/Albany_AsciiSTKMeshStruct.cpp index 69c7be1cc9..d2d31f1349 100644 --- a/src/disc/stk/Albany_AsciiSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_AsciiSTKMeshStruct.cpp @@ -399,12 +399,8 @@ AsciiSTKMeshStruct(const Teuchos::RCP& params, nsNames, ssNames, worksetSize, ebn, ebNameToIndex)); - // Create a mesh specs object for EACH side set this->initializeSideSetMeshSpecs(comm); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(comm); } AsciiSTKMeshStruct::~AsciiSTKMeshStruct() @@ -415,23 +411,8 @@ AsciiSTKMeshStruct::~AsciiSTKMeshStruct() delete [] eles; } -void -AsciiSTKMeshStruct::setFieldData( - const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - this->SetupFieldData(comm, sis, worksetSize); - this->setSideSetFieldData(comm, side_set_sis, worksetSize); -} - -void -AsciiSTKMeshStruct::setBulkData( - const Teuchos::RCP& comm, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void AsciiSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { metaData->commit(); @@ -720,8 +701,7 @@ AsciiSTKMeshStruct::setBulkData( fix_node_sharing(*bulkData); bulkData->modification_end(); - fieldAndBulkDataSet = true; - this->setSideSetFieldAndBulkData(comm, side_set_sis, worksetSize); + m_bulk_data_set = true; } Teuchos::RCP diff --git a/src/disc/stk/Albany_AsciiSTKMeshStruct.hpp b/src/disc/stk/Albany_AsciiSTKMeshStruct.hpp index 7615a4f4f7..acf9751f46 100644 --- a/src/disc/stk/Albany_AsciiSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_AsciiSTKMeshStruct.hpp @@ -25,15 +25,7 @@ class AsciiSTKMeshStruct : public GenericSTKMeshStruct ~AsciiSTKMeshStruct(); - void setFieldData (const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); - - void setBulkData (const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return false; } diff --git a/src/disc/stk/Albany_BucketArray.hpp b/src/disc/stk/Albany_BucketArray.hpp deleted file mode 100644 index b03ed865e1..0000000000 --- a/src/disc/stk/Albany_BucketArray.hpp +++ /dev/null @@ -1,181 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#ifndef ALBANY_BUCKETARRAY_HPP -#define ALBANY_BUCKETARRAY_HPP - -#include -#include -#include -#include -#include -#include - -#include - -//#define IKT_DEBUG - -namespace Albany { - -struct EntityDimension : public shards::ArrayDimTag { - - const char * name() const - { static const char n[] = "EntityDimension" ; return n ; } - - static const EntityDimension & tag() ///< Singleton - { static const EntityDimension self ; return self ; } - -private: - EntityDimension() {} - EntityDimension( const EntityDimension & ); - EntityDimension & operator = ( const EntityDimension & ); -}; - -template< class FieldType > struct BucketArray {}; - -/** \brief \ref stk::mesh::Field "Field" data \ref shards::Array "Array" - * for a given scalar field and bucket - */ -template< typename ScalarType > -struct BucketArray< stk::mesh::Field > - : public -shards::Array -{ -private: - typedef unsigned char * byte_p ; - BucketArray(); - BucketArray( const BucketArray & ); - BucketArray & operator = ( const BucketArray & ); - -public: - - typedef stk::mesh::Field field_type ; - typedef - shards::Array - array_type ; - - BucketArray( const field_type & f , const stk::mesh::Bucket & k ) - { - if (k.field_data_is_allocated(f)) { - array_type::assign( (ScalarType*)( k.field_data_location(f) ) , - k.size() ); - - } - } -}; - -template -inline size_t -get_size(stk::mesh::Bucket const&) -{ - return Tag::Size; -} - -template <> -inline size_t -get_size(stk::mesh::Bucket const&) -{ - return 0; -} - -template <> -inline size_t -get_size(stk::mesh::Bucket const& b) -{ - return b.mesh().mesh_meta_data().spatial_dimension(); -} - -//---------------------------------------------------------------------- -/** \brief \ref stk::mesh::Field "Field" data \ref shards::Array "Array" - * for a given array field and bucket - */ -template< typename ScalarType , - class Tag1 , class Tag2 , class Tag3 , class Tag4 , - class Tag5 , class Tag6 , class Tag7 > -struct BucketArray< stk::mesh::Field > - : public shards::ArrayAppend< - shards::Array , - EntityDimension >::type -{ -private: - typedef unsigned char * byte_p ; - BucketArray(); - BucketArray( const BucketArray & ); - BucketArray & operator = ( const BucketArray & ); -public: - - typedef stk::mesh::Field field_type ; - - typedef typename shards::ArrayAppend< - shards::Array , - EntityDimension >::type array_type ; - - BucketArray( const field_type & f , const stk::mesh::Bucket & b ) - { - if ( b.field_data_is_allocated(f) ) { - int stride[4]; - if (f.field_array_rank() == 1) { - stride[0] = stk::mesh::field_scalars_per_entity(f, b); - } - else if (f.field_array_rank() == 2) { - int dim0 = stk::mesh::find_restriction(f, b.entity_rank(), b.supersets()).dimension(); - stride[0] = dim0; - stride[1] = stk::mesh::field_scalars_per_entity(f, b); - } - else if (f.field_array_rank() == 3) { -#ifdef IKT_DEBUG - std::cout << "IKT field array rank = 3" << std::endl; -#endif - int dim0 = stk::mesh::find_restriction(f, b.entity_rank(), b.supersets()).dimension(); -#ifdef IKT_DEBUG - std::cout << "IKT size tag1 = " << get_size(b) << std::endl; - std::cout << "IKT size tag2 = " << get_size(b) << std::endl; - std::cout << "IKT size tag3 = " << get_size(b) << std::endl; -#endif - if (dim0 == 4) { - stride[0] = dim0; - stride[1] = get_size(b) * dim0; - stride[2] = stk::mesh::field_scalars_per_entity(f, b); - } - else { - //IKT, 12/20/18: this changes the way the qp_tensor field - //for 1D and 3D problems appears in the output exodus field. - //Fields appear like: Cauchy_Stress_1_1, ... Cauchy_Stress_8_9, - //instead of Cauchy_Stress_1_01 .. Cauchy_Stress_3_24 to make it - //more clear which entry corresponds to which component/quad point. - //I believe for 2D problems the original layout is correct, hence - //the if statement above here. - stride[0] = get_size(b); - stride[1] = get_size(b) * stride[0]; - stride[2] = stk::mesh::field_scalars_per_entity(f, b); - } -#ifdef IKT_DEBUG - std::cout << "IKT stride0, stride1, stride2 = " << stride[0] << ", " << stride[1] << ", " << stride[2] << std::endl; -#endif - } - else if (f.field_array_rank() == 4) { - int dim0 = stk::mesh::find_restriction(f, b.entity_rank(), b.supersets()).dimension(); - stride[0] = dim0; - stride[1] = get_size(b) * dim0; - stride[2] = get_size(b) * stride[1]; - stride[3] = stk::mesh::field_scalars_per_entity(f, b); - } - else { - assert(false); - } - - array_type::assign_stride( - (ScalarType*)( b.field_data_location(f) ), - stride, - (typename array_type::size_type) b.size() ); - - } - } -}; - -} //namespace Albany - -#endif // ALBANY_BUCKETARRAY_HPP diff --git a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp index 790cff5861..bf0e0b89ad 100644 --- a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.cpp @@ -117,7 +117,7 @@ ExtrudedSTKMeshStruct(const Teuchos::RCP& params, sideSetMeshStructs["basalside"] = basalMeshStruct; - const auto& basalMeshSpec = basalMeshStruct->getMeshSpecs()[0]; + const auto& basalMeshSpec = basalMeshStruct->meshSpecs[0]; std::string elem2d_name(basalMeshSpec->ctd.base->name); std::string tria = shards::getCellTopologyData >()->name; std::string quad = shards::getCellTopologyData >()->name; @@ -187,49 +187,20 @@ ExtrudedSTKMeshStruct(const Teuchos::RCP& params, // Create a mesh specs object for EACH side set this->initializeSideSetMeshSpecs(comm); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(comm); } -void ExtrudedSTKMeshStruct::setFieldData( - const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - out->setProcRankAndSize(comm->getRank(), comm->getSize()); - out->setOutputToRootOnly(0); - - // Finish to set up the basal mesh - Teuchos::RCP dummy_sis = Teuchos::rcp(new StateInfoStruct()); - auto it_sis = side_set_sis.find("basalside"); - auto& basal_sis = (it_sis==side_set_sis.end() ? dummy_sis : it_sis->second); - - this->sideSetMeshStructs.at("basalside")->setFieldData (comm, basal_sis, worksetSize); - - // Setting up the field container - this->SetupFieldData(comm, sis, worksetSize); - - this->setSideSetFieldData(comm, side_set_sis, worksetSize); -} - -void ExtrudedSTKMeshStruct::setBulkData( - const Teuchos::RCP& comm, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void ExtrudedSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { constexpr auto ELEM_RANK = stk::topology::ELEM_RANK; constexpr auto NODE_RANK = stk::topology::NODE_RANK; const auto SIDE_RANK = metaData->side_rank(); - // Finish to set up the basal mesh - Teuchos::RCP dummy_sis = Teuchos::rcp(new StateInfoStruct()); - auto it_sis = side_set_sis.find("basalside"); - auto& basal_sis = (it_sis==side_set_sis.end() ? dummy_sis : it_sis->second); - - this->sideSetMeshStructs.at("basalside")->setBulkData (comm, basal_sis, worksetSize); + // This code should always run, but just in case something changes, + // let's first check that bulk data was not set yet + if (not basalMeshStruct->isBulkDataSet()) { + basalMeshStruct->setBulkData (comm); + } constexpr auto LAYER = LayeredMeshOrdering::LAYER; @@ -592,14 +563,10 @@ void ExtrudedSTKMeshStruct::setBulkData( //fix_node_sharing(*bulkData); bulkData->modification_end(); - fieldAndBulkDataSet = true; // Check that the nodeset created from sidesets contain the right number of nodes this->checkNodeSetsFromSideSetsIntegrity (); - // We can finally extract the side set meshes and set the fields and bulk data in all of them - this->setSideSetBulkData(comm, side_set_sis, worksetSize); - if (params->get("Export 2D Data",false)) { // We export the basal mesh in GMSH format std::ofstream ofile; @@ -661,6 +628,8 @@ void ExtrudedSTKMeshStruct::setBulkData( ofile.close(); } + + m_bulk_data_set = true; } void ExtrudedSTKMeshStruct:: diff --git a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp index 01a58ae7ee..68498a0dc9 100644 --- a/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_ExtrudedSTKMeshStruct.hpp @@ -26,15 +26,7 @@ class ExtrudedSTKMeshStruct : public GenericSTKMeshStruct ~ExtrudedSTKMeshStruct() = default; - void setFieldData (const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); // empty map as default - - void setBulkData (const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); // empty map as default + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const { return false; } diff --git a/src/disc/stk/Albany_GenericSTKFieldContainer.cpp b/src/disc/stk/Albany_GenericSTKFieldContainer.cpp index 7eeba2e570..d54aab1d1b 100644 --- a/src/disc/stk/Albany_GenericSTKFieldContainer.cpp +++ b/src/disc/stk/Albany_GenericSTKFieldContainer.cpp @@ -47,7 +47,9 @@ GenericSTKFieldContainer::GenericSTKFieldContainer( params(params_), neq(neq_), numDim(numDim_), - num_params(num_params_) { + num_params(num_params_) +{ + save_solution_field = params_->get("Save Solution Field", true); } #ifdef ALBANY_SEACAS diff --git a/src/disc/stk/Albany_GenericSTKFieldContainer.hpp b/src/disc/stk/Albany_GenericSTKFieldContainer.hpp index 932092e072..a9749cc959 100644 --- a/src/disc/stk/Albany_GenericSTKFieldContainer.hpp +++ b/src/disc/stk/Albany_GenericSTKFieldContainer.hpp @@ -60,6 +60,8 @@ class GenericSTKFieldContainer : public AbstractSTKFieldContainer int neq; int numDim; int num_params{0}; + + bool save_solution_field = false; }; } // namespace Albany diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp index bebf2054c9..8d4051182e 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.cpp @@ -4,15 +4,18 @@ // in the file "license.txt" in the top-level Albany directory // //*****************************************************************// -#include -#include "Teuchos_VerboseObject.hpp" -#include "Teuchos_RCPStdSharedPtrConversions.hpp" - #include "Albany_DiscretizationFactory.hpp" #include "Albany_GenericSTKMeshStruct.hpp" #include "Albany_SideSetSTKMeshStruct.hpp" +#include +#include +#include +#include + #include "Albany_KokkosTypes.hpp" #include "Albany_Gather.hpp" +#include "Albany_CommUtils.hpp" +#include "Albany_Utils.hpp" #include "Albany_OrdinarySTKFieldContainer.hpp" #include "Albany_MultiSTKFieldContainer.hpp" @@ -21,16 +24,10 @@ #include #endif -#include "Albany_Utils.hpp" #include #include #include -#include -#include -#include -#include - // Expression reading #ifdef ALBANY_PANZER_EXPR_EVAL #include @@ -44,6 +41,11 @@ #include #endif +#include "Teuchos_VerboseObject.hpp" +#include "Teuchos_RCPStdSharedPtrConversions.hpp" + +#include + namespace Albany { @@ -55,6 +57,7 @@ GenericSTKMeshStruct::GenericSTKMeshStruct( num_params(numParams_) { metaData = Teuchos::rcp(new stk::mesh::MetaData()); + metaData->use_simple_fields(); // numDim = -1 is default flag value to postpone initialization if (numDim_>0) { @@ -64,38 +67,30 @@ GenericSTKMeshStruct::GenericSTKMeshStruct( } allElementBlocksHaveSamePhysics = true; - compositeTet = params->get("Use Composite Tet 10", false); num_time_deriv = params->get("Number Of Time Derivatives"); requiresAutomaticAura = params->get("Use Automatic Aura", false); // This is typical, can be resized for multiple material problems meshSpecs.resize(1); - - fieldAndBulkDataSet = false; - side_maps_present = false; - ignore_side_maps = false; } -void GenericSTKMeshStruct::SetupFieldData( - const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const int worksetSize) +void GenericSTKMeshStruct:: +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) { - TEUCHOS_TEST_FOR_EXCEPTION(!metaData->is_initialized(), - std::logic_error, - "LogicError: metaData->FEM_initialize(numDim) not yet called" << std::endl); - + TEUCHOS_TEST_FOR_EXCEPTION(!metaData->is_initialized(), std::logic_error, + "[GenericSTKMeshStruct::SetupFieldData] metaData->initialize(numDim) not yet called" << std::endl); if (bulkData.is_null()) { - const Teuchos::MpiComm* mpiComm = dynamic_cast* > (comm.get()); - stk::mesh::MeshBuilder meshBuilder = stk::mesh::MeshBuilder(*mpiComm->getRawMpiComm()); + auto mpiComm = getMpiCommFromTeuchosComm(comm); + stk::mesh::MeshBuilder meshBuilder = stk::mesh::MeshBuilder(mpiComm); if(requiresAutomaticAura) meshBuilder.set_aura_option(stk::mesh::BulkData::AUTO_AURA); else meshBuilder.set_aura_option(stk::mesh::BulkData::NO_AUTO_AURA); - meshBuilder.set_bucket_capacity(worksetSize); + meshBuilder.set_bucket_capacity(meshSpecs[0]->worksetSize); meshBuilder.set_add_fmwk_data(false); std::unique_ptr bulkDataPtr = meshBuilder.create(Teuchos::get_shared_ptr(metaData)); bulkData = Teuchos::rcp(bulkDataPtr.release()); @@ -133,10 +128,10 @@ void GenericSTKMeshStruct::SetupFieldData( // Build the usual Albany fields unless the user explicitly specifies the residual or solution vector layout if(user_specified_solution_components && (residual_vector.length() > 0)){ this->fieldContainer = Teuchos::rcp(new MultiSTKFieldContainer(params, - metaData, bulkData, numDim, sis, num_params)); + metaData, bulkData, numDim, num_params)); } else { this->fieldContainer = Teuchos::rcp(new OrdinarySTKFieldContainer(params, - metaData, bulkData, numDim, sis, num_params)); + metaData, bulkData, numDim, num_params)); } // Exodus is only for 2D and 3D. Have 1D version as well @@ -237,24 +232,6 @@ This function gets rid of the subset in the list. } } -Teuchos::ArrayRCP >& -GenericSTKMeshStruct::getMeshSpecs() -{ - TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs==Teuchos::null, - std::logic_error, - "meshSpecs accessed, but it has not been constructed" << std::endl); - return meshSpecs; -} - -const Teuchos::ArrayRCP >& -GenericSTKMeshStruct::getMeshSpecs() const -{ - TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs==Teuchos::null, - std::logic_error, - "meshSpecs accessed, but it has not been constructed" << std::endl); - return meshSpecs; -} - int GenericSTKMeshStruct::computeWorksetSize(const int worksetSizeMax, const int ebSizeMax) const { @@ -292,20 +269,20 @@ void GenericSTKMeshStruct::setDefaultCoordinates3d () } } -void GenericSTKMeshStruct::rebalanceInitialMeshT(const Teuchos::RCP >& comm){ +void GenericSTKMeshStruct::rebalanceInitialMesh (const Teuchos::RCP& comm){ bool rebalance = params->get("Rebalance Mesh", false); if(rebalance) { TEUCHOS_TEST_FOR_EXCEPTION (this->side_maps_present, std::runtime_error, "Error! Rebalance is not supported when side maps are present.\n"); - rebalanceAdaptedMeshT(params, comm); + rebalanceAdaptedMesh (params, comm); } } void GenericSTKMeshStruct:: -rebalanceAdaptedMeshT(const Teuchos::RCP& params_, - const Teuchos::RCP >& comm) +rebalanceAdaptedMesh (const Teuchos::RCP& params_, + const Teuchos::RCP& comm) { // Zoltan is required here #ifdef ALBANY_ZOLTAN @@ -439,7 +416,7 @@ void GenericSTKMeshStruct::checkNodeSetsFromSideSetsIntegrity () void GenericSTKMeshStruct::initializeSideSetMeshSpecs (const Teuchos::RCP& comm) { // Loop on all mesh specs - for (auto ms: this->getMeshSpecs() ) { + for (auto ms : meshSpecs ) { // Loop on all side sets of the mesh for (auto ssName : ms->ssNames) { // Get the part @@ -480,147 +457,26 @@ void GenericSTKMeshStruct::initializeSideSetMeshSpecs (const Teuchos::RCP& comm) +void GenericSTKMeshStruct::createSideMeshMaps () { - if (params->isSublist ("Side Set Discretizations")) { - const Teuchos::ParameterList& ssd_list = params->sublist("Side Set Discretizations"); - const Teuchos::Array& sideSets = ssd_list.get >("Side Sets"); - - Teuchos::RCP params_ss; - int sideDim = this->numDim - 1; - for (int i(0); i ss_mesh; - const std::string& ss_name = sideSets[i]; - params_ss = Teuchos::rcp(new Teuchos::ParameterList(ssd_list.sublist(ss_name))); - - // We must check whether a side mesh was already created elsewhere. - // If the mesh already exists, we do nothing, and we ASSUME it is a valid mesh - // This happens, for instance, for the basal mesh for extruded meshes. - if (this->sideSetMeshStructs.find(ss_name)==this->sideSetMeshStructs.end()) { - if (!params_ss->isParameter("Number Of Time Derivatives")) - params_ss->set("Number Of Time Derivatives",num_time_deriv); - - // Set sideset discretization workset size based on sideset mesh spec if a single workset is used - const auto &sideSetMeshSpecs = this->meshSpecs[0]->sideSetMeshSpecs; - auto sideSetMeshSpecIter = sideSetMeshSpecs.find(ss_name); - TEUCHOS_TEST_FOR_EXCEPTION(sideSetMeshSpecIter == sideSetMeshSpecs.end(), std::runtime_error, - "Cannot find " << ss_name << " in sideSetMeshSpecs!\n"); - - std::string method = params_ss->get("Method"); - if (method=="SideSetSTK") { - // The user said this mesh is extracted from a higher dimensional one - TEUCHOS_TEST_FOR_EXCEPTION (meshSpecs.size()!=1, std::logic_error, - "Error! So far, side set mesh extraction is allowed only from STK meshes with 1 element block.\n"); - - this->sideSetMeshStructs[ss_name] = - Teuchos::rcp(new SideSetSTKMeshStruct( - *this->meshSpecs[0], params_ss, comm, num_params)); - } else { - ss_mesh = DiscretizationFactory::createMeshStruct (params_ss,comm, num_params); - this->sideSetMeshStructs[ss_name] = Teuchos::rcp_dynamic_cast(ss_mesh,false); - TEUCHOS_TEST_FOR_EXCEPTION (this->sideSetMeshStructs[ss_name]==Teuchos::null, std::runtime_error, - "Error! Could not cast side mesh to AbstractSTKMeshStruct.\n"); - } - } - - // Checking that the side meshes have the correct dimension (in case they were loaded from file, - // and the user mistakenly gave the wrong file name) - TEUCHOS_TEST_FOR_EXCEPTION (sideDim!=this->sideSetMeshStructs[ss_name]->numDim, std::logic_error, - "Error! Mesh on side " << ss_name << " has the wrong dimension.\n"); - - // Update the side set mesh specs pointer in the mesh specs of this mesh - this->meshSpecs[0]->sideSetMeshSpecs[ss_name] = this->sideSetMeshStructs[ss_name]->getMeshSpecs(); - this->meshSpecs[0]->sideSetMeshNames.push_back(ss_name); - - // We need to create the 2D cell -> (3D cell, side_node_ids) map in the side mesh now - using ISFT = AbstractSTKFieldContainer::STKIntState; - ISFT* side_to_cell_map = &this->sideSetMeshStructs[ss_name]->metaData->declare_field (stk::topology::ELEM_RANK, "side_to_cell_map"); - stk::mesh::put_field_on_mesh(*side_to_cell_map, this->sideSetMeshStructs[ss_name]->metaData->universal_part(), 1, nullptr); + for (auto& it : sideSetMeshStructs) { + auto ss_stk_mesh = Teuchos::rcp_dynamic_cast(it.second,true); + auto ss_ms = ss_stk_mesh->meshSpecs[0]; + + // We need to create the 2D cell -> (3D cell, side_node_ids) map in the side mesh now + using ISFT = AbstractSTKFieldContainer::STKIntState; + ISFT* side_to_cell_map = &ss_stk_mesh->metaData->declare_field (stk::topology::ELEM_RANK, "side_to_cell_map"); + stk::mesh::put_field_on_mesh(*side_to_cell_map, ss_stk_mesh->metaData->universal_part(), 1, nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*side_to_cell_map, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*side_to_cell_map, Ioss::Field::TRANSIENT); #endif - // We need to create the 2D cell -> (3D cell, side_node_ids) map in the side mesh now - const int num_nodes = sideSetMeshStructs[ss_name]->getMeshSpecs()[0]->ctd.node_count; - ISFT* side_nodes_ids = &this->sideSetMeshStructs[ss_name]->metaData->declare_field (stk::topology::ELEM_RANK, "side_nodes_ids"); - stk::mesh::put_field_on_mesh(*side_nodes_ids, this->sideSetMeshStructs[ss_name]->metaData->universal_part(), num_nodes, nullptr); + // We need to create the 2D cell -> (3D cell, side_node_ids) map in the side mesh now + const int num_nodes = ss_ms->ctd.node_count; + ISFT* side_nodes_ids = &ss_stk_mesh->metaData->declare_field (stk::topology::ELEM_RANK, "side_nodes_ids"); + stk::mesh::put_field_on_mesh(*side_nodes_ids, ss_stk_mesh->metaData->universal_part(), num_nodes, nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*side_nodes_ids, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*side_nodes_ids, Ioss::Field::TRANSIENT); #endif - - // If requested, we ignore the side maps already stored in the imported side mesh (if any) - // This can be useful for side mesh of an extruded mesh, in the case it was constructed - // as side mesh of an extruded mesh with a different ordering and/or different number - // of layers. Notice that if that's the case, it probably is impossible to build a new - // set of maps, since there is no way to correctly map the side nodes to the cell nodes. - this->sideSetMeshStructs[ss_name]->ignore_side_maps = params_ss->get("Ignore Side Maps", false); - } - } -} - -void GenericSTKMeshStruct::setSideSetFieldData ( - const Teuchos::RCP& comm, - const std::map >& side_set_sis, - int worksetSize) -{ - if (this->sideSetMeshStructs.size()>0) { - // Dummy sis if not present in the maps for a given side set. - // This could happen if the side discretization has no states - Teuchos::RCP dummy_sis = Teuchos::rcp(new StateInfoStruct()); - - Teuchos::RCP params_ss; - const Teuchos::ParameterList& ssd_list = params->sublist("Side Set Discretizations"); - for (auto it : sideSetMeshStructs) { - Teuchos::RCP sideMesh; - sideMesh = Teuchos::rcp_dynamic_cast(it.second,false); - if (sideMesh!=Teuchos::null) { - // SideSetSTK mesh need to build the mesh - sideMesh->setParentMeshInfo(*this, it.first); - } - - // We check since the basal mesh for extruded stk mesh should already have it set - if (!it.second->fieldAndBulkDataSet) { - auto it_sis = side_set_sis.find(it.first); - - auto& sis = (it_sis==side_set_sis.end() ? dummy_sis : it_sis->second); - - params_ss = Teuchos::rcp(new Teuchos::ParameterList(ssd_list.sublist(it.first))); - it.second->setFieldData(comm,sis,worksetSize); // Cell equations are also defined on the side, but not vice-versa - } - } - } -} - -void GenericSTKMeshStruct::setSideSetBulkData ( - const Teuchos::RCP& comm, - const std::map >& side_set_sis, - int worksetSize) -{ - if (this->sideSetMeshStructs.size()>0) { - // Dummy sis if not present in the maps for a given side set. - // This could happen if the side discretization has no states - Teuchos::RCP dummy_sis = Teuchos::rcp(new StateInfoStruct()); - - Teuchos::RCP params_ss; - const Teuchos::ParameterList& ssd_list = params->sublist("Side Set Discretizations"); - for (auto it : sideSetMeshStructs) { - Teuchos::RCP sideMesh; - sideMesh = Teuchos::rcp_dynamic_cast(it.second,false); - if (sideMesh!=Teuchos::null) { - // SideSetSTK mesh need to build the mesh - sideMesh->setParentMeshInfo(*this, it.first); - } - - // We check since the basal mesh for extruded stk mesh should already have it set - if (!it.second->fieldAndBulkDataSet) { - auto it_sis = side_set_sis.find(it.first); - - auto& sis = (it_sis==side_set_sis.end() ? dummy_sis : it_sis->second); - - params_ss = Teuchos::rcp(new Teuchos::ParameterList(ssd_list.sublist(it.first))); - it.second->setBulkData(comm,sis,worksetSize); // Cell equations are also defined on the side, but not vice-versa - } - } } } @@ -629,10 +485,10 @@ buildCellSideNodeNumerationMap (const std::string& sideSetName, std::map& sideMap, std::map>& sideNodeMap) { - TEUCHOS_TEST_FOR_EXCEPTION (sideSetMeshStructs.find(sideSetName)==sideSetMeshStructs.end(), Teuchos::Exceptions::InvalidParameter, - "Error in 'buildSideNodeToSideSetCellNodeMap': side set " << sideSetName << " does not have a mesh.\n"); + TEUCHOS_TEST_FOR_EXCEPTION (sideSetMeshStructs.count(sideSetName)==0, Teuchos::Exceptions::InvalidParameter, + "Error in 'buildCellSideNodeNumerationMap': side set " << sideSetName << " does not have a mesh.\n"); - Teuchos::RCP side_mesh = sideSetMeshStructs.at(sideSetName); + auto side_mesh = Teuchos::rcp_dynamic_cast(sideSetMeshStructs.at(sideSetName)); // NOTE 1: the stk fields memorize maps from 2D to 3D values (since fields are defined on 2D mesh), while the input // maps map 3D values to 2D ones. For instance, if node 0 of side3D 41 is mapped to node 2 of cell2D 12, the @@ -660,8 +516,8 @@ buildCellSideNodeNumerationMap (const std::string& sideSetName, int side_lid; int num_sides; using ISFT = AbstractSTKFieldContainer::STKIntState; - ISFT* side_to_cell_map = this->sideSetMeshStructs[sideSetName]->metaData->get_field (stk::topology::ELEM_RANK, "side_to_cell_map"); - ISFT* side_nodes_ids_map = this->sideSetMeshStructs[sideSetName]->metaData->get_field (stk::topology::ELEM_RANK, "side_nodes_ids"); + ISFT* side_to_cell_map = side_mesh->metaData->get_field (stk::topology::ELEM_RANK, "side_to_cell_map"); + ISFT* side_nodes_ids_map = side_mesh->metaData->get_field (stk::topology::ELEM_RANK, "side_nodes_ids"); std::vector cell2D_nodes_ids(num_nodes), side3D_nodes_ids(num_nodes); const stk::mesh::Entity* side3D_nodes; const stk::mesh::Entity* cell2D_nodes; @@ -1441,32 +1297,13 @@ void GenericSTKMeshStruct::checkFieldIsInMesh (const std::string& fname, const s entity_rank = stk::topology::NODE_RANK; } - int dim = 1; - if (ftype.find("Vector")!=std::string::npos) { - ++dim; - } - if (ftype.find("Layered")!=std::string::npos) { - ++dim; - } - bool missing = (metaData->get_field (entity_rank, fname)==nullptr); if (missing) { - bool isFieldInMesh = false; - auto fl = metaData->get_fields(); - auto f = fl.begin(); - for (; f != fl.end(); ++f) { - isFieldInMesh = (fname == (*f)->name()); - if(isFieldInMesh) break; - } - if(isFieldInMesh) { - TEUCHOS_TEST_FOR_EXCEPTION (missing, std::runtime_error, "Error! The field '" << fname << "' in the mesh has different rank or dimensions than the ones specified\n" - << " Rank required: " << entity_rank << ", rank of field in mesh: " << (*f)->entity_rank() << "\n" - << " Dimension required: " << dim << ", dimension of field in mesh: " << (*f)->field_array_rank()+1 << "\n"); - } - else - TEUCHOS_TEST_FOR_EXCEPTION (missing, std::runtime_error, "Error! The field '" << fname << "' was not found in the mesh.\n" - << " Probably it was not registered it in the state manager (which forwards it to the mesh)\n"); + TEUCHOS_TEST_FOR_EXCEPTION ( + missing, std::runtime_error, + "Error! The field '" << fname << "' was not found in the mesh.\n" + " Probably it was not registered it in the state manager (which forwards it to the mesh)\n"); } } @@ -1516,6 +1353,7 @@ GenericSTKMeshStruct::getValidGenericSTKParameters(std::string listname) const validPL->set>("Betas BL Transform", Teuchos::tuple(0.0, 0.0, 0.0), "Beta parameters for Tanh Boundary Layer transform type"); validPL->set("Contiguous IDs", "true", "Tells Ascii mesh reader is mesh has contiguous global IDs on 1 processor."); //for LandIce problem that require transformation of STK mesh + validPL->set("Save Solution Field", "true", "Whether the solution field should be saved in the output mesh"); Teuchos::Array defaultFields; validPL->set >("Restart Fields", defaultFields, @@ -1532,7 +1370,6 @@ GenericSTKMeshStruct::getValidGenericSTKParameters(std::string listname) const validPL->set("Use Serial Mesh", false, "Read in a single mesh on PE 0 and rebalance"); validPL->set("Transfer Solution to Coordinates", false, "Copies the solution vector to the coordinates for output"); validPL->set("Set All Parts IO", false, "If true, all parts are marked as io parts"); - validPL->set("Use Composite Tet 10", false, "Flag to use the composite tet 10 basis in Intrepid"); validPL->set("Build Node Sets From Side Sets",false,"Flag to build node sets from side sets"); validPL->set("Export 3d coordinates field",false,"If true AND the mesh dimension is not already 3, export a 3d version of the coordinate field."); diff --git a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp index d93487471c..30678164ad 100644 --- a/src/disc/stk/Albany_GenericSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_GenericSTKMeshStruct.hpp @@ -19,14 +19,15 @@ class CombineAndScatterManager; class GenericSTKMeshStruct : public AbstractSTKMeshStruct { public: - Teuchos::ArrayRCP >& getMeshSpecs(); - const Teuchos::ArrayRCP >& getMeshSpecs() const; + GenericSTKMeshStruct( + const Teuchos::RCP& params, + const int numDim /*old default: -1*/, const int numParams); - //! Re-load balance adapted mesh - void rebalanceAdaptedMeshT(const Teuchos::RCP& params, - const Teuchos::RCP >& comm); + virtual ~GenericSTKMeshStruct() = default; - bool useCompositeTet(){ return compositeTet; } + //! Re-load balance adapted mesh + void rebalanceAdaptedMesh (const Teuchos::RCP& params, + const Teuchos::RCP& comm); // This routine builds two maps: side3D_id->cell2D_id, and side3D_node_lid->cell2D_node_lid. // These maps are used because the side id may differ from the cell id and the nodes order @@ -37,17 +38,9 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct std::map>& sideNodeMap); int getNumParams() const {return num_params; } -protected: - GenericSTKMeshStruct( - const Teuchos::RCP& params, - const int numDim /*old default: -1*/, const int numParams); - - virtual ~GenericSTKMeshStruct() = default; - void SetupFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const int worksetSize_); + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis); void printParts(stk::mesh::MetaData *metaData); @@ -58,7 +51,7 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct int computeWorksetSize(const int worksetSizeMax, const int ebSizeMax) const; //! Re-load balance mesh - void rebalanceInitialMeshT(const Teuchos::RCP >& comm); + void rebalanceInitialMesh (const Teuchos::RCP& comm); //! Sets all mesh parts as IO parts (will be written to file) void setAllPartsIO(); @@ -70,40 +63,19 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct void checkNodeSetsFromSideSetsIntegrity (); //! Creates empty mesh structs if required (and not already present) - void initializeSideSetMeshSpecs (const Teuchos::RCP& commT); + void initializeSideSetMeshSpecs (const Teuchos::RCP& comm); - //! Creates empty mesh structs if required (and not already present) - void initializeSideSetMeshStructs (const Teuchos::RCP& commT); - - //! Completes the creation of the side set mesh structs (if of type SideSetSTKMeshStruct) - void setSideSetFieldData( - const Teuchos::RCP& commT, - const std::map >& side_set_sis, - int worksetSize); - - void setSideSetBulkData( - const Teuchos::RCP& commT, - const std::map >& side_set_sis, - int worksetSize); - - void setSideSetFieldAndBulkData( - const Teuchos::RCP& commT, - const std::map >& side_set_sis, - int worksetSize) - { - setSideSetFieldData(commT, side_set_sis, worksetSize); - setSideSetBulkData(commT, side_set_sis, worksetSize); - } + void createSideMeshMaps (); //! Loads from file input required fields not found in the mesh - void loadRequiredInputFields (const Teuchos::RCP& commT); + void loadRequiredInputFields (const Teuchos::RCP& comm); // Routines to load, fill, or compute a field void loadField (const std::string& field_name, const Teuchos::ParameterList& params, Teuchos::RCP& field_mv, const CombineAndScatterManager& cas_manager, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, bool node, bool scalar, bool layered, const Teuchos::RCP out); void fillField (const std::string& field_name, @@ -152,12 +124,8 @@ class GenericSTKMeshStruct : public AbstractSTKMeshStruct Teuchos::RCP params; - Teuchos::ArrayRCP > meshSpecs; - bool requiresAutomaticAura; - bool compositeTet; - std::vector m_nodesets_from_sidesets; int num_params; diff --git a/src/disc/stk/Albany_GmshSTKMeshStruct.cpp b/src/disc/stk/Albany_GmshSTKMeshStruct.cpp index 1b082ba175..59a969ddd5 100644 --- a/src/disc/stk/Albany_GmshSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_GmshSTKMeshStruct.cpp @@ -30,10 +30,14 @@ #include "Albany_Utils.hpp" -Albany::GmshSTKMeshStruct::GmshSTKMeshStruct (const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams) : - GenericSTKMeshStruct (params, -1, numParams) +namespace Albany +{ + +GmshSTKMeshStruct:: +GmshSTKMeshStruct (const Teuchos::RCP& params, + const Teuchos::RCP& comm, + const int numParams) + : GenericSTKMeshStruct (params, -1, numParams) { fname = params->get("Gmsh Input Mesh File Name", "mesh.msh"); @@ -43,7 +47,7 @@ Albany::GmshSTKMeshStruct::GmshSTKMeshStruct (const Teuchos::RCPgetRank() == 0) + if (comm->getRank() == 0) { bool legacy = false; bool binary = false; @@ -51,25 +55,18 @@ Albany::GmshSTKMeshStruct::GmshSTKMeshStruct (const Teuchos::RCP nsNames; std::vector < std::string > ssNames; - set_boundaries( commT, ssNames, nsNames); + set_boundaries( comm, ssNames, nsNames); stk::topology etopology = this->get_topology(); int numEB = 0; - const stk::mesh::PartVector & all_parts = metaData->get_parts(); - for (stk::mesh::PartVector::const_iterator i = all_parts.begin(); - i != all_parts.end(); ++i) { - - stk::mesh::Part * const part = *i ; - + for (auto part : metaData->get_parts()) { if (!stk::mesh::is_auto_declared_part(*part)) { if ( part->primary_entity_rank() == stk::topology::ELEMENT_RANK) { - partVec.push_back(part); numEB++; } @@ -113,8 +104,7 @@ Albany::GmshSTKMeshStruct::GmshSTKMeshStruct (const Teuchos::RCPdeclare_part_with_topology(ebn, etopology)); shards::CellTopology shards_ctd = stk::mesh::get_cell_topology(etopology); @@ -130,18 +120,15 @@ Albany::GmshSTKMeshStruct::GmshSTKMeshStruct (const Teuchos::RCPmeshSpecs[0] = Teuchos::rcp ( - new Albany::MeshSpecsStruct (ctd, numDim, nsNames, ssNames, + new MeshSpecsStruct (ctd, numDim, nsNames, ssNames, worksetSize, partVec[0]->name(), ebNameToIndex)); // Create a mesh specs object for EACH side set - this->initializeSideSetMeshSpecs(commT); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(commT); + this->initializeSideSetMeshSpecs(comm); } -Albany::GmshSTKMeshStruct::~GmshSTKMeshStruct() +GmshSTKMeshStruct::~GmshSTKMeshStruct() { delete[] pts; @@ -198,7 +185,7 @@ Albany::GmshSTKMeshStruct::~GmshSTKMeshStruct() allowable_gmsh_versions.clear(); } -void Albany::GmshSTKMeshStruct::determine_file_type( bool& legacy, bool& binary, bool& ascii) +void GmshSTKMeshStruct::determine_file_type( bool& legacy, bool& binary, bool& ascii) { std::ifstream ifile; open_fname( ifile); @@ -222,10 +209,9 @@ void Albany::GmshSTKMeshStruct::determine_file_type( bool& legacy, bool& binary, } ifile.close(); - return; } -void Albany::GmshSTKMeshStruct::init_counters_to_zero() +void GmshSTKMeshStruct::init_counters_to_zero() { NumSides = 0; NumNodes = 0; @@ -238,11 +224,9 @@ void Albany::GmshSTKMeshStruct::init_counters_to_zero() nb_tri6 = 0; nb_lines = 0; nb_line3 = 0; - - return; } -void Albany::GmshSTKMeshStruct::init_pointers_to_null() +void GmshSTKMeshStruct::init_pointers_to_null() { pts = nullptr; tetra = nullptr; @@ -253,42 +237,26 @@ void Albany::GmshSTKMeshStruct::init_pointers_to_null() quads = nullptr; lines = nullptr; line3 = nullptr; - return; } -void Albany::GmshSTKMeshStruct::broadcast_topology( const Teuchos::RCP& commT) +void GmshSTKMeshStruct::broadcast_topology( const Teuchos::RCP& comm) { - Teuchos::broadcast(*commT, 0, 1, &this->numDim); - Teuchos::broadcast(*commT, 0, 1, &NumElemNodes); - Teuchos::broadcast(*commT, 0, 1, &NumSideNodes); - Teuchos::broadcast(*commT, 0, 1, &NumElems); - Teuchos::broadcast(*commT, 0, 1, &version_in); - - return; + Teuchos::broadcast(*comm, 0, 1, &this->numDim); + Teuchos::broadcast(*comm, 0, 1, &NumElemNodes); + Teuchos::broadcast(*comm, 0, 1, &NumSideNodes); + Teuchos::broadcast(*comm, 0, 1, &NumElems); + Teuchos::broadcast(*comm, 0, 1, &version_in); } -void Albany::GmshSTKMeshStruct::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - this->SetupFieldData(commT, sis, worksetSize); - this->setSideSetFieldData(commT, side_set_sis, worksetSize); -} - -void Albany::GmshSTKMeshStruct::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void GmshSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { metaData->commit(); bulkData->modification_begin(); // Begin modifying the mesh // Only proc 0 has loaded the file - if (commT->getRank()==0) { + if (comm->getRank()==0) { stk::mesh::PartVector singlePartVec(1); unsigned int ebNo = 0; //element block #??? @@ -322,7 +290,7 @@ void Albany::GmshSTKMeshStruct::setBulkData( if(proc_rank_field){ int* p_rank = stk::mesh::field_data(*proc_rank_field, elem); if(p_rank) - p_rank[0] = commT->getRank(); + p_rank[0] = comm->getRank(); } } @@ -375,9 +343,9 @@ void Albany::GmshSTKMeshStruct::setBulkData( params->set("Rebalance Mesh", true); // Rebalance the mesh before starting the simulation if indicated - rebalanceInitialMeshT(commT); + rebalanceInitialMesh (comm); #else - TEUCHOS_TEST_FOR_EXCEPTION(commT->getSize()>1, std::runtime_error, + TEUCHOS_TEST_FOR_EXCEPTION(comm->getSize()>1, std::runtime_error, "Error! If you use a Gmsh mesh, you need to either have ALBANY_ZOLTAN on, or run serially.\n"); #endif @@ -385,15 +353,12 @@ void Albany::GmshSTKMeshStruct::setBulkData( this->setDefaultCoordinates3d(); // Loading required input fields from file - this->loadRequiredInputFields (commT); - - // Finally, perform the setup of the (possible) side set meshes (including extraction if of type SideSetSTKMeshStruct) - this->setSideSetBulkData(commT, side_set_sis, worksetSize); + this->loadRequiredInputFields (comm); - fieldAndBulkDataSet = true; + m_bulk_data_set = true; } -Teuchos::RCP Albany::GmshSTKMeshStruct::getValidDiscretizationParameters() const +Teuchos::RCP GmshSTKMeshStruct::getValidDiscretizationParameters() const { Teuchos::RCP validPL = this->getValidGenericSTKParameters("Valid ASCII_DiscParams"); validPL->set("Gmsh Input Mesh File Name", "mesh.msh", @@ -404,7 +369,7 @@ Teuchos::RCP Albany::GmshSTKMeshStruct::getValidDi // -------------------------------- Read methods ---------------------------- // -void Albany::GmshSTKMeshStruct::loadLegacyMesh () +void GmshSTKMeshStruct::loadLegacyMesh () { std::ifstream ifile; open_fname( ifile); @@ -565,16 +530,14 @@ void Albany::GmshSTKMeshStruct::loadLegacyMesh () ifile.close(); } -void Albany::GmshSTKMeshStruct::swallow_lines_until( std::ifstream& ifile, std::string& line, std::string line_of_interest) +void GmshSTKMeshStruct::swallow_lines_until( std::ifstream& ifile, std::string& line, std::string line_of_interest) { while (std::getline (ifile, line) && line != line_of_interest) { // Keep swallowing lines... } - - return; } -void Albany::GmshSTKMeshStruct::set_NumNodes( std::ifstream& ifile) +void GmshSTKMeshStruct::set_NumNodes( std::ifstream& ifile) { std::string line; swallow_lines_until( ifile, line, "$Nodes"); @@ -597,10 +560,9 @@ void Albany::GmshSTKMeshStruct::set_NumNodes( std::ifstream& ifile) } TEUCHOS_TEST_FOR_EXCEPTION (NumNodes<=0, Teuchos::Exceptions::InvalidParameter, "Error! Invalid number of nodes.\n"); - return; } -void Albany::GmshSTKMeshStruct::load_node_data( std::ifstream& ifile) +void GmshSTKMeshStruct::load_node_data( std::ifstream& ifile) { if( version == GmshVersion::V2_2) { @@ -639,11 +601,9 @@ void Albany::GmshSTKMeshStruct::load_node_data( std::ifstream& ifile) delete[] node_IDs; } } - - return; } -void Albany::GmshSTKMeshStruct::set_num_entities( std::ifstream& ifile) +void GmshSTKMeshStruct::set_num_entities( std::ifstream& ifile) { std::string line; ifile.seekg (0, std::ios::beg); @@ -671,11 +631,9 @@ void Albany::GmshSTKMeshStruct::set_num_entities( std::ifstream& ifile) } TEUCHOS_TEST_FOR_EXCEPTION (num_entities<=0, Teuchos::Exceptions::InvalidParameter, "Error! Invalid number of mesh elements.\n"); - - return; } -void Albany::GmshSTKMeshStruct::increment_element_type( int e_type) +void GmshSTKMeshStruct::increment_element_type( int e_type) { switch (e_type) { @@ -692,11 +650,9 @@ void Albany::GmshSTKMeshStruct::increment_element_type( int e_type) TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, "Error! Element type (" << e_type << ") not supported.\n"); } - - return; } -void Albany::GmshSTKMeshStruct::set_specific_num_of_each_elements( std::ifstream& ifile) +void GmshSTKMeshStruct::set_specific_num_of_each_elements( std::ifstream& ifile) { // Gmsh lists elements and sides (and some points) all together, // and does not specify beforehand what kind of elements @@ -795,11 +751,9 @@ void Albany::GmshSTKMeshStruct::set_specific_num_of_each_elements( std::ifstream "Error! Could not determine if mesh was first or second order.\n" << "Checked for number of 2pt lines and 3pt lines, both are non-zero. \n") } - - return; } -void Albany::GmshSTKMeshStruct::size_all_element_pointers() +void GmshSTKMeshStruct::size_all_element_pointers() { // First values are the node IDs of the element, then tag lines = new int*[3]; @@ -842,14 +796,11 @@ void Albany::GmshSTKMeshStruct::size_all_element_pointers() { line3[i] = new int[nb_line3]; } - - return; } -void Albany::GmshSTKMeshStruct::set_generic_mesh_info() +void GmshSTKMeshStruct::set_generic_mesh_info() { - if (nb_tetra>0) - { + if (nb_tetra>0) { this->numDim = 3; NumElems = nb_tetra; @@ -858,9 +809,7 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 3; elems = tetra; sides = trias; - } - else if (nb_tet10>0) - { + } else if (nb_tet10>0) { this->numDim = 3; NumElems = nb_tet10; @@ -869,9 +818,7 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 6; elems = tet10; sides = tri6; - } - else if (nb_hexas>0) - { + } else if (nb_hexas>0) { this->numDim = 3; NumElems = nb_hexas; @@ -880,9 +827,7 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 4; elems = hexas; sides = quads; - } - else if (nb_trias>0) - { + } else if (nb_trias>0) { this->numDim = 2; NumElems = nb_trias; @@ -891,9 +836,7 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 2; elems = trias; sides = lines; - } - else if (nb_tri6>0) - { + } else if (nb_tri6>0) { this->numDim = 2; NumElems = nb_tri6; @@ -902,9 +845,7 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 3; elems = tri6; sides = line3; - } - else if (nb_quads>0) - { + } else if (nb_quads>0) { this->numDim = 2; NumElems = nb_quads; @@ -913,15 +854,12 @@ void Albany::GmshSTKMeshStruct::set_generic_mesh_info() NumSideNodes = 2; elems = quads; sides = lines; - } - else - { + } else { TEUCHOS_TEST_FOR_EXCEPTION (true, std::runtime_error, "Error! Invalid mesh dimension.\n"); } - return; } -void Albany::GmshSTKMeshStruct::store_element_info( +void GmshSTKMeshStruct::store_element_info( int e_type, int& iline, int& iline3, @@ -987,11 +925,9 @@ void Albany::GmshSTKMeshStruct::store_element_info( default: TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, "Error! Element type not supported; but you should have got an error before!\n"); } - - return; } -void Albany::GmshSTKMeshStruct::load_element_data( std::ifstream& ifile) +void GmshSTKMeshStruct::load_element_data( std::ifstream& ifile) { // Reset the stream to the beginning of the element section std::string line; @@ -1059,11 +995,9 @@ void Albany::GmshSTKMeshStruct::load_element_data( std::ifstream& ifile) tags.clear(); } } - - return; } -void Albany::GmshSTKMeshStruct::loadAsciiMesh () +void GmshSTKMeshStruct::loadAsciiMesh () { std::ifstream ifile; open_fname( ifile); @@ -1089,7 +1023,7 @@ void Albany::GmshSTKMeshStruct::loadAsciiMesh () } -void Albany::GmshSTKMeshStruct::loadBinaryMesh () +void GmshSTKMeshStruct::loadBinaryMesh () { std::ifstream ifile; open_fname( ifile); @@ -1346,7 +1280,7 @@ void Albany::GmshSTKMeshStruct::loadBinaryMesh () ifile.close(); } -void Albany::GmshSTKMeshStruct::set_all_nodes_boundary( std::vector& nsNames) +void GmshSTKMeshStruct::set_all_nodes_boundary( std::vector& nsNames) { std::string nsn = "Node"; nsNames.push_back(nsn); @@ -1354,11 +1288,9 @@ void Albany::GmshSTKMeshStruct::set_all_nodes_boundary( std::vector #ifdef ALBANY_SEACAS stk::io::put_io_part_attribute(*nsPartVec[nsn]); #endif - - return; } -void Albany::GmshSTKMeshStruct::set_all_sides_boundary( std::vector& ssNames) +void GmshSTKMeshStruct::set_all_sides_boundary( std::vector& ssNames) { std::string ssn = "BoundarySide"; ssNames.push_back(ssn); @@ -1367,13 +1299,11 @@ void Albany::GmshSTKMeshStruct::set_all_sides_boundary( std::vector stk::io::put_io_part_attribute(*ssPartVec[ssn]); stk::io::put_io_part_attribute(metaData->universal_part()); #endif - - return; } -void Albany::GmshSTKMeshStruct::set_boundaries( const Teuchos::RCP& commT, - std::vector& ssNames, - std::vector& nsNames) +void GmshSTKMeshStruct::set_boundaries (const Teuchos::RCP& comm, + std::vector& ssNames, + std::vector& nsNames) { set_all_nodes_boundary( nsNames); set_all_sides_boundary( ssNames); @@ -1387,14 +1317,14 @@ void Albany::GmshSTKMeshStruct::set_boundaries( const Teuchos::RCP(*commT, 0, 1, &numBdTags); + Teuchos::broadcast(*comm, 0, 1, &numBdTags); int* bdTagsArray = new int[numBdTags]; std::set::iterator it=bdTags.begin(); for (int k=0; it!=bdTags.end(); ++it,++k) { bdTagsArray[k] = *it; } - Teuchos::broadcast(*commT, 0, numBdTags, bdTagsArray); + Teuchos::broadcast(*comm, 0, numBdTags, bdTagsArray); // Adding boundary nodesets and sidesets separating different labels for (int k=0; k physical_surface_names; std::map physical_volume_names; - get_physical_names( physical_surface_names, commT, 2); - get_physical_names( physical_volume_names, commT, 3); + get_physical_names( physical_surface_names, comm, 2); + get_physical_names( physical_volume_names, comm, 3); std::map< std::string, int>::iterator name_it; for( name_it = physical_surface_names.begin(); name_it != physical_surface_names.end(); name_it++) @@ -1437,11 +1367,9 @@ void Albany::GmshSTKMeshStruct::set_boundaries( const Teuchos::RCPget_topology(); metaData->declare_part_with_topology(volume_i.str(), etopology); - - return; } -void Albany::GmshSTKMeshStruct::add_sideset( std::string sideset_name, int tag, std::vector& ssNames) +void GmshSTKMeshStruct::add_sideset( std::string sideset_name, int tag, std::vector& ssNames) { std::stringstream ssn_i; ssn_i << "BoundarySideSet" << sideset_name; @@ -1465,12 +1391,9 @@ void Albany::GmshSTKMeshStruct::add_sideset( std::string sideset_name, int tag, #ifdef ALBANY_SEACAS stk::io::put_io_part_attribute(*ssPartVec[ssn_i.str()]); #endif - - return; } - -void Albany::GmshSTKMeshStruct::add_nodeset( std::string nodeset_name, int tag, std::vector& nsNames) +void GmshSTKMeshStruct::add_nodeset( std::string nodeset_name, int tag, std::vector& nsNames) { std::stringstream nsn_i; nsn_i << "BoundaryNodeSet" << nodeset_name; @@ -1482,11 +1405,9 @@ void Albany::GmshSTKMeshStruct::add_nodeset( std::string nodeset_name, int tag, #ifdef ALBANY_SEACAS stk::io::put_io_part_attribute(*nsPartVec[nsn_i.str()]); #endif - - return; } -stk::topology Albany::GmshSTKMeshStruct::get_topology() +stk::topology GmshSTKMeshStruct::get_topology() { stk::topology etopology; switch (this->numDim) { @@ -1536,15 +1457,13 @@ stk::topology Albany::GmshSTKMeshStruct::get_topology() return etopology; } -void Albany::GmshSTKMeshStruct::set_allowable_gmsh_versions() +void GmshSTKMeshStruct::set_allowable_gmsh_versions() { allowable_gmsh_versions.insert( 2.2); allowable_gmsh_versions.insert( 4.1); - - return; } -bool Albany::GmshSTKMeshStruct::set_version_enum_from_float() +bool GmshSTKMeshStruct::set_version_enum_from_float() { bool can_read = true; if( version_in == (float)2.2 ) @@ -1563,7 +1482,7 @@ bool Albany::GmshSTKMeshStruct::set_version_enum_from_float() return can_read; } -void Albany::GmshSTKMeshStruct::check_version () +void GmshSTKMeshStruct::check_version () { // Tell user what gmsh version we're reading Teuchos::RCP out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); @@ -1584,22 +1503,18 @@ void Albany::GmshSTKMeshStruct::check_version () TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cannot read this version of gmsh *.msh file!"); } - - return; } -void Albany::GmshSTKMeshStruct::open_fname( std::ifstream& ifile) +void GmshSTKMeshStruct::open_fname( std::ifstream& ifile) { ifile.open(fname.c_str()); if (!ifile.is_open()) { TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Error! Cannot open mesh file '" << fname << "'.\n"); } - - return; } -void Albany::GmshSTKMeshStruct::get_name_for_physical_names( std::string& name, int& id, std::ifstream& ifile, int& dim) +void GmshSTKMeshStruct::get_name_for_physical_names( std::string& name, int& id, std::ifstream& ifile, int& dim) { std::string line; @@ -1619,11 +1534,9 @@ void Albany::GmshSTKMeshStruct::get_name_for_physical_names( std::string& name, name.erase( std::remove(name.begin(), name.end(), '"'), name.end()); name = "_" + name; } - - return; } -void Albany::GmshSTKMeshStruct::get_physical_tag_to_surface_tag_map( +void GmshSTKMeshStruct::get_physical_tag_to_surface_tag_map( std::ifstream& ifile, std::map& physical_surface_tags, int num_surfaces) @@ -1668,11 +1581,9 @@ void Albany::GmshSTKMeshStruct::get_physical_tag_to_surface_tag_map( physical_surface_tags.insert( std::make_pair( physical_tag, surface_tag)); } } - - return; } -void Albany::GmshSTKMeshStruct::get_physical_tag_to_volume_tag_map( +void GmshSTKMeshStruct::get_physical_tag_to_volume_tag_map( std::ifstream& ifile, std::map& physical_volume_tags, int num_volumes) @@ -1717,11 +1628,9 @@ void Albany::GmshSTKMeshStruct::get_physical_tag_to_volume_tag_map( physical_volume_tags.insert( std::make_pair( physical_tag, volume_tag)); } } - - return; } -void Albany::GmshSTKMeshStruct::read_physical_names_from_file( std::map& physical_names, int dim_) +void GmshSTKMeshStruct::read_physical_names_from_file( std::map& physical_names, int dim_) { std::ifstream ifile; open_fname( ifile); @@ -1837,50 +1746,44 @@ void Albany::GmshSTKMeshStruct::read_physical_names_from_file( std::map names, - int* tags_array, - int pair_number, - const Teuchos::RCP& commT, - std::map< std::string, int>& physical_names) +void GmshSTKMeshStruct:: +broadcast_name_tag_pair (std::vector< std::string> names, + int* tags_array, + int pair_number, + const Teuchos::RCP& comm, + std::map< std::string, int>& physical_names) { std::string name; - if( commT->getRank() == 0) - { + if( comm->getRank() == 0) { name = names[pair_number]; int strsize = name.size(); - Teuchos::broadcast( *commT, 0, &strsize); + Teuchos::broadcast( *comm, 0, &strsize); char* ptr = (strsize) ? (&name[0]) : 0; - Teuchos::broadcast( *commT, 0, strsize, ptr); - } - else - { + Teuchos::broadcast( *comm, 0, strsize, ptr); + } else { int strsize; - Teuchos::broadcast( *commT, 0, &strsize); + Teuchos::broadcast( *comm, 0, &strsize); name.resize( strsize); char* ptr = (strsize) ? (&name[0]) : 0; - Teuchos::broadcast( *commT, 0, strsize, ptr); + Teuchos::broadcast( *comm, 0, strsize, ptr); } int tag = tags_array[pair_number]; physical_names.insert( std::make_pair( name, tag)); - - return; } - -void Albany::GmshSTKMeshStruct::broadcast_physical_names( std::map& physical_names, - const Teuchos::RCP& commT) +void GmshSTKMeshStruct:: +broadcast_physical_names (std::map& physical_names, + const Teuchos::RCP& comm) { // Broadcast the number of name-tag pairs int num_pairs = physical_names.size(); - Teuchos::broadcast(*commT, 0, 1, &num_pairs); + Teuchos::broadcast(*comm, 0, 1, &num_pairs); // First unpack the names and tags from the map. // Only proc 0 will be doing anything here. @@ -1901,25 +1804,24 @@ void Albany::GmshSTKMeshStruct::broadcast_physical_names( std::map(*commT, 0, num_pairs, tags_array); + Teuchos::broadcast(*comm, 0, num_pairs, tags_array); for( int i = 0; i < num_pairs; i++) { - broadcast_name_tag_pair( names, tags_array, i, commT, physical_names); + broadcast_name_tag_pair( names, tags_array, i, comm, physical_names); } delete[] tags_array; - return; } -void Albany::GmshSTKMeshStruct::get_physical_names( std::map& physical_names, - const Teuchos::RCP& commT, +void GmshSTKMeshStruct::get_physical_names( std::map& physical_names, + const Teuchos::RCP& comm, int dim) { - if( commT->getRank() == 0 ) + if( comm->getRank() == 0 ) { read_physical_names_from_file( physical_names, dim); } - broadcast_physical_names( physical_names, commT); - - return; + broadcast_physical_names( physical_names, comm); } + +} // Albany diff --git a/src/disc/stk/Albany_GmshSTKMeshStruct.hpp b/src/disc/stk/Albany_GmshSTKMeshStruct.hpp index 78349b82a2..d42c6353f9 100644 --- a/src/disc/stk/Albany_GmshSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_GmshSTKMeshStruct.hpp @@ -23,20 +23,12 @@ class GmshSTKMeshStruct : public GenericSTKMeshStruct public: GmshSTKMeshStruct (const Teuchos::RCP& params, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, const int numParams); ~GmshSTKMeshStruct(); - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); - - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return false; } @@ -48,7 +40,7 @@ class GmshSTKMeshStruct : public GenericSTKMeshStruct // Set boundary information. // Includes sideset and nodeset names and counts. - void set_boundaries( const Teuchos::RCP& commT, + void set_boundaries( const Teuchos::RCP& comm, std::vector& ssNames, std::vector& nsNames); @@ -73,12 +65,12 @@ class GmshSTKMeshStruct : public GenericSTKMeshStruct // Gets the physical name-tag pairs for version 4.1 meshes void get_physical_names( std::map& physical_names, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, int dim); // Share physical_names map with all other processes void broadcast_physical_names( std::map& physical_names, - const Teuchos::RCP& commT); + const Teuchos::RCP& comm); // Read the physical names for Gmsh V 4.1 // to populate the physical_names map @@ -93,7 +85,7 @@ class GmshSTKMeshStruct : public GenericSTKMeshStruct void determine_file_type( bool& legacy, bool& binary, bool& ascii); // Broadcast topology of the mesh from 0 to all over procs - void broadcast_topology( const Teuchos::RCP& commT); + void broadcast_topology( const Teuchos::RCP& comm); // Sets NumNodes for ascii msh files void set_NumNodes( std::ifstream& ifile); @@ -144,7 +136,7 @@ class GmshSTKMeshStruct : public GenericSTKMeshStruct void broadcast_name_tag_pair( std::vector< std::string> names, int* tags_array, int pair_number, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, std::map< std::string, int>& physical_names); // Reads a single physical name from ifile. diff --git a/src/disc/stk/Albany_IossSTKMeshStruct.cpp b/src/disc/stk/Albany_IossSTKMeshStruct.cpp index c74d3b3995..9ab0f73b71 100644 --- a/src/disc/stk/Albany_IossSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_IossSTKMeshStruct.cpp @@ -8,8 +8,10 @@ #ifdef ALBANY_SEACAS -#include +#include "Albany_Utils.hpp" +#include "Albany_CommUtils.hpp" +#include #include "Teuchos_VerboseObject.hpp" #include @@ -24,9 +26,8 @@ #include -#include "Albany_Utils.hpp" +#include -#include namespace { @@ -45,33 +46,33 @@ void get_element_block_sizes(stk::io::StkMeshIoBroker &mesh_data, } // Anonymous namespace -Albany::IossSTKMeshStruct::IossSTKMeshStruct( - const Teuchos::RCP& params_, - const Teuchos::RCP& commT, - const int numParams_) : - GenericSTKMeshStruct(params_, -1, numParams_), - out(Teuchos::VerboseObjectBase::getDefaultOStream()), - useSerialMesh(false), - periodic(params->get("Periodic BC", false)), - m_hasRestartSolution(false), - m_restartDataTime(-1.0), - m_solutionFieldHistoryDepth(0) +namespace Albany +{ + +IossSTKMeshStruct:: +IossSTKMeshStruct(const Teuchos::RCP& params_, + const Teuchos::RCP& comm, + const int numParams_) + : GenericSTKMeshStruct(params_, -1, numParams_) + , out(Teuchos::VerboseObjectBase::getDefaultOStream()) + , useSerialMesh(false) + , periodic(params->get("Periodic BC", false)) { params->validateParameters(*getValidDiscretizationParameters(),0); usePamgen = (params->get("Method","Exodus") == "Pamgen"); - std::vector entity_rank_names = stk::mesh::entity_rank_names(); - const Teuchos::MpiComm* theComm = dynamic_cast* > (commT.get()); + auto mpiComm = getMpiCommFromTeuchosComm(comm); - mesh_data = Teuchos::rcp(new stk::io::StkMeshIoBroker(*theComm->getRawMpiComm())); + mesh_data = Teuchos::rcp(new stk::io::StkMeshIoBroker(mpiComm)); + mesh_data->use_simple_fields(); // Use Greg Sjaardema's capability to repartition on the fly. // Several partitioning choices: rcb, rib, hsfc, kway, kway-gemo, linear, random // linear does not require Zoltan or metis - if (params->get("Use Serial Mesh", false) && commT->getSize() > 1){ + if (params->get("Use Serial Mesh", false) && comm->getSize() > 1){ // Option external reads the nemesis files, and must be the default #ifdef ALBANY_ZOLTAN mesh_data->property_add(Ioss::Property("DECOMPOSITION_METHOD", "rib")); @@ -111,8 +112,8 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( typedef Teuchos::Array StringArray; const StringArray additionalNodeSets = params->get("Additional Node Sets", StringArray()); - for (StringArray::const_iterator it = additionalNodeSets.begin(), it_end = additionalNodeSets.end(); it != it_end; ++it) { - stk::mesh::Part &newNodeSet = metaData->declare_part(*it, stk::topology::NODE_RANK); + for (const auto& nsn : additionalNodeSets) { + stk::mesh::Part &newNodeSet = metaData->declare_part(nsn, stk::topology::NODE_RANK); if (!stk::io::is_part_io_part(newNodeSet)) { stk::mesh::Field * const distrFactorfield = metaData->get_field(stk::topology::NODE_RANK, "distribution_factors"); if (distrFactorfield != NULL){ @@ -132,11 +133,7 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( std::vector nsNames; int numEB = 0; - for (stk::mesh::PartVector::const_iterator i = all_parts.begin(); - i != all_parts.end(); ++i) { - - stk::mesh::Part * const part = *i ; - + for (const auto& part : all_parts) { if (!stk::mesh::is_auto_declared_part(*part)) { if ( part->primary_entity_rank() == stk::topology::ELEMENT_RANK) { @@ -157,6 +154,9 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( } } + *out << "IOSS-STK: number of node sets = " << nsPartVec.size() << std::endl; + *out << "IOSS-STK: number of side sets = " << ssPartVec.size() << std::endl; + // the method 'initializesidesetmeshspecs' requires that ss parts store a valid stk topology. // therefore, we try to retrieve the topology of this part using stk stuff. auto r = mesh_data->get_input_ioss_region(); @@ -177,19 +177,6 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( cullSubsetParts(ssNames, ssPartVec); // Eliminate sidesets that are subsets of other sidesets -#if 0 - // for debugging, print out the parts now - std::map::iterator it; - - for(it = ssPartVec.begin(); it != ssPartVec.end(); ++it){ // loop over the parts in the map - - // for each part in turn, get the name of parts that are a subset of it - - print(*out, "Found \n", *it->second); - } - // end debugging -#endif - int worksetSizeMax = params->get("Workset Size", -1); // Get number of elements per element block using Ioss for use @@ -210,17 +197,13 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( } // Construct MeshSpecsStruct - { - const CellTopologyData& ctd = *elementBlockTopologies_[0].getCellTopologyData(); - this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct( - ctd, numDim, nsNames, ssNames, worksetSize, partVec[0]->name(), - ebNameToIndex)); - } + const CellTopologyData& ctd = *elementBlockTopologies_[0].getCellTopologyData(); + this->meshSpecs[0] = Teuchos::rcp(new MeshSpecsStruct( + ctd, numDim, nsNames, ssNames, worksetSize, partVec[0]->name(), + ebNameToIndex)); - { - const Ioss::Region& inputRegion = *(mesh_data->get_input_ioss_region()); - m_solutionFieldHistoryDepth = inputRegion.get_property("state_count").get_int(); - } + const Ioss::Region& inputRegion = *(mesh_data->get_input_ioss_region()); + m_solutionFieldHistoryDepth = inputRegion.get_property("state_count").get_int(); // Upon request, add a nodeset for each sideset if (params->get("Build Node Sets From Side Sets",false)) @@ -233,7 +216,7 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( this->setAllPartsIO(); // Create a mesh specs object for EACH side set - this->initializeSideSetMeshSpecs(commT); + this->initializeSideSetMeshSpecs(comm); // Get upper bound on sideset workset sizes by using Ioss element counts on side blocks if (worksetSize == ebSizeMax) { @@ -271,12 +254,9 @@ Albany::IossSTKMeshStruct::IossSTKMeshStruct( } } } - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(commT); } -Albany::IossSTKMeshStruct::~IossSTKMeshStruct() +IossSTKMeshStruct::~IossSTKMeshStruct() { // Explicitly delete these in exactly this order. There are three // reasons. First, StkMeshIoBroker does not have a MetaData getter that @@ -295,143 +275,96 @@ Albany::IossSTKMeshStruct::~IossSTKMeshStruct() mesh_data = Teuchos::null; } -void -Albany::IossSTKMeshStruct::setFieldData ( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void IossSTKMeshStruct:: +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) { - this->SetupFieldData(commT, sis, worksetSize); + GenericSTKMeshStruct::setFieldData(comm, sis); if(mesh_data->is_bulk_data_null()) mesh_data->set_bulk_data(*bulkData); - *out << "IOSS-STK: number of node sets = " << nsPartVec.size() << std::endl; - *out << "IOSS-STK: number of side sets = " << ssPartVec.size() << std::endl; - - std::vector missing; // Restart index to read solution from exodus file. - int index = params->get("Restart Index",-1); // Default to no restart - double res_time = params->get("Restart Time",-1.0); // Default to no restart - Ioss::Region& region = *(mesh_data->get_input_ioss_region()); - /* - * The following code block reads a single mesh on PE 0, then distributes the mesh across - * the other processors. stk_rebalance is used, which requires Zoltan - * - * This code is only compiled if ALBANY_ZOLTAN is true - */ - -#ifdef ALBANY_ZOLTAN // rebalance needs Zoltan - - if(useSerialMesh){ - - // trick to avoid hanging + if (params->isParameter("Restart Index")) { + TEUCHOS_TEST_FOR_EXCEPTION ( + params->isParameter("Restart Time"), std::logic_error, + "Error! Do not provide both 'Restart Index' and 'Restart Time'.\n"); + + // User has specified a time step to restart at + int index = params->get("Restart Index"); // Default to no restart + + const auto& region = *mesh_data->get_input_ioss_region(); + m_restartDataTime = region.get_state_time(index); + m_hasRestartSolution = true; + + *out << "Restart Index set, reading solution index : " << index << std::endl; + } else if (params->isParameter("Restart Time")) { + m_restartDataTime = params->get("Restart Time"); + // User has specified a time to restart at + m_hasRestartSolution = true; + *out << "Restart solution time set, reading solution time : " << m_restartDataTime << std::endl; + } else { + *out << "Neither restart index or time are set. Not reading solution data from exodus file"<< std::endl; + } - if(commT->getRank() == 0){ // read in the mesh on PE 0 - mesh_data->populate_bulk_data(); + if (m_hasRestartSolution) { + mesh_data->add_all_mesh_fields_as_input_fields(); // KL: this adds "solution field" + Teuchos::Array default_field = {{"solution", "solution_dot", "solution_dotdot"}}; + const auto& restart_fields = params->get >("Restart Fields", default_field); - if (this->numDim!=3) - { - // Try to load 3d coordinates (if present in the input file) - loadOrSetCoordinates3d(index); - } + // Check if states are available in the mesh + const auto& region = *mesh_data->get_input_ioss_region(); + const auto& node_blocks = region.get_node_blocks(); - // Read solution from exodus file. - if (index >= 0) { // User has specified a time step to restart at - m_restartDataTime = region.get_state_time(index); - m_hasRestartSolution = true; - } - else if (res_time >= 0) { // User has specified a time to restart at - m_restartDataTime = res_time; - m_hasRestartSolution = true; + for (const auto& st_ptr : *sis) { + auto& st = *st_ptr; + if (std::find(restart_fields.begin(),restart_fields.end(),st.name)==restart_fields.end()) { + continue; } - else { - *out << "Neither restart index or time are set. Not reading solution data from exodus file"<< std::endl; + st.restartDataAvailable = node_blocks.size()>0; + for (const auto& block : node_blocks) { + st.restartDataAvailable &= block->field_exists(st.name); } } - } // End UseSerialMesh - reading mesh on PE 0 - - else -#endif - - /* - * The following code block reads a single mesh when Albany is compiled serially, or a - * Nemspread fileset if ALBANY_MPI is true. - * - */ - - { // running in Serial or Parallel read from Nemspread files - if (this->numDim!=3) - { - // Try to load 3d coordinates (if present in the input file) - loadOrSetCoordinates3d(index); - } + } - if (!usePamgen) - { - // Read solution from exodus file. - if (index >= 0) - { // User has specified a time step to restart at - m_restartDataTime = region.get_state_time(index); - m_hasRestartSolution = true; - } - else if (res_time >= 0) - { // User has specified a time to restart at - m_restartDataTime = res_time; - m_hasRestartSolution = true; - } - else - { - *out << "Restart Index not set. Not reading solution from exodus (" << index << ")"<< std::endl; - } - } +} - } // End Parallel Read - or running in serial +void IossSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) +{ + metaData->commit(); - if(m_hasRestartSolution){ +#ifdef ALBANY_ZOLTAN // rebalance needs Zoltan + // The following code block reads a single mesh on PE 0, then distributes the mesh across + // the other processors. stk_rebalance is used, which requires Zoltan + if (useSerialMesh){ - Teuchos::Array default_field = {{"solution", "solution_dot", "solution_dotdot"}}; - auto& restart_fields = params->get >("Restart Fields", default_field); - - // Get the fields to be used for restart - - // See what state data was initialized from the stk::io request - // This should be propagated into stk::io - const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks(); - - // Uncomment to print what fields are in the exodus file - // const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks(); - // for (const auto& block : elem_blocks) { - // Ioss::NameList exo_eb_fld_names; - // block->field_describe(&exo_eb_fld_names); - // for (const auto& name : exo_eb_fld_names) { - // *out << "Found field \"" << name << "\" in elem blocks of exodus file\n"; - // } - // } - // for (const auto& block : node_blocks) { - // Ioss::NameList exo_nb_fld_names; - // block->field_describe(&exo_nb_fld_names); - // for (const auto& name : exo_nb_fld_names) { - // *out << "Found field \"" << name << "\" in node blocks of exodus file\n"; - // } - // } + // trick to avoid hanging + bulkData->modification_begin(); - for (const auto& st_ptr : *sis) { - auto& st = *st_ptr; - for (const auto& block : node_blocks) { - if (block->field_exists(st.name)) { - for (const auto& restart_field : restart_fields) { - if (st.name==restart_field) { - *out << "Restarting from field \"" << st.name << "\" found in exodus file.\n"; - st.restartDataAvailable = true; - break; - } - } - } - } + if(comm->getRank() == 0){ + mesh_data->populate_bulk_data(); + } else { + // trick to avoid hanging + bulkData->modification_begin(); bulkData->modification_begin(); } + bulkData->modification_end(); + } else +#endif + { + // The following code block reads a single mesh when Albany is compiled serially, or a + // Nemspread fileset if ALBANY_MPI is true. + bulkData->modification_begin(); + mesh_data->populate_bulk_data(); + bulkData->modification_end(); + } + + // Note: we cannot load fields/coords during setFieldData, since we need a valid bulkData + std::vector missing; + if (m_hasRestartSolution) { + mesh_data->read_defined_input_fields(m_restartDataTime, &missing); // Read global mesh variables. Should we emit warnings at all? for (auto& it : fieldContainer->getMeshVectorStates()) { bool found = mesh_data->get_global (it.first, it.second, false); // Last variable is abort_if_not_found. We don't want that. @@ -452,25 +385,14 @@ Albany::IossSTKMeshStruct::setFieldData ( if (!found) *out << " *** WARNING *** Mesh scalar integer state '" << it.first << "' was not found in the mesh database.\n"; } - } else { - // We put all the fields as 'missing' - const stk::mesh::FieldVector& fields = metaData->get_fields(); - for (decltype(fields.size()) i=0; iname()); - missing.push_back(stk::io::MeshField(fields[i],fields[i]->name())); - } } // If this is a boundary mesh, the side_map/side_node_map may already be present, so we check - side_maps_present = true; - bool coherence = true; - for (const auto& it : missing) { - if (it.field()->name()=="side_to_cell_map" || it.field()->name()=="side_nodes_ids") - { - side_maps_present = false; - coherence = !coherence; // Both fields should be present or absent so coherence should change exactly twice - } - } + auto region = mesh_data->get_input_ioss_region(); + const auto& elem_blocks = region->get_element_blocks(); + auto *eb = elem_blocks[0]; + side_maps_present = eb->field_exists("side_to_cell_map") and eb->field_exists("side_nodes_ids"); + const bool coherence = eb->field_exists("side_to_cell_map") == eb->field_exists("side_nodes_ids"); TEUCHOS_TEST_FOR_EXCEPTION (!coherence, std::runtime_error, "Error! The maps 'side_to_cell_map' and 'side_nodes_ids' should either both be present or both missing, but only one of them was found in the mesh file.\n"); if (useSerialMesh) @@ -478,140 +400,13 @@ Albany::IossSTKMeshStruct::setFieldData ( // Only proc 0 actually read the mesh, and can confirm whether or not the side maps were present // Unfortunately, Teuchos does not have a specialization for type bool when it comes to communicators, so we need ints int bool_to_int = side_maps_present ? 1 : 0; - Teuchos::broadcast(*commT,0,1,&bool_to_int); + Teuchos::broadcast(*comm,0,1,&bool_to_int); side_maps_present = bool_to_int == 1 ? true : false; } - // Loading required input fields from file - //this->loadRequiredInputFields (commT); - - // Rebalance the mesh before starting the simulation if indicated - rebalanceInitialMeshT(commT); - - // Check that the nodeset created from sidesets contain the right number of nodes - this->checkNodeSetsFromSideSetsIntegrity (); - - this->setSideSetFieldData(commT, side_set_sis, worksetSize); -} - -void -Albany::IossSTKMeshStruct::setBulkData ( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - mesh_data->add_all_mesh_fields_as_input_fields(); // KL: this adds "solution field" - std::vector missing; - - metaData->commit(); - - // Restart index to read solution from exodus file. - int index = params->get("Restart Index",-1); // Default to no restart - double res_time = params->get("Restart Time",-1.0); // Default to no restart - Ioss::Region& region = *(mesh_data->get_input_ioss_region()); - /* - * The following code block reads a single mesh on PE 0, then distributes the mesh across - * the other processors. stk_rebalance is used, which requires Zoltan - * - * This code is only compiled if ALBANY_ZOLTAN is true - */ - -#ifdef ALBANY_ZOLTAN // rebalance needs Zoltan - - if(useSerialMesh){ - - // trick to avoid hanging - bulkData->modification_begin(); - - if(commT->getRank() == 0){ // read in the mesh on PE 0 - - - //stk::io::process_mesh_bulk_data(region, *bulkData); - mesh_data->populate_bulk_data(); - - if (this->numDim!=3) - { - // Try to load 3d coordinates (if present in the input file) - loadOrSetCoordinates3d(index); - } - - //bulkData = &mesh_data->bulk_data(); - - // Read solution from exodus file. - if (index >= 0) { // User has specified a time step to restart at - *out << "Restart Index set, reading solution index : " << index << std::endl; - mesh_data->read_defined_input_fields(index, &missing); - m_restartDataTime = region.get_state_time(index); - m_hasRestartSolution = true; - } - else if (res_time >= 0) { // User has specified a time to restart at - *out << "Restart solution time set, reading solution time : " << res_time << std::endl; - mesh_data->read_defined_input_fields(res_time, &missing); - m_restartDataTime = res_time; - m_hasRestartSolution = true; - } - else { - *out << "Neither restart index or time are set. Not reading solution data from exodus file"<< std::endl; - } - } - else { - // trick to avoid hanging - bulkData->modification_begin(); bulkData->modification_begin(); - } - - bulkData->modification_end(); - - } // End UseSerialMesh - reading mesh on PE 0 - - else -#endif - - /* - * The following code block reads a single mesh when Albany is compiled serially, or a - * Nemspread fileset if ALBANY_MPI is true. - * - */ - - { // running in Serial or Parallel read from Nemspread files - bulkData->modification_begin(); - mesh_data->populate_bulk_data(); - if (this->numDim!=3) - { - // Try to load 3d coordinates (if present in the input file) - loadOrSetCoordinates3d(index); - } - - if (!usePamgen) - { - // Read solution from exodus file. - if (index >= 0) - { // User has specified a time step to restart at - *out << "Restart Index set, reading solution index : " << index << std::endl; - mesh_data->read_defined_input_fields(index, &missing); - m_restartDataTime = region.get_state_time(index); - m_hasRestartSolution = true; - } - else if (res_time >= 0) - { // User has specified a time to restart at - *out << "Restart solution time set, reading solution time : " << res_time << std::endl; - mesh_data->read_defined_input_fields(res_time, &missing); - m_restartDataTime = res_time; - m_hasRestartSolution = true; - } - else - { - *out << "Restart Index not set. Not reading solution from exodus (" << index << ")"<< std::endl; - } - } - - bulkData->modification_end(); - - } // End Parallel Read - or running in serial - // Check if the input mesh is layered (i.e., if it stores layers info) std::string state_name = "layer_thickness_ratio"; - if(mesh_data->has_input_global(state_name)) { + if (mesh_data->has_input_global(state_name)) { // layer ratios std::vector ltr; mesh_data->get_global (state_name, ltr, true); @@ -667,119 +462,23 @@ Albany::IossSTKMeshStruct::setBulkData ( this->local_cell_layers_data->bot_side_pos = this->meshSpecs[0]->ctd.side_count - 2; } - if(m_hasRestartSolution){ - - Teuchos::Array default_field = {{"solution", "solution_dot", "solution_dotdot"}}; - auto& restart_fields = params->get >("Restart Fields", default_field); - - // Get the fields to be used for restart - - // See what state data was initialized from the stk::io request - // This should be propagated into stk::io - const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks(); - - // Uncomment to print what fields are in the exodus file - // const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks(); - // for (const auto& block : elem_blocks) { - // Ioss::NameList exo_eb_fld_names; - // block->field_describe(&exo_eb_fld_names); - // for (const auto& name : exo_eb_fld_names) { - // *out << "Found field \"" << name << "\" in elem blocks of exodus file\n"; - // } - // } - // for (const auto& block : node_blocks) { - // Ioss::NameList exo_nb_fld_names; - // block->field_describe(&exo_nb_fld_names); - // for (const auto& name : exo_nb_fld_names) { - // *out << "Found field \"" << name << "\" in node blocks of exodus file\n"; - // } - // } - - for (const auto& st_ptr : *sis) { - auto& st = *st_ptr; - for (const auto& block : node_blocks) { - if (block->field_exists(st.name)) { - for (const auto& restart_field : restart_fields) { - if (st.name==restart_field) { - *out << "Restarting from field \"" << st.name << "\" found in exodus file.\n"; - st.restartDataAvailable = true; - break; - } - } - } - } - } - - // Read global mesh variables. Should we emit warnings at all? - for (auto& it : fieldContainer->getMeshVectorStates()) { - bool found = mesh_data->get_global (it.first, it.second, false); // Last variable is abort_if_not_found. We don't want that. - if (!found) - *out << " *** WARNING *** Mesh vector state '" << it.first << "' was not found in the mesh database.\n"; - } - - for (auto& it : fieldContainer->getMeshScalarIntegerStates()) { - bool found = mesh_data->get_global (it.first, it.second, false); // Last variable is abort_if_not_found. We don't want that. - if (!found) - *out << " *** WARNING *** Mesh scalar integer state '" << it.first << "' was not found in the mesh database.\n"; - } - } else { - // We put all the fields as 'missing' - const stk::mesh::FieldVector& fields = metaData->get_fields(); - for (decltype(fields.size()) i=0; iname()); - missing.push_back(stk::io::MeshField(fields[i],fields[i]->name())); - } - } - - // If this is a boundary mesh, the side_map/side_node_map may already be present, so we check - side_maps_present = true; - bool coherence = true; - for (const auto& it : missing) - { - if (it.field()->name()=="side_to_cell_map" || it.field()->name()=="side_nodes_ids") - { - side_maps_present = false; - coherence = !coherence; // Both fields should be present or absent so coherence should change exactly twice - } - } - TEUCHOS_TEST_FOR_EXCEPTION (!coherence, std::runtime_error, "Error! The maps 'side_to_cell_map' and 'side_nodes_ids' should either both be present or both missing, but only one of them was found in the mesh file.\n"); - - if (useSerialMesh) - { - // Only proc 0 actually read the mesh, and can confirm whether or not the side maps were present - // Unfortunately, Teuchos does not have a specialization for type bool when it comes to communicators, so we need ints - int bool_to_int = side_maps_present ? 1 : 0; - Teuchos::broadcast(*commT,0,1,&bool_to_int); - side_maps_present = bool_to_int == 1 ? true : false; - } + // Set 3d coords *before* loading fields, or else you won't have x/y/z available + loadOrSetCoordinates3d(); - // Loading required input fields from file - this->loadRequiredInputFields (commT); + // Load any required field from file + this->loadRequiredInputFields (comm); // Rebalance the mesh before starting the simulation if indicated - rebalanceInitialMeshT(commT); + rebalanceInitialMesh(comm); // Check that the nodeset created from sidesets contain the right number of nodes this->checkNodeSetsFromSideSetsIntegrity (); - // Finally, perform the setup of the (possible) side set meshes (including extraction if of type SideSetSTKMeshStruct) - this->setSideSetBulkData(commT, side_set_sis, worksetSize); - - fieldAndBulkDataSet = true; -} - -double -Albany::IossSTKMeshStruct::getSolutionFieldHistoryStamp(int step) const -{ - TEUCHOS_ASSERT(step >= 0 && step < m_solutionFieldHistoryDepth); - - const int index = step + 1; // 1-based step indexing - const Ioss::Region & inputRegion = *(mesh_data->get_input_ioss_region()); - return inputRegion.get_state_time(index); + m_bulk_data_set = true; } void -Albany::IossSTKMeshStruct::loadOrSetCoordinates3d(int index) +IossSTKMeshStruct::loadOrSetCoordinates3d() { const std::string coords3d_name = "coordinates3d"; @@ -787,8 +486,12 @@ Albany::IossSTKMeshStruct::loadOrSetCoordinates3d(int index) const Ioss::NodeBlockContainer& node_blocks = region->get_node_blocks(); Ioss::NodeBlock *nb = node_blocks[0]; - if (nb->field_exists(coords3d_name) && index>0) - { + const int current_step = region->get_current_state(); + const int index = params->isParameter("Restart Index") + ? params->get("Restart Index") + : current_step; + + if (nb->field_exists(coords3d_name) and index!=-1) { // The field "coordinates3d" exists in the input mesh // (which must then come from a previous Albany run), so load it. std::vector nodes; @@ -801,40 +504,41 @@ Albany::IossSTKMeshStruct::loadOrSetCoordinates3d(int index) // or different to the one desired, we set it to the restart // index state. If the state was already set, we also take // care of resetting the index back to the original configuration. - const int current_step = region->get_current_state(); - if (current_step != index) { - if (current_step != -1) { - region->end_state(current_step); - } - region->begin_state(index); + + if (current_step != -1) { + region->end_state(current_step); } + region->begin_state(index); stk::io::field_data_from_ioss(*bulkData, this->getCoordinatesField3d(), nodes, nb, coords3d_name); + region->end_state(index); + if (current_step != -1) { - region->begin_state(index); + region->begin_state(current_step); } - } - else - { + } else { if (nb->field_exists(coords3d_name)) { // The 3d coords field exists in the input mesh, but the restart index was // not set in the input file. We issue a warning, and load default coordinates *out << "WARNING! The field 'coordinates3d' was found in the input mesh, but no restart index was specified.\n" << " Albany will set the 3d coordinates to the 'default' ones (filling native coordinates with trailing zeros).\n"; } - // The input mesh does not store the 'coordinates3d' field - // (perhaps the mesh does not come from a previous Albany run). - // Hence, we initialize coordinates3d with coordinates, - // and we fill 'extra' dimensions with 0's. Hopefully, this is ok. - - // Use GenericSTKMeshStruct functionality - this->setDefaultCoordinates3d(); + setDefaultCoordinates3d(); } } +double +IossSTKMeshStruct::getSolutionFieldHistoryStamp(int step) const +{ + TEUCHOS_ASSERT(step >= 0 && step < m_solutionFieldHistoryDepth); + + const int index = step + 1; // 1-based step indexing + const Ioss::Region & inputRegion = *(mesh_data->get_input_ioss_region()); + return inputRegion.get_state_time(index); +} void -Albany::IossSTKMeshStruct::loadSolutionFieldHistory(int step) +IossSTKMeshStruct::loadSolutionFieldHistory(int step) { TEUCHOS_ASSERT(step >= 0 && step < m_solutionFieldHistoryDepth); @@ -843,7 +547,7 @@ Albany::IossSTKMeshStruct::loadSolutionFieldHistory(int step) } Teuchos::RCP -Albany::IossSTKMeshStruct::getValidDiscretizationParameters() const +IossSTKMeshStruct::getValidDiscretizationParameters() const { Teuchos::RCP validPL = this->getValidGenericSTKParameters("Valid IOSS_DiscParams"); @@ -861,4 +565,6 @@ Albany::IossSTKMeshStruct::getValidDiscretizationParameters() const return validPL; } +} // namespace Albany + #endif // ALBANY_SEACAS diff --git a/src/disc/stk/Albany_IossSTKMeshStruct.hpp b/src/disc/stk/Albany_IossSTKMeshStruct.hpp index 414debaa05..cedb4511df 100644 --- a/src/disc/stk/Albany_IossSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_IossSTKMeshStruct.hpp @@ -21,22 +21,17 @@ namespace Albany { class IossSTKMeshStruct : public GenericSTKMeshStruct { - public: + public: IossSTKMeshStruct (const Teuchos::RCP& params, const Teuchos::RCP& commT, const int numParams); ~IossSTKMeshStruct(); - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis); - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); int getSolutionFieldHistoryDepth() const {return m_solutionFieldHistoryDepth;} double getSolutionFieldHistoryStamp(int step) const; @@ -48,11 +43,11 @@ namespace Albany { //! If restarting, convenience function to return restart data time double restartDataTime() const {return m_restartDataTime;} - private: + private: - Ioss::Init::Initializer ioInit; + void loadOrSetCoordinates3d (); - void loadOrSetCoordinates3d (int index); + Ioss::Init::Initializer ioInit; Teuchos::RCP getValidDiscretizationParameters() const; @@ -62,9 +57,9 @@ namespace Albany { bool periodic; Teuchos::RCP mesh_data; - bool m_hasRestartSolution; + bool m_hasRestartSolution = false; double m_restartDataTime; - int m_solutionFieldHistoryDepth; + int m_solutionFieldHistoryDepth = 0; }; diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp index 02d3d8f81e..208c3a84dd 100644 --- a/src/disc/stk/Albany_MultiSTKFieldContainer.cpp +++ b/src/disc/stk/Albany_MultiSTKFieldContainer.cpp @@ -40,7 +40,6 @@ MultiSTKFieldContainer::MultiSTKFieldContainer( const Teuchos::RCP& metaData_, const Teuchos::RCP& bulkData_, const int numDim_, - const Teuchos::RCP& sis, const int num_params_) : GenericSTKFieldContainer( params_, @@ -80,8 +79,6 @@ MultiSTKFieldContainer::MultiSTKFieldContainer( #endif } - this->addStateStructs(sis); - initializeProcRankField(); } @@ -104,137 +101,112 @@ MultiSTKFieldContainer::MultiSTKFieldContainer( { typedef typename AbstractSTKFieldContainer::STKFieldType SFT; - sol_vector_name.resize(solution_vector.size()); - sol_index.resize(solution_vector.size()); + if (save_solution_field) { + sol_vector_name.resize(solution_vector.size()); + sol_index.resize(solution_vector.size()); - // Check the input + // Check the input - auto const num_derivs = solution_vector[0].size(); - for (auto i = 1; i < solution_vector.size(); ++i) { - ALBANY_ASSERT( - solution_vector[i].size() == num_derivs, - "\n*** ERROR ***\n" - "Number of derivatives for each variable is different.\n" - "Check definition of solution vector and its derivatives.\n"); - } + auto const num_derivs = solution_vector[0].size(); + for (auto i = 1; i < solution_vector.size(); ++i) { + ALBANY_ASSERT( + solution_vector[i].size() == num_derivs, + "\n*** ERROR ***\n" + "Number of derivatives for each variable is different.\n" + "Check definition of solution vector and its derivatives.\n"); + } - for (int vec_num = 0; vec_num < solution_vector.size(); vec_num++) { - if (solution_vector[vec_num].size() == - 0) { // Do the default solution vector + for (int vec_num = 0; vec_num < solution_vector.size(); vec_num++) { + if (solution_vector[vec_num].size() == + 0) { // Do the default solution vector - std::string name = params_->get( - sol_tag_name[vec_num], sol_id_name[vec_num]); - SFT* solution = - &metaData_->declare_field(stk::topology::NODE_RANK, name); - stk::mesh::put_field_on_mesh( - *solution, metaData_->universal_part(), neq_, nullptr); + std::string name = params_->get( + sol_tag_name[vec_num], sol_id_name[vec_num]); + SFT* solution = + &metaData_->declare_field(stk::topology::NODE_RANK, name); + stk::mesh::put_field_on_mesh( + *solution, metaData_->universal_part(), neq_, nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); #endif - sol_vector_name[vec_num].push_back(name); - sol_index[vec_num].push_back(this->neq); - } else if (solution_vector[vec_num].size() == 1) { // User is just renaming - // the entire solution - // vector + sol_vector_name[vec_num].push_back(name); + sol_index[vec_num].push_back(this->neq); + } else if (solution_vector[vec_num].size() == 1) { // User is just renaming + // the entire solution + // vector - SFT* solution = &metaData_->declare_field( - stk::topology::NODE_RANK, solution_vector[vec_num][0]); - stk::mesh::put_field_on_mesh( - *solution, metaData_->universal_part(), neq_, nullptr); + SFT* solution = &metaData_->declare_field( + stk::topology::NODE_RANK, solution_vector[vec_num][0]); + stk::mesh::put_field_on_mesh( + *solution, metaData_->universal_part(), neq_, nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); #endif - sol_vector_name[vec_num].push_back(solution_vector[vec_num][0]); - sol_index[vec_num].push_back(neq_); + sol_vector_name[vec_num].push_back(solution_vector[vec_num][0]); + sol_index[vec_num].push_back(neq_); - } else { // user is breaking up the solution into multiple fields + } else { // user is breaking up the solution into multiple fields - // make sure the number of entries is even + // make sure the number of entries is even - TEUCHOS_TEST_FOR_EXCEPTION( - (solution_vector[vec_num].size() % 2), - std::logic_error, - "Error in input file: specification of solution vector layout is " - "incorrect." - << std::endl); + TEUCHOS_TEST_FOR_EXCEPTION( + (solution_vector[vec_num].size() % 2), + std::logic_error, + "Error in input file: specification of solution vector layout is " + "incorrect." + << std::endl); - int len, accum = 0; + int len, accum = 0; - for (int i = 0; i < solution_vector[vec_num].size(); i += 2) { - if (solution_vector[vec_num][i + 1] == "V") { - len = numDim_; // vector - accum += len; - SFT* solution = &metaData_->declare_field( - stk::topology::NODE_RANK, solution_vector[vec_num][i]); - stk::mesh::put_field_on_mesh( - *solution, metaData_->universal_part(), len, nullptr); + for (int i = 0; i < solution_vector[vec_num].size(); i += 2) { + if (solution_vector[vec_num][i + 1] == "V") { + len = numDim_; // vector + accum += len; + SFT* solution = &metaData_->declare_field( + stk::topology::NODE_RANK, solution_vector[vec_num][i]); + stk::mesh::put_field_on_mesh( + *solution, metaData_->universal_part(), len, nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); #endif - sol_vector_name[vec_num].push_back(solution_vector[vec_num][i]); - sol_index[vec_num].push_back(len); - - } else if (solution_vector[vec_num][i + 1] == "S") { - len = 1; // scalar - accum += len; - SFT* solution = &metaData_->declare_field( - stk::topology::NODE_RANK, solution_vector[vec_num][i]); - stk::mesh::put_field_on_mesh( - *solution, metaData_->universal_part(), nullptr); + sol_vector_name[vec_num].push_back(solution_vector[vec_num][i]); + sol_index[vec_num].push_back(len); + + } else if (solution_vector[vec_num][i + 1] == "S") { + len = 1; // scalar + accum += len; + SFT* solution = &metaData_->declare_field( + stk::topology::NODE_RANK, solution_vector[vec_num][i]); + stk::mesh::put_field_on_mesh( + *solution, metaData_->universal_part(), nullptr); #ifdef ALBANY_SEACAS - stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); + stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT); #endif - sol_vector_name[vec_num].push_back(solution_vector[vec_num][i]); - sol_index[vec_num].push_back(len); - - } else { - TEUCHOS_TEST_FOR_EXCEPTION( - true, - std::logic_error, - "Error in input file: specification of solution vector layout is " - "incorrect." - << std::endl); + sol_vector_name[vec_num].push_back(solution_vector[vec_num][i]); + sol_index[vec_num].push_back(len); + + } else { + TEUCHOS_TEST_FOR_EXCEPTION( + true, + std::logic_error, + "Error in input file: specification of solution vector layout is " + "incorrect." + << std::endl); + } } + TEUCHOS_TEST_FOR_EXCEPTION( + accum != neq_, + std::logic_error, + "Error in input file: specification of solution vector layout is " + "incorrect." + << std::endl); } - TEUCHOS_TEST_FOR_EXCEPTION( - accum != neq_, - std::logic_error, - "Error in input file: specification of solution vector layout is " - "incorrect." - << std::endl); } } - // Do the coordinates - this->coordinates_field = - &metaData_->declare_field(stk::topology::NODE_RANK, "coordinates"); - stk::mesh::put_field_on_mesh( - *this->coordinates_field, metaData_->universal_part(), numDim_, nullptr); -#ifdef ALBANY_SEACAS - stk::io::set_field_role(*this->coordinates_field, Ioss::Field::MESH); -#endif - - if (numDim_ == 3) { - this->coordinates_field3d = this->coordinates_field; - } else { - this->coordinates_field3d = &metaData_->declare_field( - stk::topology::NODE_RANK, "coordinates3d"); - stk::mesh::put_field_on_mesh( - *this->coordinates_field3d, metaData_->universal_part(), 3, nullptr); -#ifdef ALBANY_SEACAS - if (params_->get("Export 3d coordinates field", false)) { - stk::io::set_field_role( - *this->coordinates_field3d, Ioss::Field::TRANSIENT); - } -#endif - } - this->addStateStructs(sis); - - //initializeProcRankField(); - } void MultiSTKFieldContainer:: @@ -451,25 +423,9 @@ fillVectorImpl (Thyra_Vector& field_vector, auto* raw_field = metaData->get_field(field_entity_rank, field_name); ALBANY_EXPECT (raw_field != nullptr, "Error! Something went wrong while retrieving a field.\n"); - const int rank = raw_field->field_array_rank(); - - using SFT = typename AbstractSTKFieldContainer::STKFieldType; - using Helper = STKFieldContainerHelper; - if (rank == 0) { - - const SFT* field = this->metaData->template get_field( - field_entity_rank, field_name); - Helper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); - } else if (rank == 1) { - const SFT* field = this->metaData->template get_field( - field_entity_rank, field_name); - Helper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); - } else { - TEUCHOS_TEST_FOR_EXCEPTION( - true, - std::runtime_error, - "Error! Only scalar and vector fields supported so far.\n"); - } + + const auto* field = this->metaData->template get_field(field_entity_rank, field_name); + STKFieldContainerHelper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); } void MultiSTKFieldContainer:: @@ -501,22 +457,9 @@ saveVectorImpl (const Thyra_Vector& field_vector, auto* raw_field = metaData->get_field(field_entity_rank, field_name); ALBANY_EXPECT (raw_field != nullptr, "Error! Something went wrong while retrieving a field.\n"); - const int rank = raw_field->field_array_rank(); - - using SFT = typename AbstractSTKFieldContainer::STKFieldType; - using Helper = STKFieldContainerHelper; - if (rank == 0) { - SFT* field = this->metaData->template get_field( - stk::topology::NODE_RANK, field_name); - Helper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); - } else if (rank == 1) { - SFT* field = this->metaData->template get_field( - stk::topology::NODE_RANK, field_name); - Helper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); - } else { - TEUCHOS_TEST_FOR_EXCEPTION (true, std::runtime_error, - "Error! Only scalar and vector fields supported so far.\n"); - } + + auto* field = this->metaData->template get_field(field_entity_rank, field_name); + STKFieldContainerHelper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped, components); } } // namespace Albany diff --git a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp index 1a3f3adc8a..63591ecc16 100644 --- a/src/disc/stk/Albany_MultiSTKFieldContainer.hpp +++ b/src/disc/stk/Albany_MultiSTKFieldContainer.hpp @@ -20,7 +20,6 @@ class MultiSTKFieldContainer : public GenericSTKFieldContainer const Teuchos::RCP& metaData_, const Teuchos::RCP& bulkData_, const int numDim_, - const Teuchos::RCP& sis, const int num_params); MultiSTKFieldContainer( diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp index df36e3aafe..1043cb711a 100644 --- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp +++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp @@ -45,7 +45,6 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( const Teuchos::RCP& metaData_, const Teuchos::RCP& bulkData_, const int numDim_, - const Teuchos::RCP& sis, const int num_params_) : GenericSTKFieldContainer( params_, @@ -60,13 +59,11 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( #endif // Start STK stuff - this->coordinates_field = - metaData_->get_field(stk::topology::NODE_RANK, "coordinates"); - + this->coordinates_field = metaData_->get_field(stk::topology::NODE_RANK, "coordinates"); + //STK throws when declaring a field that has been already declared if(this->coordinates_field == nullptr) { - this->coordinates_field = - &metaData_->declare_field(stk::topology::NODE_RANK, "coordinates"); + this->coordinates_field = &metaData_->declare_field(stk::topology::NODE_RANK, "coordinates"); } stk::mesh::put_field_on_mesh( *this->coordinates_field, metaData_->universal_part(), numDim_, nullptr); @@ -88,10 +85,7 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( #endif } - this->addStateStructs(sis); - initializeProcRankField(); - } OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( @@ -117,10 +111,9 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( #endif //IKT FIXME? - currently won't write dxdp to output file if problem is steady, //as this output doesn't work in same way. May want to change in the future. - bool output_sens_field = false; const auto& sens_method = params_->get("Sensitivity Method","NONE"); - if (this->num_params > 0 && num_time_deriv > 0 && sens_method != "None") output_sens_field = true; + output_sens_field = this->num_params > 0 && num_time_deriv > 0 && sens_method != "None"; //Create tag and id arrays for sensitivity field (dxdp or dgdp) std::vector sol_sens_tag_name_vec; @@ -145,40 +138,41 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( sol_sens_id_name_vec[0] = "sensitivity dg" + std::to_string(resp_fn_index) + "dp" + std::to_string(param_sens_index); } - solution_field.resize(num_time_deriv + 1); - solution_field_dtk.resize(num_time_deriv + 1); - solution_field_dxdp.resize(this->num_params); - - for (int num_vecs = 0; num_vecs <= num_time_deriv; num_vecs++) { - solution_field[num_vecs] = &metaData_->declare_field( - stk::topology::NODE_RANK, - params_->get( - sol_tag_name[num_vecs], sol_id_name[num_vecs])); - stk::mesh::put_field_on_mesh( - *solution_field[num_vecs], metaData_->universal_part(), neq_, nullptr); // KL: this + if (save_solution_field) { + solution_field.resize(num_time_deriv + 1); + solution_field_dtk.resize(num_time_deriv + 1); + solution_field_dxdp.resize(this->num_params); -#if defined(ALBANY_DTK) - if (output_dtk_field == true) { - solution_field_dtk[num_vecs] = &metaData_->declare_field( + for (int num_vecs = 0; num_vecs <= num_time_deriv; num_vecs++) { + solution_field[num_vecs] = &metaData_->declare_field( stk::topology::NODE_RANK, params_->get( - sol_dtk_tag_name[num_vecs], sol_dtk_id_name[num_vecs])); + sol_tag_name[num_vecs], sol_id_name[num_vecs])); stk::mesh::put_field_on_mesh( - *solution_field_dtk[num_vecs], - metaData_->universal_part(), - neq_, - nullptr); - } + *solution_field[num_vecs], metaData_->universal_part(), neq_, nullptr); // KL: this +#if defined(ALBANY_DTK) + if (output_dtk_field == true) { + solution_field_dtk[num_vecs] = &metaData_->declare_field( + stk::topology::NODE_RANK, + params_->get( + sol_dtk_tag_name[num_vecs], sol_dtk_id_name[num_vecs])); + stk::mesh::put_field_on_mesh( + *solution_field_dtk[num_vecs], + metaData_->universal_part(), + neq_, + nullptr); + } #endif #ifdef ALBANY_SEACAS - stk::io::set_field_role(*solution_field[num_vecs], Ioss::Field::TRANSIENT); + stk::io::set_field_role(*solution_field[num_vecs], Ioss::Field::TRANSIENT); #if defined(ALBANY_DTK) - if (output_dtk_field == true) - stk::io::set_field_role( - *solution_field_dtk[num_vecs], Ioss::Field::TRANSIENT); + if (output_dtk_field == true) + stk::io::set_field_role( + *solution_field_dtk[num_vecs], Ioss::Field::TRANSIENT); #endif #endif + } } //Transient sensitivities output to Exodus @@ -203,9 +197,6 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer( } this->addStateStructs(sis); - - //initializeProcRankField(); - } void @@ -270,9 +261,11 @@ saveSolnVector (const Thyra_Vector& soln, const dof_mgr_ptr_t& sol_dof_mgr, const bool overlapped) { - saveVectorImpl (soln, solution_field[0]->name(), sol_dof_mgr, overlapped); + if (save_solution_field) { + saveVectorImpl (soln, solution_field[0]->name(), sol_dof_mgr, overlapped); + } - if (soln_dxdp != Teuchos::null) { + if (soln_dxdp != Teuchos::null and output_sens_field) { TEUCHOS_TEST_FOR_EXCEPTION( soln_dxdp->domain()->dim() != this->num_params, std::runtime_error, "Error in saveSolnVector! Wrong number of vectors in soln_dxdp.\n" @@ -293,7 +286,9 @@ saveSolnVector (const Thyra_Vector& soln, const bool overlapped) { saveSolnVector(soln,soln_dxdp, sol_dof_mgr, overlapped); - saveVectorImpl (soln_dot, solution_field[1]->name(), sol_dof_mgr, overlapped); + if (save_solution_field) { + saveVectorImpl (soln_dot, solution_field[1]->name(), sol_dof_mgr, overlapped); + } } void OrdinarySTKFieldContainer:: @@ -305,7 +300,9 @@ saveSolnVector (const Thyra_Vector& soln, const bool overlapped) { saveSolnVector(soln, soln_dxdp, soln_dot, sol_dof_mgr, overlapped); - saveVectorImpl (soln_dotdot, solution_field[2]->name(), sol_dof_mgr, overlapped); + if (save_solution_field) { + saveVectorImpl (soln_dotdot, solution_field[2]->name(), sol_dof_mgr, overlapped); + } } void OrdinarySTKFieldContainer:: @@ -318,7 +315,7 @@ saveSolnMultiVector (const Thyra_MultiVector& soln, saveVectorImpl (*soln.col(icomp), solution_field[icomp]->name(), sol_dof_mgr, overlapped); } - if (soln_dxdp != Teuchos::null) { + if (soln_dxdp != Teuchos::null and output_sens_field) { TEUCHOS_TEST_FOR_EXCEPTION( soln_dxdp->domain()->dim() != this->num_params, std::runtime_error, "Error in saveSolnVector! Wrong number of vectors in soln_dxdp.\n" @@ -345,9 +342,7 @@ transferSolutionToCoords() TEUCHOS_TEST_FOR_EXCEPTION (!this->solutionFieldContainer, std::logic_error, "Error OrdinarySTKFieldContainer::transferSolutionToCoords not called from a solution field container.\n"); - using SFT = typename AbstractSTKFieldContainer::STKFieldType; - using Helper = STKFieldContainerHelper; - Helper::copySTKField(*solution_field[0], *this->coordinates_field); + STKFieldContainerHelper::copySTKField(*solution_field[0], *this->coordinates_field); } void OrdinarySTKFieldContainer:: @@ -378,25 +373,9 @@ fillVectorImpl(Thyra_Vector& field_vector, auto* raw_field = metaData->get_field(field_entity_rank, field_name); ALBANY_EXPECT (raw_field != nullptr, "Error! Something went wrong while retrieving a field.\n"); - const int rank = raw_field->field_array_rank(); - - using SFT = typename AbstractSTKFieldContainer::STKFieldType; - using Helper = STKFieldContainerHelper; - - if (rank == 0) { - const SFT* field = this->metaData->template get_field( - field_entity_rank, field_name); - Helper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); - } else if (rank == 1) { - const SFT* field = this->metaData->template get_field( - field_entity_rank, field_name); - Helper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); - } else { - TEUCHOS_TEST_FOR_EXCEPTION( - true, - std::runtime_error, - "Error! Only scalar and vector fields supported so far.\n"); - } + + const auto* field = this->metaData->template get_field(field_entity_rank, field_name); + STKFieldContainerHelper::fillVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); } void OrdinarySTKFieldContainer:: @@ -427,25 +406,9 @@ saveVectorImpl (const Thyra_Vector& field_vector, auto* raw_field = metaData->get_field(field_entity_rank, field_name); ALBANY_EXPECT (raw_field != nullptr, "Error! Something went wrong while retrieving a field.\n"); - const int rank = raw_field->field_array_rank(); - - using SFT = typename AbstractSTKFieldContainer::STKFieldType; - using Helper = STKFieldContainerHelper; - - if (rank == 0) { - SFT* field = this->metaData->template get_field( - stk::topology::NODE_RANK, field_name); - Helper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); - } else if (rank == 1) { - SFT* field = this->metaData->template get_field( - stk::topology::NODE_RANK, field_name); - Helper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); - } else { - TEUCHOS_TEST_FOR_EXCEPTION( - true, - std::runtime_error, - "Error! Only scalar and vector fields supported so far.\n"); - } + + auto* field = this->metaData->template get_field(field_entity_rank, field_name); + STKFieldContainerHelper::saveVector(field_vector, *field, *this->bulkData, field_dof_mgr, overlapped); } } // namespace Albany diff --git a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp index d9a32f8bb6..5db1ae9d3c 100644 --- a/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp +++ b/src/disc/stk/Albany_OrdinarySTKFieldContainer.hpp @@ -19,7 +19,6 @@ class OrdinarySTKFieldContainer : public GenericSTKFieldContainer const Teuchos::RCP& metaData_, const Teuchos::RCP& bulkData_, const int numDim_, - const Teuchos::RCP& sis, const int num_params); OrdinarySTKFieldContainer( @@ -132,6 +131,8 @@ class OrdinarySTKFieldContainer : public GenericSTKFieldContainer Teuchos::Array solution_field_dxdp; AbstractSTKFieldContainer::STKFieldType* residual_field; + + bool output_sens_field = false; }; } // namespace Albany diff --git a/src/disc/stk/Albany_STK3DPointStruct.cpp b/src/disc/stk/Albany_STK3DPointStruct.cpp index fd814c2631..af6d3d3f68 100644 --- a/src/disc/stk/Albany_STK3DPointStruct.cpp +++ b/src/disc/stk/Albany_STK3DPointStruct.cpp @@ -13,10 +13,11 @@ namespace Albany { //Constructor for meshes read from ASCII file -STK3DPointStruct::STK3DPointStruct(const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams) : - GenericSTKMeshStruct(params, 3, numParams) +STK3DPointStruct:: +STK3DPointStruct(const Teuchos::RCP& params, + const Teuchos::RCP& comm, + const int numParams) + : GenericSTKMeshStruct(params, 3, numParams) { std::cout << "---3DPoint constructor---" << std::endl; partVec.push_back(&metaData->declare_part_with_topology("Block0", stk::topology::PARTICLE)); @@ -43,37 +44,16 @@ STK3DPointStruct::STK3DPointStruct(const Teuchos::RCP& p std::cout << "---3DPoint constructor done---" << std::endl; // Create a mesh specs object for EACH side set - this->initializeSideSetMeshSpecs(commT); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(commT); + this->initializeSideSetMeshSpecs(comm); } -STK3DPointStruct::~STK3DPointStruct() {}; - -void -STK3DPointStruct::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - std::cout << "---3DPoint::setFieldData---" << std::endl; - SetupFieldData(commT, sis, worksetSize); - this->setSideSetFieldData(commT, side_set_sis, worksetSize); -} - -void -STK3DPointStruct::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void STK3DPointStruct:: +setBulkData (const Teuchos::RCP& comm) { std::cout << "---3DPoint::setBulkData---" << std::endl; metaData->commit(); bulkData->modification_begin(); // Begin modifying the mesh - //TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::buildMesh(commT); + //TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::buildMesh(comm); stk::mesh::PartVector nodePartVec; stk::mesh::PartVector singlePartVec(1); singlePartVec[0] = partVec[0]; // Get the element block part to put the element in. @@ -87,12 +67,11 @@ STK3DPointStruct::setBulkData( bulkData->modification_end(); - fieldAndBulkDataSet = true; - this->setSideSetBulkData(commT, side_set_sis, worksetSize); + m_bulk_data_set = true; } void -STK3DPointStruct::buildMesh(const Teuchos::RCP& /* commT */) +STK3DPointStruct::buildMesh(const Teuchos::RCP& /* comm */) { std::cout << "---3DPoint::buildMesh---" << std::endl; } diff --git a/src/disc/stk/Albany_STK3DPointStruct.hpp b/src/disc/stk/Albany_STK3DPointStruct.hpp index 9c10a85803..c7dae838af 100644 --- a/src/disc/stk/Albany_STK3DPointStruct.hpp +++ b/src/disc/stk/Albany_STK3DPointStruct.hpp @@ -22,21 +22,12 @@ namespace Albany { //! Default constructor STK3DPointStruct(const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams); + const Teuchos::RCP& comm, + const int numParams); - ~STK3DPointStruct(); + ~STK3DPointStruct() = default; - //! Sets mesh generation parameters - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); - - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return false; } @@ -47,7 +38,7 @@ namespace Albany { private: //! Build the mesh - void buildMesh(const Teuchos::RCP& commT); + void buildMesh(const Teuchos::RCP& comm); //! Build a parameter list that contains valid input parameters Teuchos::RCP diff --git a/src/disc/stk/Albany_STKDiscretization.cpp b/src/disc/stk/Albany_STKDiscretization.cpp index ae0626a522..1a76d8fff4 100644 --- a/src/disc/stk/Albany_STKDiscretization.cpp +++ b/src/disc/stk/Albany_STKDiscretization.cpp @@ -8,7 +8,6 @@ #include #include -#include "Albany_BucketArray.hpp" #include "Albany_Macros.hpp" #include "Albany_NodalGraphUtils.hpp" #include "Albany_STKDiscretization.hpp" @@ -77,7 +76,8 @@ STKDiscretization::STKDiscretization( { if (stkMeshStruct->sideSetMeshStructs.size() > 0) { for (auto it : stkMeshStruct->sideSetMeshStructs) { - auto side_disc = Teuchos::rcp(new STKDiscretization(discParams, neq, it.second, comm)); + auto stk_mesh = Teuchos::rcp_dynamic_cast(it.second,true); + auto side_disc = Teuchos::rcp(new STKDiscretization(discParams, neq, stk_mesh, comm)); sideSetDiscretizations.insert(std::make_pair(it.first, side_disc)); sideSetDiscretizationsSTK.insert(std::make_pair(it.first, side_disc)); } @@ -348,8 +348,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -379,8 +378,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -411,8 +409,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -431,8 +428,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -451,8 +447,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -476,8 +471,7 @@ STKDiscretization::transformMesh() stkMeshStruct->PBCStruct.scale[0] *= L; stkMeshStruct->PBCStruct.scale[1] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -507,8 +501,7 @@ STKDiscretization::transformMesh() #endif stkMeshStruct->PBCStruct.scale[0] *= L; stk::mesh::Field* surfaceHeight_field = - metaData->get_field>( - stk::topology::NODE_RANK, "surface_height"); + metaData->get_field(stk::topology::NODE_RANK, "surface_height"); const auto numOverlapNodes = getLocalSubdim(getOverlapNodeVectorSpace()); for (int i = 0; i < numOverlapNodes; i++) { double* x = stk::mesh::field_data(*coordinates_field, overlapnodes[i]); @@ -1319,7 +1312,7 @@ STKDiscretization::computeWorksetInfo() } else { for (int i = 0; i < numBuckets; ++i) { m_wsPhysIndex[i] = - stkMeshStruct->getMeshSpecs()[0]->ebNameToIndex[m_wsEBNames[i]]; + stkMeshStruct->meshSpecs[0]->ebNameToIndex[m_wsEBNames[i]]; } } @@ -1333,10 +1326,6 @@ STKDiscretization::computeWorksetInfo() // Clear map if remeshing if (!elemGIDws.empty()) { elemGIDws.clear(); } - typedef stk::mesh::Cartesian NodeTag; - typedef stk::mesh::Cartesian ElemTag; - typedef stk::mesh::Cartesian CompTag; - for (int b = 0; b < numBuckets; b++) { stk::mesh::Bucket& buck = *buckets[b]; m_ws_elem_coords[b].resize(buck.size()); @@ -1635,7 +1624,7 @@ STKDiscretization::computeSideSets() // Save the index of the element block that this elem lives in sStruct.elem_ebIndex = - stkMeshStruct->getMeshSpecs()[0]->ebNameToIndex[m_wsEBNames[workset]]; + 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]; @@ -1800,7 +1789,7 @@ STKDiscretization::computeSideSets() 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->getMeshSpecs()[0]->ctd)); + 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(); @@ -1837,7 +1826,7 @@ STKDiscretization::computeSideSets() // 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->getMeshSpecs()[0]->ctd; + auto ctd = stkMeshStruct->meshSpecs[0]->ctd; // Ensure we have ONE cell per layer. const auto topo_hexa = shards::getCellTopologyData>(); diff --git a/src/disc/stk/Albany_STKFieldContainerHelper.cpp b/src/disc/stk/Albany_STKFieldContainerHelper.cpp index 7402b47b94..b463bdd9fa 100644 --- a/src/disc/stk/Albany_STKFieldContainerHelper.cpp +++ b/src/disc/stk/Albany_STKFieldContainerHelper.cpp @@ -5,11 +5,227 @@ //*****************************************************************// #include "Albany_STKFieldContainerHelper.hpp" -#include "Albany_STKFieldContainerHelper_Def.hpp" -#include "Albany_AbstractSTKFieldContainer.hpp" +#include "Albany_STKUtils.hpp" +#include "Albany_GlobalLocalIndexer.hpp" +#include "Albany_ThyraUtils.hpp" -namespace Albany { +#include +#include +#include + +namespace Albany +{ + +// Fill the result vector +// Create a multidimensional array view of the +// solution field data for this bucket of nodes. +// The array is two dimensional (NumberNodes, FieldDim) +// where FieldDim=1 for scalar field, and 1+ for vector fields +void STKFieldContainerHelper:: +fillVector ( Thyra_Vector& field_thyra, + const FieldType& field_stk, + const stk::mesh::BulkData& bulkData, + const Teuchos::RCP& dof_mgr, + const bool overlapped, + const std::vector& components) +{ + constexpr int rank = FieldRank::n; + static_assert(rank==0 || rank==1, + "Error! Can only handle Scalar and Vector fields for now.\n"); + + using ScalarT = typename FieldScalar::type; + using ViewT = Kokkos::View; + + const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host(); + const bool restricted = dof_mgr->part_name()!=dof_mgr->elem_block_name(); + + auto data = getNonconstLocalData(field_thyra); + constexpr auto ELEM_RANK = stk::topology::ELEM_RANK; + const auto& elems = dof_mgr->getAlbanyConnManager()->getElementsInBlock(); + const int num_elems = elems.size(); + const auto indexer = dof_mgr->indexer(); + +#ifdef ALBANY_DEBUG + // Safety check + if (elems.size()>0) { + const auto& e0 = bulkData.get_entity(ELEM_RANK,elems[0]+1); + const auto& n0 = *bulkData.begin_nodes(e0); + TEUCHOS_TEST_FOR_EXCEPTION (stk::mesh::field_scalars_per_entity(field_stk,n0) const std::vector& + { + return dof_mgr->getGIDFieldOffsets(eq); + }; + for (int ielem=0; ielemgetElementGIDs(ielem); + for (int i=0; igetLocalElement(gids[get_offsets(0)[i]]); + if (owned_lid<0) { + continue; + } + } + auto stk_data = stk::mesh::field_data(field_stk,nodes[i]); + for (auto fid : components) { + const auto& offsets = get_offsets(fid); + const auto lid = elem_dof_lids(ielem,offsets[i]); + if (!restricted ||lid>=0) { + data[lid] = stk_data[fid]; + } + } + } + } +} + +// Shortcut for 'get all fields' +void STKFieldContainerHelper:: +fillVector ( Thyra_Vector& field_thyra, + const FieldType& field_stk, + const stk::mesh::BulkData& bulkData, + const Teuchos::RCP& dof_mgr, + const bool overlapped) +{ + std::vector components(dof_mgr->getNumFields()); + std::iota(components.begin(),components.end(),0); + fillVector(field_thyra,field_stk,bulkData,dof_mgr,overlapped,components); +} + +void STKFieldContainerHelper:: +saveVector(const Thyra_Vector& field_thyra, + FieldType& field_stk, + const stk::mesh::BulkData& bulkData, + const Teuchos::RCP& dof_mgr, + const bool overlapped, + const std::vector& components) +{ + constexpr int rank = FieldRank::n; + static_assert(rank==0 || rank==1, + "Error! Can only handle scalar and vector fields for now.\n"); + + using ScalarT = typename FieldScalar::type; + using ViewT = Kokkos::View; + + const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host(); + const bool restricted = dof_mgr->part_name()!=dof_mgr->elem_block_name(); + + auto data = getLocalData(field_thyra); + constexpr auto ELEM_RANK = stk::topology::ELEM_RANK; + const auto& elems = dof_mgr->getAlbanyConnManager()->getElementsInBlock(); + const int num_elems = elems.size(); + const auto indexer = dof_mgr->indexer(); + +#ifdef ALBANY_DEBUG + // Safety check + if (elems.size()>0) { + const auto& e0 = bulkData.get_entity(ELEM_RANK,elems[0]+1); + const auto& n0 = *bulkData.begin_nodes(e0); + TEUCHOS_TEST_FOR_EXCEPTION (stk::mesh::field_scalars_per_entity(field_stk,n0) const std::vector& + { + return dof_mgr->getGIDFieldOffsets(eq); + }; + for (int ielem=0; ielemgetElementGIDs(ielem); + for (int i=0; igetLocalElement(gids[get_offsets(0)[i]]); + if (owned_lid<0) { + continue; + } + } + auto stk_data = stk::mesh::field_data(field_stk,nodes[i]); + for (auto fid : components) { + const auto& offsets = get_offsets(fid); + const auto lid = elem_dof_lids(ielem,offsets[i]); + if (!restricted or lid>=0) { + stk_data[fid] = data[lid]; + } + } + } + } +} + +void STKFieldContainerHelper:: +copySTKField(const FieldType& source, + FieldType& target) +{ + constexpr int rank = FieldRank::n; + static_assert(rank==0 || rank==1, + "Error! Can only handle scalar and vector fields for now.\n"); + + using ScalarT = typename FieldScalar::type; + using SrcViewT = Kokkos::View; + using TgtViewT = Kokkos::View< ScalarT**,Kokkos::LayoutRight,Kokkos::HostSpace>; + + const stk::mesh::BulkData& mesh = source.get_mesh(); + const stk::mesh::BucketVector& bv = mesh.buckets(stk::topology::NODE_RANK); + + for (auto it : bv) { + const stk::mesh::Bucket& bucket = *it; + const int num_nodes_in_bucket = bucket.size(); + const int num_source_components = stk::mesh::field_scalars_per_entity(target, bucket); + const int num_target_components = stk::mesh::field_scalars_per_entity(target, bucket); + + const int uneven_downsampling = num_source_components % num_target_components; + TEUCHOS_TEST_FOR_EXCEPTION(uneven_downsampling, std::logic_error, + "Error in stk fields: specification of coordinate vector vs. solution layout is incorrect.\n"); + + auto src_stk_data = stk::mesh::field_data(source,bucket); + auto tgt_stk_data = stk::mesh::field_data(target,bucket); + + SrcViewT src_view (src_stk_data,num_nodes_in_bucket,num_source_components); + TgtViewT tgt_view (tgt_stk_data,num_nodes_in_bucket,num_target_components); + + for(int i=0; i& dof_mgr, + const bool overlapped) +{ + std::vector components(dof_mgr->getNumFields()); + std::iota(components.begin(),components.end(),0); + saveVector(field_thyra,field_stk,bulkData,dof_mgr,overlapped,components); +} -template struct STKFieldContainerHelper; } // namespace Albany diff --git a/src/disc/stk/Albany_STKFieldContainerHelper.hpp b/src/disc/stk/Albany_STKFieldContainerHelper.hpp index e621c15964..8ec7cb5a9d 100644 --- a/src/disc/stk/Albany_STKFieldContainerHelper.hpp +++ b/src/disc/stk/Albany_STKFieldContainerHelper.hpp @@ -17,9 +17,10 @@ namespace Albany { class GlobalLocalIndexer; // FieldType can be either scalar or vector, the code is the same. -template struct STKFieldContainerHelper { + using FieldType = stk::mesh::Field; + // Fill (aka get) and save (aka set) methods // If passing no components, we do all fields in the dof mgr diff --git a/src/disc/stk/Albany_STKFieldContainerHelper_Def.hpp b/src/disc/stk/Albany_STKFieldContainerHelper_Def.hpp deleted file mode 100644 index c757cdcf94..0000000000 --- a/src/disc/stk/Albany_STKFieldContainerHelper_Def.hpp +++ /dev/null @@ -1,236 +0,0 @@ -//*****************************************************************// -// Albany 3.0: Copyright 2016 Sandia Corporation // -// This Software is released under the BSD license detailed // -// in the file "license.txt" in the top-level Albany directory // -//*****************************************************************// - -#include "Albany_STKFieldContainerHelper.hpp" -#include "Albany_STKUtils.hpp" -#include "Albany_GlobalLocalIndexer.hpp" -#include "Albany_ThyraUtils.hpp" - -#include -#include -#include - -namespace Albany -{ - -// Fill the result vector -// Create a multidimensional array view of the -// solution field data for this bucket of nodes. -// The array is two dimensional (NumberNodes, FieldDim) -// where FieldDim=1 for scalar field, and 1+ for vector fields -template -void STKFieldContainerHelper:: -fillVector ( Thyra_Vector& field_thyra, - const FieldType& field_stk, - const stk::mesh::BulkData& bulkData, - const Teuchos::RCP& dof_mgr, - const bool overlapped, - const std::vector& components) -{ - constexpr int rank = FieldRank::n; - static_assert(rank==0 || rank==1, - "Error! Can only handle Scalar and Vector fields for now.\n"); - - using ScalarT = typename FieldScalar::type; - using ViewT = Kokkos::View; - - const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host(); - const bool restricted = dof_mgr->part_name()!=dof_mgr->elem_block_name(); - - auto data = getNonconstLocalData(field_thyra); - constexpr auto ELEM_RANK = stk::topology::ELEM_RANK; - const auto& elems = dof_mgr->getAlbanyConnManager()->getElementsInBlock(); - const int num_elems = elems.size(); - const auto indexer = dof_mgr->indexer(); - -#ifdef ALBANY_DEBUG - // Safety check - if (elems.size()>0) { - const auto& e0 = bulkData.get_entity(ELEM_RANK,elems[0]+1); - const auto& n0 = *bulkData.begin_nodes(e0); - TEUCHOS_TEST_FOR_EXCEPTION (stk::mesh::field_scalars_per_entity(field_stk,n0) const std::vector& - { - return dof_mgr->getGIDFieldOffsets(eq); - }; - for (int ielem=0; ielemgetElementGIDs(ielem); - for (int i=0; igetLocalElement(gids[get_offsets(0)[i]]); - if (owned_lid<0) { - continue; - } - } - auto stk_data = stk::mesh::field_data(field_stk,nodes[i]); - for (auto fid : components) { - const auto& offsets = get_offsets(fid); - const auto lid = elem_dof_lids(ielem,offsets[i]); - if (!restricted ||lid>=0) { - data[lid] = stk_data[fid]; - } - } - } - } -} - -// Shortcut for 'get all fields' -template -void STKFieldContainerHelper:: -fillVector ( Thyra_Vector& field_thyra, - const FieldType& field_stk, - const stk::mesh::BulkData& bulkData, - const Teuchos::RCP& dof_mgr, - const bool overlapped) -{ - std::vector components(dof_mgr->getNumFields()); - std::iota(components.begin(),components.end(),0); - fillVector(field_thyra,field_stk,bulkData,dof_mgr,overlapped,components); -} - -template -void STKFieldContainerHelper:: -saveVector(const Thyra_Vector& field_thyra, - FieldType& field_stk, - const stk::mesh::BulkData& bulkData, - const Teuchos::RCP& dof_mgr, - const bool overlapped, - const std::vector& components) -{ - constexpr int rank = FieldRank::n; - static_assert(rank==0 || rank==1, - "Error! Can only handle scalar and vector fields for now.\n"); - - using ScalarT = typename FieldScalar::type; - using ViewT = Kokkos::View; - - const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host(); - const bool restricted = dof_mgr->part_name()!=dof_mgr->elem_block_name(); - - auto data = getLocalData(field_thyra); - constexpr auto ELEM_RANK = stk::topology::ELEM_RANK; - const auto& elems = dof_mgr->getAlbanyConnManager()->getElementsInBlock(); - const int num_elems = elems.size(); - const auto indexer = dof_mgr->indexer(); - -#ifdef ALBANY_DEBUG - // Safety check - if (elems.size()>0) { - const auto& e0 = bulkData.get_entity(ELEM_RANK,elems[0]+1); - const auto& n0 = *bulkData.begin_nodes(e0); - TEUCHOS_TEST_FOR_EXCEPTION (stk::mesh::field_scalars_per_entity(field_stk,n0) const std::vector& - { - return dof_mgr->getGIDFieldOffsets(eq); - }; - for (int ielem=0; ielemgetElementGIDs(ielem); - for (int i=0; igetLocalElement(gids[get_offsets(0)[i]]); - if (owned_lid<0) { - continue; - } - } - auto stk_data = stk::mesh::field_data(field_stk,nodes[i]); - for (auto fid : components) { - const auto& offsets = get_offsets(fid); - const auto lid = elem_dof_lids(ielem,offsets[i]); - if (!restricted or lid>=0) { - stk_data[fid] = data[lid]; - } - } - } - } -} - -template -void STKFieldContainerHelper:: -copySTKField(const FieldType& source, - FieldType& target) -{ - constexpr int rank = FieldRank::n; - static_assert(rank==0 || rank==1, - "Error! Can only handle scalar and vector fields for now.\n"); - - using ScalarT = typename FieldScalar::type; - using SrcViewT = Kokkos::View; - using TgtViewT = Kokkos::View< ScalarT**,Kokkos::LayoutRight,Kokkos::HostSpace>; - - const stk::mesh::BulkData& mesh = source.get_mesh(); - const stk::mesh::BucketVector& bv = mesh.buckets(stk::topology::NODE_RANK); - - for (auto it : bv) { - const stk::mesh::Bucket& bucket = *it; - const int num_nodes_in_bucket = bucket.size(); - const int num_source_components = stk::mesh::field_scalars_per_entity(target, bucket); - const int num_target_components = stk::mesh::field_scalars_per_entity(target, bucket); - - const int uneven_downsampling = num_source_components % num_target_components; - TEUCHOS_TEST_FOR_EXCEPTION(uneven_downsampling, std::logic_error, - "Error in stk fields: specification of coordinate vector vs. solution layout is incorrect.\n"); - - auto src_stk_data = stk::mesh::field_data(source,bucket); - auto tgt_stk_data = stk::mesh::field_data(target,bucket); - - SrcViewT src_view (src_stk_data,num_nodes_in_bucket,num_source_components); - TgtViewT tgt_view (tgt_stk_data,num_nodes_in_bucket,num_target_components); - - for(int i=0; i -void STKFieldContainerHelper:: -saveVector (const Thyra_Vector& field_thyra, - FieldType& field_stk, - const stk::mesh::BulkData& bulkData, - const Teuchos::RCP& dof_mgr, - const bool overlapped) -{ - std::vector components(dof_mgr->getNumFields()); - std::iota(components.begin(),components.end(),0); - saveVector(field_thyra,field_stk,bulkData,dof_mgr,overlapped,components); -} - - -} // namespace Albany diff --git a/src/disc/stk/Albany_STKNodeFieldContainer.hpp b/src/disc/stk/Albany_STKNodeFieldContainer.hpp index 67b63ce54a..b120573d4c 100644 --- a/src/disc/stk/Albany_STKNodeFieldContainer.hpp +++ b/src/disc/stk/Albany_STKNodeFieldContainer.hpp @@ -8,7 +8,6 @@ #define ALBANY_STK_NODE_FIELD_CONTAINER_HPP #include "Albany_AbstractNodeFieldContainer.hpp" -#include "Albany_BucketArray.hpp" #include "Albany_ThyraUtils.hpp" #include "Albany_GlobalLocalIndexer.hpp" #include "Albany_StateInfoStruct.hpp" // For StateView diff --git a/src/disc/stk/Albany_SideSetSTKMeshStruct.cpp b/src/disc/stk/Albany_SideSetSTKMeshStruct.cpp index ff3393951b..b7519a19b2 100644 --- a/src/disc/stk/Albany_SideSetSTKMeshStruct.cpp +++ b/src/disc/stk/Albany_SideSetSTKMeshStruct.cpp @@ -7,6 +7,7 @@ #include #include "Albany_SideSetSTKMeshStruct.hpp" +#include "Albany_CommUtils.hpp" #include "Teuchos_RCPStdSharedPtrConversions.hpp" @@ -27,10 +28,11 @@ namespace Albany { -SideSetSTKMeshStruct::SideSetSTKMeshStruct (const MeshSpecsStruct& inputMeshSpecs, - const Teuchos::RCP& params, - const Teuchos::RCP& commT, - const int numParams) : +SideSetSTKMeshStruct:: +SideSetSTKMeshStruct (const MeshSpecsStruct& inputMeshSpecs, + const Teuchos::RCP& params, + const Teuchos::RCP& comm, + const int numParams) : GenericSTKMeshStruct(params, -1, numParams) { @@ -103,8 +105,8 @@ SideSetSTKMeshStruct::SideSetSTKMeshStruct (const MeshSpecsStruct& inputMeshSpec this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct(ctd, this->numDim, nsNames, ssNames, worksetSize, ebn, ebNameToIndex)); - const Teuchos::MpiComm* mpiComm = dynamic_cast* > (commT.get()); - stk::mesh::MeshBuilder meshBuilder = stk::mesh::MeshBuilder(*mpiComm->getRawMpiComm()); + auto mpiComm = getMpiCommFromTeuchosComm(comm); + stk::mesh::MeshBuilder meshBuilder = stk::mesh::MeshBuilder(mpiComm); meshBuilder.set_aura_option(stk::mesh::BulkData::NO_AUTO_AURA); meshBuilder.set_bucket_capacity(worksetSize); meshBuilder.set_add_fmwk_data(false); @@ -112,35 +114,22 @@ SideSetSTKMeshStruct::SideSetSTKMeshStruct (const MeshSpecsStruct& inputMeshSpec bulkData = Teuchos::rcp(bulkDataPtr.release()); } -SideSetSTKMeshStruct::~SideSetSTKMeshStruct() -{ - // Nothing to be done here -} - -void SideSetSTKMeshStruct::setParentMeshInfo (const AbstractSTKMeshStruct& parentMeshStruct_, - const std::string& sideSetName) +void SideSetSTKMeshStruct:: +setParentMeshInfo (const AbstractSTKMeshStruct& parentMeshStruct_, + const std::string& sideSetName) { + TEUCHOS_TEST_FOR_EXCEPTION (not parentMeshStruct.is_null(), std::logic_error, + "[SideSetSTKMeshStruct::setParentMeshInfo] Parent mesh was already set.\n"); parentMeshStruct = Teuchos::rcpFromRef(parentMeshStruct_); parentMeshSideSetName = sideSetName; } -void SideSetSTKMeshStruct::setFieldData ( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& /*side_set_sis*/) -{ - this->SetupFieldData(commT, sis, worksetSize); -} - -void SideSetSTKMeshStruct::setBulkData ( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /*side_set_sis*/) +void SideSetSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { - TEUCHOS_TEST_FOR_EXCEPTION (parentMeshStruct->ssPartVec.find(parentMeshSideSetName)==parentMeshStruct->ssPartVec.end(), std::logic_error, - "Error! The side set " << parentMeshSideSetName << " is not present in the input mesh.\n"); + TEUCHOS_TEST_FOR_EXCEPTION ( + parentMeshStruct->ssPartVec.count(parentMeshSideSetName)==0, std::logic_error, + "Error! The side set " << parentMeshSideSetName << " is not present in the input mesh.\n"); // Extracting the side part and updating the selector const stk::mesh::Part& ss_part = *parentMeshStruct->ssPartVec.find(parentMeshSideSetName)->second; @@ -149,11 +138,10 @@ void SideSetSTKMeshStruct::setBulkData ( const stk::mesh::MetaData& inputMetaData = *parentMeshStruct->metaData; const stk::mesh::BulkData& inputBulkData = *parentMeshStruct->bulkData; - typedef AbstractSTKFieldContainer::STKFieldType STKFieldType; - const STKFieldType& parent_coordinates_field = *parentMeshStruct->getCoordinatesField(); - const STKFieldType& parent_coordinates_field3d = *parentMeshStruct->getCoordinatesField3d(); - STKFieldType& coordinates_field = *fieldContainer->getCoordinatesField(); - STKFieldType& coordinates_field3d = *fieldContainer->getCoordinatesField3d(); + const auto& parent_coordinates_field = *parentMeshStruct->getCoordinatesField(); + const auto& parent_coordinates_field3d = *parentMeshStruct->getCoordinatesField3d(); + auto& coordinates_field = *fieldContainer->getCoordinatesField(); + auto& coordinates_field3d = *fieldContainer->getCoordinatesField3d(); // Now we can extract the entities std::vector sides, nodes; @@ -211,11 +199,13 @@ void SideSetSTKMeshStruct::setBulkData ( } } - // Loading the fields from file - this->loadRequiredInputFields (commT); - // Insertion of entities end bulkData->modification_end(); + + // Loading the fields from file + this->loadRequiredInputFields (comm); + + m_bulk_data_set = true; } Teuchos::RCP SideSetSTKMeshStruct::getValidDiscretizationParameters() const diff --git a/src/disc/stk/Albany_SideSetSTKMeshStruct.hpp b/src/disc/stk/Albany_SideSetSTKMeshStruct.hpp index 57fa9cdc7f..44e826b3eb 100644 --- a/src/disc/stk/Albany_SideSetSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_SideSetSTKMeshStruct.hpp @@ -21,17 +21,14 @@ class SideSetSTKMeshStruct : public GenericSTKMeshStruct const Teuchos::RCP& commT, const int numParams); - virtual ~SideSetSTKMeshStruct(); + virtual ~SideSetSTKMeshStruct() = default; void setFieldData (const Teuchos::RCP& comm, const Teuchos::RCP& sis, const unsigned int worksetSize, const std::map >& side_set_sis = {}); - void setBulkData (const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); void setParentMeshInfo (const AbstractSTKMeshStruct& parentMeshStruct_, const std::string& sideSetName); diff --git a/src/disc/stk/Albany_TmplSTKMeshStruct.hpp b/src/disc/stk/Albany_TmplSTKMeshStruct.hpp index bc0e4fd727..17458cf169 100644 --- a/src/disc/stk/Albany_TmplSTKMeshStruct.hpp +++ b/src/disc/stk/Albany_TmplSTKMeshStruct.hpp @@ -89,18 +89,13 @@ class TmplSTKMeshStruct : public GenericSTKMeshStruct { TmplSTKMeshStruct(const Teuchos::RCP& params, const Teuchos::RCP& commT, const int numParams); - ~TmplSTKMeshStruct() {}; + ~TmplSTKMeshStruct() = default; //! Sets mesh generation parameters - void setFieldData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis); - void setBulkData (const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis = {}); + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return false; } @@ -183,17 +178,9 @@ template<> void TmplSTKMeshStruct<1>::buildMesh(const Teuchos::RCP void TmplSTKMeshStruct<2>::buildMesh(const Teuchos::RCP& commT); template<> void TmplSTKMeshStruct<3>::buildMesh(const Teuchos::RCP& commT); -template<> void TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis); - -template<> void TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis); +template<> +void TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >:: +setBulkData (const Teuchos::RCP& comm); template<> Teuchos::RCP TmplSTKMeshStruct<0>::getValidDiscretizationParameters() const; template<> Teuchos::RCP TmplSTKMeshStruct<1>::getValidDiscretizationParameters() const; diff --git a/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp b/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp index e63f08f807..6ba8827e2d 100644 --- a/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp +++ b/src/disc/stk/Albany_TmplSTKMeshStruct_Def.hpp @@ -33,7 +33,7 @@ namespace Albany { template TmplSTKMeshStruct::TmplSTKMeshStruct( const Teuchos::RCP& params, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, const int numParams) : GenericSTKMeshStruct(params, traits_type::size, numParams), periodic_x(params->get("Periodic_x BC", false)), @@ -129,7 +129,7 @@ TmplSTKMeshStruct::TmplSTKMeshStruct( // Format and print out information about the mesh that is being generated if constexpr (Dim>0) // Avoid warning for pointless comparison in for loop - if (commT->getRank()==0){ // Not reached for 0D problems + if (comm->getRank()==0){ // Not reached for 0D problems std::cout <<"TmplSTKMeshStruct:: Creating " << Dim << "D mesh of size "; @@ -265,7 +265,7 @@ TmplSTKMeshStruct::TmplSTKMeshStruct( // so that the problem setup can know the worksetSize // Distribute the elems equally. Build total_elems elements, with nodeIDs starting at StartIndex - elem_map = Teuchos::rcp(new Tpetra_Map(total_elems, StartIndex, commT, Tpetra::GloballyDistributed)); + elem_map = Teuchos::rcp(new Tpetra_Map(total_elems, StartIndex, comm, Tpetra::GloballyDistributed)); int worksetSize = this->computeWorksetSize(worksetSizeMax, elem_map->getLocalNumElements() * (triangles ? 2 : 1)); @@ -279,13 +279,11 @@ TmplSTKMeshStruct::TmplSTKMeshStruct( } // Construct MeshSpecsStruct - { - const CellTopologyData& ctd = *elementBlockTopologies_[0].getCellTopologyData(); + const CellTopologyData& ctd = *elementBlockTopologies_[0].getCellTopologyData(); - this->meshSpecs[0] = Teuchos::rcp(new MeshSpecsStruct(ctd, numDim, - nsNames, ssNames, worksetSize, EBSpecs[0].name, - ebNameToIndex)); - } + this->meshSpecs[0] = Teuchos::rcp(new MeshSpecsStruct(ctd, numDim, + nsNames, ssNames, worksetSize, EBSpecs[0].name, + ebNameToIndex)); // Upon request, add a nodeset for each sideset if (params->get("Build Node Sets From Side Sets",false)) { @@ -293,19 +291,14 @@ TmplSTKMeshStruct::TmplSTKMeshStruct( } // Create a mesh specs object for EACH side set - this->initializeSideSetMeshSpecs(commT); - - // Initialize the requested sideset mesh struct in the mesh - this->initializeSideSetMeshStructs(commT); + this->initializeSideSetMeshSpecs(comm); } template void -TmplSTKMeshStruct::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) +TmplSTKMeshStruct:: +setFieldData (const Teuchos::RCP& comm, + const Teuchos::RCP& sis) { // Create global mesh: Dim-D structured, rectangular std::vector h_dim[traits_type::size]; @@ -325,39 +318,29 @@ TmplSTKMeshStruct::setFieldData( x[idx][i] = x[idx][i - 1] + h_dim[idx][i - 1]; // place the coordinates of the element nodes } - SetupFieldData(commT, sis, worksetSize); - this->setSideSetFieldData(commT, side_set_sis, worksetSize); + GenericSTKMeshStruct::setFieldData(comm, sis); } template -void -TmplSTKMeshStruct::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void TmplSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { metaData->commit(); - // STK bulkData->modification_begin(); // Begin modifying the mesh - buildMesh(commT); + buildMesh(comm); - // STK fix_node_sharing(*bulkData); bulkData->modification_end(); this->setDefaultCoordinates3d(); - this->loadRequiredInputFields (commT); + this->loadRequiredInputFields (comm); // Rebalance the mesh before starting the simulation if indicated - rebalanceInitialMeshT(commT); - - fieldAndBulkDataSet = true; + rebalanceInitialMesh(comm); - // Finally, setup the side set meshes (if any) - this->setSideSetBulkData(commT, side_set_sis, worksetSize); + m_bulk_data_set = true; } template @@ -419,6 +402,8 @@ EBSpecsStruct<0>::numElems(int /* i */){ template<> void EBSpecsStruct<0>::calcElemSizes(std::vector h[]){ + // For Dim=0, the elem sizes vector was not sized yet + h[0].resize(1); h[0][0] = 1.0; } @@ -518,7 +503,7 @@ EBSpecsStruct<3>::Initialize(int i, const Teuchos::RCP& // Specializations to build the mesh for each dimension template<> void -TmplSTKMeshStruct<0>::buildMesh(const Teuchos::RCP& /* commT */) +TmplSTKMeshStruct<0>::buildMesh(const Teuchos::RCP& /* comm */) { stk::mesh::PartVector nodePartVec; stk::mesh::PartVector singlePartVec(1); @@ -541,38 +526,26 @@ TmplSTKMeshStruct<0>::buildMesh(const Teuchos::RCP& /* commT template<> void -TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::setFieldData( - const Teuchos::RCP& commT, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& /*side_set_sis*/) -{ - SetupFieldData(commT, sis, worksetSize); -} - -template<> -void -TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::setBulkData( - const Teuchos::RCP& commT, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /*side_set_sis*/) +TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >:: +setBulkData(const Teuchos::RCP& comm) { metaData->commit(); // STK bulkData->modification_begin(); // Begin modifying the mesh - TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::buildMesh(commT); + TmplSTKMeshStruct<0, albany_stk_mesh_traits<0> >::buildMesh(comm); // STK fix_node_sharing(*bulkData); bulkData->modification_end(); + + m_bulk_data_set = true; } template<> void -TmplSTKMeshStruct<1>::buildMesh(const Teuchos::RCP& commT) +TmplSTKMeshStruct<1>::buildMesh(const Teuchos::RCP& comm) { stk::mesh::PartVector nodePartVec; stk::mesh::PartVector singlePartVec(1); @@ -641,7 +614,7 @@ TmplSTKMeshStruct<1>::buildMesh(const Teuchos::RCP& commT) if(proc_rank_field){ int* p_rank = stk::mesh::field_data(*proc_rank_field, elem); if(p_rank) - p_rank[0] = commT->getRank(); + p_rank[0] = comm->getRank(); } // Set node sets. There are no side sets currently with 1D problems (only 2D and 3D) @@ -669,7 +642,7 @@ TmplSTKMeshStruct<1>::buildMesh(const Teuchos::RCP& commT) template<> void -TmplSTKMeshStruct<2>::buildMesh(const Teuchos::RCP& /* commT */) +TmplSTKMeshStruct<2>::buildMesh(const Teuchos::RCP& /* comm */) { stk::mesh::PartVector nodePartVec; stk::mesh::PartVector singlePartVec(1); @@ -910,7 +883,7 @@ TmplSTKMeshStruct<2>::buildMesh(const Teuchos::RCP& /* commT template<> void -TmplSTKMeshStruct<3>::buildMesh(const Teuchos::RCP& commT) +TmplSTKMeshStruct<3>::buildMesh(const Teuchos::RCP& comm) { stk::mesh::PartVector nodePartVec; stk::mesh::PartVector singlePartVec(1); @@ -1012,7 +985,7 @@ TmplSTKMeshStruct<3>::buildMesh(const Teuchos::RCP& commT) if(proc_rank_field){ int* p_rank = stk::mesh::field_data(*proc_rank_field, elem); if(p_rank) - p_rank[0] = commT->getRank(); + p_rank[0] = comm->getRank(); } double* coord; diff --git a/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.cpp b/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.cpp old mode 100755 new mode 100644 index 91f2daa81d..5b661328c2 --- a/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.cpp +++ b/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.cpp @@ -25,7 +25,6 @@ #include #endif -//#include #include #include @@ -332,7 +331,7 @@ constructMesh(const Teuchos::RCP& comm, const Teuchos::RCP& sis, const unsigned int worksetSize) { - this->SetupFieldData(comm, sis, worksetSize); + this->setFieldData(comm, sis); metaData->commit(); diff --git a/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.hpp b/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.hpp index c9cb98a34d..05441b669a 100644 --- a/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.hpp +++ b/src/landIce/interfaceWithCISM/Albany_CismSTKMeshStruct.hpp @@ -4,14 +4,11 @@ // in the file "license.txt" in the top-level Albany directory // //*****************************************************************// - #ifndef ALBANY_CISM_STKMESHSTRUCT_HPP #define ALBANY_CISM_STKMESHSTRUCT_HPP #include "Albany_GenericSTKMeshStruct.hpp" -//#include - namespace Albany { class CismSTKMeshStruct : public GenericSTKMeshStruct { @@ -57,20 +54,7 @@ namespace Albany { ~CismSTKMeshStruct() = default; - void setFieldData( - const Teuchos::RCP& /* comm */, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /* side_set_sis */ = {}) - { - // Nothing to do here - } - - void setBulkData( - const Teuchos::RCP& /* comm */, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /* side_set_sis */ = {}) + void setBulkData (const Teuchos::RCP& /* comm */) { // Nothing to do here } @@ -93,7 +77,6 @@ namespace Albany { void setRestartDataTime(double restartT) {restartTime = restartT; } private: - //Ioss::Init::Initializer ioInit; Teuchos::RCP getValidDiscretizationParameters() const; @@ -146,9 +129,8 @@ namespace Albany { int debug_output_verbosity; void resizeVec(std::vector > &vec , const unsigned int rows , const unsigned int columns); void resizeVec(std::vector > &vec , const unsigned int rows , const unsigned int columns); - - protected: }; -} -#endif +} // namespace Albany + +#endif // ALBANY_CISM_STKMESHSTRUCT_HPP diff --git a/src/landIce/interfaceWithCISM/ali_driver.cpp b/src/landIce/interfaceWithCISM/ali_driver.cpp index 08c233ac13..2ef4f39552 100644 --- a/src/landIce/interfaceWithCISM/ali_driver.cpp +++ b/src/landIce/interfaceWithCISM/ali_driver.cpp @@ -520,7 +520,7 @@ void ali_driver_init(int /* argc */, int /* exec_mode */, AliToGlimmer * ftg_ptr albanyApp->createMeshSpecs(meshStruct); albanyApp->buildProblem(); - meshStruct->constructMesh(reducedMpiCommT, discParams, albanyApp->getStateMgr().getStateInfoStruct(), meshStruct->getMeshSpecs()[0]->worksetSize); + meshStruct->constructMesh(reducedMpiCommT, discParams, albanyApp->getStateMgr().getStateInfoStruct(), meshStruct->meshSpecs[0]->worksetSize); //Create nodeVS //global_node_id_owned_map_Ptr is 1-based, so nodeVS is 1-based diff --git a/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.cpp b/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.cpp index 56c802d0a8..5a6a6876db 100644 --- a/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.cpp +++ b/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.cpp @@ -176,24 +176,10 @@ MpasSTKMeshStruct(const Teuchos::RCP& params, ebn, ebNameToIndex)); this->initializeSideSetMeshSpecs(comm); - this->initializeSideSetMeshStructs(comm); } -void MpasSTKMeshStruct::setFieldData( - const Teuchos::RCP& comm, - const Teuchos::RCP& sis, - const unsigned int worksetSize, - const std::map >& side_set_sis) -{ - this->SetupFieldData(comm, sis, worksetSize); - this->setSideSetFieldData(comm, side_set_sis, worksetSize); -} - -void MpasSTKMeshStruct::setBulkData( - const Teuchos::RCP& comm, - const Teuchos::RCP& /* sis */, - const unsigned int worksetSize, - const std::map >& side_set_sis) +void MpasSTKMeshStruct:: +setBulkData (const Teuchos::RCP& comm) { constexpr auto LAYER = LayeredMeshOrdering::LAYER; constexpr auto COLUMN = LayeredMeshOrdering::COLUMN; @@ -406,7 +392,6 @@ void MpasSTKMeshStruct::setBulkData( this->loadRequiredInputFields (comm); - this->setSideSetBulkData(comm, side_set_sis, worksetSize); } Teuchos::RCP @@ -415,16 +400,4 @@ MpasSTKMeshStruct::getValidDiscretizationParameters() const return this->getValidGenericSTKParameters("Valid MpasSTKMeshStructParams"); } -int -MpasSTKMeshStruct::prismType(int const* prismVertexIds, int& minIndex) -{ - int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}}; - minIndex = std::min_element (prismVertexIds, prismVertexIds + 3) - prismVertexIds; - - int v1 (prismVertexIds[PrismVerticesMap[minIndex][1]]); - int v2 (prismVertexIds[PrismVerticesMap[minIndex][2]]); - - return v1 > v2; -} - } // namespace Albany diff --git a/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.hpp b/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.hpp index 1767293543..cb297e39b9 100644 --- a/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.hpp +++ b/src/landIce/interfaceWithMPAS/Albany_MpasSTKMeshStruct.hpp @@ -18,7 +18,7 @@ class MpasSTKMeshStruct : public GenericSTKMeshStruct public: MpasSTKMeshStruct(const Teuchos::RCP& params, - const Teuchos::RCP& commT, + const Teuchos::RCP& comm, const std::vector& indexToVertexID, const std::vector& vertexProcIDs, const std::vector& verticesCoords, @@ -39,17 +39,7 @@ class MpasSTKMeshStruct : public GenericSTKMeshStruct ~MpasSTKMeshStruct( ) = default; - void setFieldData( - const Teuchos::RCP& /* comm */, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /* side_set_sis */ = {}); - - void setBulkData( - const Teuchos::RCP& /* comm */, - const Teuchos::RCP& /* sis */, - const unsigned int /* worksetSize */, - const std::map >& /* side_set_sis */ = {}); + void setBulkData (const Teuchos::RCP& comm); //! Flag if solution has a restart values -- used in Init Cond bool hasRestartSolution() const {return hasRestartSol; } @@ -91,8 +81,6 @@ class MpasSTKMeshStruct : public GenericSTKMeshStruct const std::vector& iceMarginEdgesIds; int numLayers; LayeredMeshOrdering Ordering; - - int prismType(int const* prismVertexIds, int& minIndex); }; } // namespace Albany diff --git a/src/landIce/interfaceWithMPAS/Interface.cpp b/src/landIce/interfaceWithMPAS/Interface.cpp index 8537d55fdd..7c2f0dde2b 100644 --- a/src/landIce/interfaceWithMPAS/Interface.cpp +++ b/src/landIce/interfaceWithMPAS/Interface.cpp @@ -21,26 +21,29 @@ #include "LandIce_ProblemFactory.hpp" #include "LandIce_StokesFO.hpp" +#include "Albany_SideSetSTKMeshStruct.hpp" #include "Albany_MpasSTKMeshStruct.hpp" -#include "Teuchos_ParameterList.hpp" -#include "Teuchos_RCP.hpp" #include "Albany_Utils.hpp" #include "Albany_SolverFactory.hpp" -#include "Teuchos_XMLParameterListHelpers.hpp" -#include "Teuchos_StandardCatchMacros.hpp" -#include -#include -#include "Piro_PerformSolve.hpp" -#include "Albany_OrdinarySTKFieldContainer.hpp" -#include "Albany_STKDiscretization.hpp" -#include "Thyra_DetachedVectorView.hpp" -#include "Teuchos_YamlParameterListHelpers.hpp" -#include "Teuchos_StackedTimer.hpp" -#include "Teuchos_TimeMonitor.hpp" #include "Albany_GlobalLocalIndexer.hpp" #include "Albany_StringUtils.hpp" // for 'upper_case' +#include "Albany_OrdinarySTKFieldContainer.hpp" +#include "Albany_STKDiscretization.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include #ifdef ALBANY_SEACAS #include @@ -159,9 +162,9 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, const auto& basalFrictionParams = basalParams.sublist("Basal Friction Coefficient"); const auto betaType = util::upper_case(basalFrictionParams.get("Type")); - Teuchos::RCP ss_ms; - ss_ms = meshStruct->sideSetMeshStructs.at("basalside"); - betaField = ss_ms->metaData->get_field (stk::topology::NODE_RANK, "beta"); + auto ss_mesh = meshStruct->sideSetMeshStructs.at("basalside"); + auto ss_mesh_stk = Teuchos::rcp_dynamic_cast(ss_mesh,true); + betaField = ss_mesh_stk->metaData->get_field (stk::topology::NODE_RANK, "beta"); for (int j = 0; j < numVertices3D; ++j) { int ib = (ordering == 0) * (j % lVertexColumnShift) @@ -384,9 +387,9 @@ void velocity_solver_solve_fo(int nLayers, int globalVerticesStride, } } - if (Teuchos::nonnull(ss_ms) && !betaData.empty() && (betaField!=nullptr)) { + if (Teuchos::nonnull(ss_mesh_stk) && !betaData.empty() && (betaField!=nullptr)) { for(int ib = 0; ib < (int) indexToVertexID.size(); ++ib) { - stk::mesh::Entity node = ss_ms->bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[ib]); + stk::mesh::Entity node = ss_mesh_stk->bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[ib]); const double* betaVal = stk::mesh::field_data(*betaField,node); betaData[ib] = betaVal[0]; } @@ -812,6 +815,20 @@ void velocity_solver_extrude_3d_grid(int nLayers, int globalTrianglesStride, nLayers, num_params, Ordering)); } + // Create side meshes + for (const auto& ss : ss_pl.get>("Side Sets")) { + const auto ms = meshStruct->meshSpecs[0]; + auto params_ss = Teuchos::rcpFromRef(ss_pl.sublist(ss)); + params_ss->set("Number Of Time Derivatives",discretizationList->get("Number Of Time Derivatives")); + + auto ss_mesh = Teuchos::rcp(new Albany::SideSetSTKMeshStruct(*ms, params_ss, mpiComm, num_params)); + + ss_mesh->setParentMeshInfo(*meshStruct, ss); + + meshStruct->sideSetMeshStructs[ss] = ss_mesh; + } + meshStruct->createSideMeshMaps(); + albanyApp->createMeshSpecs(meshStruct); albanyApp->buildProblem(); diff --git a/tests/landIce/AsciiMeshes/Humboldt/CMakeLists.txt b/tests/landIce/AsciiMeshes/Humboldt/CMakeLists.txt index 2502877bbb..fc0b0f9d4a 100644 --- a/tests/landIce/AsciiMeshes/Humboldt/CMakeLists.txt +++ b/tests/landIce/AsciiMeshes/Humboldt/CMakeLists.txt @@ -2,24 +2,8 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.msh ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.msh COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.exo ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.exo COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.exo.4.0 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.exo.4.0 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.exo.4.1 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.exo.4.1 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.exo.4.2 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.exo.4.2 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_2d.exo.4.3 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_2d.exo.4.3 COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_contiguous_2d.exo ${CMAKE_CURRENT_BINARY_DIR}/humboldt_contiguous_2d.exo COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_contiguous_2d.exo.4.0 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_contiguous_2d.exo.4.0 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_contiguous_2d.exo.4.1 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_contiguous_2d.exo.4.1 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_contiguous_2d.exo.4.2 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_contiguous_2d.exo.4.2 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/humboldt_contiguous_2d.exo.4.3 - ${CMAKE_CURRENT_BINARY_DIR}/humboldt_contiguous_2d.exo.4.3 COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/basal_friction.ascii ${CMAKE_CURRENT_BINARY_DIR}/basal_friction.ascii COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/basal_friction_zeroed.ascii diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo index f711eb5ee6..54cd77b3e3 100644 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo and b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.0 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.0 deleted file mode 100644 index 40e50d53df..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.0 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.1 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.1 deleted file mode 100644 index b6d2d90764..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.1 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.2 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.2 deleted file mode 100644 index 820fe2c616..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.2 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.3 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.3 deleted file mode 100644 index 63a8953697..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_2d.exo.4.3 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo index 04cc157625..e11b53570d 100644 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo and b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.0 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.0 deleted file mode 100644 index aca89de216..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.0 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.1 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.1 deleted file mode 100644 index 5e27a3d664..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.1 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.2 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.2 deleted file mode 100644 index df34f951ec..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.2 and /dev/null differ diff --git a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.3 b/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.3 deleted file mode 100644 index e46306d5ae..0000000000 Binary files a/tests/landIce/AsciiMeshes/Humboldt/humboldt_contiguous_2d.exo.4.3 and /dev/null differ diff --git a/tests/landIce/Enthalpy/input_kleiner_B.yaml b/tests/landIce/Enthalpy/input_kleiner_B.yaml index f0c919e0cd..35b9d0cf27 100644 --- a/tests/landIce/Enthalpy/input_kleiner_B.yaml +++ b/tests/landIce/Enthalpy/input_kleiner_B.yaml @@ -124,7 +124,7 @@ ANONYMOUS: Field Origin: File Field Value: 0.0 Field 4: - Number Of Layers: 41 + Number Of Layers: 81 Vector Dim: 2 Field Name: velocity Field Type: Node Layered Vector diff --git a/tests/landIce/ExoMeshes/CMakeLists.txt b/tests/landIce/ExoMeshes/CMakeLists.txt index b5245d1792..5e7251a144 100644 --- a/tests/landIce/ExoMeshes/CMakeLists.txt +++ b/tests/landIce/ExoMeshes/CMakeLists.txt @@ -44,8 +44,6 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/gis20km_in.exo.4.2 ${CMAKE_CURRENT_BINARY_DIR}/gis20km_in.exo.4.2 COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/gis20km_in.exo.4.3 ${CMAKE_CURRENT_BINARY_DIR}/gis20km_in.exo.4.3 COPYONLY) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/square_quad4.jou - ${CMAKE_CURRENT_BINARY_DIR}/square_quad4.jou COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/gis_unstruct_2d.exo ${CMAKE_CURRENT_BINARY_DIR}/gis_unstruct_2d.exo COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/gis_unstruct_2d.exo.4.0 diff --git a/tests/landIce/ExoMeshes/square_quad.exo b/tests/landIce/ExoMeshes/square_quad.exo index 3bb580ac5f..ac1bdd00a8 100644 Binary files a/tests/landIce/ExoMeshes/square_quad.exo and b/tests/landIce/ExoMeshes/square_quad.exo differ diff --git a/tests/landIce/ExoMeshes/square_quad4.jou b/tests/landIce/ExoMeshes/square_quad4.jou deleted file mode 100644 index eabf8d8139..0000000000 --- a/tests/landIce/ExoMeshes/square_quad4.jou +++ /dev/null @@ -1,41 +0,0 @@ -## cubitx -## Cubit Version 12.1 -## Cubit Build 44018 -## Revised 2010-08-12 10:36:21 -0600 (Thu, 12 Aug 2010) -## Running 01/14/2013 09:47:25 AM -## Command Options: -## -warning = On -## -information = On -journal off -brick x 1 y 1 z 1 -move body 1 x 0.5 y 0.5 z 0.5 -curve 1 to 8 interval 5 -curve 9 to 12 interval 1 -block 1 surf 1 -block 1 element quad4 -mesh surf 1 - -sideset 11 curve 3 -sideset 11 name "sideset0" -sideset 12 curve 1 -sideset 12 name "sideset1" -sideset 13 curve 4 -sideset 13 name "sideset2" -sideset 14 curve 2 -sideset 14 name "sideset3" - - -nodeset 21 curve 3 -nodeset 21 name "nodeset0" -nodeset 22 curve 1 -nodeset 22 name "nodeset1" -nodeset 23 curve 4 -nodeset 23 name "nodeset2" -nodeset 24 curve 2 -nodeset 24 name "nodeset3" -nodeset 99 vertex in surf 1 with x_coord == 0 and y_coord == 0 -nodeset 99 name "nodeset99" - - - -export genesis "square_quad4.exo" overwrite diff --git a/tests/landIce/ExoMeshes/square_quad9.jou b/tests/landIce/ExoMeshes/square_quad9.jou deleted file mode 100644 index 7f53a9a83b..0000000000 --- a/tests/landIce/ExoMeshes/square_quad9.jou +++ /dev/null @@ -1,43 +0,0 @@ -## cubitx -## Cubit Version 12.1 -## Cubit Build 44018 -## Revised 2010-08-12 10:36:21 -0600 (Thu, 12 Aug 2010) -## Running 01/14/2013 09:47:25 AM -## Command Options: -## -warning = On -## -information = On -journal off -brick x 1 y 1 z 1 -move body 1 x 0.5 y 0.5 z 0.5 -curve 1 to 8 interval 160 -curve 9 to 12 interval 1 -block 1 surf 1 -block 1 element quad9 -mesh surf 1 - -sideset 11 curve 3 -sideset 11 name "sideset0" -sideset 12 curve 1 -sideset 12 name "sideset1" -sideset 13 curve 4 -sideset 13 name "sideset2" -sideset 14 curve 2 -sideset 14 name "sideset3" - - -nodeset 21 curve 3 -nodeset 21 name "nodeset0" -nodeset 22 curve 1 -nodeset 22 name "nodeset1" -nodeset 23 curve 4 -nodeset 23 name "nodeset2" -nodeset 24 curve 2 -nodeset 24 name "nodeset3" -nodeset 99 vertex in surf 1 with x_coord == 0 and y_coord == 0 -nodeset 99 name "nodeset99" - - - -export genesis "square_quad9.exo" overwrite - - diff --git a/tests/landIce/FO_GIS/CMakeLists.txt b/tests/landIce/FO_GIS/CMakeLists.txt index 32358ae465..bc0abbdc53 100644 --- a/tests/landIce/FO_GIS/CMakeLists.txt +++ b/tests/landIce/FO_GIS/CMakeLists.txt @@ -133,11 +133,17 @@ if (ALBANY_IFPACK2) endif() if (NOT ALBANY_PARALELL_EXODUS) - set (testNameHumboldtDecompMesh ${testNameRoot}_Humboldt_decompMesh) + set (testNameHumboldtDecompMesh ${testNameRoot}_Humboldt_2d_decompMesh) add_test (NAME ${testNameHumboldtDecompMesh} COMMAND ${SerialSeacasDecomp.exe} -processors ${MPIMNP} humboldt_2d.exo WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../AsciiMeshes/Humboldt) - set_tests_properties (${testNameHumboldtDecompMesh} PROPERTIES FIXTURES_SETUP humboldtMeshSetup) + set_tests_properties (${testNameHumboldtDecompMesh} PROPERTIES FIXTURES_SETUP humboldtMeshSetup2d) + + set (testNameHumboldtDecompMesh ${testNameRoot}_Humboldt_contiguous_2d_decompMesh) + add_test (NAME ${testNameHumboldtDecompMesh} + COMMAND ${SerialSeacasDecomp.exe} -processors ${MPIMNP} humboldt_contiguous_2d.exo + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../AsciiMeshes/Humboldt) + set_tests_properties (${testNameHumboldtDecompMesh} PROPERTIES FIXTURES_SETUP humboldtMeshSetupContiguous2d) endif() if (ALBANY_FROSCH) @@ -153,47 +159,44 @@ if (ALBANY_FROSCH) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_frosch_fluxdiv.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_frosch_fluxdiv.yaml) add_test(${testNameRoot}_Humboldt_FluxDiv_Tpetra_FROSch ${Albany.exe} input_fo_humboldt_frosch_fluxdiv.yaml) + set_tests_properties(${testNameRoot}_Humboldt_FluxDiv_Tpetra_FROSch + PROPERTIES + LABELS "LandIce;Tpetra;Forward") configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_frosch_power_law.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_frosch_power_law.yaml) add_test(${testNameRoot}_Humboldt_powerLaw_Tpetra_FROSch ${Albany.exe} input_fo_humboldt_frosch_power_law.yaml) + set_tests_properties(${testNameRoot}_Humboldt_powerLaw_Tpetra_FROSch + PROPERTIES + LABELS "LandIce;Tpetra;Forward") configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_frosch_effect_press.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_frosch_effect_press.yaml) add_test(${testNameRoot}_Humboldt_effectPress_Tpetra_FROSch ${Albany.exe} input_fo_humboldt_frosch_effect_press.yaml) + set_tests_properties(${testNameRoot}_Humboldt_effectPress_Tpetra_FROSch + PROPERTIES + LABELS "LandIce;Tpetra;Forward") configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_frosch_pressurized_bed.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_frosch_pressurized_bed.yaml) add_test(${testNameRoot}_Humboldt_pressurizedBed_Tpetra_FROSch ${Albany.exe} input_fo_humboldt_frosch_pressurized_bed.yaml) - + set_tests_properties(${testNameRoot}_Humboldt_pressurizedBed_Tpetra_FROSch + PROPERTIES + LABELS "LandIce;Tpetra;Forward") + if (NOT ALBANY_PARALELL_EXODUS) set_tests_properties(${testNameRoot}_Humboldt_FluxDiv_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward" - FIXTURES_REQUIRED humboldtMeshSetup) + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetupContiguous2d) + + set_tests_properties(${testNameRoot}_Humboldt_pressurizedBed_Tpetra_FROSch + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetup2d) set_tests_properties(${testNameRoot}_Humboldt_powerLaw_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward" - FIXTURES_REQUIRED humboldtMeshSetup) + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetup2d) set_tests_properties(${testNameRoot}_Humboldt_effectPress_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward" - FIXTURES_REQUIRED humboldtMeshSetup) + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetup2d) - else () - set_tests_properties(${testNameRoot}_Humboldt_FluxDiv_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward") - - set_tests_properties(${testNameRoot}_Humboldt_powerLaw_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward") - - set_tests_properties(${testNameRoot}_Humboldt_effectPress_Tpetra_FROSch - PROPERTIES - LABELS "LandIce;Tpetra;Forward") endif() endif() @@ -342,12 +345,19 @@ if (ALBANY_IFPACK2) set (testName ${testNameRoot}_Humboldt_Analysis) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_analysis.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_analysis.yaml) - add_test(${testNameRoot}_Humboldt_Analysis ${AlbanyAnalysis.exe} input_fo_humboldt_analysis.yaml) + add_test(${testName} ${AlbanyAnalysis.exe} input_fo_humboldt_analysis.yaml) set (testName ${testNameRoot}_Humboldt_Analysis_Mat_Free_Reg) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_fo_humboldt_analysis_mat_free_reg.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_fo_humboldt_analysis_mat_free_reg.yaml) - add_test(${testNameRoot}_Humboldt_Analysis_Analysis_Mat_Free_Reg ${AlbanyAnalysis.exe} input_fo_humboldt_analysis_mat_free_reg.yaml) + add_test(${testName} ${AlbanyAnalysis.exe} input_fo_humboldt_analysis_mat_free_reg.yaml) + + if (NOT ALBANY_PARALELL_EXODUS) + set_tests_properties(${testNameRoot}_Humboldt_Analysis + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetupContiguous2d) + set_tests_properties(${testNameRoot}_Humboldt_Analysis_Mat_Free_Reg + PROPERTIES FIXTURES_REQUIRED humboldtMeshSetupContiguous2d) + endif() endif() endif() diff --git a/tests/landIce/FO_GIS/input_fo_gis_populate_meshes.yaml b/tests/landIce/FO_GIS/input_fo_gis_populate_meshes.yaml index d7a81667d7..5adaa68770 100644 --- a/tests/landIce/FO_GIS/input_fo_gis_populate_meshes.yaml +++ b/tests/landIce/FO_GIS/input_fo_gis_populate_meshes.yaml @@ -20,6 +20,7 @@ ANONYMOUS: Method: Ioss Number Of Time Derivatives: 0 Use Serial Mesh: false + Save Solution Field: false Exodus Input File Name: ../ExoMeshes/gis_unstruct_2d.exo Exodus Output File Name: gis_unstruct_basal_populated.exo Required Fields Info: @@ -48,6 +49,7 @@ ANONYMOUS: upperside: Method: SideSetSTK Number Of Time Derivatives: 0 + Save Solution Field: false Exodus Output File Name: gis_unstruct_surface_populated.exo Required Fields Info: Number Of Fields: 2 diff --git a/tests/unit/Albany_UnitTestSetupHelpers.hpp b/tests/unit/Albany_UnitTestSetupHelpers.hpp index 107097a412..0ad44f56a6 100644 --- a/tests/unit/Albany_UnitTestSetupHelpers.hpp +++ b/tests/unit/Albany_UnitTestSetupHelpers.hpp @@ -65,7 +65,7 @@ createTestDisc (const Teuchos::RCP& comm, layout->dimensions(info->dim); sis->push_back(info); } - return factory.createDiscretization(neq,sis); + return factory.createDiscretization(neq,{},sis,{}); } // Helper to get topology/FEBasis/cubature given mesh dim