From f3c00a7cfa3d4608345836936cb03b4d32098f82 Mon Sep 17 00:00:00 2001 From: Alan Williams Date: Wed, 1 Apr 2020 12:23:19 -0600 Subject: [PATCH 1/2] Snapshot of stk as of 4/01/2020 --- packages/stk/CMakeLists.txt | 2 +- packages/stk/cmake/STK_Trilinos_config.h.in | 2 +- .../stk/stk_balance/stk_balance/balance.cpp | 6 +- .../stk_balance/stk_balance/balanceUtils.cpp | 15 + .../stk_balance/stk_balance/balanceUtils.hpp | 7 +- .../internal/balanceCommandLine.hpp | 11 +- .../internal/privateDeclarations.cpp | 36 +- packages/stk/stk_io/stk_io/IossBridge.cpp | 8 +- packages/stk/stk_io/stk_io/MeshField.cpp | 2 + packages/stk/stk_io/stk_io/StkIoUtils.cpp | 4 +- packages/stk/stk_io/stk_io/StkIoUtils.hpp | 4 +- .../stk/stk_mesh/stk_mesh/base/BulkData.cpp | 15 +- .../stk/stk_mesh/stk_mesh/base/BulkData.hpp | 3 - .../stk_mesh/base/DeviceMeshManager.cpp | 5 +- .../stk/stk_mesh/stk_mesh/base/EntityKey.hpp | 5 +- .../stk_mesh/stk_mesh/base/FieldParallel.cpp | 4 + .../stk_mesh/stk_mesh/base/GetEntities.cpp | 25 +- .../stk_mesh/base/HostMeshManager.cpp | 5 +- packages/stk/stk_mesh/stk_mesh/base/Ngp.hpp | 3 +- .../stk/stk_mesh/stk_mesh/base/NgpAtomics.hpp | 2 +- .../stk_mesh/base/NgpMultistateField.hpp | 58 +- packages/stk/stk_mesh/stk_mesh/base/Types.hpp | 2 + .../stk_mesh/baseImpl/EntityRepository.cpp | 255 +++++- .../stk_mesh/baseImpl/EntityRepository.hpp | 35 +- .../stk_mesh/baseImpl/MeshModification.cpp | 2 + .../stk_mesh/NgpFieldAccess.cpp | 8 +- .../stk_mesh/calculate_centroid.hpp | 2 +- .../stk/stk_simd/stk_simd/disimd/DISimd.hpp | 2 +- .../stk_simd/disimd/DISimdDoubleLoadStore.hpp | 49 +- .../stk_simd/disimd/DISimdDoubleMath.hpp | 18 +- packages/stk/stk_simd/stk_simd/disimd/avx.hpp | 52 +- .../stk/stk_simd/stk_simd/disimd/avx512.hpp | 60 +- .../stk/stk_simd/stk_simd/disimd/neon.hpp | 16 + .../stk/stk_simd/stk_simd/disimd/scalar.hpp | 7 + .../stk/stk_simd/stk_simd/disimd/simd.hpp | 11 +- .../stk_simd/stk_simd/disimd/simd_common.hpp | 33 +- packages/stk/stk_simd/stk_simd/disimd/sse.hpp | 49 +- packages/stk/stk_simd/stk_simd/disimd/vsx.hpp | 16 + .../topology_detail/meta_functions.hpp | 13 +- .../TransferCopyByIdStkMeshAdapter.hpp | 3 +- .../stk_balance/UnitTestTransientFields.cpp | 729 ++++++++++-------- .../stk_mesh/ngp/NgpMeshManagerTest.cpp | 15 +- .../stk_mesh/ngp/NgpMultistateFieldTest.cpp | 47 +- .../stk_mesh/ngp/NgpParallelSumTest.cpp | 8 +- .../stk_mesh/ngp/TestNgpMeshUpdate.cpp | 8 +- .../stk_unit_tests/stk_mesh/ngp/howToNgp.cpp | 62 +- .../stk_mesh/ngp/ngpFieldTest.cpp | 17 +- .../stk_ngp/NgpFieldPerformanceTest.cpp | 3 +- .../stk_util/diag/UnitTestTimer.cpp | 20 +- .../stk/stk_util/stk_util/stk_kokkos_macros.h | 26 +- 50 files changed, 1197 insertions(+), 593 deletions(-) diff --git a/packages/stk/CMakeLists.txt b/packages/stk/CMakeLists.txt index 51e4c3e44566..af12202feb4d 100644 --- a/packages/stk/CMakeLists.txt +++ b/packages/stk/CMakeLists.txt @@ -66,7 +66,7 @@ IF (TPL_ENABLE_CUDA) SET(HAVE_STK_CUDA ON) ENDIF() -#SET(STK_KOKKOS_SIMD ON) +SET(STK_KOKKOS_SIMD ON) SET(STK_VOLATILE_SIMD ON) IF ( ${PACKAGE_NAME}_DISABLE_MPI_NEIGHBOR_COMM ) diff --git a/packages/stk/cmake/STK_Trilinos_config.h.in b/packages/stk/cmake/STK_Trilinos_config.h.in index bb654ed71c12..838489459240 100644 --- a/packages/stk/cmake/STK_Trilinos_config.h.in +++ b/packages/stk/cmake/STK_Trilinos_config.h.in @@ -52,7 +52,7 @@ #cmakedefine HAVE_STK_CUDA -//#cmakedefine STK_KOKKOS_SIMD +#cmakedefine STK_KOKKOS_SIMD #cmakedefine STK_VOLATILE_SIMD diff --git a/packages/stk/stk_balance/stk_balance/balance.cpp b/packages/stk/stk_balance/stk_balance/balance.cpp index 4c802f75d9c2..73fb969216fc 100644 --- a/packages/stk/stk_balance/stk_balance/balance.cpp +++ b/packages/stk/stk_balance/stk_balance/balance.cpp @@ -477,7 +477,7 @@ void run_stk_balance_with_settings(const std::string& outputFilename, const std: InternalMesh mesh(comm, balanceSettings.getCoordinateFieldName()); - initial_decomp_and_balance(mesh, balanceSettings, exodusFilename, outputFilename, initialDecompMethod); + initial_decomp_and_balance(mesh, balanceSettings, exodusFilename, outputFilename, balanceSettings.getInitialDecompMethod()); } StkBalanceSettings create_balance_settings(const stk::balance::ParsedOptions & options) @@ -525,6 +525,10 @@ StkBalanceSettings create_balance_settings(const stk::balance::ParsedOptions & o balanceSettings.setDecompMethod(options.decompMethod); } + if (options.is_option_provided(stk::balance::ParsedOptions::INITIAL_DECOMP_METHOD)) { + balanceSettings.setInitialDecompMethod(options.initialDecompMethod); + } + return balanceSettings; } diff --git a/packages/stk/stk_balance/stk_balance/balanceUtils.cpp b/packages/stk/stk_balance/stk_balance/balanceUtils.cpp index ef62ca38cb1f..b7c3cbd9e874 100644 --- a/packages/stk/stk_balance/stk_balance/balanceUtils.cpp +++ b/packages/stk/stk_balance/stk_balance/balanceUtils.cpp @@ -16,6 +16,10 @@ namespace balance ////////////////////////////////////////////////////////////////////////// +BalanceSettings::BalanceSettings() + : initialDecompMethod("RIB") +{} + size_t BalanceSettings::getNumNodesRequiredForConnection(stk::topology element1Topology, stk::topology element2Topology) const { return 1; @@ -131,6 +135,16 @@ std::string BalanceSettings::getDecompMethod() const return std::string("parmetis"); } +void BalanceSettings::setInitialDecompMethod(const std::string& method) +{ + initialDecompMethod = method; +} + +std::string BalanceSettings::getInitialDecompMethod() const +{ + return initialDecompMethod; +} + std::string BalanceSettings::getCoordinateFieldName() const { return std::string("coordinates"); @@ -379,6 +393,7 @@ void GraphCreationSettings::setDecompMethod(const std::string& input_method) { method = input_method; } + void GraphCreationSettings::setToleranceForFaceSearch(double tol) { m_UseConstantToleranceForFaceSearch = true; diff --git a/packages/stk/stk_balance/stk_balance/balanceUtils.hpp b/packages/stk/stk_balance/stk_balance/balanceUtils.hpp index 15df94bf222a..f182e6c96032 100644 --- a/packages/stk/stk_balance/stk_balance/balanceUtils.hpp +++ b/packages/stk/stk_balance/stk_balance/balanceUtils.hpp @@ -77,7 +77,7 @@ class DecompositionChangeList class BalanceSettings { public: - BalanceSettings() {} + BalanceSettings(); virtual ~BalanceSettings() {} enum GraphOption @@ -123,6 +123,9 @@ class BalanceSettings virtual void setDecompMethod(const std::string& method) ; virtual std::string getDecompMethod() const ; + virtual void setInitialDecompMethod(const std::string& method) ; + virtual std::string getInitialDecompMethod() const ; + virtual std::string getCoordinateFieldName() const ; virtual bool shouldPrintMetrics() const; @@ -150,6 +153,8 @@ class BalanceSettings virtual double getNodeBalancerTargetLoadBalance() const; virtual unsigned getNodeBalancerMaxIterations() const; +private: + std::string initialDecompMethod; }; class BasicGeometricSettings : public BalanceSettings diff --git a/packages/stk/stk_balance/stk_balance/internal/balanceCommandLine.hpp b/packages/stk/stk_balance/stk_balance/internal/balanceCommandLine.hpp index b1d6624fa41b..5d04bac2c801 100644 --- a/packages/stk/stk_balance/stk_balance/internal/balanceCommandLine.hpp +++ b/packages/stk/stk_balance/stk_balance/internal/balanceCommandLine.hpp @@ -56,13 +56,15 @@ struct ParsedOptions double faceSearchRelTol; bool useContactSearch; std::string decompMethod; + std::string initialDecompMethod; enum ParserFlags { APP_TYPE = 1, FACE_SEARCH_ABS_TOL = 2, FACE_SEARCH_REL_TOL = 4, CONTACT_SEARCH = 8, - DECOMP_METHOD = 16 + DECOMP_METHOD = 16, + INITIAL_DECOMP_METHOD = 32 }; ParsedOptions() @@ -78,6 +80,7 @@ struct ParsedOptions faceSearchRelTol(0.15), useContactSearch(true), decompMethod("parmetis"), + initialDecompMethod("RIB"), m_parsedState(0) { } @@ -112,6 +115,12 @@ struct ParsedOptions m_parsedState |= DECOMP_METHOD; } + void set_initial_decomp_method(std::string method) + { + initialDecompMethod = method; + m_parsedState |= INITIAL_DECOMP_METHOD; + } + bool is_option_provided(ParserFlags parserFlags) const { return (m_parsedState & parserFlags); diff --git a/packages/stk/stk_balance/stk_balance/internal/privateDeclarations.cpp b/packages/stk/stk_balance/stk_balance/internal/privateDeclarations.cpp index 8609bef43c97..854b6b5e871a 100644 --- a/packages/stk/stk_balance/stk_balance/internal/privateDeclarations.cpp +++ b/packages/stk/stk_balance/stk_balance/internal/privateDeclarations.cpp @@ -875,43 +875,49 @@ bool shouldOmitSpiderElement(const stk::mesh::BulkData & stkMeshBulkData, return omitConnection; } -void fix_spider_elements(const BalanceSettings & balanceSettings, stk::mesh::BulkData & stkMeshBulkData) +bool found_element_at_end_of_leg(int newOwner, const stk::mesh::BulkData & bulk) { - stk::mesh::Ghosting * customAura = stk::tools::create_custom_aura(stkMeshBulkData, stkMeshBulkData.mesh_meta_data().globally_shared_part(), "customAura"); + return (newOwner < bulk.parallel_size()); +} - stk::mesh::MetaData & meta = stkMeshBulkData.mesh_meta_data(); - const stk::mesh::Field & beamConnectivityCountField = *balanceSettings.getSpiderConnectivityCountField(stkMeshBulkData); +void fix_spider_elements(const BalanceSettings & balanceSettings, stk::mesh::BulkData & bulk) +{ + stk::mesh::Ghosting * customAura = stk::tools::create_custom_aura(bulk, bulk.mesh_meta_data().globally_shared_part(), "customAura"); + + stk::mesh::MetaData & meta = bulk.mesh_meta_data(); + const stk::mesh::Field & beamConnectivityCountField = *balanceSettings.getSpiderConnectivityCountField(bulk); stk::mesh::EntityVector beams; stk::mesh::Part & beamPart = meta.get_topology_root_part(stk::topology::BEAM_2); - stk::mesh::get_selected_entities(beamPart & meta.locally_owned_part(), stkMeshBulkData.buckets(stk::topology::ELEM_RANK), beams); + stk::mesh::get_selected_entities(beamPart & meta.locally_owned_part(), bulk.buckets(stk::topology::ELEM_RANK), beams); stk::mesh::EntityProcVec beamsToMove; for (stk::mesh::Entity beam : beams) { - if (isElementPartOfSpider(stkMeshBulkData, beamConnectivityCountField, beam)) { - const stk::mesh::Entity* nodes = stkMeshBulkData.begin_nodes(beam); + if (isElementPartOfSpider(bulk, beamConnectivityCountField, beam)) { + const stk::mesh::Entity* nodes = bulk.begin_nodes(beam); const int node1ConnectivityCount = *stk::mesh::field_data(beamConnectivityCountField, nodes[0]); const int node2ConnectivityCount = *stk::mesh::field_data(beamConnectivityCountField, nodes[1]); const stk::mesh::Entity endNode = (node1ConnectivityCount < node2ConnectivityCount) ? nodes[0] : nodes[1]; - const stk::mesh::Entity* elements = stkMeshBulkData.begin_elements(endNode); - const unsigned numElements = stkMeshBulkData.num_elements(endNode); - int newOwner = stkMeshBulkData.parallel_size() - 1; + const stk::mesh::Entity* elements = bulk.begin_elements(endNode); + const unsigned numElements = bulk.num_elements(endNode); + int newOwner = std::numeric_limits::max(); + for (unsigned i = 0; i < numElements; ++i) { - if (stkMeshBulkData.bucket(elements[i]).topology() != stk::topology::BEAM_2) { - newOwner = std::min(newOwner, stkMeshBulkData.parallel_owner_rank(elements[i])); + if (bulk.bucket(elements[i]).topology() != stk::topology::BEAM_2) { + newOwner = std::min(newOwner, bulk.parallel_owner_rank(elements[i])); } } - if (newOwner != stkMeshBulkData.parallel_rank()) { + if (found_element_at_end_of_leg(newOwner, bulk) && (newOwner != bulk.parallel_rank())) { beamsToMove.push_back(std::make_pair(beam, newOwner)); } } } - stk::tools::destroy_custom_aura(stkMeshBulkData, customAura); + stk::tools::destroy_custom_aura(bulk, customAura); - stkMeshBulkData.change_entity_owner(beamsToMove); + bulk.change_entity_owner(beamsToMove); } void keep_spiders_on_original_proc(stk::mesh::BulkData &bulk, const stk::balance::BalanceSettings & balanceSettings, DecompositionChangeList &changeList) diff --git a/packages/stk/stk_io/stk_io/IossBridge.cpp b/packages/stk/stk_io/stk_io/IossBridge.cpp index 6b225446da81..f715e6ff1d5d 100644 --- a/packages/stk/stk_io/stk_io/IossBridge.cpp +++ b/packages/stk/stk_io/stk_io/IossBridge.cpp @@ -326,6 +326,8 @@ namespace { throw std::runtime_error(errmsg.str()); } + field->sync_to_host(); + field->modify_on_host(); for (size_t i=0; i < entity_count; ++i) { if (mesh.is_valid(entities[i])) { T *fld_data = static_cast(stk::mesh::field_data(*field, entities[i])); @@ -365,6 +367,8 @@ namespace { stk::mesh::MetaData &meta = stk::mesh::MetaData::get(*stk_part); stk::mesh::Selector selector = (meta.globally_shared_part() | meta.locally_owned_part()) & *stk_part; + field->sync_to_host(); + field->modify_on_host(); for (size_t i=0; i < entity_count; ++i) { if (mesh.is_valid(entities[i])) { const stk::mesh::Bucket &bucket = mesh.bucket(entities[i]); @@ -392,9 +396,10 @@ namespace { std::vector io_field_data(entity_count*field_component_count); + field->sync_to_host(); for (size_t i=0; i < entity_count; ++i) { if (mesh.is_valid(entities[i]) && mesh.entity_rank(entities[i]) == field->entity_rank()) { - T *fld_data = static_cast(stk::mesh::field_data(*field, entities[i])); + const T *fld_data = static_cast(stk::mesh::field_data(*field, entities[i])); if (fld_data != nullptr) { for(size_t j=0; jentity_rank() == stk::topology::NODE_RANK)) { + nodeFactorVar->sync_to_host(); for(auto& node : nodes) { df.push_back(*(double*) (stk::mesh::field_data(*nodeFactorVar, node))); } diff --git a/packages/stk/stk_io/stk_io/MeshField.cpp b/packages/stk/stk_io/stk_io/MeshField.cpp index b23ad836f585..1a75adcd5a23 100644 --- a/packages/stk/stk_io/stk_io/MeshField.cpp +++ b/packages/stk/stk_io/stk_io/MeshField.cpp @@ -250,6 +250,8 @@ double MeshField::restore_field_data(stk::mesh::BulkData &bulk, const stk::mesh::EntityRank rank = field_part.get_entity_rank(); stk::io::get_input_entity_list(io_entity, rank, bulk, entity_list); + m_field->sync_to_host(); + m_field->modify_on_host(); for (size_t i=0; i < entity_list.size(); ++i) { if (bulk.is_valid(entity_list[i])) { double *fld_data = static_cast(stk::mesh::field_data(*m_field, entity_list[i])); diff --git a/packages/stk/stk_io/stk_io/StkIoUtils.cpp b/packages/stk/stk_io/stk_io/StkIoUtils.cpp index 01fe2f7fb664..5a43613164e8 100644 --- a/packages/stk/stk_io/stk_io/StkIoUtils.cpp +++ b/packages/stk/stk_io/stk_io/StkIoUtils.cpp @@ -442,7 +442,7 @@ bool isSidesetSupported(const stk::mesh::BulkData &bulk, const stk::mesh::Entity } -stk::mesh::FieldVector get_transient_fields(stk::mesh::MetaData &meta) +stk::mesh::FieldVector get_transient_fields(const stk::mesh::MetaData &meta) { stk::mesh::FieldVector fields; @@ -456,7 +456,7 @@ stk::mesh::FieldVector get_transient_fields(stk::mesh::MetaData &meta) return fields; } -stk::mesh::FieldVector get_transient_fields(stk::mesh::MetaData &meta, const stk::mesh::EntityRank rank) +stk::mesh::FieldVector get_transient_fields(const stk::mesh::MetaData &meta, const stk::mesh::EntityRank rank) { stk::mesh::FieldVector fields; diff --git a/packages/stk/stk_io/stk_io/StkIoUtils.hpp b/packages/stk/stk_io/stk_io/StkIoUtils.hpp index 90017486acea..c86e9fdd6142 100644 --- a/packages/stk/stk_io/stk_io/StkIoUtils.hpp +++ b/packages/stk/stk_io/stk_io/StkIoUtils.hpp @@ -95,8 +95,8 @@ bool should_reconstruct_sideset(const stk::mesh::BulkData& bulkData, const stk:: bool isSidesetSupported(const stk::mesh::BulkData &bulk, const stk::mesh::EntityVector &sides, const stk::mesh::impl::ParallelPartInfo ¶llelPartInfo); -stk::mesh::FieldVector get_transient_fields(stk::mesh::MetaData &meta); -stk::mesh::FieldVector get_transient_fields(stk::mesh::MetaData &meta, const stk::mesh::EntityRank rank); +stk::mesh::FieldVector get_transient_fields(const stk::mesh::MetaData &meta); +stk::mesh::FieldVector get_transient_fields(const stk::mesh::MetaData &meta, const stk::mesh::EntityRank rank); const stk::mesh::Part& get_sideset_parent(const stk::mesh::Part& sidesetPart); diff --git a/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp b/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp index 98b756549c50..a6915eedad70 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp @@ -593,7 +593,7 @@ void BulkData::update_ngp_mesh() const { if (m_ngpMeshManager == nullptr) { -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH m_ngpMeshManager = new DeviceMeshManager(*this); #else m_ngpMeshManager = new HostMeshManager(*this); @@ -1446,7 +1446,7 @@ void BulkData::record_entity_deletion(Entity entity) { const EntityKey key = entity_key(entity); set_mesh_index(entity, 0, 0); - m_entity_repo->destroy_entity(key); + m_entity_repo->destroy_entity(key, entity); notifier.notify_local_entities_created_or_deleted(key.rank()); notifier.notify_local_buckets_changed(key.rank()); m_meshModification.mark_entity_as_deleted(entity.local_offset()); @@ -1651,6 +1651,8 @@ void BulkData::declare_entities(stk::topology::rank_t rank, const IDVECTOR& newI this->internal_set_owner(new_entity, parallel_rank()); } + //now clear the created-entity cache... + begin_entities(rank); } template @@ -2180,6 +2182,8 @@ void BulkData::update_field_data_states(FieldBase* field) for(int state=0; state(state); field_set[state] = field->field_state(state_identifier); + field_set[state]->sync_to_host(); + field_set[state]->modify_on_host(); } int outer_idx = 0; @@ -2205,8 +2209,11 @@ void BulkData::update_field_data_states() const int num_state = field.number_of_states(); i += num_state ; - - + for(int s=0; s < num_state; ++s) + { + field_set[outer_idx+s]->sync_to_host(); + field_set[outer_idx+s]->modify_on_host(); + } if (num_state > 1) { for ( int b = 0, be = field_set[outer_idx]->get_meta_data_for_field().size(); b < be; ++b) { diff --git a/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp b/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp index 2cb3a3fcfb7d..9fb08c9d799f 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp @@ -109,9 +109,6 @@ enum class FaceCreationBehavior; namespace stk { namespace mesh { -typedef std::unordered_map::const_iterator const_entity_iterator; -typedef std::unordered_map::iterator entity_iterator; - void communicate_field_data(const Ghosting & ghosts, const std::vector & fields); void communicate_field_data(const BulkData & mesh, const std::vector & fields); void parallel_sum_including_ghosts(const BulkData & mesh, const std::vector & fields); diff --git a/packages/stk/stk_mesh/stk_mesh/base/DeviceMeshManager.cpp b/packages/stk/stk_mesh/stk_mesh/base/DeviceMeshManager.cpp index 25653c20bb9c..8b6d246849e0 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/DeviceMeshManager.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/DeviceMeshManager.cpp @@ -1,6 +1,7 @@ #include "DeviceMeshManager.hpp" #include "BulkData.hpp" -#include "stk_mesh/base/NgpMesh.hpp" +#include "NgpMesh.hpp" +#include "stk_util/stk_kokkos_macros.h" namespace stk { namespace mesh { @@ -18,7 +19,7 @@ DeviceMeshManager::~DeviceMeshManager() stk::mesh::NgpMesh & DeviceMeshManager::get_mesh() { -#ifndef KOKKOS_ENABLE_CUDA +#ifndef STK_USE_DEVICE_MESH ThrowErrorMsg("DeviceMeshManager does not have host mesh"); static HostMesh hostMesh; return hostMesh; // Never used; Avoids build problems on the CPU diff --git a/packages/stk/stk_mesh/stk_mesh/base/EntityKey.hpp b/packages/stk/stk_mesh/stk_mesh/base/EntityKey.hpp index def10900cd7f..9ee23f6c9ea8 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/EntityKey.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/EntityKey.hpp @@ -38,7 +38,6 @@ #include // for size_t #include // for uint64_t #include // for ostream -#include #include // for EntityId, EntityRank #include // for ThrowAssertMsg @@ -47,7 +46,7 @@ namespace mesh { struct EntityKey { - enum entity_key_t : uint64_t { + enum entity_key_t { RANK_SHIFT = 56ULL , MIN_ID = 0ULL , MAX_ID = (1ULL << RANK_SHIFT) - 1ULL @@ -101,7 +100,7 @@ struct HashValueForEntityKey STK_FUNCTION size_t operator()(EntityKey k) const { - return std::hash{}(static_cast(k)); + return static_cast(k.m_value); } }; diff --git a/packages/stk/stk_mesh/stk_mesh/base/FieldParallel.cpp b/packages/stk/stk_mesh/stk_mesh/base/FieldParallel.cpp index f7cc44500b96..016af3d1d0e7 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/FieldParallel.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/FieldParallel.cpp @@ -238,6 +238,7 @@ void parallel_op_impl(const BulkData& mesh, std::vector fields ThrowRequireMsg(f.type_is(), "Please don't mix fields with different primitive types in the same parallel assemble operation"); + f.sync_to_host(); const BucketIndices& bktIndices = mesh.volatile_fast_shared_comm_map(f.entity_rank())[proc]; for(size_t i=0; i fields const BucketIndices& bktIndices = mesh.volatile_fast_shared_comm_map(f.entity_rank())[iproc]; const std::vector& ords = bktIndices.ords; + f.sync_to_host(); + f.modify_on_host(); + size_t ords_offset = 0; for(size_t i=0; isize(); } + const BucketVector::const_iterator ie = input_buckets.end(); + BucketVector::const_iterator ik = input_buckets.begin(); + + for ( ; ik != ie ; ++ik ) { + const Bucket & k = ** ik ; + if ( selector( k ) ) { count += k.size(); } } return count ; @@ -96,14 +100,17 @@ void get_selected_entities( const Selector & selector , { size_t count = count_selected_entities(selector,input_buckets); - entities.clear(); - entities.reserve(count); + entities.resize(count); + + const BucketVector::const_iterator ie = input_buckets.end(); + BucketVector::const_iterator ik = input_buckets.begin(); - for (const Bucket* bptr : input_buckets) { - const Bucket& bkt = *bptr; - if ( selector( bkt ) ) { - for(Entity e : bkt) { - entities.push_back(e); + for ( size_t j = 0 ; ik != ie ; ++ik ) { + const Bucket & k = ** ik ; + if ( selector( k ) ) { + const size_t n = k.size(); + for ( size_t i = 0; i < n; ++i, ++j ) { + entities[j] = k[i] ; } } } diff --git a/packages/stk/stk_mesh/stk_mesh/base/HostMeshManager.cpp b/packages/stk/stk_mesh/stk_mesh/base/HostMeshManager.cpp index 74975c6ab0bf..c27b5a95f4dd 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/HostMeshManager.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/HostMeshManager.cpp @@ -1,6 +1,7 @@ #include "HostMeshManager.hpp" #include "BulkData.hpp" -#include "stk_mesh/base/NgpMesh.hpp" +#include "NgpMesh.hpp" +#include "stk_util/stk_kokkos_macros.h" namespace stk { namespace mesh { @@ -18,7 +19,7 @@ HostMeshManager::~HostMeshManager() stk::mesh::NgpMesh & HostMeshManager::get_mesh() { -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH ThrowErrorMsg("HostMeshManager does not have device mesh"); static DeviceMesh deviceMesh; return deviceMesh; // Never used; Avoids build problems on the GPU diff --git a/packages/stk/stk_mesh/stk_mesh/base/Ngp.hpp b/packages/stk/stk_mesh/stk_mesh/base/Ngp.hpp index 34a2b5f8665a..3fa5ac2bd901 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Ngp.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Ngp.hpp @@ -35,6 +35,7 @@ #define STK_MESH_NGP_HPP #include "Kokkos_Core.hpp" +#include "stk_util/stk_kokkos_macros.h" namespace stk { namespace mesh { @@ -46,7 +47,7 @@ template class DeviceField; template class ConstHostField; template class ConstDeviceField; -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH using NgpMesh = stk::mesh::DeviceMesh; template using NgpField = stk::mesh::DeviceField; template using NgpConstField = stk::mesh::ConstDeviceField; diff --git a/packages/stk/stk_mesh/stk_mesh/base/NgpAtomics.hpp b/packages/stk/stk_mesh/stk_mesh/base/NgpAtomics.hpp index acb8ac541af1..bf700e149daa 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/NgpAtomics.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/NgpAtomics.hpp @@ -34,7 +34,7 @@ #ifndef STK_MESH_NGPATOMICS_HPP #define STK_MESH_NGPATOMICS_HPP -#include "stk_util/stk_config.h" +#include "stk_util/stk_kokkos_macros.h" #include "Kokkos_Core.hpp" namespace stk { diff --git a/packages/stk/stk_mesh/stk_mesh/base/NgpMultistateField.hpp b/packages/stk/stk_mesh/stk_mesh/base/NgpMultistateField.hpp index 3d995de1c9b2..a92f5d0a50ed 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/NgpMultistateField.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/NgpMultistateField.hpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include @@ -120,6 +121,30 @@ class NgpMultistateField fieldNewState.modify_on_device(); } + STK_FUNCTION + T& get_new(const stk::mesh::NgpMesh& ngpMesh, stk::mesh::Entity entity, int component) const + { + return get_field_new_state().get(ngpMesh, entity, component); + } + + STK_FUNCTION + T& get_new(const stk::mesh::FastMeshIndex& entity, int component) const + { + return get_field_new_state().get(entity, component); + } + + STK_FUNCTION + T get_old(stk::mesh::FieldState state, const stk::mesh::NgpMesh& ngpMesh, stk::mesh::Entity entity, int component) const + { + return get_field_old_state(state).get(ngpMesh, entity, component); + } + + STK_FUNCTION + T get_old(stk::mesh::FieldState state, const stk::mesh::FastMeshIndex& entity, int component) const + { + return get_field_old_state(state).get(entity, component); + } + private: void clear_sync_state() { @@ -158,39 +183,6 @@ class NgpMultistateField stk::mesh::NgpConstField fieldOldStates[MAX_NUM_FIELD_STATES]; }; -template -class ConvenientNgpMultistateField : public NgpMultistateField -{ -public: - ConvenientNgpMultistateField(stk::mesh::BulkData& bulk, const stk::mesh::FieldBase &stkField) : - NgpMultistateField(bulk, stkField) - {} - - STK_FUNCTION - T& get_new(const stk::mesh::NgpMesh& ngpMesh, stk::mesh::Entity entity, int component) const - { - return NgpMultistateField::get_field_new_state().get(ngpMesh, entity, component); - } - - STK_FUNCTION - T get_old(stk::mesh::FieldState state, const stk::mesh::NgpMesh& ngpMesh, stk::mesh::Entity entity, int component) const - { - return NgpMultistateField::get_field_old_state(state).get(ngpMesh, entity, component); - } - - STK_FUNCTION - T& get_new(const stk::mesh::FastMeshIndex& entity, int component) const - { - return NgpMultistateField::get_field_new_state().get(entity, component); - } - - STK_FUNCTION - T get_old(stk::mesh::FieldState state, const stk::mesh::FastMeshIndex& entity, int component) const - { - return NgpMultistateField::get_field_old_state(state).get(entity, component); - } -}; - } } diff --git a/packages/stk/stk_mesh/stk_mesh/base/Types.hpp b/packages/stk/stk_mesh/stk_mesh/base/Types.hpp index cbfd48a325e5..30fa77dbd210 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Types.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Types.hpp @@ -157,6 +157,8 @@ typedef std::map, BucketVector> SelectorBucketMa typedef std::vector VolatileFastSharedCommMap; typedef std::map > EntityToDependentProcessorsMap; +typedef std::vector >::const_iterator const_entity_iterator; +typedef std::vector >::iterator entity_iterator; typedef unsigned Ordinal; static const Ordinal InvalidOrdinal = static_cast(-1); // std::numeric_limits::max(); diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.cpp index 99f7e0ef6dae..1e1ae1e862b9 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.cpp @@ -44,8 +44,38 @@ namespace stk { namespace mesh { namespace impl { +struct EntityKeyEntityLess { +inline bool operator()(const std::pair& key_ent_pair, const EntityKey& key) const +{ + return key_ent_pair.first.m_value < key.m_value; +} +inline bool operator()(const EntityKey& key, const std::pair& key_ent_pair) const +{ + return key.m_value < key_ent_pair.first.m_value; +} +}; + +struct match_EntityKey { + match_EntityKey(const EntityKey& key) + : m_key(key) + {} + bool operator()(const std::pair& item) const + { return item.first == m_key; } + + bool operator==(const std::pair& item) const + { return item.first == m_key; } + +private: + const EntityKey& m_key; +}; + EntityRepository::EntityRepository() - : m_entities(stk::topology::NUM_RANKS) + : m_entities(stk::topology::NUM_RANKS), + m_create_cache(stk::topology::NUM_RANKS), + m_update_cache(stk::topology::NUM_RANKS), + m_destroy_cache(stk::topology::NUM_RANKS), + m_maxCreateCacheSize(512), + m_maxUpdateCacheSize(4096) { } @@ -62,40 +92,237 @@ size_t capacity_in_bytes(const VecType& v) size_t EntityRepository::heap_memory_in_bytes() const { size_t bytes = 0; - for(auto entityMap : m_entities) { bytes += entityMap.size()*sizeof(EntityKeyEntity); } + for(auto vec : m_entities) { bytes += capacity_in_bytes(vec); } + for(auto vec : m_create_cache) { bytes += capacity_in_bytes(vec); } + for(auto vec : m_update_cache) { bytes += capacity_in_bytes(vec); } + for(auto vec : m_destroy_cache) { bytes += capacity_in_bytes(vec); } return bytes; } -std::pair +void EntityRepository::clear_all_cache() +{ + for(EntityRank rank=stk::topology::BEGIN_RANK; rank& destroy = m_destroy_cache[rank]; + std::sort(destroy.begin(), destroy.end()); + EntityKeyEntityVector& entities = m_entities[rank]; + size_t destroyIdx = 0; + EntityKeyEntityVector::iterator start = std::lower_bound(entities.begin(), entities.end(), destroy[0], EntityKeyEntityLess()); + EntityKeyEntityVector::iterator end = std::upper_bound(entities.begin(), entities.end(), destroy.back(), EntityKeyEntityLess()); + size_t startIdx = std::distance(entities.begin(), start); + size_t endIdx = std::distance(entities.begin(), end); + size_t keep = startIdx; + for(size_t i=startIdx; i keep) { + entities[keep] = entities[i]; + } + ++keep; + } + if (endIdx < entities.size()) { + size_t len = entities.size() - endIdx; + std::memmove(&entities[keep], &entities[endIdx], len*sizeof(EntityKeyEntity)); + keep += len; + } + entities.resize(keep); + + destroy.clear(); + size_t num = std::max(m_entities[stk::topology::NODE_RANK].size(), + m_entities[stk::topology::ELEM_RANK].size()); + unsigned possibleCacheSize = num/1000; + m_maxUpdateCacheSize = std::max(m_maxUpdateCacheSize, possibleCacheSize); + } +} + +void EntityRepository::clear_updated_entity_cache(EntityRank rank) const +{ + if (!m_update_cache[rank].empty()) { + std::vector >& update = m_update_cache[rank]; + std::sort(update.begin(), update.end()); + EntityKeyEntityVector::iterator iter = m_entities[rank].begin(); + EntityKeyEntityVector::iterator end = m_entities[rank].end(); + for(const std::pair& oldnew : update) { + EntityKeyEntityVector::iterator thisIter = std::lower_bound(iter, end, oldnew.first, EntityKeyEntityLess()); + if (thisIter != end && thisIter->first == oldnew.first) { + thisIter->first = oldnew.second; + iter = thisIter+1; + } + } + m_update_cache[rank].clear(); + std::sort(m_entities[rank].begin(), m_entities[rank].end()); + } +} + +void EntityRepository::clear_created_entity_cache(EntityRank rank) const +{ + if (!m_create_cache[rank].empty()) { + std::sort(m_create_cache[rank].begin(), m_create_cache[rank].end()); + unsigned numOld = m_entities[rank].size(); + m_entities[rank].insert(m_entities[rank].end(), m_create_cache[rank].begin(), m_create_cache[rank].end()); + if (numOld > 0) { + const EntityKey& firstNewKey = m_create_cache[rank][0].first; + const EntityKey& lastOldKey = m_entities[rank][numOld-1].first; + const bool isOverlap = (firstNewKey < lastOldKey); + if (isOverlap) { + EntityKeyEntityVector::iterator oldEnd = m_entities[rank].begin()+numOld; + EntityKeyEntityVector::iterator startOfOverlap = std::lower_bound(m_entities[rank].begin(), oldEnd, firstNewKey, EntityKeyEntityLess()); + EntityKeyEntityVector::iterator endOfOverlap = std::lower_bound(oldEnd, m_entities[rank].end(), lastOldKey, EntityKeyEntityLess()); + std::inplace_merge(startOfOverlap, oldEnd, endOfOverlap); + } + } + m_create_cache[rank].clear(); + size_t num = std::max(m_entities[stk::topology::NODE_RANK].size(), + m_entities[stk::topology::ELEM_RANK].size()); + unsigned possibleCacheSize = num/1000; + m_maxCreateCacheSize = std::max(m_maxCreateCacheSize, possibleCacheSize); + } +} + +void EntityRepository::clear_cache(EntityRank rank) const +{ + clear_created_entity_cache(rank); + + clear_updated_entity_cache(rank); + + clear_destroyed_entity_cache(rank); +} + +std::pair +EntityRepository::add_to_cache(const EntityKey& key) +{ + bool inserted_new_entity = false; + EntityRank rank = key.rank(); + EntityKeyEntityVector& cache = m_create_cache[rank]; + + if (cache.size() >= m_maxCreateCacheSize) { + clear_cache(rank); + } + + EntityKeyEntityVector& entities = m_entities[rank]; + EntityKeyEntityVector::iterator iter = std::lower_bound(entities.begin(), entities.end(), + key, EntityKeyEntityLess()); + Entity entity; + if (iter == entities.end() || iter->first != key) { + cache.emplace_back(key, Entity()); + iter = cache.begin()+(cache.size()-1); + inserted_new_entity = true; + } + else { + inserted_new_entity = false; + } + + if (cache.size() >= m_maxCreateCacheSize) { + clear_cache(rank); + iter = std::lower_bound(entities.begin(), entities.end(), key, EntityKeyEntityLess()); + } + + return std::make_pair(iter, inserted_new_entity); +} + +stk::mesh::entity_iterator EntityRepository::get_from_cache(const EntityKey& key) const +{ + if (!m_create_cache[key.rank()].empty()) { + EntityKeyEntityVector& cache = m_create_cache[key.rank()]; + EntityKeyEntityVector::iterator iter = + std::find_if(cache.begin(), cache.end(), match_EntityKey(key)); + if (iter != cache.end()) { + return iter; + } + } + return m_create_cache[key.rank()].end(); +} + +std::pair EntityRepository::internal_create_entity( const EntityKey & key) { - EntityKeyEntityMap& entMap = m_entities[key.rank()]; - return entMap.insert(EntityKeyEntity(key,Entity())); + if (key.rank() > m_entities.size()) { + m_entities.resize(key.rank()); + m_create_cache.resize(key.rank()); + m_update_cache.resize(key.rank()); + m_destroy_cache.resize(key.rank()); + } + + clear_updated_entity_cache(key.rank()); + clear_destroyed_entity_cache(key.rank()); + + entity_iterator ent = get_from_cache(key); + if (ent != m_create_cache[key.rank()].end()) { + return std::make_pair(ent, false); + } + + return add_to_cache(key); } Entity EntityRepository::get_entity(const EntityKey &key) const { - const EntityRank rank = key.rank(); + EntityRank rank = key.rank(); + if (!m_destroy_cache[rank].empty()) { + const std::vector& destroyed = m_destroy_cache[rank]; + if (destroyed.size() < 64) { + std::vector::const_iterator iter = std::find(destroyed.begin(), destroyed.end(), key); + if (iter != destroyed.end()) { + return Entity(); + } + } + else { + clear_destroyed_entity_cache(rank); + } + } + + clear_updated_entity_cache(rank); + + ThrowErrorMsgIf( ! key.is_valid(), + "Invalid key: " << key.rank() << " " << key.id()); + + if (rank >= m_entities.size()) { + return Entity(); + } - entity_iterator ent = m_entities[rank].find(key); - if (ent != m_entities[rank].end()) { + entity_iterator ent = get_from_cache(key); + if (ent != m_create_cache[rank].end()) { return ent->second; } - return Entity(); + const EntityKeyEntityVector& entities = m_entities[rank]; + + const EntityKeyEntityVector::const_iterator iter = std::lower_bound(entities.begin(), entities.end(), key, EntityKeyEntityLess()); + + return (iter != entities.end() && (iter->first==key)) ? iter->second : Entity() ; } void EntityRepository::update_entity_key(EntityKey new_key, EntityKey old_key, Entity entity) { - const EntityRank rank = new_key.rank(); + EntityRank rank = new_key.rank(); + clear_created_entity_cache(rank); + clear_destroyed_entity_cache(rank); + + if (m_update_cache[rank].size() >= m_maxUpdateCacheSize) { + clear_cache(rank); + } - m_entities[rank].erase(old_key); - m_entities[rank].insert(EntityKeyEntity(new_key, entity)); + m_update_cache[rank].emplace_back(old_key, new_key); } -void EntityRepository::destroy_entity(EntityKey key) +void EntityRepository::destroy_entity(EntityKey key, Entity entity) { - m_entities[key.rank()].erase(key); + EntityRank rank = key.rank(); + clear_created_entity_cache(rank); + clear_updated_entity_cache(rank); + + if (m_destroy_cache[rank].size() >= m_maxUpdateCacheSize) { + clear_cache(rank); + } + + m_destroy_cache[rank].push_back(key); } } // namespace impl diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.hpp index 87350ac5c2bb..db8e6c3deb6e 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityRepository.hpp @@ -36,7 +36,7 @@ #define stk_mesh_baseImpl_EntityRepository_hpp #include // for size_t -#include +#include // for map, map<>::value_compare #include // for Entity #include // for pair #include "stk_mesh/base/EntityKey.hpp" // for EntityKey @@ -51,11 +51,9 @@ class EntityRepository { public: typedef std::pair EntityKeyEntity; - typedef std::unordered_map EntityKeyEntityMap; - typedef EntityKeyEntityMap::const_iterator const_iterator; - typedef EntityKeyEntityMap::iterator iterator; - typedef EntityKeyEntityMap::const_iterator const_entity_iterator; - typedef EntityKeyEntityMap::iterator entity_iterator; + typedef std::vector EntityKeyEntityVector; + typedef EntityKeyEntityVector::const_iterator const_iterator; + typedef EntityKeyEntityVector::iterator iterator; EntityRepository(); @@ -65,10 +63,12 @@ class EntityRepository { const_entity_iterator begin_rank(EntityRank ent_rank) const { + clear_cache(ent_rank); return m_entities[ent_rank].begin(); } const_entity_iterator end_rank(EntityRank ent_rank) const { + clear_cache(ent_rank); return m_entities[ent_rank].end(); } @@ -82,17 +82,34 @@ class EntityRepository { void update_entity_key(EntityKey new_key, EntityKey old_key, Entity entity); - void destroy_entity(EntityKey key); + void destroy_entity(EntityKey key, Entity entity); void update_num_ranks(unsigned numRanks) { m_entities.resize(numRanks); + m_create_cache.resize(numRanks); + m_update_cache.resize(numRanks); + m_destroy_cache.resize(numRanks); } + void clear_all_cache(); + void clear_cache(EntityRank rank) const; + size_t heap_memory_in_bytes() const; private: - - mutable std::vector m_entities; + void clear_destroyed_entity_cache(EntityRank rank) const; + void clear_updated_entity_cache(EntityRank rank) const; + void clear_created_entity_cache(EntityRank rank) const; + entity_iterator get_from_cache(const EntityKey& key) const; + std::pair add_to_cache(const EntityKey& key); + + mutable std::vector m_entities; + + mutable std::vector m_create_cache; + mutable std::vector > > m_update_cache; + mutable std::vector > m_destroy_cache; + mutable unsigned m_maxCreateCacheSize; + mutable unsigned m_maxUpdateCacheSize; //disable copy constructor and assignment operator EntityRepository(const EntityRepository &); diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModification.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModification.cpp index c10108220e7c..1dc902cde250 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModification.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModification.cpp @@ -73,6 +73,8 @@ bool MeshModification::internal_modification_end(modification_optimization opt) ThrowAssertMsg(m_bulkData.add_fmwk_data() || impl::check_no_shared_elements_or_higher(m_bulkData)==0, "BulkData::modification_end ERROR, Sharing of entities with rank ELEMENT_RANK or higher is not allowed."); + m_bulkData.m_entity_repo->clear_all_cache(); + if(m_bulkData.parallel_size() > 1) { // Resolve modification or deletion of shared entities diff --git a/packages/stk/stk_performance_tests/stk_mesh/NgpFieldAccess.cpp b/packages/stk/stk_performance_tests/stk_mesh/NgpFieldAccess.cpp index 0f42021661a5..e178b7db9d72 100644 --- a/packages/stk/stk_performance_tests/stk_mesh/NgpFieldAccess.cpp +++ b/packages/stk/stk_performance_tests/stk_mesh/NgpFieldAccess.cpp @@ -102,7 +102,7 @@ TEST_F(NgpFieldAccess, Centroid) { if (get_parallel_size() != 1) return; - const int NUM_RUNS = 2000; + const int NUM_RUNS = 400; const int ELEMS_PER_DIM = 100; declare_centroid_field(); @@ -121,7 +121,7 @@ TEST_F(NgpFieldAccess, ConstCentroid) { if (get_parallel_size() != 1) return; - const int NUM_RUNS = 2000; + const int NUM_RUNS = 400; const int ELEMS_PER_DIM = 100; declare_centroid_field(); @@ -140,7 +140,7 @@ TEST_F(NgpFieldAccess, CentroidMultiBlock) { if (get_parallel_size() != 1) return; - const int NUM_RUNS = 2000; + const int NUM_RUNS = 5; const int ELEMS_PER_DIM = 100; const int NUM_BLOCKS = 100; @@ -164,7 +164,7 @@ TEST_F(NgpFieldAccess, ConstCentroidMultiBlock) { if (get_parallel_size() != 1) return; - const int NUM_RUNS = 2000; + const int NUM_RUNS = 5; const int ELEMS_PER_DIM = 100; const int NUM_BLOCKS = 100; diff --git a/packages/stk/stk_performance_tests/stk_mesh/calculate_centroid.hpp b/packages/stk/stk_performance_tests/stk_mesh/calculate_centroid.hpp index a40b90858c22..0163e4ce9c51 100644 --- a/packages/stk/stk_performance_tests/stk_mesh/calculate_centroid.hpp +++ b/packages/stk/stk_performance_tests/stk_mesh/calculate_centroid.hpp @@ -130,7 +130,7 @@ void calculate_centroid_using_coord_field(const stk::mesh::BulkData &bulk, const const stk::mesh::FieldBase& coords = *bulk.mesh_meta_data().coordinate_field(); CoordFieldType ngpCoords(bulk, coords); stk::mesh::NgpField ngpCentroid(bulk, centroid); - stk::mesh::NgpMesh ngpMesh(bulk); + stk::mesh::NgpMesh& ngpMesh = bulk.get_updated_ngp_mesh(); calculate_centroid(ngpMesh, ngpCoords, selector, ngpCentroid); diff --git a/packages/stk/stk_simd/stk_simd/disimd/DISimd.hpp b/packages/stk/stk_simd/stk_simd/disimd/DISimd.hpp index e50ec3a196af..a65631da5c96 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/DISimd.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/DISimd.hpp @@ -37,7 +37,7 @@ #ifndef SIMD_ALWAYS_INLINE //currently necessary to avoid the 'always_inline' defined in simd.hpp -#define SIMD_ALWAYS_INLINE +#define SIMD_ALWAYS_INLINE __attribute__((always_inline)) #endif #include "./simd.hpp" diff --git a/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleLoadStore.hpp b/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleLoadStore.hpp index f3545890bbc4..6b52567190bf 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleLoadStore.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleLoadStore.hpp @@ -37,21 +37,52 @@ namespace stk { namespace simd { +using native_simd = SIMD_NAMESPACE::simd; + +namespace impl { + +template +STK_MATH_FORCE_INLINE +void store_strided(double* x, const simd::Double& z, const int offset) +{ + for(int i=0; i +STK_MATH_FORCE_INLINE +void store_strided<2>(double* x, const simd::Double& z, const int offset) +{ + x[0] = z[0]; x[offset] = z[1]; +} + +template<> +STK_MATH_FORCE_INLINE +void store_strided<4>(double* x, const simd::Double& z, const int offset) +{ + x[0] = z[0]; x[offset] = z[1]; x[2*offset] = z[2]; x[3*offset] = z[3]; +} + +template<> +STK_MATH_FORCE_INLINE +void store_strided<8>(double* x, const simd::Double& z, const int offset) +{ + x[0] = z[0]; x[offset] = z[1]; x[2*offset] = z[2]; x[3*offset] = z[3]; + x[4*offset] = z[4]; x[5*offset] = z[5]; x[6*offset] = z[6]; x[7*offset] = z[7]; +} + +} // namespace impl + STK_MATH_FORCE_INLINE simd::Double load_aligned(const double* x) { - return simd::Double(SIMD_NAMESPACE::simd(x, SIMD_NAMESPACE::element_aligned_tag())); + return simd::Double(native_simd(x, SIMD_NAMESPACE::element_aligned_tag())); } STK_MATH_FORCE_INLINE simd::Double load(const double* x) { //FIXME: supposed to be un-aligned load... - return simd::Double(SIMD_NAMESPACE::simd(x, SIMD_NAMESPACE::element_aligned_tag())); + return simd::Double(native_simd(x, SIMD_NAMESPACE::element_aligned_tag())); } STK_MATH_FORCE_INLINE simd::Double load(const double* x, const int offset) { - simd::Double result; - for(int i=0; i(x, z, offset); } } // namespace simd diff --git a/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleMath.hpp b/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleMath.hpp index 348d022c0b4f..5439d3272eee 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleMath.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/DISimdDoubleMath.hpp @@ -297,27 +297,15 @@ STK_MATH_FORCE_INLINE simd::Double erf(const simd::Double& a) { } STK_MATH_FORCE_INLINE simd::Double multiplysign(const simd::Double& x, const simd::Double& y) { // return x times sign of y - simd::Double tmp; - for (int i=0; i < simd::ndoubles; ++i) { - tmp[i] = x[i]*std::copysign(1.0, y[i]); - } - return tmp; + return simd::Double(SIMD_NAMESPACE::multiplysign(x._data, y._data)); } STK_MATH_FORCE_INLINE simd::Double copysign(const simd::Double& x, const simd::Double& y) { // return abs(x) times sign of y - simd::Double tmp; - for (int i=0; i < simd::ndoubles; ++i) { - tmp[i] = std::copysign(x[i], y[i]); - } - return tmp; + return simd::Double(SIMD_NAMESPACE::copysign(x._data, y._data)); } STK_MATH_FORCE_INLINE simd::Double abs(const simd::Double& x) { - simd::Double tmp; - for (int i=0; i < simd::ndoubles; ++i) { - tmp[i] = std::abs(x[i]); - } - return tmp; + return simd::Double(SIMD_NAMESPACE::abs(x._data)); } STK_MATH_FORCE_INLINE simd::Double min(const simd::Double& x, const simd::Double& y) { diff --git a/packages/stk/stk_simd/stk_simd/disimd/avx.hpp b/packages/stk/stk_simd/stk_simd/disimd/avx.hpp index 0ec77d0402ff..d1b4ba6edd03 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/avx.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/avx.hpp @@ -105,6 +105,11 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float value) :m_value(_mm256_set1_ps(value)) {} + SIMD_ALWAYS_INLINE inline simd( + float a, float b, float c, float d, + float e, float f, float g, float h) + :m_value(_mm256_setr_ps(a, b, c, d, e, f, g, h)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -115,9 +120,13 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags /*flags*/) + :m_value(_mm256_loadu_ps(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(float const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride], + ptr[4*stride], ptr[5*stride], ptr[6*stride], ptr[7*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m256 const& value_in) :m_value(value_in) {} @@ -151,6 +160,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + __m256 const sign_mask = _mm256_set1_ps(-0.f); + return simd(_mm256_xor_ps(a.get(), _mm256_and_ps(sign_mask, b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + __m256 const sign_mask = _mm256_set1_ps(-0.); + return simd(_mm256_xor_ps(_mm256_andnot_ps(sign_mask, a.get()), _mm256_and_ps(sign_mask, b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m256 sign_mask = _mm256_set1_ps(-0.f); // -0.f = 1 << 31 return simd(_mm256_andnot_ps(sign_mask, a.get())); @@ -244,10 +263,18 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return 4; } SIMD_ALWAYS_INLINE inline simd(double value) :m_value(_mm256_set1_pd(value)) {} + SIMD_ALWAYS_INLINE inline simd( + double a, double b, double c, double d) + :m_value(_mm256_setr_pd(a, b, c, d)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -264,9 +291,12 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags /*flags*/) + :m_value(_mm256_loadu_pd(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(double const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m256d const& value_in) :m_value(value_in) {} @@ -305,6 +335,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + __m256d const sign_mask = _mm256_set1_pd(-0.f); + return simd(_mm256_xor_pd(a.get(), _mm256_and_pd(sign_mask, b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + __m256d const sign_mask = _mm256_set1_pd(-0.f); + return simd(_mm256_xor_pd(_mm256_andnot_pd(sign_mask, a.get()), _mm256_and_pd(sign_mask, b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m256d const sign_mask = _mm256_set1_pd(-0.f); // -0.f = 1 << 31 return simd(_mm256_andnot_pd(sign_mask, a.get())); diff --git a/packages/stk/stk_simd/stk_simd/disimd/avx512.hpp b/packages/stk/stk_simd/stk_simd/disimd/avx512.hpp index 21e4fdf814ce..f592b9821be2 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/avx512.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/avx512.hpp @@ -107,6 +107,15 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float value) :m_value(_mm512_set1_ps(value)) {} + SIMD_ALWAYS_INLINE inline simd( + float a, float b, float c, float d, + float e, float f, float g, float h, + float i, float j, float k, float l, + float m, float n, float o, float p) + :m_value(_mm512_setr_ps( + a, b, c, d, e, f, g, h, + i, j, k, l, m, n, o, p)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -117,9 +126,15 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags /*flags*/) + :m_value(_mm512_loadu_ps(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(float const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride], + ptr[4*stride], ptr[5*stride], ptr[6*stride], ptr[7*stride], + ptr[8*stride], ptr[9*stride], ptr[10*stride], ptr[11*stride], + ptr[12*stride], ptr[13*stride], ptr[14*stride], ptr[15*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m512 const& value_in) :m_value(value_in) {} @@ -153,6 +168,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + static simd sign_mask(-0.f); + return simd(_mm512_xor_ps(a.get(), _mm512_and_ps(sign_mask.get(), b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + static simd sign_mask(-0.f); + return simd(_mm512_xor_ps(_mm512_andnot_ps(sign_mask.get(), a.get()) , _mm512_and_ps(sign_mask.get(), b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m512 const rhs = a.get(); return reinterpret_cast<__m512>(_mm512_and_epi32(reinterpret_cast<__m512i>(rhs), _mm512_set1_epi32(0x7fffffff))); @@ -242,10 +267,19 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return 8; } SIMD_ALWAYS_INLINE inline simd(double value) :m_value(_mm512_set1_pd(value)) {} + SIMD_ALWAYS_INLINE inline simd( + double a, double b, double c, double d, + double e, double f, double g, double h) + :m_value(_mm512_setr_pd(a, b, c, d, e, f, g, h)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -262,9 +296,13 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags /*flags*/) + :m_value(_mm512_loadu_pd(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(double const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride], + ptr[4*stride], ptr[5*stride], ptr[6*stride], ptr[7*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m512d const& value_in) :m_value(value_in) {} @@ -303,6 +341,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + static simd sign_mask(-0.0); + return simd(_mm512_xor_pd(a.get(), _mm512_and_pd(sign_mask.get(), b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + static simd sign_mask(-0.0); + return simd(_mm512_xor_pd(_mm512_andnot_pd(sign_mask.get(), a.get()) , _mm512_and_pd(sign_mask.get(), b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m512d const rhs = a.get(); return reinterpret_cast<__m512d>(_mm512_and_epi64(_mm512_set1_epi64(0x7FFFFFFFFFFFFFFF), diff --git a/packages/stk/stk_simd/stk_simd/disimd/neon.hpp b/packages/stk/stk_simd/stk_simd/disimd/neon.hpp index 821e39ca4545..ed23ae071140 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/neon.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/neon.hpp @@ -105,6 +105,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float value) :m_value(vdupq_n_f32(value)) {} + SIMD_ALWAYS_INLINE inline simd(float a, float b, float c, float d) + :m_value((float32x4_t){a, b, c, d}) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -118,6 +121,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags flags) { copy_from(ptr, flags); } + SIMD_ALWAYS_INLINE inline simd(float const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(float32x4_t const& value_in) :m_value(value_in) {} @@ -230,10 +236,17 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return 2; } SIMD_ALWAYS_INLINE inline simd(double value) :m_value(vdupq_n_f64(value)) {} + SIMD_ALWAYS_INLINE inline simd(double a, double b) + :m_value((float64x2_t){a, b}) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -253,6 +266,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags flags) { copy_from(ptr, flags); } + SIMD_ALWAYS_INLINE inline simd(double const* ptr, int stride) + :simd(ptr[0], ptr[stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(float64x2_t const& value_in) :m_value(value_in) {} diff --git a/packages/stk/stk_simd/stk_simd/disimd/scalar.hpp b/packages/stk/stk_simd/stk_simd/disimd/scalar.hpp index 794fc05ce2be..9d6daf2ede00 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/scalar.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/scalar.hpp @@ -94,6 +94,10 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE static constexpr int size() { return 1; } SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd(T value) :m_value(value) @@ -117,6 +121,9 @@ class simd { SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd(T const* ptr, Flags flags) { copy_from(ptr, flags); } + SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd(T const* ptr, int /*stride*/) + : m_value(ptr[0]) + {} SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd operator*(simd const& other) const { return simd(m_value * other.m_value); } diff --git a/packages/stk/stk_simd/stk_simd/disimd/simd.hpp b/packages/stk/stk_simd/stk_simd/disimd/simd.hpp index 6fa8a0525c90..d4985cdb900d 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/simd.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/simd.hpp @@ -53,7 +53,12 @@ #include "vector_size.hpp" #endif -#if defined(__SSE__) +#ifdef __CUDACC__ +#include "cuda_warp.hpp" + +#else + +#ifdef __SSE__ #include "sse.hpp" #endif @@ -69,12 +74,10 @@ #include "neon.hpp" #endif -#if defined(__VSX__) && (!defined(__CUDACC__)) +#ifdef __VSX__ #include "vsx.hpp" #endif -#ifdef __CUDACC__ -#include "cuda_warp.hpp" #endif namespace SIMD_NAMESPACE { diff --git a/packages/stk/stk_simd/stk_simd/disimd/simd_common.hpp b/packages/stk/stk_simd/stk_simd/disimd/simd_common.hpp index ee4e7e7cfb21..e3cd8aa7bae7 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/simd_common.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/simd_common.hpp @@ -164,6 +164,37 @@ SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd operator/(simd c return a / simd(b); } +template +SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd multiplysign(simd a, simd b) { + T tmp_a[simd::size()]; + T tmp_b[simd::size()]; + a.copy_to(tmp_a, element_aligned_tag()); + b.copy_to(tmp_b, element_aligned_tag()); + for (int i = 0; i < simd::size(); ++i) tmp_a[i] = tmp_a[i]*std::copysign(1.0, tmp_b[i]); + a.copy_from(tmp_a, element_aligned_tag()); + return a; +} + +template +SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd copysign(simd a, simd b) { + T tmp_a[simd::size()]; + T tmp_b[simd::size()]; + a.copy_to(tmp_a, element_aligned_tag()); + b.copy_to(tmp_b, element_aligned_tag()); + for (int i = 0; i < simd::size(); ++i) tmp_a[i] = std::copysign(tmp_a[i], tmp_b[i]); + a.copy_from(tmp_a, element_aligned_tag()); + return a; +} + +template +SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd abs(simd a) { + T tmp[simd::size()]; + a.copy_to(tmp, element_aligned_tag()); + for (int i = 0; i < simd::size(); ++i) tmp[i] = std::abs(tmp[i]); + a.copy_from(tmp, element_aligned_tag()); + return a; +} + template SIMD_ALWAYS_INLINE SIMD_HOST_DEVICE inline simd cbrt(simd a) { T tmp[simd::size()]; @@ -207,7 +238,7 @@ class simd_storage { SIMD_ALWAYS_INLINE inline simd_storage() = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return simd::size(); } - SIMD_ALWAYS_INLINE inline + SIMD_ALWAYS_INLINE explicit inline simd_storage(simd const& value) { value.copy_to(m_value, element_aligned_tag()); } diff --git a/packages/stk/stk_simd/stk_simd/disimd/sse.hpp b/packages/stk/stk_simd/stk_simd/disimd/sse.hpp index d8f6ebab3894..ed7420786fb7 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/sse.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/sse.hpp @@ -125,6 +125,10 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float value) :m_value(_mm_set1_ps(value)) {} + SIMD_ALWAYS_INLINE inline simd( + float a, float b, float c, float d) + :m_value(_mm_setr_ps(a, b, c, d)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -135,9 +139,12 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags /*flags*/) + :m_value(_mm_loadu_ps(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(float const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m128 const& value_in) :m_value(value_in) {} @@ -171,6 +178,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + __m128 const sign_mask = _mm_set1_ps(-0.); + return simd(_mm_xor_ps(a.get(), _mm_and_ps(sign_mask, b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + __m128 const sign_mask = _mm_set1_ps(-0.); + return simd(_mm_xor_ps(_mm_andnot_ps(sign_mask, a.get()), _mm_and_ps(sign_mask, b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m128 const sign_mask = _mm_set1_ps(-0.f); // -0.f = 1 << 31 return simd(_mm_andnot_ps(sign_mask, a.get())); @@ -266,10 +283,17 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return 2; } SIMD_ALWAYS_INLINE inline simd(double value) :m_value(_mm_set1_pd(value)) {} + SIMD_ALWAYS_INLINE inline simd(double a, double b) + :m_value(_mm_setr_pd(a, b)) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -286,9 +310,12 @@ class simd { return *this; } template - SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags flags) { - copy_from(ptr, flags); - } + SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags /*flags*/) + :m_value(_mm_loadu_pd(ptr)) + {} + SIMD_ALWAYS_INLINE inline simd(double const* ptr, int stride) + :simd(ptr[0], ptr[stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__m128d const& value_in) :m_value(value_in) {} @@ -327,6 +354,16 @@ class simd { } }; +SIMD_ALWAYS_INLINE inline simd multiplysign(simd const& a, simd const& b) { + __m128d const sign_mask = _mm_set1_pd(-0.); + return simd(_mm_xor_pd(a.get(), _mm_and_pd(sign_mask, b.get()))); +} + +SIMD_ALWAYS_INLINE inline simd copysign(simd const& a, simd const& b) { + __m128d const sign_mask = _mm_set1_pd(-0.); + return simd(_mm_xor_pd(_mm_andnot_pd(sign_mask, a.get()), _mm_and_pd(sign_mask, b.get()))); +} + SIMD_ALWAYS_INLINE inline simd abs(simd const& a) { __m128d const sign_mask = _mm_set1_pd(-0.); // -0. = 1 << 63 return simd(_mm_andnot_pd(sign_mask, a.get())); diff --git a/packages/stk/stk_simd/stk_simd/disimd/vsx.hpp b/packages/stk/stk_simd/stk_simd/disimd/vsx.hpp index 0782f872f88a..1c4c98913296 100644 --- a/packages/stk/stk_simd/stk_simd/disimd/vsx.hpp +++ b/packages/stk/stk_simd/stk_simd/disimd/vsx.hpp @@ -114,6 +114,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float value) :m_value(vec_splats(value)) {} + SIMD_ALWAYS_INLINE inline simd(float a, float b, float c, float d) + :m_value((__vector float){a, b, c, d}) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -127,6 +130,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(float const* ptr, Flags flags) { copy_from(ptr, flags); } + SIMD_ALWAYS_INLINE inline simd(float const* ptr, int stride) + :simd(ptr[0], ptr[stride], ptr[2*stride], ptr[3*stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__vector float const& value_in) :m_value(value_in) {} @@ -247,10 +253,17 @@ class simd { using mask_type = simd_mask; using storage_type = simd_storage; SIMD_ALWAYS_INLINE inline simd() = default; + SIMD_ALWAYS_INLINE inline simd(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd(simd&&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd const&) = default; + SIMD_ALWAYS_INLINE inline simd& operator=(simd&&) = default; SIMD_ALWAYS_INLINE inline static constexpr int size() { return 2; } SIMD_ALWAYS_INLINE inline simd(double value) :m_value(vec_splats(value)) {} + SIMD_ALWAYS_INLINE inline simd(double a, double b) + :m_value((__vector double){a, b}) + {} SIMD_ALWAYS_INLINE inline simd(storage_type const& value) { copy_from(value.data(), element_aligned_tag()); @@ -270,6 +283,9 @@ class simd { SIMD_ALWAYS_INLINE inline simd(double const* ptr, Flags flags) { copy_from(ptr, flags); } + SIMD_ALWAYS_INLINE inline simd(double const* ptr, int stride) + :simd(ptr[0], ptr[stride]) + {} SIMD_ALWAYS_INLINE inline constexpr simd(__vector double const& value_in) :m_value(value_in) {} diff --git a/packages/stk/stk_topology/stk_topology/topology_detail/meta_functions.hpp b/packages/stk/stk_topology/stk_topology/topology_detail/meta_functions.hpp index c62d126e131b..9dd1830efcd2 100644 --- a/packages/stk/stk_topology/stk_topology/topology_detail/meta_functions.hpp +++ b/packages/stk/stk_topology/stk_topology/topology_detail/meta_functions.hpp @@ -85,13 +85,22 @@ constexpr bool defined_on_spatial_dimension_() //------------------------------------------------------------------------------ + template STK_INLINE_FUNCTION -constexpr topology::topology_t face_topology_() +constexpr typename std::enable_if::type +face_topology_() { - return (FaceOrdinal < Topology::num_faces) ? Topology::face_topology_vector[FaceOrdinal] : topology::INVALID_TOPOLOGY; + return Topology::face_topology_vector[FaceOrdinal]; } +template +STK_INLINE_FUNCTION +constexpr typename std::enable_if= Topology::num_faces, topology::topology_t>::type +face_topology_() +{ + return topology::INVALID_TOPOLOGY; +} //------------------------------------------------------------------------------ template diff --git a/packages/stk/stk_transfer/stk_transfer/copy_by_id/TransferCopyByIdStkMeshAdapter.hpp b/packages/stk/stk_transfer/stk_transfer/copy_by_id/TransferCopyByIdStkMeshAdapter.hpp index ca67f6c7597d..5f8f92ed1d74 100644 --- a/packages/stk/stk_transfer/stk_transfer/copy_by_id/TransferCopyByIdStkMeshAdapter.hpp +++ b/packages/stk/stk_transfer/stk_transfer/copy_by_id/TransferCopyByIdStkMeshAdapter.hpp @@ -207,7 +207,8 @@ public : stk::mesh::FieldBase* field=m_transfer_fields[fieldIndex]; - DataTypeKey::data_t dataType; + DataTypeKey::data_t dataType( DataTypeKey::data_t::INVALID_TYPE ); + if(field->type_is()) { dataType = DataTypeKey::data_t::UNSIGNED_INTEGER; diff --git a/packages/stk/stk_unit_tests/stk_balance/UnitTestTransientFields.cpp b/packages/stk/stk_unit_tests/stk_balance/UnitTestTransientFields.cpp index 069da80a0c38..595f3d673b2c 100644 --- a/packages/stk/stk_unit_tests/stk_balance/UnitTestTransientFields.cpp +++ b/packages/stk/stk_unit_tests/stk_balance/UnitTestTransientFields.cpp @@ -3,8 +3,6 @@ #include #include #include -#include -#include #include #include @@ -21,406 +19,495 @@ #include #include -namespace +class StkBalanceRunner { +public: + StkBalanceRunner(MPI_Comm c) + : m_comm(c) + { } -void unlink_serial_file(const std::string &baseName) -{ - unlink(baseName.c_str()); -} + void set_filename(const std::string& name) + { + m_options.m_inFile = name; + } -void unlink_parallel_file(const std::string &baseName) -{ - std::string file = stk::balance::internal::get_parallel_filename(stk::parallel_machine_rank(MPI_COMM_WORLD), - stk::parallel_machine_size(MPI_COMM_WORLD), - baseName); + void set_output_dir(const std::string& name) + { + m_options.outputDirectory = name; + } - unlink(file.c_str()); -} + void set_app_type_defaults(stk::balance::AppTypeDefaults defaults) + { + m_options.set_app_type_default(defaults); + } + void set_initial_decomp_method(const std::string& method) + { + m_options.set_initial_decomp_method(method); + } + + void set_decomp_method(const std::string& method) + { + m_options.set_decomp_method(method); + } + + void run_end_to_end() const + { + stk::balance::run_stk_rebalance(m_options, m_comm); + } + +private: + MPI_Comm m_comm; + stk::balance::ParsedOptions m_options; +}; -void verify_expected_time_steps_in_parallel_file(const std::string ¶llelMeshName, const std::vector& expectedTimeSteps) +class MeshFromFile { - stk::io::StkMeshIoBroker broker; - stk::mesh::MetaData metaData; - stk::mesh::BulkData bulkData(metaData, MPI_COMM_WORLD); - stk::io::fill_mesh_preexisting(broker, parallelMeshName, bulkData); - std::vector inputTimeSteps = broker.get_time_steps(); - EXPECT_EQ(expectedTimeSteps, inputTimeSteps); -} +public: + MeshFromFile(const MPI_Comm& c) + : m_comm(c), + m_empty(true), + bulk(meta, m_comm) + { } + + void fill_from_serial(const std::string& fileName) + { + stk::io::fill_mesh_with_auto_decomp(fileName, bulk, broker); + m_empty = false; + } + + void fill_from_parallel(const std::string& baseName) + { + stk::io::fill_mesh_preexisting(broker, baseName, bulk); + m_empty = false; + } + + bool is_empty() const { return m_empty; } + +private: + MPI_Comm m_comm; + bool m_empty; -void verify_expected_global_variable_in_parallel_file(const std::string ¶llelMeshName, const std::vector& expectedTimeSteps, - const std::string& globalVariableName) +public: + stk::mesh::MetaData meta; + stk::mesh::BulkData bulk; + stk::io::StkMeshIoBroker broker; +}; + +class TransientWriter { - stk::io::StkMeshIoBroker broker; - stk::mesh::MetaData metaData; - stk::mesh::BulkData bulkData(metaData, MPI_COMM_WORLD); - stk::io::fill_mesh_preexisting(broker, parallelMeshName, bulkData); - std::vector names; - broker.get_global_variable_names(names); - ASSERT_EQ(3u, names.size()); - EXPECT_EQ(globalVariableName+"_double", names[0]); - EXPECT_EQ(globalVariableName+"_int", names[1]); - EXPECT_EQ(globalVariableName+"_real_vec", names[2]); - - std::vector inputTimeSteps = broker.get_time_steps(); - EXPECT_EQ(expectedTimeSteps, inputTimeSteps); - double testTimeValue = 0; - const double epsilon = std::numeric_limits::epsilon(); - for(size_t i=0; i goldVecValue(3, expectedTimeSteps[i]); - std::vector testVecValue; - success = broker.get_global(globalVariableName+"_real_vec", testVecValue); - EXPECT_EQ(goldVecValue, testVecValue); +public: + TransientWriter(MPI_Comm c, const std::string& name) + : m_comm(c), + m_fileBaseName(name), + m_fieldName("myFieldName"), + m_varName("TestTime") + { } + + void set_field_name(const std::string& name) + { + m_fieldName = name; + } + + void set_global_variable_name(const std::string& name) + { + m_varName = name; + } + + void set_time_steps(const std::vector& steps) + { + m_timeSteps = steps; + } + + void write_static_mesh(const std::string& meshDesc) const + { + stk::unit_test_util::generated_mesh_to_file_in_serial(meshDesc, m_fileBaseName); + } + + void write_transient_mesh(const std::string& meshDesc) const + { + stk::unit_test_util::generated_mesh_with_transient_data_to_file_in_serial( + meshDesc, m_fileBaseName, m_fieldName, m_varName, m_timeSteps, m_fieldSetter); + } + + void write_two_element_mesh_with_sideset(stk::unit_test_util::ElementOrdering elemOrdering) const + { + if (stk::parallel_machine_rank(m_comm) == 0) { + stk::mesh::MetaData meta(3); + stk::mesh::BulkData bulk(meta, MPI_COMM_SELF); + + stk::unit_test_util::create_AB_mesh_with_sideset_and_field( + bulk, stk::unit_test_util::LEFT, elemOrdering, "dummyField"); + stk::io::write_mesh_with_fields(m_fileBaseName, bulk, 1, 1.0); } -} + } -void verify_input_transient_fields(stk::io::StkMeshIoBroker &stkIo, const std::vector& expectedTimeSteps) +private: + MPI_Comm m_comm; + std::string m_fileBaseName; + std::string m_fieldName; + std::string m_varName; + std::vector m_timeSteps; + stk::unit_test_util::IdAndTimeFieldValueSetter m_fieldSetter; +}; + +class TransientVerifier { - EXPECT_EQ(stkIo.get_num_time_steps(), (int)expectedTimeSteps.size()); - EXPECT_EQ(expectedTimeSteps, stkIo.get_time_steps()); +public: + TransientVerifier(const MPI_Comm& c) + : m_comm(c), + m_epsilon(std::numeric_limits::epsilon()) + { } - stk::io::FieldNameToPartVector fieldNamePartVector = stkIo.get_nodal_var_names(); + void verify_num_transient_fields(const MeshFromFile& mesh, unsigned expectedNumFields) const + { + stk::io::FieldNameToPartVector fieldNamePartVector = mesh.broker.get_nodal_var_names(); + EXPECT_EQ(fieldNamePartVector.size(), expectedNumFields); + } - EXPECT_EQ(fieldNamePartVector.size(), 2u); -} + void verify_time_steps(const MeshFromFile& mesh, const std::vector& expectedTimeSteps) const + { + EXPECT_EQ(expectedTimeSteps, mesh.broker.get_time_steps()); + } -void verify_input_static_fields(stk::io::StkMeshIoBroker &stkIo) -{ - EXPECT_EQ(0, stkIo.get_num_time_steps()); + void verify_global_variables_at_each_time_step(MeshFromFile& mesh, + const std::string& globalVariableName, + const std::vector& expectedTimeSteps) const + { + verify_time_steps(mesh, expectedTimeSteps); + verify_global_variable_names(mesh, globalVariableName); - stk::io::FieldNameToPartVector fieldNamePartVector = stkIo.get_nodal_var_names(); + int index = 0; + for(double timeStep : expectedTimeSteps) { + index++; + mesh.broker.read_defined_input_fields(index); - EXPECT_EQ(fieldNamePartVector.size(), 0u); -} + verify_global_int(mesh, globalVariableName, index); -void verify_transient_field_values(stk::mesh::BulkData& bulk, stk::mesh::FieldBase* field, double timeStep) -{ - const stk::mesh::BucketVector & entityBuckets = bulk.get_buckets(field->entity_rank(),bulk.mesh_meta_data().locally_owned_part()); - for (size_t bucketIndex = 0; bucketIndex < entityBuckets.size(); ++bucketIndex) { - stk::mesh::Bucket & entityBucket = * entityBuckets[bucketIndex]; - for (size_t entityIndex = 0; entityIndex < entityBucket.size(); ++entityIndex) { - stk::mesh::Entity entity = entityBucket[entityIndex]; - double * data = static_cast (stk::mesh::field_data(*field, entity)); - unsigned numEntriesPerEntity = stk::mesh::field_scalars_per_entity(*field, entity); - for(unsigned i=0; i(bulk.identifier(entity)), data[i]); - } + verify_global_double(mesh, globalVariableName, timeStep); + verify_global_real_vec(mesh, globalVariableName, timeStep); } -} + } -void verify_transient_data_from_parallel_file(const std::string &inputParallelFile, const std::vector& expectedTimeSteps) -{ - stk::mesh::MetaData meta; - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker broker; - stk::io::fill_mesh_preexisting(broker, inputParallelFile, bulk); + void verify_sideset_orientation(const MeshFromFile& mesh, + const stk::mesh::EntityId expectedId, + const stk::mesh::ConnectivityOrdinal expectedOrdinal) const + { + std::vector sidesets = mesh.bulk.get_sidesets(); + + if (stk::parallel_machine_rank(m_comm) == 0) { + ASSERT_EQ(sidesets.size(), 1u); + ASSERT_EQ(sidesets[0]->size(), 1u); + + const stk::mesh::SideSetEntry sideSetEntry = (*sidesets[0])[0]; + EXPECT_EQ(mesh.bulk.identifier(sideSetEntry.element), expectedId); + EXPECT_EQ(sideSetEntry.side, expectedOrdinal); + } + else { + ASSERT_EQ(sidesets.size(), 1u); + EXPECT_EQ(sidesets[0]->size(), 0u); + } + } + + void compare_entity_rank_names(const MeshFromFile& meshA, const MeshFromFile& meshB) const + { + EXPECT_EQ(meshA.meta.entity_rank_names(), meshB.meta.entity_rank_names()); + } - stk::mesh::FieldVector transientFields = stk::io::get_transient_fields(meta); + void verify_transient_field_names(const MeshFromFile& mesh, const std::string& fieldBaseName) const + { + const stk::mesh::FieldVector transientFields = stk::io::get_transient_fields(mesh.meta, stk::topology::NODE_RANK); + verify_transient_field_name(transientFields[0], fieldBaseName+"_scalar"); + verify_transient_field_name(transientFields[1], fieldBaseName+"_vector"); + } + + void verify_transient_fields(MeshFromFile& mesh) const + { + stk::io::StkMeshIoBroker& broker = mesh.broker; + + const stk::mesh::FieldVector transientFields = stk::io::get_transient_fields(mesh.meta); + const std::vector timeSteps = broker.get_time_steps(); - verify_input_transient_fields(broker, expectedTimeSteps); - std::vector timeSteps = broker.get_time_steps(); for(int iStep=0; iStep goldNames = {baseName+"_double", + baseName+"_int", + baseName+"_real_vec"}; + std::vector names; + mesh.broker.get_global_variable_names(names); + EXPECT_EQ(names, goldNames); + } - verify_input_static_fields(broker); -} + void verify_global_double(const MeshFromFile& mesh, const std::string& variable, double goldValue) const + { + double value; + EXPECT_TRUE(mesh.broker.get_global(variable+"_double", value)); + EXPECT_NEAR(value, goldValue, m_epsilon); + } -void verify_static_data_from_serial_file(const std::string &inputSerialFile) -{ - stk::mesh::MetaData meta; - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker broker; - broker.property_add(Ioss::Property("DECOMPOSITION_METHOD", "LINEAR")); - stk::io::fill_mesh_preexisting(broker, inputSerialFile, bulk); + void verify_global_int(const MeshFromFile& mesh, const std::string& variable, int goldValue) const + { + int value; + EXPECT_TRUE(mesh.broker.get_global(variable+"_int", value)); + EXPECT_EQ(value, goldValue); + } - verify_input_static_fields(broker); -} + void verify_global_real_vec(const MeshFromFile& mesh, const std::string& variable, double goldValue) const + { + std::vector values; + EXPECT_TRUE(mesh.broker.get_global(variable+"_real_vec", values)); -void verify_transient_field_name(stk::io::StkMeshIoBroker &brokerA, stk::io::StkMeshIoBroker &brokerB, const std::string &fieldBaseName) -{ - EXPECT_EQ(brokerA.meta_data().entity_rank_names(), brokerB.meta_data().entity_rank_names()); + EXPECT_EQ(values.size(), 3u); + for (double value : values) { + EXPECT_NEAR(value, goldValue, m_epsilon); + } + } - stk::mesh::FieldVector transientFieldsA = stk::io::get_transient_fields(brokerA.meta_data(), stk::topology::NODE_RANK); - stk::mesh::FieldVector transientFieldsB = stk::io::get_transient_fields(brokerB.meta_data(), stk::topology::NODE_RANK); + void verify_transient_field_name(stk::mesh::FieldBase* field, const std::string& fieldName) const + { + EXPECT_EQ(0, strcasecmp(fieldName.c_str(), field->name().c_str())) + << fieldName << " " << field->name(); + } - EXPECT_EQ(2u, transientFieldsA.size()); - EXPECT_EQ(2u, transientFieldsB.size()); + void verify_transient_field_values(const stk::mesh::BulkData& bulk, stk::mesh::FieldBase* field, double timeStep) const + { + const stk::mesh::BucketVector & entityBuckets = bulk.get_buckets(field->entity_rank(),bulk.mesh_meta_data().locally_owned_part()); - std::string scalarFieldName = fieldBaseName + "_scalar"; - std::string vectorFieldName = fieldBaseName + "_vector"; + for (size_t bucketIndex = 0; bucketIndex < entityBuckets.size(); ++bucketIndex) { + stk::mesh::Bucket & entityBucket = * entityBuckets[bucketIndex]; - EXPECT_EQ(0, strcasecmp(scalarFieldName.c_str(), transientFieldsA[0]->name().c_str())) << scalarFieldName << " " << transientFieldsA[0]->name(); - EXPECT_EQ(0, strcasecmp(vectorFieldName.c_str(), transientFieldsA[1]->name().c_str())) << vectorFieldName << " " << transientFieldsA[1]->name(); + for (size_t entityIndex = 0; entityIndex < entityBucket.size(); ++entityIndex) { + stk::mesh::Entity entity = entityBucket[entityIndex]; - EXPECT_EQ(0, strcasecmp(scalarFieldName.c_str(), transientFieldsB[0]->name().c_str())) << scalarFieldName << " " << transientFieldsB[0]->name(); - EXPECT_EQ(0, strcasecmp(vectorFieldName.c_str(), transientFieldsB[1]->name().c_str())) << vectorFieldName << " " << transientFieldsB[1]->name(); -} + double * data = static_cast (stk::mesh::field_data(*field, entity)); + unsigned numEntriesPerEntity = stk::mesh::field_scalars_per_entity(*field, entity); -void verify_transient_transfer(const std::string& serialOutputMeshName, - const std::string& parallelOutputMeshName, - const std::string& fieldBaseName, - const std::vector& expectedTimeSteps) -{ - stk::mesh::MetaData metaA; - stk::mesh::BulkData bulkA(metaA, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker brokerA; - brokerA.property_add(Ioss::Property("DECOMPOSITION_METHOD", "LINEAR")); - stk::io::fill_mesh_preexisting(brokerA, serialOutputMeshName, bulkA); - - stk::mesh::MetaData metaB; - stk::mesh::BulkData bulkB(metaB, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker brokerB; - stk::io::fill_mesh_preexisting(brokerB, parallelOutputMeshName, bulkB); - - verify_input_transient_fields(brokerA, expectedTimeSteps); - verify_transient_field_name(brokerA, brokerB, fieldBaseName); - verify_transient_data_from_parallel_file(parallelOutputMeshName, expectedTimeSteps); -} + for(unsigned i=0; i(bulk.identifier(entity)), data[i]); + } + } + } + } -void verify_static_transfer(const std::string& serialOutputMeshName, - const std::string& parallelOutputMeshName) -{ - verify_static_data_from_serial_file(serialOutputMeshName); - verify_static_data_from_parallel_file(parallelOutputMeshName); -} + MPI_Comm m_comm; + const double m_epsilon; +}; -TEST(TestTransientFieldBalance, verifyNumberOfSteps) +class TransientFieldBalance : public stk::unit_test_util::MeshFixture { - if(stk::parallel_machine_size(MPI_COMM_WORLD) == 2) - { - std::string outputMeshName = "twenty_hex_transient.e"; - std::string fieldBaseName = "myField"; - std::string globalVariableName = "TestTime"; - const std::vector expectedTimeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; +public: + TransientFieldBalance() + : fileBaseName("transient.e"), + balanceRunner(get_comm()), + writer(get_comm(), fileBaseName), + verifier(get_comm()), + m_initialMeshPtr(new MeshFromFile(get_comm())), + m_balancedMeshPtr(new MeshFromFile(get_comm())) + { + balanceRunner.set_filename(fileBaseName); + balanceRunner.set_output_dir("."); + balanceRunner.set_app_type_defaults(stk::balance::SD_DEFAULTS); + balanceRunner.set_decomp_method("rcb"); + } - stk::unit_test_util::IdAndTimeFieldValueSetter fieldSetter; - stk::unit_test_util::generated_mesh_with_transient_data_to_file_in_serial("1x1x20", outputMeshName, fieldBaseName, globalVariableName, expectedTimeSteps, fieldSetter); + MeshFromFile& get_initial_mesh() + { + if (m_initialMeshPtr->is_empty()) read_initial_mesh(); + return *m_initialMeshPtr; + } - stk::mesh::MetaData meta; - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker broker; - stk::io::fill_mesh_with_auto_decomp(outputMeshName, bulk, broker); - verify_input_transient_fields(broker, expectedTimeSteps); + MeshFromFile& get_balanced_mesh() + { + if (m_balancedMeshPtr->is_empty()) read_balanced_mesh(); + return *m_balancedMeshPtr; + } - stk::balance::ParsedOptions options(outputMeshName, "."); - options.set_app_type_default(stk::balance::SD_DEFAULTS); - stk::balance::run_stk_rebalance(options, MPI_COMM_WORLD); + void cleanup_files() + { + unlink_serial_file(fileBaseName); + unlink_parallel_file(fileBaseName); + } - verify_expected_time_steps_in_parallel_file(outputMeshName, expectedTimeSteps); +private: + void read_initial_mesh() + { + m_initialMeshPtr->fill_from_serial(fileBaseName); + } - unlink_serial_file(outputMeshName.c_str()); - unlink_parallel_file(outputMeshName.c_str()); - } -} + void read_balanced_mesh() + { + m_balancedMeshPtr->fill_from_parallel(fileBaseName); + } -void check_sideset_orientation(const stk::mesh::BulkData& bulk, - const std::vector & sidesets, - const stk::mesh::EntityId expectedId, - const stk::mesh::ConnectivityOrdinal expectedOrdinal) + void unlink_parallel_file(const std::string& baseName) + { + std::string parallelName = stk::io::construct_parallel_filename(baseName, get_parallel_size(), get_parallel_rank()); + unlink_serial_file(parallelName); + } + + void unlink_serial_file(const std::string& fileName) + { + unlink(fileName.c_str()); + } + +protected: + const std::string fileBaseName; + StkBalanceRunner balanceRunner; + TransientWriter writer; + const TransientVerifier verifier; + +private: + const std::unique_ptr m_initialMeshPtr; + const std::unique_ptr m_balancedMeshPtr; +}; + +TEST_F(TransientFieldBalance, verifyStaticDataTransfer) { - const stk::mesh::SideSetEntry sideSetEntry = (*sidesets[0])[0]; - EXPECT_EQ(expectedId, bulk.identifier(sideSetEntry.element)); - EXPECT_EQ(expectedOrdinal, sideSetEntry.side); + if (get_parallel_size() == 2 || + get_parallel_size() == 4 || + get_parallel_size() == 8 || + get_parallel_size() == 16) { + + writer.write_static_mesh("1x4x4"); + + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_time_steps(initialMesh, {}); + verifier.verify_num_transient_fields(initialMesh, 0u); + + balanceRunner.run_end_to_end(); + + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_time_steps(balancedMesh, {}); + verifier.verify_num_transient_fields(balancedMesh, 0u); + + cleanup_files(); + } } -TEST(TestTransientFieldBalance, verifyTransientDataTransferWithSidesets) +TEST_F(TransientFieldBalance, verifyNumberOfSteps) { - std::string serialOutputMeshName = "two_hex_transient.e"; - std::string parallelOutputMeshName = "two_hex_transient_balanced.e"; - const std::vector expectedTimeSteps = {0.0}; - const stk::mesh::EntityId expectedId = 1; - const stk::mesh::ConnectivityOrdinal expectedOrdinal = 5; + if (get_parallel_size() != 2) return; - if (stk::parallel_machine_size(MPI_COMM_WORLD) == 2) - { - // Build the target serial mesh and write it to a file - if (stk::parallel_machine_rank(MPI_COMM_WORLD) == 0) - { - stk::mesh::MetaData meta(3); - stk::mesh::BulkData bulk(meta, MPI_COMM_SELF); + const std::vector timeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; - stk::unit_test_util::create_AB_mesh_with_sideset_and_field(bulk, stk::unit_test_util::LEFT, stk::unit_test_util::INCREASING, "dummyField"); - stk::io::write_mesh_with_fields(serialOutputMeshName, bulk, 1, 1.0); + writer.set_time_steps(timeSteps); + writer.write_transient_mesh("1x1x20"); - std::vector sidesets = bulk.get_sidesets(); - ASSERT_EQ(1u, sidesets.size()); - EXPECT_EQ(1u, sidesets[0]->size()); + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_time_steps(initialMesh, timeSteps); + verifier.verify_num_transient_fields(initialMesh, 2u); - check_sideset_orientation(bulk, sidesets, expectedId, expectedOrdinal); - } + balanceRunner.run_end_to_end(); + + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_time_steps(balancedMesh, timeSteps); + + cleanup_files(); +} - // Read, balance, and write the mesh on two processors - stk::balance::BasicZoltan2Settings rcbOptions; - stk::balance::run_stk_balance_with_settings(parallelOutputMeshName, serialOutputMeshName, MPI_COMM_WORLD, rcbOptions); +TEST_F(TransientFieldBalance, verifyGlobalVariable) +{ + if (get_parallel_size() != 2) return; - stk::mesh::MetaData meta(3); - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); + const std::vector timeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; + const std::string globalVariableName = "test_time"; - // Read the decomposed mesh to make sure the sideset maintains its polarity - stk::io::fill_mesh(parallelOutputMeshName, bulk); + writer.set_time_steps(timeSteps); + writer.set_global_variable_name(globalVariableName); + writer.write_transient_mesh("1x1x20"); - std::vector sidesets = bulk.get_sidesets(); - if ((sidesets.size() == 1u) && (sidesets[0]->size() == 1u)) - { - check_sideset_orientation(bulk, sidesets, expectedId, expectedOrdinal); - } - } - unlink_serial_file(serialOutputMeshName); - unlink_parallel_file(parallelOutputMeshName); + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_time_steps(initialMesh, timeSteps); + verifier.verify_num_transient_fields(initialMesh, 2u); + + balanceRunner.run_end_to_end(); + + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_global_variables_at_each_time_step(balancedMesh, globalVariableName, timeSteps); + + cleanup_files(); } -TEST(TestTransientFieldBalance, verifyTransientDataTransferWithSidesetsOnMovedElements) +TEST_F(TransientFieldBalance, verifyTransientDataTransferOnFourProcessors) { - std::string serialOutputMeshName = "two_hex_transient.e"; - std::string parallelOutputMeshName = "two_hex_transient_balanced.e"; - const std::vector expectedTimeSteps = {0.0}; - const stk::mesh::EntityId expectedId = 2; - const stk::mesh::ConnectivityOrdinal expectedOrdinal = 5; + if (get_parallel_size() != 4) return; - if (stk::parallel_machine_size(MPI_COMM_WORLD) == 2) - { - // Build the target serial mesh and write it to a file - if (stk::parallel_machine_rank(MPI_COMM_WORLD) == 0) - { - stk::mesh::MetaData meta(3); - stk::mesh::BulkData bulk(meta, MPI_COMM_SELF); + const std::vector timeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; + const std::string fieldName = "myField"; - stk::unit_test_util::create_AB_mesh_with_sideset_and_field(bulk, stk::unit_test_util::LEFT, stk::unit_test_util::DECREASING, "dummyField"); - stk::io::write_mesh_with_fields(serialOutputMeshName, bulk, 1, 1.0); + writer.set_time_steps(timeSteps); + writer.set_field_name(fieldName); + writer.write_transient_mesh("1x4x4"); - std::vector sidesets = bulk.get_sidesets(); - ASSERT_EQ(1u, sidesets.size()); - ASSERT_EQ(1u, sidesets[0]->size()); + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_time_steps(initialMesh, timeSteps); + verifier.verify_num_transient_fields(initialMesh, 2u); + verifier.verify_transient_field_names(initialMesh, fieldName); - check_sideset_orientation(bulk, sidesets, expectedId, expectedOrdinal); - } + balanceRunner.run_end_to_end(); - { - // Read, balance, and write the mesh on two processors - stk::balance::BasicZoltan2Settings rcbOptions; - const std::string initialDecompMethod = "BLOCK"; - stk::mesh::MetaData meta; - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); - stk::balance::initial_decomp_and_balance(bulk, rcbOptions, serialOutputMeshName, parallelOutputMeshName, initialDecompMethod); - } + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_time_steps(balancedMesh, timeSteps); + verifier.verify_num_transient_fields(balancedMesh, 2u); + verifier.verify_transient_field_names(balancedMesh, fieldName); - { - stk::mesh::MetaData meta(3); - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); + verifier.compare_entity_rank_names(initialMesh, balancedMesh); - // Read the decomposed mesh to make sure the sideset maintains its polarity - stk::io::fill_mesh(parallelOutputMeshName, bulk); + verifier.verify_transient_fields(balancedMesh); - std::vector sidesets = bulk.get_sidesets(); - if (stk::parallel_machine_rank(MPI_COMM_WORLD) == 0) { - ASSERT_EQ(1u, sidesets.size()); - ASSERT_EQ(1u, sidesets[0]->size()); - check_sideset_orientation(bulk, sidesets, expectedId, expectedOrdinal); - } - } - } - unlink_serial_file(serialOutputMeshName); - unlink_parallel_file(parallelOutputMeshName); + cleanup_files(); } -TEST(TestTransientFieldBalance, verifyTransientDataTransferOnFourProcessors) +TEST_F(TransientFieldBalance, verifyTransientDataTransferWithSidesets) { - std::string serialOutputMeshName = "sixteen_hex_transient.e"; - std::string parallelOutputMeshName = "sixteen_hex_transient_balanced.e"; - std::string fieldBaseName = "myField"; - std::string globalVariableName = "TestTime"; - const std::vector expectedTimeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; + if (get_parallel_size() != 2) return; - if(stk::parallel_machine_size(MPI_COMM_WORLD) == 4) - { - stk::unit_test_util::IdAndTimeFieldValueSetter fieldSetter; - stk::unit_test_util::generated_mesh_with_transient_data_to_file_in_serial("1x4x4", serialOutputMeshName, fieldBaseName, globalVariableName, expectedTimeSteps, fieldSetter); + const stk::mesh::EntityId expectedId = 1; + const stk::mesh::ConnectivityOrdinal expectedOrdinal = 5; - stk::balance::BasicZoltan2Settings rcbOptions; - stk::balance::run_stk_balance_with_settings(parallelOutputMeshName, serialOutputMeshName, MPI_COMM_WORLD, rcbOptions); + writer.write_two_element_mesh_with_sideset(stk::unit_test_util::INCREASING); - verify_transient_transfer(serialOutputMeshName, parallelOutputMeshName, fieldBaseName, expectedTimeSteps); + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_sideset_orientation(initialMesh, expectedId, expectedOrdinal); - unlink_serial_file(serialOutputMeshName); - unlink_parallel_file(parallelOutputMeshName); - } + balanceRunner.run_end_to_end(); + + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_sideset_orientation(balancedMesh, expectedId, expectedOrdinal); + + cleanup_files(); } -TEST(TestTransientFieldBalance, verifyStaticDataTransfer) +TEST_F(TransientFieldBalance, verifyTransientDataTransferWithSidesetsOnMovedElements) { - std::string serialOutputMeshName = "sixteen_hex_static.e"; - std::string parallelOutputMeshName = "sixteen_hex_static_balanced.e"; + if (get_parallel_size() != 2) return; - if((stk::parallel_machine_size(MPI_COMM_WORLD) == 2) || - (stk::parallel_machine_size(MPI_COMM_WORLD) == 4) || - (stk::parallel_machine_size(MPI_COMM_WORLD) == 8) || - (stk::parallel_machine_size(MPI_COMM_WORLD) == 16)) - { - stk::unit_test_util::generated_mesh_to_file_in_serial("1x4x4", serialOutputMeshName); + const stk::mesh::EntityId expectedId = 2; + const stk::mesh::ConnectivityOrdinal expectedOrdinal = 5; - stk::balance::BasicZoltan2Settings rcbOptions; - stk::balance::run_stk_balance_with_settings(parallelOutputMeshName, serialOutputMeshName, MPI_COMM_WORLD, rcbOptions); + writer.write_two_element_mesh_with_sideset(stk::unit_test_util::DECREASING); - verify_static_transfer(serialOutputMeshName, parallelOutputMeshName); + MeshFromFile& initialMesh = get_initial_mesh(); + verifier.verify_sideset_orientation(initialMesh, expectedId, expectedOrdinal); - unlink_serial_file(serialOutputMeshName); - unlink_parallel_file(parallelOutputMeshName); - } -} + balanceRunner.set_initial_decomp_method("BLOCK"); + balanceRunner.run_end_to_end(); -TEST(TestTransientFieldBalance, verifyGlobalVariable) -{ - if(stk::parallel_machine_size(MPI_COMM_WORLD) == 2) - { - std::string outputMeshName = "twenty_hex_transient.e"; - std::string fieldBaseName = "myField"; - const std::vector expectedTimeSteps = {0.0, 1.0, 2.0, 3.0, 4.0}; - std::string globalVariableName = "test_time"; - stk::unit_test_util::IdAndTimeFieldValueSetter fieldSetter; - stk::unit_test_util::generated_mesh_with_transient_data_to_file_in_serial("1x1x20", outputMeshName, fieldBaseName, globalVariableName, expectedTimeSteps, fieldSetter); - - stk::mesh::MetaData meta; - stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); - stk::io::StkMeshIoBroker broker; - stk::io::fill_mesh_with_auto_decomp(outputMeshName, bulk, broker); - verify_input_transient_fields(broker, expectedTimeSteps); - - stk::balance::ParsedOptions options(outputMeshName, "."); - options.set_app_type_default(stk::balance::SD_DEFAULTS); - stk::balance::run_stk_rebalance(options, MPI_COMM_WORLD); - - verify_expected_global_variable_in_parallel_file(outputMeshName, expectedTimeSteps, globalVariableName); - - unlink_serial_file(outputMeshName.c_str()); - unlink_parallel_file(outputMeshName.c_str()); - } + MeshFromFile& balancedMesh = get_balanced_mesh(); + verifier.verify_sideset_orientation(balancedMesh, expectedId, expectedOrdinal); + + cleanup_files(); } TEST(TestTransientFieldBalance, verifyThrowIfInputFileEqualsOutputFile) @@ -434,7 +521,7 @@ TEST(TestTransientFieldBalance, verifyThrowIfInputFileEqualsOutputFile) stk::balance::BasicZoltan2Settings rcbOptions; EXPECT_THROW(stk::balance::run_stk_balance_with_settings(parallelOutputMeshName, serialMeshName, MPI_COMM_WORLD, rcbOptions), std::logic_error); - unlink_serial_file(serialMeshName); + unlink(serialMeshName.c_str()); } } @@ -449,8 +536,6 @@ TEST(TestTransientFieldBalance, verifyThrowIfInputFileEqualsDotSlashOutputFile) stk::balance::BasicZoltan2Settings rcbOptions; EXPECT_THROW(stk::balance::run_stk_balance_with_settings(parallelOutputMeshName, serialMeshName, MPI_COMM_WORLD, rcbOptions), std::logic_error); - unlink_serial_file(serialMeshName); + unlink(serialMeshName.c_str()); } } - -} diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMeshManagerTest.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMeshManagerTest.cpp index f259da6a502f..3cc73318adca 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMeshManagerTest.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMeshManagerTest.cpp @@ -1,5 +1,4 @@ #include -//#include #include #include #include @@ -35,10 +34,10 @@ TEST(NgpMeshManager, NgpMesh_FromBulkData) stk::mesh::NgpMesh & ngpMesh = bulk.get_updated_ngp_mesh(); -#ifndef KOKKOS_ENABLE_CUDA - EXPECT_TRUE((std::is_same::value)); -#else +#ifdef STK_USE_DEVICE_MESH EXPECT_TRUE((std::is_same::value)); +#else + EXPECT_TRUE((std::is_same::value)); #endif } @@ -48,10 +47,10 @@ TEST(NgpMeshManager, NgpMesh_Update) stk::mesh::BulkData bulk(meta, MPI_COMM_WORLD); stk::io::fill_mesh("generated:1x1x4", bulk); -#ifndef KOKKOS_ENABLE_CUDA - stk::mesh::NgpMeshManager * ngpMeshManager = new stk::mesh::HostMeshManager(bulk); -#else +#ifdef STK_USE_DEVICE_MESH stk::mesh::NgpMeshManager * ngpMeshManager = new stk::mesh::DeviceMeshManager(bulk); +#else + stk::mesh::NgpMeshManager * ngpMeshManager = new stk::mesh::HostMeshManager(bulk); #endif bulk.modification_begin(); @@ -59,7 +58,7 @@ TEST(NgpMeshManager, NgpMesh_Update) stk::mesh::NgpMesh& ngpMesh = ngpMeshManager->get_mesh(); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_FALSE(ngpMesh.is_up_to_date()); #else EXPECT_TRUE(ngpMesh.is_up_to_date()); diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMultistateFieldTest.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMultistateFieldTest.cpp index bf4337e3ab67..895be934e22f 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMultistateFieldTest.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpMultistateFieldTest.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -24,7 +25,7 @@ void assign_value_to_field(stk::mesh::BulkData &bulk, Field ngpField, double val }); ngpField.modify_on_device(); } -void assign_value_to_statedfield(stk::mesh::BulkData &bulk, stk::mesh::ConvenientNgpMultistateField ngpStatedField, double value) +void assign_value_to_statedfield(stk::mesh::BulkData &bulk, stk::mesh::NgpMultistateField ngpStatedField, double value) { stk::mesh::NgpMesh & ngpMesh = bulk.get_updated_ngp_mesh(); stk::mesh::for_each_entity_run(ngpMesh, stk::topology::ELEM_RANK, bulk.mesh_meta_data().universal_part(), @@ -51,11 +52,31 @@ class StatedFields : public stk::unit_test_util::MeshFixture }); } + template + void test_field_state_has_value_on_device(stk::mesh::BulkData& bulk, Field ngpMultiStateField, stk::mesh::FieldState state, double expectedValue) + { + stk::mesh::NgpMesh& ngpMesh = bulk.get_updated_ngp_mesh(); + + if(state == stk::mesh::StateNew) { + stk::mesh::for_each_entity_run(ngpMesh, stk::topology::ELEM_RANK, bulk.mesh_meta_data().universal_part(), + KOKKOS_LAMBDA(const stk::mesh::FastMeshIndex& entity) + { + NGP_ThrowRequire(expectedValue == ngpMultiStateField.get_new(entity, 0)); + }); + } else { + stk::mesh::for_each_entity_run(ngpMesh, stk::topology::ELEM_RANK, bulk.mesh_meta_data().universal_part(), + KOKKOS_LAMBDA(const stk::mesh::FastMeshIndex& entity) + { + NGP_ThrowRequire(expectedValue == ngpMultiStateField.get_old(state, entity, 0)); + }); + } + } + protected: - stk::mesh::ConvenientNgpMultistateField create_field_with_num_states(unsigned numStates) + stk::mesh::NgpMultistateField create_field_with_num_states(unsigned numStates) { make_mesh_with_field(numStates); - return stk::mesh::ConvenientNgpMultistateField(get_bulk(), *stkField); + return stk::mesh::NgpMultistateField(get_bulk(), *stkField); } void make_mesh_with_field(unsigned numStates) { @@ -193,9 +214,25 @@ TEST_F(StatedFields, field_swap) test_field_has_value_on_device(get_bulk(), ngpFieldNew, 1); test_field_has_value_on_device(get_bulk(), ngpFieldOld, 0); } + +TEST_F(StatedFields, new_field_swap) +{ + stk::mesh::NgpMultistateField ngpStatedField = create_field_with_num_states(2); + verify_can_assign_values_per_state(ngpStatedField); + + test_field_state_has_value_on_device(get_bulk(), ngpStatedField, stk::mesh::StateNew, 0); + test_field_state_has_value_on_device(get_bulk(), ngpStatedField, stk::mesh::StateOld, 1); + + get_bulk().update_field_data_states(); + ngpStatedField.increment_state(); + + test_field_state_has_value_on_device(get_bulk(), ngpStatedField, stk::mesh::StateNew, 1); + test_field_state_has_value_on_device(get_bulk(), ngpStatedField, stk::mesh::StateOld, 0); +} + TEST_F(StatedFields, gettingFieldData_direct) { - stk::mesh::ConvenientNgpMultistateField ngpStatedField = create_field_with_num_states(3); + stk::mesh::NgpMultistateField ngpStatedField = create_field_with_num_states(3); verify_initial_value_set_per_state(ngpStatedField); assign_value_to_statedfield(get_bulk(), ngpStatedField, 2); @@ -216,7 +253,7 @@ TEST_F(StatedFields, gettingFieldData_direct) } TEST_F(StatedFields, updateFromHostStkField) { - stk::mesh::ConvenientNgpMultistateField ngpStatedField = create_field_with_num_states(1); + stk::mesh::NgpMultistateField ngpStatedField = create_field_with_num_states(1); verify_initial_value_set_per_state(ngpStatedField); stk::mesh::NgpField ngpField = ngpStatedField.get_field_new_state(); diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpParallelSumTest.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpParallelSumTest.cpp index 5e96989e93a4..b93778119f25 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpParallelSumTest.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpParallelSumTest.cpp @@ -138,7 +138,7 @@ NGP_TEST_F(NgpParallelSum, simpleVersion) check_field_on_device(ngpMesh, deviceUserField, deviceGoldValues); } -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH NGP_TEST_F(NgpParallelSum, simpleVersion_noSyncToDeviceAfterwards) { const double initValue = 0.0; @@ -186,7 +186,7 @@ NGP_TEST_F(NgpCopyOwnedToShared, simpleVersion) check_field_on_device(ngpMesh, deviceUserField, deviceGoldValues); } -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH NGP_TEST_F(NgpCopyOwnedToShared, simpleVersion_noSyncToDeviceAfterwards) { const double initValue = 0.0; @@ -234,7 +234,7 @@ NGP_TEST_F(NgpCommunicateFieldData, simpleVersion_takesGhosting) check_field_on_device(ngpMesh, deviceUserField, deviceGoldValues); } -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH NGP_TEST_F(NgpCommunicateFieldData, simpleVersion_takesGhosting_noSyncToDeviceAfterwards) { const double initValue = 0.0; @@ -282,7 +282,7 @@ NGP_TEST_F(NgpCommunicateFieldData, simpleVersion_takesBulkData) check_field_on_device(ngpMesh, deviceUserField, deviceGoldValues); } -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH NGP_TEST_F(NgpCommunicateFieldData, simpleVersion_takesBulkData_noSyncToDeviceAfterwards) { const double initValue = 0.0; diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/TestNgpMeshUpdate.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/TestNgpMeshUpdate.cpp index 50a745492af8..15764d71a2af 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/TestNgpMeshUpdate.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/TestNgpMeshUpdate.cpp @@ -34,7 +34,7 @@ TEST_F(UpdateNgpMesh, lazyAutoUpdate) get_bulk().modification_begin(); get_bulk().modification_end(); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_FALSE(ngpMesh->is_up_to_date()); ngpMesh = &get_bulk().get_updated_ngp_mesh(); // Trigger update EXPECT_TRUE(ngpMesh->is_up_to_date()); @@ -56,7 +56,7 @@ TEST_F(UpdateNgpMesh, manualUpdate) get_bulk().modification_begin(); get_bulk().modification_end(); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_FALSE(ngpMesh.is_up_to_date()); get_bulk().update_ngp_mesh(); EXPECT_TRUE(ngpMesh.is_up_to_date()); @@ -75,7 +75,7 @@ TEST_F(UpdateNgpMesh, OnlyOneDeviceMesh_InternalAndExternal) // Create first NgpMesh inside BulkData get_bulk().get_updated_ngp_mesh(); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_THROW(stk::mesh::NgpMesh secondNgpMesh(get_bulk()), std::logic_error); #else EXPECT_NO_THROW(stk::mesh::NgpMesh secondNgpMesh(get_bulk())); @@ -88,7 +88,7 @@ TEST_F(UpdateNgpMesh, OnlyOneDeviceMesh_TwoExternal) stk::mesh::NgpMesh firstNgpMesh(get_bulk()); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_THROW(stk::mesh::NgpMesh secondNgpMesh(get_bulk()), std::logic_error); #else EXPECT_NO_THROW(stk::mesh::NgpMesh secondNgpMesh(get_bulk())); diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/howToNgp.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/howToNgp.cpp index 63ba101dc791..9bdeca25b859 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/howToNgp.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/howToNgp.cpp @@ -708,19 +708,14 @@ unsigned count_num_elems(stk::mesh::NgpMesh ngpMesh, stk::mesh::EntityRank rank, stk::mesh::Part &part) { -#ifdef KOKKOS_ENABLE_CUDA - using DeviceSpace = Kokkos::CudaSpace; -#else - using DeviceSpace = Kokkos::HostSpace; -#endif - Kokkos::View numElems("numElems", 1); + Kokkos::View numElems("numElems", 1); stk::mesh::for_each_entity_run(ngpMesh, rank, part, KOKKOS_LAMBDA(const stk::mesh::FastMeshIndex& entity) { unsigned fieldValue = static_cast(ngpField(entity, 0)); Kokkos::atomic_add(&numElems(0), fieldValue); }); - Kokkos::View::HostMirror numElemsHost = + Kokkos::View::HostMirror numElemsHost = Kokkos::create_mirror_view(numElems); Kokkos::deep_copy(numElemsHost, numElems); return numElemsHost(0); @@ -809,49 +804,10 @@ void verify_state_new_has_value(stk::mesh::BulkData &bulk, } } -void set_new_as_old_plus_one_on_device(stk::mesh::NgpMesh &ngpMesh, - stk::mesh::EntityRank rank, - stk::mesh::Selector sel, - stk::mesh::NgpMultistateField &ngpMultistateField) -{ - ngpMultistateField.sync_to_device(); - int component = 0; - stk::mesh::for_each_entity_run(ngpMesh, rank, sel, - KOKKOS_LAMBDA(const stk::mesh::FastMeshIndex& entity) - { - stk::mesh::NgpConstField stateNField = ngpMultistateField.get_field_old_state(stk::mesh::StateN); - stk::mesh::NgpField stateNp1Field = ngpMultistateField.get_field_new_state(); - stateNp1Field(entity, component) = stateNField(entity, component) + 1; - }); - ngpMultistateField.modify_on_device(); -} - -TEST_F(NgpHowTo, useMultistateFields) -{ - int numStates = 2; - int initialValue = 0; - stk::mesh::Field &stkField = create_field_with_num_states_and_init(get_meta(), "myField", numStates, initialValue); - setup_mesh("generated:1x1x4", stk::mesh::BulkData::AUTO_AURA); - - stk::mesh::NgpMultistateField ngpMultistateField(get_bulk(), stkField); - stk::mesh::NgpMesh & ngpMesh = get_bulk().get_updated_ngp_mesh(); - - set_new_as_old_plus_one_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); - int newStateValue = 1; - verify_state_new_has_value(get_bulk(), ngpMesh, stkField, ngpMultistateField, newStateValue); - - get_bulk().update_field_data_states(); - ngpMultistateField.increment_state(); - - set_new_as_old_plus_one_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); - newStateValue = 2; - verify_state_new_has_value(get_bulk(), ngpMesh, stkField, ngpMultistateField, newStateValue); -} - -void set_new_as_old_plus_one_in_convenient_field_on_device(stk::mesh::NgpMesh &ngpMesh, +void set_new_as_old_plus_one_in_field_on_device(stk::mesh::NgpMesh &ngpMesh, stk::mesh::EntityRank rank, stk::mesh::Selector sel, - stk::mesh::ConvenientNgpMultistateField &ngpMultistateField) + stk::mesh::NgpMultistateField &ngpMultistateField) { int component = 0; stk::mesh::for_each_entity_run(ngpMesh, rank, sel, @@ -862,24 +818,24 @@ void set_new_as_old_plus_one_in_convenient_field_on_device(stk::mesh::NgpMesh &n ngpMultistateField.modify_on_device(); } -TEST_F(NgpHowTo, useConvenientMultistateFields) +TEST_F(NgpHowTo, useMultistateFields) { int numStates = 2; int initialValue = 0; stk::mesh::Field &stkField = create_field_with_num_states_and_init(get_meta(), "myField", numStates, initialValue); setup_mesh("generated:1x1x4", stk::mesh::BulkData::AUTO_AURA); - stk::mesh::ConvenientNgpMultistateField ngpMultistateField(get_bulk(), stkField); - stk::mesh::NgpMesh & ngpMesh = get_bulk().get_updated_ngp_mesh(); + stk::mesh::NgpMultistateField ngpMultistateField(get_bulk(), stkField); + stk::mesh::NgpMesh ngpMesh(get_bulk()); - set_new_as_old_plus_one_in_convenient_field_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); + set_new_as_old_plus_one_in_field_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); int newStateValue = 1; verify_state_new_has_value(get_bulk(), ngpMesh, stkField, ngpMultistateField, newStateValue); get_bulk().update_field_data_states(); ngpMultistateField.increment_state(); - set_new_as_old_plus_one_in_convenient_field_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); + set_new_as_old_plus_one_in_field_on_device(ngpMesh, stk::topology::ELEM_RANK, get_meta().universal_part(), ngpMultistateField); newStateValue = 2; verify_state_new_has_value(get_bulk(), ngpMesh, stkField, ngpMultistateField, newStateValue); } diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/ngpFieldTest.cpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/ngpFieldTest.cpp index 887285b0981d..0165da453a55 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/ngpFieldTest.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/ngpFieldTest.cpp @@ -202,7 +202,7 @@ TEST_F(NgpFieldFixture, GetNgpField) stk::mesh::NgpField & ngpIntField = stk::mesh::get_updated_ngp_field(stkIntField); EXPECT_EQ(ngpIntField.get_ordinal(), stkIntField.mesh_meta_data_ordinal()); -#ifdef KOKKOS_ENABLE_CUDA +#ifdef STK_USE_DEVICE_MESH EXPECT_TRUE((std::is_same&>::value)); #else EXPECT_TRUE((std::is_same&>::value)); @@ -250,12 +250,14 @@ TEST_F(NgpFieldFixture, ModifyAndSync) sync_field_to_host(stkIntField); check_field_on_host(get_bulk(), stkIntField, multiplier*multiplier); +#ifdef STK_USE_DEVICE_MESH size_t expectedSyncsToDevice = 2; size_t expectedSyncsToHost = 1; -#ifndef KOKKOS_ENABLE_CUDA - expectedSyncsToDevice = 0; - expectedSyncsToHost = 0; +#else + size_t expectedSyncsToDevice = 0; + size_t expectedSyncsToHost = 0; #endif + EXPECT_EQ(expectedSyncsToDevice, stkIntField.num_syncs_to_device()); EXPECT_EQ(expectedSyncsToHost, stkIntField.num_syncs_to_host()); @@ -289,11 +291,12 @@ TEST_F(NgpFieldFixture, UpdateNgpFieldAfterMeshMod_WithMostCurrentDataOnHost) sync_field_to_host(stkIntField); check_field_on_host(get_bulk(), stkIntField, multiplier*multiplier); +#ifdef STK_USE_DEVICE_MESH size_t expectedSyncsToDevice = 2; size_t expectedSyncsToHost = 1; -#ifndef KOKKOS_ENABLE_CUDA - expectedSyncsToDevice = 0; - expectedSyncsToHost = 0; +#else + size_t expectedSyncsToDevice = 0; + size_t expectedSyncsToHost = 0; #endif EXPECT_EQ(expectedSyncsToDevice, stkIntField.num_syncs_to_device()); diff --git a/packages/stk/stk_unit_tests/stk_ngp/NgpFieldPerformanceTest.cpp b/packages/stk/stk_unit_tests/stk_ngp/NgpFieldPerformanceTest.cpp index 0b0cb36f2282..cc9e74cabd1e 100644 --- a/packages/stk/stk_unit_tests/stk_ngp/NgpFieldPerformanceTest.cpp +++ b/packages/stk/stk_unit_tests/stk_ngp/NgpFieldPerformanceTest.cpp @@ -121,7 +121,8 @@ TEST_F(NgpFieldPerf, constFieldDataAccessIsFasterThanFieldDataAccess) warm_up_gpu(); double constTime = time_field_data_access>(); double nonConstTime = time_field_data_access>(); - std::cerr << "non-const time: " << nonConstTime << ", const time: " << constTime << '\n'; + double diffPerf = (nonConstTime - constTime) * 100 / nonConstTime; + std::cerr << "non-const time: " << nonConstTime << ", const time: " << constTime << " percentage gain: " << diffPerf << "%%\n"; //only expect nonConstTime to be faster in release-mode, on cude, for large problem sizes if (stk::unit_test_util::get_command_line_option("-dim",20) >= 60) diff --git a/packages/stk/stk_unit_tests/stk_util/diag/UnitTestTimer.cpp b/packages/stk/stk_unit_tests/stk_util/diag/UnitTestTimer.cpp index 9f267c368f5b..c11ff5593ef2 100644 --- a/packages/stk/stk_unit_tests/stk_util/diag/UnitTestTimer.cpp +++ b/packages/stk/stk_unit_tests/stk_util/diag/UnitTestTimer.cpp @@ -33,7 +33,8 @@ // #include // for size_t -#include // for sleep +#include +#include #include // for sin #include // for ostringstream, etc #include // for printTimersTable @@ -77,7 +78,7 @@ enum { namespace { double -work(int repetitions = 100000) +work(int repetitions = 1000) { double x = 1.0; @@ -174,21 +175,21 @@ TEST(UnitTestTimer, UnitTest) std::ostringstream oss; oss << x << std::endl; - ::sleep(1); + std::this_thread::sleep_for(std::chrono::milliseconds(10)); lap_timer.lap(); stk::diag::MetricTraits::Type lap_time = lap_timer.getMetric().getLap(); - ASSERT_TRUE(lap_time >= 1.0); + EXPECT_NEAR(0.01, lap_time, 1.0e-3); - ::sleep(1); + std::this_thread::sleep_for(std::chrono::milliseconds(10)); lap_timer.stop(); lap_time = lap_timer.getMetric().getLap(); - ASSERT_TRUE(lap_time >= 2.0); + EXPECT_NEAR(0.02, lap_time, 1.0e-3); } { @@ -216,7 +217,7 @@ TEST(UnitTestTimer, UnitTest) stk::diag::TimeBlock _time2(second_timer_on); stk::diag::TimeBlock _time3(second_timer_off); - ::sleep(1); + std::this_thread::sleep_for(std::chrono::milliseconds(10)); } // Grab previous subtimer and run 100 laps @@ -311,10 +312,6 @@ TEST(UnitTestTimer, UnitTest) object_vector[j].run(); stk::diag::printTimersTable(strout, unitTestTimer(), stk::diag::METRICS_ALL, true); - - std::cout << strout.str() << std::endl; - -// dw().m(LOG_TIMER) << strout.str() << stk::diag::dendl; } } @@ -348,5 +345,4 @@ TEST(UnitTestTimer, YuugeNumberOfTimers) } stk::diag::printTimersTable(strout, rootTimer, stk::diag::METRICS_ALL, false, MPI_COMM_WORLD); - std::cout << strout.str() << std::endl; } diff --git a/packages/stk/stk_util/stk_util/stk_kokkos_macros.h b/packages/stk/stk_util/stk_util/stk_kokkos_macros.h index 6ef567e0f8d0..70e1e95573a8 100644 --- a/packages/stk/stk_util/stk_util/stk_kokkos_macros.h +++ b/packages/stk/stk_util/stk_util/stk_kokkos_macros.h @@ -1,21 +1,19 @@ #ifndef STK_KOKKOS_MACROS_H #define STK_KOKKOS_MACROS_H -//When we are ready to depend on kokkos we will #include Kokkos_Macros.hpp -//and define our STK_FUNCTION macros in terms of the corresponding kokkos macros. -//Until then, do our own (duplicated) implementation. - -#if defined(__CUDA_ARCH__) - -#define STK_INLINE_FUNCTION __device__ __host__ inline -#define STK_FUNCTION __device__ __host__ - -#else - -#define STK_INLINE_FUNCTION inline -#define STK_FUNCTION - +#include + +// This should eventually need to be supplemented with checks for ROCM and +// other accelerator platforms +// +#ifdef KOKKOS_ENABLE_CUDA + #ifndef STK_USE_DEVICE_MESH + #define STK_USE_DEVICE_MESH + #endif #endif +#define STK_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION +#define STK_FUNCTION KOKKOS_FUNCTION + #endif From 3ea6bf2027bc54872a26d4eed9ac89ee2fe08f8b Mon Sep 17 00:00:00 2001 From: Alan Williams Date: Wed, 1 Apr 2020 14:46:35 -0600 Subject: [PATCH 2/2] Percept snapshot as of 4/01/2020 --- packages/percept/src/adapt/NodeRegistry.cpp | 71 -- packages/percept/src/adapt/NodeRegistry.hpp | 5 - .../percept/src/adapt/NodeRegistry_KOKKOS.cpp | 914 ------------------ .../percept/src/adapt/NodeRegistry_KOKKOS.hpp | 448 --------- packages/percept/src/adapt/Refiner.cpp | 2 +- .../src/adapt/UniformRefinerPattern.cpp | 4 +- packages/percept/src/adapt/main/MeshAdapt.cpp | 13 +- .../src/adapt/sierra_element/CellTopology.cpp | 321 ------ .../src/adapt/sierra_element/CellTopology.hpp | 37 - .../adapt/sierra_element/RefinementKey.cpp | 99 -- .../adapt/sierra_element/RefinementKey.hpp | 79 -- .../sierra_element/RefinementTopology.cpp | 699 -------------- .../sierra_element/RefinementTopology.hpp | 82 -- .../sierra_element/StdMeshObjTopologies.cpp | 206 ---- packages/percept/src/percept/FieldBLAS.hpp | 4 +- .../percept/src/percept/GeometryVerifier.cpp | 6 +- packages/percept/src/percept/PerceptMesh.cpp | 71 +- packages/percept/src/percept/PerceptMesh.hpp | 13 - .../percept/src/percept/TopologyVerifier.cpp | 2 +- .../src/percept/eigen_verify/EigenVerify.cpp | 2 +- .../src/percept/function/FieldFunction.hpp | 5 +- .../src/percept/function/FunctionOperator.hpp | 4 +- .../function/internal/IntegratedOp.hpp | 10 +- .../recovery/GeometryRecoverySplineFit.cpp | 2 +- .../mesh/geometry/volume/VolumeUtil.cpp | 1 - .../mesh/geometry/volume/VolumeUtil.hpp | 66 -- .../mesh/mod/smoother/JacobianUtil.cpp | 55 -- .../mesh/mod/smoother/JacobianUtil.hpp | 67 -- .../mesh/mod/smoother/SGridJacobianUtil.hpp | 3 - .../mesh/mod/smoother/SmootherMetric.hpp | 97 -- .../mesh_transfer/RotationTranslation.hpp | 4 +- packages/percept/src/percept/norm/H1Norm.hpp | 5 +- .../src/percept/norm/IntrepidManager.cpp | 2 +- packages/percept/src/percept/norm/Norm.hpp | 14 +- .../stk_rebalance_utils/RebalanceUtils.cpp | 2 +- 35 files changed, 47 insertions(+), 3368 deletions(-) delete mode 100644 packages/percept/src/adapt/NodeRegistry_KOKKOS.cpp delete mode 100644 packages/percept/src/adapt/NodeRegistry_KOKKOS.hpp delete mode 100644 packages/percept/src/adapt/sierra_element/RefinementKey.cpp delete mode 100644 packages/percept/src/adapt/sierra_element/RefinementKey.hpp diff --git a/packages/percept/src/adapt/NodeRegistry.cpp b/packages/percept/src/adapt/NodeRegistry.cpp index d4fd35cf9312..60cea730238d 100644 --- a/packages/percept/src/adapt/NodeRegistry.cpp +++ b/packages/percept/src/adapt/NodeRegistry.cpp @@ -3185,75 +3185,4 @@ } } } - - bool NodeRegistry::verifyAllKeysInKokkosNR(NodeRegistry_KOKKOS * nrk, SetOfEntities& nodesMappedTo, unsigned& noKeysNotInCommon) - { //checks if all the key value pairs of the NodeRegistry being called on match key value pairs in the NodeRegistry_KOKKOS being passed as an argument - //madbrew: not sure if this is really a fair comparison since the node id's in NodeIdsOnSubDimEntityType can differ due to the fact that map loops differ between the registry's when adding id's - //this causes the registry's to have the same number of id's and even the same sets of id's but a different mapping between parent (i.e. edge/face) and children node(s) - nodesMappedTo.clear(); - noKeysNotInCommon = 0; - if(!nrk) - { - std::cout << "NodeRegistry::compare_to_kokkos_NR : Invalid NodeRegistry_KOKKOS pointer, aborting NR diff\n"; - return false; - } - - bool kUOmap_contains_bUOmap = true; - - SubDimCellToDataMap * map = &m_cell_2_data_map; - for (SubDimCellToDataMap::iterator iter = map->begin(); iter != map->end(); ++iter) - { - const SubDimCell_SDCEntityType& subDimEntity_BOOST = (*iter).first; - SubDimCellData& data_BOOST = (*iter).second; - - SubDimCellData * data_KOKKOS = nrk->getFromMapPtr(subDimEntity_BOOST); - - if(!data_KOKKOS){ - kUOmap_contains_bUOmap = false; - if(kUOmap_contains_bUOmap) - std::cout << "NodeRegistry::compare_to_kokkos_NR : key with m_HashCode = " << subDimEntity_BOOST.getHash() << " NOT FOUND kokkos map\n"; - noKeysNotInCommon++; - } - else - { - if ( !( std::get<0>(*data_KOKKOS) == std::get<0>(data_BOOST) )) - { -// kUOmap_contains_bUOmap = false; - std::cout << "NodeRegistry::compare_to_kokkos_NR : NodeIdsOnSubDimEntityType on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap NodeIdsOnSubDimEntityType = "<< std::get<0>(*data_KOKKOS) << " bOUmap NodeIdsOnSubDimEntityType = " << std::get<0>(data_BOOST) << "\n"; - } - if(!( std::get<1>(*data_KOKKOS) == std::get<1>(data_BOOST) )) - { -// kUOmap_contains_bUOmap = false; - std::cout << "NodeRegistry::compare_to_kokkos_NR : EntityKey on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap EntityKey = "<< std::get<1>(*data_KOKKOS) << " bOUmap EntityKey = " << std::get<1>(data_BOOST) << "\n"; - } - if(!( std::get<2>(*data_KOKKOS) == std::get<2>(data_BOOST) ) ) - { -// kUOmap_contains_bUOmap = false; - std::cout << "NodeRegistry::compare_to_kokkos_NR : unsigned char on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap unsigned char = "<< std::get<2>(*data_KOKKOS) << " bOUmap unsigned char = " << std::get<2>(data_BOOST) << "\n"; - } - if(!( std::get<3>(*data_KOKKOS) == std::get<3>(data_BOOST) ) ) - { -// kUOmap_contains_bUOmap = false; - std::cout << "NodeRegistry::compare_to_kokkos_NR : unsigned char on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap unsigned char = "<< std::get<3>(*data_KOKKOS) << " bOUmap unsigned char = " << std::get<3>(data_BOOST) << "\n"; - } - if(!( std::get<4>(*data_KOKKOS) == std::get<4>(data_BOOST) ) ) - { -// kUOmap_contains_bUOmap = false; - std::cout << "NodeRegistry::compare_to_kokkos_NR : Double2 on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap Double2 = "<< std::get<4>(*data_KOKKOS) << " bOUmap Double2 = " << std::get<4>(data_BOOST) << "\n"; - } - } - if(!kUOmap_contains_bUOmap) - std::cout << "NodeRegistry::compare_to_kokkos_NR : key with m_HashCode = " << subDimEntity_BOOST.getHash() << " doesn't match anything in kokkos map\n\n"; - for(unsigned iEnt=0;iEnt(data_BOOST).m_entity_id_vector.size();iEnt++){ - unsigned ent_val = std::get<0>(data_BOOST).m_entity_id_vector[iEnt]; - nodesMappedTo.insert(stk::mesh::Entity(ent_val)); - } - } - return noKeysNotInCommon==0; - } }//percept diff --git a/packages/percept/src/adapt/NodeRegistry.hpp b/packages/percept/src/adapt/NodeRegistry.hpp index 0a3df61ddc28..ff9ccc217909 100644 --- a/packages/percept/src/adapt/NodeRegistry.hpp +++ b/packages/percept/src/adapt/NodeRegistry.hpp @@ -56,15 +56,12 @@ #define DEBUG_NR_DEEP 0 #include -#include // use old PerceptMesh/BulkData create entities if set to 1 - if 0, use PerceptMesh ID server which is much faster (doesn't use DistributedIndex) #define USE_CREATE_ENTITIES 0 namespace percept { - class NodeRegistry_KOKKOS; - using std::vector; using std::map; using std::set; @@ -350,8 +347,6 @@ void mod_begin(); void mod_end(const std::string& msg=""); - bool verifyAllKeysInKokkosNR(NodeRegistry_KOKKOS * nrk, SetOfEntities& nodesMappedTo, unsigned& noKeysNotInCommon); - private: percept::PerceptMesh& m_eMesh; Refiner *m_refiner; diff --git a/packages/percept/src/adapt/NodeRegistry_KOKKOS.cpp b/packages/percept/src/adapt/NodeRegistry_KOKKOS.cpp deleted file mode 100644 index 51c0ea41eed2..000000000000 --- a/packages/percept/src/adapt/NodeRegistry_KOKKOS.cpp +++ /dev/null @@ -1,914 +0,0 @@ -// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering -// Solutions of Sandia, LLC (NTESS). Under the terms of Contract -// DE-NA0003525 with NTESS, the U.S. Government retains certain rights -// in this software. -// -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - namespace percept { - - void NodeRegistry_KOKKOS::initialize() - { - m_cell_2_data_map.clear(); - } - - void NodeRegistry_KOKKOS:: - beginRegistration(int ireg, int nreg) - { -// if (CHECK_DB) checkDB("beginRegistration"); - - m_nodes_to_ghost.resize(0); - m_pseudo_entities.clear(); - m_state = NRS_START_REGISTER_NODE; - if (m_debug) - std::cout << "P[" << m_eMesh.get_rank() << "] tmp NodeRegistry_KOKKOS::beginRegistration" << std::endl; - } - - void NodeRegistry_KOKKOS::resizeMap() - { //attempts to resize the map to fit only it's current key/val pairs. The kokkos rehash doesn't allow for downsizing the map which is problematic for our use case of overestimating the number of nodes needed to avoid resizing the map - //in the node registration process - //madbrew : Note that this calls the size() function which is NOT a device function and cannot be called from a parallel kernel. Further, the Kokkos documentation says that this function has undefined behavior when erasable() is true - //It appears to be able to reliably tell me how many entries are in the map in a unit test, so I'm guessing the undefined behavior arises when the size function is called simultaneously with an insert/removal - unsigned current_cap = m_cell_2_data_map.capacity(); - unsigned current_size = m_cell_2_data_map.size(); - if(current_size>current_cap) - throw std::runtime_error("Something went horribly wrong, map size is greater than capacity"); - if( ( (current_cap - current_size) < m_waste_tolerance ) || current_size == 128) //since the map always sizes itself to the closest multiple of 128, 128 is the minimum size. No need to resize - return; //if the map is already full or an acceptable amount of space is being wasted, then do not resize - - m_cell_2_data_map.rehash(current_size); - } - - void NodeRegistry_KOKKOS:: - endRegistration(int ireg, int nreg) - { - if (m_debug) - std::cout << "P[" << m_eMesh.get_rank() << "] tmp NodeRegistry_KOKKOS::endRegistration start" << std::endl; - - - communicate_marks(); - - removeUnmarkedSubDimEntities(); - - std::cout << "map's size = " << m_cell_2_data_map.size() << std::endl; - resizeMap(); - std::cout << "map's capacity = " << m_cell_2_data_map.capacity() << std::endl; - - mod_begin(); - this->createNewNodesInParallel(); - - { - m_nodes_to_ghost.resize(0); - - if (m_debug) - std::cout << "P[" << m_eMesh.get_rank() << "] tmp NodeRegistry_KOKKOS::endRegistration end" << std::endl; - -// if (CHECK_DB) checkDB("endRegistration - after rehash"); - } - m_state = NRS_END_REGISTER_NODE; - } - - /// when a sub-dim entity is visited during node registration but is flagged as not being marked, and thus not requiring - /// any new nodes, we flag it with NR_MARK_NONE, then remove it here - void NodeRegistry_KOKKOS::removeUnmarkedSubDimEntities() - { - EXCEPTWATCH; -// SubDimCellToDataMap_KOKKOS::iterator iter; - bool debug = false; - -// for (iter = m_cell_2_data_map.begin(); iter != m_cell_2_data_map.end(); ++iter) -// { -// SubDimCellData& nodeId_elementOwnerId = (*iter).second; - for(unsigned iMapIndx = 0; iMapIndx(nodeId_elementOwnerId); - if (debug) std::cout << "tmp SRK nodeIds_onSE.size= " << nodeIds_onSE.size() << std::endl; - if (nodeIds_onSE.size()) - { - unsigned mark = nodeIds_onSE.m_mark; - unsigned is_marked = mark & NR_MARK; - unsigned is_not_marked = mark & NR_MARK_NONE; - - if (debug) - std::cout << "tmp SRK is_marked= " << is_marked << " is_not_marked= " << is_not_marked << std::endl; - if (!is_marked && is_not_marked) - { - // check if the node is a "hanging node" in which case it has relations - bool found = false; - for (unsigned in=0; in < nodeIds_onSE.size(); ++in) - { - if(!m_eMesh.is_valid(nodeIds_onSE[in])) continue; - if (debug) std::cout << "is valid = " << m_eMesh.is_valid(nodeIds_onSE[in]) << std::endl; - size_t num_rels = m_eMesh.get_bulk_data()->count_relations(nodeIds_onSE[in]); - if (num_rels) - { - if (debug) std::cout << "tmp SRK num_rels is non-zero in removeUnmarkedSubDimEntities, id= " << m_eMesh.identifier(nodeIds_onSE[in]) << std::endl; - found = true; - break; - } - } - - if (debug) - std::cout << "P[" << m_eMesh.get_rank() << " removeUnmarkedSubDimEntities:: tmp SRK for node= " << m_eMesh.identifier(nodeIds_onSE[0]) << " FOUND mark = " << mark - << " NR_MARK =" << NR_MARK << " NR_MARK_NONE= " << NR_MARK_NONE - << " resize to 0 (delete unmarked entity) = " << (!found) - << std::endl; - if (!found) - { - nodeIds_onSE.resize(0); - } - } - } - } - } - - /// Register the need for a new node on the sub-dimensional entity @param subDimEntity on element @param element. - /// If the element is a ghost element, the entity is still registered: the locality/ownership of the new entity - /// can be determined by the locality of the element (ghost or not). - bool NodeRegistry_KOKKOS::registerNeedNewNode(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes,const CellTopologyData * const bucket_topo_data) - { - bool ret_val = false; - SubDimCell_SDCEntityType subDimEntity(&m_eMesh); - bool foundGhostNode = getSubDimEntity(subDimEntity, element, needed_entity_rank.first, iSubDimOrd, bucket_topo_data); - if (foundGhostNode) - return ret_val; - - if (s_use_new_ownership_check && m_eMesh.isGhostElement(element)) - return false; - - static SubDimCellData empty_SubDimCellData; - - SubDimCellData* nodeId_elementOwnerId_ptr = getFromMapPtr(subDimEntity); - SubDimCellData& nodeId_elementOwnerId = (nodeId_elementOwnerId_ptr ? *nodeId_elementOwnerId_ptr : empty_SubDimCellData); - bool is_empty = nodeId_elementOwnerId_ptr == 0; - bool is_not_empty_but_data_cleared = (!is_empty && std::get(nodeId_elementOwnerId).size() == 0); - - // if empty or if my id is the smallest, make this element the owner - stk::mesh::EntityId db_id = std::get(nodeId_elementOwnerId).id(); - stk::mesh::EntityRank db_rank = std::get(nodeId_elementOwnerId).rank(); - - stk::mesh::EntityId element_id = m_eMesh.identifier(element); - stk::mesh::EntityId element_rank = m_eMesh.entity_rank(element); - bool should_put_in_id = (element_id < db_id); - bool should_put_in_rank_gt = (element_rank > db_rank); - bool should_put_in_rank_gte = (element_rank >= db_rank); - bool should_put_in = should_put_in_rank_gt || (should_put_in_id && should_put_in_rank_gte); - - unsigned smark=0; - if (!is_empty) - { - unsigned& mark = std::get(nodeId_elementOwnerId).m_mark; - if (needNodes) - mark |= NR_MARK; - else - mark |= NR_MARK_NONE; - smark = mark; - } - - /// once it's in, the assertion should be: - /// owning_elementId < non_owning_elementId && owning_elementRank >= non_owning_elementRank - /// - if (is_empty || is_not_empty_but_data_cleared || should_put_in) - { - // create SubDimCellData for the map rhs - // add one to iSubDimOrd for error checks later - VERIFY_OP_ON(m_eMesh.identifier(element), !=, 0, "hmmm registerNeedNewNode #1"); - VERIFY_OP_ON(m_eMesh.is_valid(element), ==, true, "hmmm registerNeedNewNode #2"); - if (is_empty || is_not_empty_but_data_cleared) - { - unsigned numNewNodes = needed_entity_rank.second; - - SubDimCellData data = std::forward_as_tuple( NodeIdsOnSubDimEntityType(numNewNodes, stk::mesh::Entity(), smark), stk::mesh::EntityKey(m_eMesh.entity_rank(element), m_eMesh.identifier(element)), (unsigned char)needed_entity_rank.first, (unsigned char)(iSubDimOrd+1), Double2() ); - NodeIdsOnSubDimEntityType& nid_new = std::get(data); - if (needNodes) - nid_new.m_mark |= NR_MARK; - else - nid_new.m_mark |= NR_MARK_NONE; - - smark = nid_new.m_mark; - putInMap(subDimEntity, data); - } - else - { - stk::mesh::Entity owning_element = m_eMesh.get_entity(db_rank, db_id); - VERIFY_OP_ON(m_eMesh.is_valid(owning_element), ==, true, "hmmm"); - NodeIdsOnSubDimEntityType& nid = std::get(nodeId_elementOwnerId); - SubDimCellData data = std::forward_as_tuple(nid, stk::mesh::EntityKey(m_eMesh.entity_rank(element), m_eMesh.identifier(element)), (unsigned char)needed_entity_rank.first, (unsigned char)(iSubDimOrd+1), Double2() ); - putInMap(subDimEntity,data); - } - ret_val = true; - } - - // all debug code from here on - - bool debug = false; - if (debug && !is_empty) - { - unsigned originalMark = std::get(nodeId_elementOwnerId).m_mark; - std::cout << m_eMesh.rank() << " registerNeedNewNode element= " << m_eMesh.print_entity_compact(element) << "\n needNodes= " << needNodes << " iSubDimOrd= " << iSubDimOrd - << " is_empty= " << is_empty << " is_not_empty_but_data_cleared= " << is_not_empty_but_data_cleared << " should_put_in= " << should_put_in - << " originalMark= " << originalMark << " currentMark= " << smark << " NR_MARK= " << NR_MARK - << m_eMesh.demangled_stacktrace(20) - << std::endl; - } - if (debug) - { - SubDimCellData* a_nodeId_elementOwnerId_ptr = getFromMapPtr(subDimEntity); - SubDimCellData& a_nodeId_elementOwnerId = (a_nodeId_elementOwnerId_ptr ? *a_nodeId_elementOwnerId_ptr : empty_SubDimCellData); - bool a_is_empty = a_nodeId_elementOwnerId_ptr == 0; - unsigned& gotMark = std::get(a_nodeId_elementOwnerId).m_mark; - unsigned nidSize = std::get(a_nodeId_elementOwnerId).size(); - - std::ostringstream sout; - sout << "P[" << m_eMesh.get_rank() << "] registerNeedNewNode:: element= " << m_eMesh.identifier(element) << " nidSize= " << nidSize - << " nid= " << (nidSize ? (int)m_eMesh.identifier(std::get(nodeId_elementOwnerId)[0]) : -1); - - sout << " smark= " << smark << " gotMark= " << gotMark << " needNodes= " << needNodes << " isG= " << m_eMesh.isGhostElement(element) - << " is_empty= " << is_empty - << " a_is_empty= " << a_is_empty - << " should_put_in= " << should_put_in - << " should_put_in_id= " << should_put_in_id - << " should_put_in_rank_gt= " << should_put_in_rank_gt - << " should_put_in_rank_gte= " << should_put_in_rank_gte - << " needed_entity_rank= " - << needed_entity_rank.first << " subDimEntity= "; - - for (unsigned k=0; k < subDimEntity.size(); k++) - { - sout << " " << m_eMesh.identifier(subDimEntity[k]) << " "; - } - std::cout << sout.str() << std::endl; - } - - //deletethis -/* SubDimCellData * tester = getFromMapPtr(subDimEntity); - if(!tester) - std::cout << " Cannot access the entity just registered in the map\n\n"; - else - std::cout << " Can access the entity just registered in the map\n\n";*/ - //deletethis - return ret_val; - } - - - /** - * For any given subDimEntity (edge, tri or quad face), we can define an "owner" for it as the processor - * that owns the node with the minimum ID. This is implemented below through use of the ordered std::map. - * However, a proc can "own" a subDimEntity, but not have any elements use that subDimEntity. - * - * Thus, we use the more global rule that the for all subDimEntity's that are used by an element, the - * proc with minimum rank is the one that actually "owns" the subDimEntity - this is implemented in the - * routines that unpack the node id's that get sent to all sharing procs. - * - * So, the main function of this routine is to tell the caller all the sharing procs of the subDimEntity - * and thus all procs that have an element that use this subDimEntity will send their info to all other - * sharing procs, but only the one with the min proc rank will not unpack and override its node ID. - * Other procs will unpack and override if the proc rank is smaller than their rank. - * - * Note: for consistency, and for the unpack to work correctly and assign ownership of the data, we - * return the current proc in the list of sharing procs @param other_procs - * - */ - - int NodeRegistry_KOKKOS::proc_owning_subdim_entity(const SubDimCell_SDCEntityType& subDimEntity, std::vector& other_procs, bool& all_shared) - { - other_procs.resize(0); - all_shared = false; - - VERIFY_OP_ON(subDimEntity.size(), >=, 1, "bad size"); - - int my_proc = m_eMesh.get_rank(); - - bool debug = false; - - std::vector sharing_procs; - std::map count_per_proc; - all_shared = true; - for (unsigned i = 0; i < subDimEntity.size(); ++i) - { - if (m_eMesh.shared(subDimEntity[i])) - { - m_eMesh.get_bulk_data()->comm_shared_procs(m_eMesh.entity_key(subDimEntity[i]), sharing_procs); - - sharing_procs.push_back(my_proc); - - for (unsigned j=0; j < sharing_procs.size(); ++j) - { - ++count_per_proc[sharing_procs[j]]; - } - } - else - { - all_shared = false; - } - } - if (debug) - { - std::cerr << m_eMesh.rank() << "tmp srk all_shared= " << all_shared - << "\n n0= " << m_eMesh.print_entity_compact(subDimEntity[0]) - << "\n n1= " << m_eMesh.print_entity_compact(subDimEntity[1]) - << std::endl; - } - if (!all_shared) - return -1; - - VERIFY_OP_ON(count_per_proc.size(), >=, 1, "bad size"); - int proc_owner = -1; - for (auto& counts : count_per_proc) - { - if (counts.second == subDimEntity.size()) - { - if (proc_owner < 0) - proc_owner = counts.first; - other_procs.push_back(counts.first); - } - } - if (debug) { - std::cerr << m_eMesh.rank() << "tmp srk proc_owner= " << proc_owner << std::endl; - } - VERIFY_OP_ON(proc_owner, >=, 0, "bad proc_owner"); - return proc_owner; - } - - - /// check the newly registered node from the registry, which does one of three things, depending on what mode we are in: - /// 1. counts buffer in prep for sending (just does a pack) - /// 2. packs the buffer (after buffers are alloc'd) - /// 3. returns the new node after all communications are done - - - void NodeRegistry_KOKKOS:: - doForAllSubEntities(ElementFunctionPrototype function, const stk::mesh::Entity element, vector& needed_entity_ranks,const CellTopologyData * const bucket_topo_data) - { - const CellTopologyData * const cell_topo_data = (bucket_topo_data ? bucket_topo_data : m_eMesh.get_cell_topology(element)); - - shards::CellTopology cell_topo(cell_topo_data); - - for (unsigned ineed_ent=0; ineed_ent < needed_entity_ranks.size(); ineed_ent++) - { - unsigned numSubDimNeededEntities = 0; - stk::mesh::EntityRank needed_entity_rank = needed_entity_ranks[ineed_ent].first; - - if (needed_entity_rank == m_eMesh.edge_rank()) - { - numSubDimNeededEntities = cell_topo_data->edge_count; - } - else if (needed_entity_rank == m_eMesh.face_rank()) - { - numSubDimNeededEntities = cell_topo_data->side_count; - } - else if (needed_entity_rank == stk::topology::ELEMENT_RANK) - { - numSubDimNeededEntities = 1; - } - - curr_funct = function; - doForAllSubEntitiesFunctor for_each(this,element,needed_entity_ranks,bucket_topo_data,numSubDimNeededEntities,ineed_ent); - for_each.run(); - -/* - for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++) - { - (this ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, true, bucket_topo_data); - - } // iSubDimOrd -*/ - } // ineed_ent - } - - unsigned NodeRegistry_KOKKOS::local_size() - { - unsigned sz=0; -// for (SubDimCellToDataMap_KOKKOS::iterator cell_iter = m_cell_2_data_map.begin(); cell_iter != m_cell_2_data_map.end(); ++cell_iter) -// { -// SubDimCellData& data = (*cell_iter).second; - for(unsigned iMapIndx = 0; iMapIndx(nodeId_elementOwnerId).id(); - - NodeIdsOnSubDimEntityType& nodeIds_onSE = std::get(nodeId_elementOwnerId); - if (nodeIds_onSE.size()) - { - stk::mesh::EntityRank erank = std::get(nodeId_elementOwnerId).rank(); - stk::mesh::Entity owning_element = get_entity_element(*m_eMesh.get_bulk_data(), erank, owning_elementId); - - if (!m_eMesh.is_valid(owning_element)) - { - if (!s_use_new_ownership_check) - { - std::cout << "tmp owning_element = null, owning_elementId= " << owning_elementId - << std::endl; - throw std::logic_error("logic: hmmm #5.2"); - } - else - continue; - } - if (!m_eMesh.isGhostElement(owning_element)) - { - sz += nodeIds_onSE.size(); - } - } - } - return sz; - } - - /// after registering all needed nodes, this method is used to request new nodes on this processor - void NodeRegistry_KOKKOS::createNewNodesInParallel() - { - stk::mesh::Part* new_nodes_part = m_eMesh.get_non_const_part("refine_new_nodes_part"); - VERIFY_OP_ON(new_nodes_part, !=, 0, "new_nodes_part"); - unsigned num_nodes_needed = local_size(); - - // assert( bulk data is in modifiable mode) - // create new entities on this proc - vector new_nodes; - - if (m_useAddNodeSharing) - { - stk::diag::Timer *timerCE_ = 0; - if (m_refiner) - { - timerCE_ = new stk::diag::Timer("NR_CreateEnt", m_refiner->rootTimer()); - timerCE_->start(); - } - -#if USE_CREATE_ENTITIES - m_eMesh.createEntities( stk::topology::NODE_RANK, num_nodes_needed, new_nodes); -#else - m_eMesh.initializeIdServer(); - stk::mesh::Part& nodePart = m_eMesh.get_fem_meta_data()->get_topology_root_part(stk::topology::NODE); - stk::mesh::PartVector nodeParts(1, &nodePart); - m_eMesh.getEntitiesUsingIdServer( m_eMesh.node_rank(), num_nodes_needed, new_nodes, nodeParts); -#endif - if (timerCE_) - { - timerCE_->stop(); - delete timerCE_; - } - } - std::vector ids(num_nodes_needed); - - if (new_nodes_part) - { - stk::mesh::Selector selector(m_eMesh.get_fem_meta_data()->locally_owned_part() ); - std::vector add_parts(1, new_nodes_part); - std::vector remove_parts; - for (unsigned ind = 0; ind < new_nodes.size(); ind++) - { - if (m_eMesh.m_new_nodes_field) - { - NewNodesType_type *ndata = stk::mesh::field_data(*m_eMesh.m_new_nodes_field, new_nodes[ind]); - if (ndata) - { - ndata[0] = static_cast(1); - } - } - m_eMesh.get_bulk_data()->change_entity_parts( new_nodes[ind], add_parts, remove_parts ); - } - } - // set map values to new node id's -// unsigned inode=0; - - Kokkos::ViewthreadSafeiNode("threadSafeiNode",1); - Kokkos::View::HostMirror threadSafeiNode_mir("threadSafeiNode_mir",1); - threadSafeiNode_mir(0) = 0; - Kokkos::deep_copy(threadSafeiNode,threadSafeiNode_mir); - - // Malachi Phillips Jun 21, 2018 - // - // This needs a capture by *this in order to work on __device__ - // - // Capture by *this cannot work on __host__, unless using a complaint c++17 compiler - // - Kokkos::parallel_for(Kokkos::RangePolicy(0,m_cell_2_data_map.capacity()), [&] (const unsigned int iMapIndx) - { - const SubDimCell_SDCEntityType& subDimEntity = m_cell_2_data_map.key_at(iMapIndx); - unsigned iData = m_cell_2_data_map.find(subDimEntity); - if(iData == Kokkos::UnorderedMapInvalidIndex) - //continue; //continue doesn't work the same way in pf - (void)iData; - else{ - SubDimCellData& data = m_cell_2_data_map.value_at(iData); - NodeIdsOnSubDimEntityType& nodeIds_onSE = std::get(data); - if (!nodeIds_onSE.size()) - // continue; //continue doesn't work the same way in pf - (void)iData; - else{ - stk::mesh::EntityId owning_elementId = std::get(data).id(); - -#ifndef __CUDACC__ - if (!owning_elementId) - { - throw std::logic_error("logic: hmmm #5.4.0"); - } -#endif - - stk::mesh::EntityRank erank = std::get(data).rank(); - stk::mesh::Entity owning_element = get_entity_element(*m_eMesh.get_bulk_data(), erank, owning_elementId); - - if (!m_eMesh.is_valid(owning_element)) - { -#ifndef __CUDACC__ - if (!s_use_new_ownership_check) - throw std::logic_error("logic: hmmm #5.4"); -#endif - } - - else if (s_use_new_ownership_check && m_eMesh.isGhostElement(owning_element)) - { - //std::cerr << m_eMesh.rank() << " owning_element= " << m_eMesh.print_entity_compact(owning_element) << std::endl; - // FIXME - //VERIFY_MSG("found ghost"); - // continue; //continue doesn't work the same way in pf - } - - else if (!m_eMesh.isGhostElement(owning_element)) - { - if (nodeIds_onSE.m_entity_id_vector.size() != nodeIds_onSE.size()) - { -#ifndef __CUDACC__ - throw std::logic_error("NodeRegistry_KOKKOS:: createNewNodesInParallel logic err #0.0"); -#endif - } - - for (unsigned ii = 0; ii < nodeIds_onSE.size(); ii++) - { - unsigned localiNode = Kokkos::atomic_fetch_add(&threadSafeiNode(0),1); //fetches old value then increments the index in a thread safe manner - - VERIFY_OP(localiNode, < , num_nodes_needed, "UniformRefiner::doBreak() too many nodes"); -#ifndef __CUDACC__ - if ( DEBUG_NR_UNREF) - { - std::cout << "tmp createNewNodesInParallel: old node id= " << (m_eMesh.is_valid(nodeIds_onSE[ii]) ? toString(m_eMesh.identifier(nodeIds_onSE[ii])) : std::string("null")) << std::endl; - std::cout << "tmp createNewNodesInParallel: new node="; - m_eMesh.print_entity(std::cout, new_nodes[localiNode]); - } -#endif - - // if already exists from a previous iteration/call to doBreak, don't reset it and just use the old node - if (m_eMesh.is_valid(nodeIds_onSE[ii])) - { -#ifndef __CUDACC__ - if (DEBUG_NR_UNREF) - { - std::cout << "tmp createNewNodesInParallel: old node id is no-null, re-using it= " << (m_eMesh.is_valid(nodeIds_onSE[ii]) ? toString(m_eMesh.identifier(nodeIds_onSE[ii])) : std::string("null")) << std::endl; - std::cout << "tmp createNewNodesInParallel: new node="; - m_eMesh.print_entity(std::cout, new_nodes[localiNode]); - } -#endif - } - else - { - nodeIds_onSE[ii] = new_nodes[localiNode]; - nodeIds_onSE.m_entity_id_vector[ii] = m_eMesh.identifier(new_nodes[localiNode]); - } - - // inode++; - } - } - } - } - });//iMapIndx - } - - void NodeRegistry_KOKKOS::mod_begin() - { - if (m_refiner) - { - m_refiner->mod_begin(); - } - else - { - m_eMesh.get_bulk_data()->modification_begin(); - } - } - - void NodeRegistry_KOKKOS::mod_end(const std::string& msg) - { - if (m_refiner) - { - m_refiner->mod_end(0,"NReg:"+msg); - } - else - { - m_eMesh.get_bulk_data()->modification_end(); - } - } - - void NodeRegistry_KOKKOS::communicate_marks() - { - - // only one stage now - each proc sends to other sharing procs and each can |= and accumulate locally - - for (int stage = 0; stage < 1; ++stage) - { - stk::CommSparse commAll (m_eMesh.parallel()); - - communicate_marks_pack(commAll, stage); - - commAll.allocate_buffers(); - - communicate_marks_pack(commAll, stage); - commAll.communicate(); - communicate_marks_unpack(commAll); - } - } - - void NodeRegistry_KOKKOS::communicate_marks_pack(stk::CommSparse& commAll, int stage) - { - CommDataType buffer_entry; - - SubDimCellToDataMap_KOKKOS& map = m_cell_2_data_map; - -// for (SubDimCellToDataMap_KOKKOS::iterator iter = map.begin(); iter != map.end(); ++iter) -// { -// const SubDimCell_SDCEntityType& subDimEntity = (*iter).first; -// SubDimCellData& nodeId_elementOwnerId = (*iter).second; - for(unsigned iMapIndx = 0; iMapIndx(std::get(nodeId_elementOwnerId)); - - static std::vector procs_to_send_to; - - bool all_shared = false; - int new_owning_proc = proc_owning_subdim_entity(subDimEntity, procs_to_send_to, all_shared); - (void)new_owning_proc; - - bool need_to_send = false; - if (stage == 0) - { - need_to_send = all_shared && (procs_to_send_to.size() != 0); - } - - if (!need_to_send) - continue; - - NodeIdsOnSubDimEntityType& nodeIds_onSE = std::get(nodeId_elementOwnerId); - - if (nodeIds_onSE.size()) - { - unsigned mark = nodeIds_onSE.m_mark; - unsigned is_marked = mark & NR_MARK; - unsigned is_not_marked = mark & NR_MARK_NONE; - // only need to send if it's got non-zero mark info - if (is_marked || is_not_marked) - { - std::get(buffer_entry) = subDimEntity.size(); - std::get(buffer_entry) = owning_element_subDim_rank; - - for (unsigned inode=0; inode < subDimEntity.size(); ++inode) - std::get(buffer_entry)[inode] = m_eMesh.entity_key(subDimEntity[inode]); - - for (unsigned jprocs = 0; jprocs < procs_to_send_to.size(); ++jprocs) - { - if (procs_to_send_to[jprocs] != m_eMesh.get_rank()) - { - commAll.send_buffer( procs_to_send_to[jprocs] ).pack< CommDataType > (buffer_entry); - commAll.send_buffer( procs_to_send_to[jprocs] ).pack(mark); - } - } - } - } - } - } - - void NodeRegistry_KOKKOS::communicate_marks_unpack(stk::CommSparse& commAll) - { - unsigned proc_size = m_eMesh.get_parallel_size(); - - CommDataType buffer_entry; - - //try - { - for(unsigned from_proc = 0; from_proc < proc_size; ++from_proc ) - { - stk::CommBuffer & recv_buffer = commAll.recv_buffer( from_proc ); - - while ( recv_buffer.remaining() ) - { - recv_buffer.unpack< CommDataType >( buffer_entry ); - unsigned mark=0; - recv_buffer.unpack< unsigned > (mark); - - { - unsigned subDimEntitySize = std::get(buffer_entry); - - SubDimCell_SDCEntityType subDimEntity(&m_eMesh); - //getSubDimEntity(subDimEntity, owning_element, needed_entity_rank, iSubDimOrd); - for (unsigned inode = 0; inode < subDimEntitySize; ++inode) - { - stk::mesh::Entity node = m_eMesh.get_entity(std::get(buffer_entry)[inode]); - VERIFY_OP_ON(m_eMesh.is_valid(node), ==, true, "bad node"); - subDimEntity.insert(node); - } - - static SubDimCellData empty_SubDimCellData; - - SubDimCellData* nodeId_elementOwnerId_ptr = getFromMapPtr(subDimEntity); - SubDimCellData& nodeId_elementOwnerId = (nodeId_elementOwnerId_ptr ? *nodeId_elementOwnerId_ptr : empty_SubDimCellData); - bool is_empty = nodeId_elementOwnerId_ptr == 0; - - //VERIFY_OP_ON(is_empty, !=, true, "hmmm"); - if (!is_empty) - { - NodeIdsOnSubDimEntityType& nodeIds_onSE = std::get(nodeId_elementOwnerId); - - // accumulation from all other procs - nodeIds_onSE.m_mark |= mark; - } - } - } - } - } - } - - bool NodeRegistry_KOKKOS::verifyAllKeysInBoostNR(NodeRegistry * nr, SetOfEntities& nodesMappedTo, unsigned& noKeysNotInCommon) - { //checks to if all the nodes within the kokkos NR have the same mapping and nodes within the boost NR - //madbrew: not sure if this is really a fair comparison since the node id's in NodeIdsOnSubDimEntityType can differ due to the fact that map loops differ between the registry's when adding id's - //this causes the registry's to have the same number of id's and even the same sets of id's but a different mapping between parent (i.e. edge/face) and children node(s) - nodesMappedTo.clear(); - noKeysNotInCommon = 0; - if(!nr) - { - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : Invalid NodeRegistry_KOKKOS pointer, aborting NR diff\n"; - return false; - } - - bool bUOmap_contains_kUOmap = true; - - SubDimCellToDataMap_KOKKOS * map = &m_cell_2_data_map; - for (unsigned iMapIndx = 0; iMapIndx < map->capacity(); iMapIndx++) - { - const SubDimCell_SDCEntityType& subDimEntity_KOKKOS = map->key_at(iMapIndx); - unsigned iVal = map->find(subDimEntity_KOKKOS); - if(iVal==Kokkos::UnorderedMapInvalidIndex) - continue; - SubDimCellData& data_KOKKOS = map->value_at(iVal); - - SubDimCellData * data_BOOST = nr->getFromMapPtr(subDimEntity_KOKKOS); - - if(!data_BOOST){ - bUOmap_contains_kUOmap = false; - if(bUOmap_contains_kUOmap) - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : key with m_HashCode = " << subDimEntity_KOKKOS.getHash() << " NOT FOUND kokkos map\n"; - noKeysNotInCommon++; - } - else - { - if ( !( std::get<0>(data_KOKKOS) == std::get<0>(*data_BOOST) )) - { -// bUOmap_contains_kUOmap = false; - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : NodeIdsOnSubDimEntityType on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap NodeIdsOnSubDimEntityType = "<< std::get<0>(data_KOKKOS) << " bOUmap NodeIdsOnSubDimEntityType = " << std::get<0>(*data_BOOST) << "\n"; - } - if(!( std::get<1>(data_KOKKOS) == std::get<1>(*data_BOOST) )) - { -// bUOmap_contains_kUOmap = false; - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : EntityKey on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap EntityKey = "<< std::get<1>(data_KOKKOS) << " bOUmap EntityKey = " << std::get<1>(*data_BOOST) << "\n"; - } - if(!( std::get<2>(data_KOKKOS) == std::get<2>(*data_BOOST) ) ) - { -// bUOmap_contains_kUOmap = false; - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : unsigned char on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap unsigned char = "<< std::get<2>(data_KOKKOS) << " bOUmap unsigned char = " << std::get<2>(*data_BOOST) << "\n"; - } - if(!( std::get<3>(data_KOKKOS) == std::get<3>(*data_BOOST) ) ) - { -// bUOmap_contains_kUOmap = false; - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : unsigned char on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap unsigned char = "<< std::get<3>(data_KOKKOS) << " bOUmap unsigned char = " << std::get<3>(*data_BOOST) << "\n"; - } - if(!( std::get<4>(data_KOKKOS) == std::get<4>(*data_BOOST) ) ) - { -// bUOmap_contains_kUOmap = false; - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : Double2 on tuple's do not match between NodeRegistry and KokkosNodeRegistry\n"; - std::cout << " kOUmap Double2 = "<< std::get<4>(data_KOKKOS) << " bOUmap Double2 = " << std::get<4>(*data_BOOST) << "\n"; - } - } - - if(!bUOmap_contains_kUOmap) - std::cout << "NodeRegistry_KOKKOS::compare_to_boost_NR : key with m_HashCode = " << subDimEntity_KOKKOS.getHash() << " doesn't match anything in kokkos map\n\n"; - for(unsigned iEnt=0;iEnt(data_KOKKOS).m_entity_id_vector.size();iEnt++){ - unsigned ent_val = std::get<0>(data_KOKKOS).m_entity_id_vector[iEnt]; - nodesMappedTo.insert(stk::mesh::Entity(ent_val)); - } - } - return noKeysNotInCommon==0; - } - - bool NodeRegistry_KOKKOS:: - getSubDimEntity(SubDimCell_SDCEntityType& subDimEntity, const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd, - const CellTopologyData * const bucket_topo_data - ) - { - subDimEntity.clear(); - // in the case of elements, we don't share any nodes so we just make a map of element id to node - if (needed_entity_rank == stk::topology::ELEMENT_RANK) - { - subDimEntity.resize(1); - subDimEntity[0] = element ; - subDimEntity.updateHashCode(); - return false; - } - - const CellTopologyData * const cell_topo_data = (bucket_topo_data ? bucket_topo_data : m_eMesh.get_cell_topology(element) ); - - stk::mesh::Entity const * const elem_nodes = m_eMesh.get_bulk_data()->begin_nodes(element); - - const unsigned * inodes = 0; - unsigned nSubDimNodes = 0; - static const unsigned edge_nodes_2[2] = {0,1}; - static const unsigned face_nodes_3[3] = {0,1,2}; - static const unsigned face_nodes_4[4] = {0,1,2,3}; - - // special case for faces in 3D - if (needed_entity_rank == m_eMesh.face_rank() && needed_entity_rank == m_eMesh.entity_rank(element)) - { - nSubDimNodes = cell_topo_data->vertex_count; - - // note, some cells have sides with both 3 and 4 nodes (pyramid, prism) - if (nSubDimNodes ==3 ) - inodes = face_nodes_3; - else - inodes = face_nodes_4; - - } - // special case for edges in 2D - else if (needed_entity_rank == m_eMesh.edge_rank() && needed_entity_rank == m_eMesh.entity_rank(element)) - { - nSubDimNodes = cell_topo_data->vertex_count; - - if (nSubDimNodes == 2 ) - { - inodes = edge_nodes_2; - } - else - { - throw std::runtime_error("NodeRegistry_KOKKOS bad for edges"); - } - } - else if (needed_entity_rank == m_eMesh.edge_rank()) - { - inodes = cell_topo_data->edge[iSubDimOrd].node; - nSubDimNodes = 2; - } - else if (needed_entity_rank == m_eMesh.face_rank()) - { - nSubDimNodes = cell_topo_data->side[iSubDimOrd].topology->vertex_count; - // note, some cells have sides with both 3 and 4 nodes (pyramid, prism) - inodes = cell_topo_data->side[iSubDimOrd].node; - } - - subDimEntity.resize(nSubDimNodes); - for (unsigned jnode = 0; jnode < nSubDimNodes; jnode++) - { - subDimEntity[jnode] = elem_nodes[inodes[jnode]] ; - } - bool foundGhostNode = false; - if (m_checkForGhostedNodes) - { - for (unsigned jnode = 0; jnode < nSubDimNodes; jnode++) - { - stk::mesh::Bucket& bucket = m_eMesh.bucket(subDimEntity[jnode]); - if (!bucket.owned() && !bucket.shared()) - foundGhostNode = true; - } - } - - subDimEntity.sort(); - subDimEntity.updateHashCode(); - return foundGhostNode; - } - - }//percept diff --git a/packages/percept/src/adapt/NodeRegistry_KOKKOS.hpp b/packages/percept/src/adapt/NodeRegistry_KOKKOS.hpp deleted file mode 100644 index 8f0d31ac44d7..000000000000 --- a/packages/percept/src/adapt/NodeRegistry_KOKKOS.hpp +++ /dev/null @@ -1,448 +0,0 @@ -// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering -// Solutions of Sandia, LLC (NTESS). Under the terms of Contract -// DE-NA0003525 with NTESS, the U.S. Government retains certain rights -// in this software. -// -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -/** - * NodeRegistry_KOKKOS: class that defines a map of key->value where key is a SubDimCell - * (e.g. an edge or face bewtween two elements), and value contains various data - * about the SubDimCell, such as which element owns the SubDimCell, an array of - * new nodes needed on it for refinement, etc. - */ - -#ifndef adapt_NodeRegistry_KOKKOS_hpp -#define adapt_NodeRegistry_KOKKOS_hpp - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include - -#include -#include -#include - -#include -#include -#include - -#include - -#include -#include - -#define DEBUG_PRINT_11 0 -#define NR_PRINT(a) do { if (DEBUG_PRINT_11) std::cout << #a << " = " << a ; } while(0) -#define NR_PRINT_OUT(a,out) do { if (DEBUG_PRINT_11) out << #a << " = " << a << std::endl; } while(0) - -#define STK_ADAPT_NODEREGISTRY_USE_ENTITY_REPO 0 -#define STK_ADAPT_NODEREGISTRY_DO_REHASH 1 - -#define DEBUG_NR_UNREF 0 -#define DEBUG_NR_DEEP 0 - -#include -#include -#include -#include -// use old PerceptMesh/BulkData create entities if set to 1 - if 0, use PerceptMesh ID server which is much faster (doesn't use DistributedIndex) -#define USE_CREATE_ENTITIES 0 - - namespace percept { - - class NodeRegistry; - - using std::vector; - using std::map; - using std::set; - - class Refiner; - - /// map of the node ids on a sub-dim entity to the data on the sub-dim entity - typedef Kokkos::UnorderedMap, my_fast_equal_to > SubDimCellToDataMap_KOKKOS; - - //======================================================================================================================== - //======================================================================================================================== - //======================================================================================================================== - class NodeRegistry_KOKKOS - { - public: - - // this switches new method of determining ownership of subDimEntitys - to be removed - static const bool s_use_new_ownership_check = true; - - static const unsigned NR_MARK_NONE = 1u << 0; - static const unsigned NR_MARK = 1u << 1; - - //======================================================================================================================== - // high-level interface - - NodeRegistry_KOKKOS(percept::PerceptMesh& eMesh, - Refiner *refiner=0, - bool useCustomGhosting = true,unsigned size_hint=100, unsigned waste_tol=0) : m_eMesh(eMesh), - m_refiner(refiner), - m_comm_all( new stk::CommSparse(eMesh.parallel()) ), - m_cell_2_data_map(size_hint), - m_pseudo_entities(*eMesh.get_bulk_data()), - m_useCustomGhosting(useCustomGhosting), - m_useAddNodeSharing(false), - m_checkForGhostedNodes(false), - m_gee_cnt(0), m_gen_cnt(0), - m_debug(false), - m_state(NRS_NONE), - m_waste_tolerance(waste_tol) - { - m_useCustomGhosting = false; - m_useAddNodeSharing = !m_useCustomGhosting; - } - - ~NodeRegistry_KOKKOS() { - if (m_comm_all) - delete m_comm_all; - } - - void init_comm_all(); - void clear_dangling_nodes(SetOfEntities* nodes_to_be_deleted); - void initialize(); - - // avoid adding subDimEntity if it has a ghosted node - void setCheckForGhostedNodes(bool val) { m_checkForGhostedNodes = val; } - bool getCheckForGhostedNodes() { return m_checkForGhostedNodes; } - -#define CHECK_DB 0 - - void beginRegistration(int ireg=0, int nreg=1); - void endRegistration(int ireg=0, int nreg=1); - - void beginCheckForRemote(int ireg=0, int nreg=1); - void endCheckForRemote(int ireg=0, int nreg=1); - - void beginGetFromRemote(int ireg=0, int nreg=1); - void endGetFromRemote(int ireg=0, int nreg=1); - - - /// Register the need for a new node on the sub-dimensional entity @param subDimEntity on element @param element. - /// If the element is a ghost element, the entity is still registered: the locality/ownership of the new entity - /// can be determined by the locality of the element (ghost or not). - bool registerNeedNewNode(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes,const CellTopologyData * const bucket_topo_data); - - /// check the newly registered node from the registry, which does one of three things, depending on what mode we are in: - /// 1. counts buffer in prep for sending (just does a pack) - /// 2. packs the buffer (after buffers are alloc'd) - /// 3. returns the new node after all communications are done - bool checkForRemote(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes_notUsed,const CellTopologyData * const bucket_topo_data); - bool getFromRemote(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes_notUsed,const CellTopologyData * const bucket_topo_data); - - - /// util - /// used to find which proc "owns" a subDimEntity - see more comments in implementation - int proc_owning_subdim_entity(const SubDimCell_SDCEntityType& subDimEntity, std::vector& other_procs, bool& all_shared); - - /// makes coordinates of this new node be the centroid of its sub entity - void prolongateCoords(const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd); - void prolongateCoordsAllSubDims(stk::mesh::Entity element); - void prolongateField(const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd, stk::mesh::FieldBase *field); - - /// do interpolation for all fields, for the given nodes - void prolongateFieldNodeVector(std::vector & nodes); - - double spacing_edge(std::vector& nodes, - unsigned iv0, unsigned iv1, unsigned nsz, unsigned nsp, double lspc[8][3], double den_xyz[3], double *coord[8]); - - void normalize_spacing(stk::mesh::Entity element, std::vector &nodes, - unsigned nsz, unsigned nsp, double spc[8][3], double den_xyz[3], double *coord[8]); - - void prolongate(stk::mesh::FieldBase *field, unsigned *subDimSize_in=0, bool useFindValidCentroid = true); - - /// do interpolation for all fields - void prolongateFields(const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd); - - /// do interpolation for all fields - void prolongateFields(); - - /// check for adding new nodes to existing parts based on sub-entity part ownership - void addToExistingParts(const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd); - - void add_rbars(std::vector >& rbar_types ); - - /// Check for adding new nodes to existing parts based on sub-entity part ownership. - /// This version does it in bulk and thus avoids repeats on shared sub-dim entities. - void addToExistingPartsNew(); - - inline SubDimCellData& getNewNodeAndOwningElement(SubDimCell_SDCEntityType& subDimEntity) - { -// return m_cell_2_data_map[subDimEntity]; - unsigned mapIndx = m_cell_2_data_map.find(subDimEntity); - return m_cell_2_data_map.value_at(mapIndx); - } - - inline SubDimCellData * getFromMapPtr(const SubDimCell_SDCEntityType& subDimEntity) const - { -// const SubDimCellToDataMap_KOKKOS::const_iterator i = m_cell_2_data_map.find( subDimEntity ); -// return i != m_cell_2_data_map.end() ? (const_cast(&(i->second))) : 0 ; - unsigned mapIndx = m_cell_2_data_map.find(subDimEntity); - if(mapIndx == Kokkos::UnorderedMapInvalidIndex){ -// std::cout << "Entity not found. Map's capacity = " << m_cell_2_data_map.capacity() << std::endl; - return 0; - } - return &m_cell_2_data_map.value_at(mapIndx); - } - - inline SubDimCellData& getFromMap(const SubDimCell_SDCEntityType& subDimEntity) const - { -// const SubDimCellToDataMap_KOKKOS::const_iterator i = m_cell_2_data_map.find( subDimEntity ); -// return * (const_cast(&(i->second))); - unsigned mapIndx = m_cell_2_data_map.find(subDimEntity); - return m_cell_2_data_map.value_at(mapIndx); - } - - // allow client code to insert into map (e.g. for nodes created on-the-fly - see quad transition code) - void forceInMap(std::vector& vecNodes, unsigned mark, stk::mesh::Entity element, stk::mesh::EntityRank needed_rank, unsigned iSubDimOrd ); - - inline void putInMap(SubDimCell_SDCEntityType& subDimEntity, SubDimCellData& data) - { - bool debug = false; - if (debug) - { - SubDimCellData& nodeId_elementOwnderId = data; - - NodeIdsOnSubDimEntityType& nodeIds_onSE = std::get(nodeId_elementOwnderId); - stk::mesh::EntityId owning_elementId = std::get(nodeId_elementOwnderId).id(); - - if (1) - std::cout << "put in map: nodeIds_onSE.size= " << (nodeIds_onSE.size()) - << " owning_elementId= " << owning_elementId - << " subDimEntity.size= " << subDimEntity.size() - << "\n" << m_eMesh.demangled_stacktrace(20) - << std::endl; - } - Kokkos::UnorderedMapInsertResult insert_result = m_cell_2_data_map.insert(subDimEntity,data); - if(insert_result.failed()){ - std::cout << " the entity failed to be inserted. Rehashing with an increased size limit\n"; - unsigned current_cap = m_cell_2_data_map.capacity(); - unsigned size_delta = 128; //arbitrarily chosen. Need to find a smarter way to size the map - m_cell_2_data_map.rehash(current_cap + size_delta); //NOTE: THIS IS ~NOT~ A DEVICE FUNCTION, IT ~CANNOT~ BE CALLED FROM A PARALLEL KERNEL. - insert_result = m_cell_2_data_map.insert(subDimEntity,data); - if(insert_result.failed()) - std::cout << "Failed second insert attempt\n"; - } - } - - /// @param needNodes should be true in general; it's used by registerNeedNewNode to generate actual data or not on the subDimEntity - /// For local refinement, subDimEntity's needs are not always known uniquely by the pair {elementId, iSubDimOrd}; for example, in - /// an element-based marking scheme, the shared face between two elements may be viewed differently. So, we need the ability to - /// override the default behavior of always creating new nodes on the subDimEntity, but still allow the entity to be created in - /// the NodeRegistry_KOKKOS databse. - - typedef bool (NodeRegistry_KOKKOS::*ElementFunctionPrototype)( const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes, - const CellTopologyData * const bucket_topo_data - ); - - /// this is a helper method that loops over all sub-dimensional entities whose rank matches on of those in @param needed_entity_ranks - /// and registers that sub-dimensional entity as needing a new node. - /// @param isGhost should be true if this element is a ghost, in which case this will call the appropriate method to set up for - // communications - - void doForAllSubEntities(ElementFunctionPrototype function, const stk::mesh::Entity element, vector& needed_entity_ranks, const CellTopologyData * const bucket_topo_data); - - /// fill - /// @param subDimEntity with the stk::mesh::EntityId's of - /// the ordinal @param iSubDimOrd sub-dimensional entity of - /// @param element of rank - /// @param needed_entity_rank - /// - bool getSubDimEntity(SubDimCell_SDCEntityType& subDimEntity, const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd, - const CellTopologyData * const bucket_topo_data_element_0 =0 - ); - - //======================================================================================================================== - // low-level interface - - unsigned total_size(); - unsigned local_size(); - - /// Replace element ownership - /// When remeshing during unrefinement, replace ownership of sub-dim entities by non-deleted elements - bool replaceElementOwnership(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes, - const CellTopologyData * const bucket_topo_data - ); - - /// communicates pending deletes and removes subDimEntity's as requested - bool replaceElementOwnership(); - - bool initializeEmpty(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes_notUsed, - const CellTopologyData * const bucket_topo_data); - - void setAllReceivedNodeData(); - - /// when a sub-dim entity is visited during node registration but is flagged as not being marked, and thus not requiring - /// any new nodes, we flag it with NR_MARK_NONE, then remove it here - void removeUnmarkedSubDimEntities(); - - void communicate_marks(); - - private: - - void communicate_marks_pack(stk::CommSparse& commAll, int stage); - void communicate_marks_unpack(stk::CommSparse& commAll); - - void do_add_node_sharing_comm(); - - public: - - bool is_empty( const stk::mesh::Entity element, stk::mesh::EntityRank needed_entity_rank, unsigned iSubDimOrd); - void query(std::ostream& out, stk::mesh::EntityId elementId, int rank, unsigned iSubDimOrd, std::string msg="", stk::mesh::EntityId *edge=0, SubDimCell_SDCEntityType* subDimEntityIn=0); - void query(stk::mesh::EntityId elementId, int rank, unsigned iSubDimOrd, std::string msg="", stk::mesh::EntityId *edge=0, SubDimCell_SDCEntityType* subDimEntityIn=0) - { - query(std::cout, elementId, rank, iSubDimOrd, msg, edge, subDimEntityIn); - } - - static std::string print_string(PerceptMesh& eMesh, const SubDimCell_SDCEntityType& subDimEntity) - { - std::ostringstream ost; - ost << "SubDimCell:"; - for (unsigned ii=0; ii < subDimEntity.size(); ++ii) - { - ost << " " << eMesh.identifier(subDimEntity[ii]); - } - return ost.str(); - } - - static void print(PerceptMesh& eMesh, const SubDimCell_SDCEntityType& subDimEntity) - { - std::cout << print_string(eMesh, subDimEntity); - } - - inline stk::mesh::Entity get_entity(stk::mesh::BulkData& bulk, stk::mesh::EntityRank rank, stk::mesh::EntityId id) const - { - return bulk.get_entity(rank, id); - } - - inline stk::mesh::Entity get_entity_element(stk::mesh::BulkData& bulk, stk::mesh::EntityRank rank, stk::mesh::EntityId id) const - { - return get_entity(bulk, rank, id); - } - - NodeIdsOnSubDimEntityType* getNewNodesOnSubDimEntity(const stk::mesh::Entity element, stk::mesh::EntityRank& needed_entity_rank, unsigned iSubDimOrd); - -// void checkDB(std::string msg=""); - - /// allocate the send/recv buffers for all-to-all communication - bool allocateBuffers(); - void communicate(); - void unpack(); - - /// after registering all needed nodes, this method is used to request new nodes on this processor - void createNewNodesInParallel(); - - /// unpacks the incoming information in @param buffer_entry and adds that information to my local node registry - /// (i.e. the map of sub-dimensional entity to global node id is updated) - void createNodeAndConnect(CommDataType& buffer_entry, NodeIdsOnSubDimEntityType& nodeIds_onSE, unsigned from_proc, vector& nodes_to_ghost); - - public: - inline SubDimCellToDataMap_KOKKOS& getMap() { return m_cell_2_data_map; } - inline PerceptMesh& getMesh() { return m_eMesh; } - inline bool getUseCustomGhosting() { return m_useCustomGhosting; } - - /// remove any sub-dim entities from the map that have a node in deleted_nodes - it is assumed - /// that the list of deleted_nodes is parallel-consistent (a node being deleted on proc A, owned - /// by proc B, must also be in proc B's deleted_nodes list) - void cleanDeletedNodes(SetOfEntities& deleted_nodes, - SetOfEntities& kept_nodes_orig_minus_kept_nodes, - SubDimCellToDataMap_KOKKOS& to_save, - bool debug=false); - - void cleanInvalidNodes(bool debug=false); - - // further cleanup of the NodeRegistry_KOKKOS db - some elements get deleted on some procs but the ghost elements - // are still in the db - the easiest way to detect this is as used here: during modification_begin(), - // cleanup of cached transactions are performed, then we find these by seeing if our element id's are - // no longer in the stk_mesh db. - void clear_element_owner_data_phase_2(bool resize_nodeId_data=true, bool mod_begin_end=true, SetOfEntities * elemsToBeDeleted=0); - void dumpDB(std::string msg=""); - - // estimate of memory used by this object - unsigned get_memory_usage(); - - void mod_begin(); - void mod_end(const std::string& msg=""); - - bool verifyAllKeysInBoostNR(NodeRegistry * nr, SetOfEntities& nodesMappedTo, unsigned& noKeysNotInCommon); - - private: - percept::PerceptMesh& m_eMesh; - Refiner *m_refiner; - stk::CommSparse * m_comm_all; - SubDimCellToDataMap_KOKKOS m_cell_2_data_map; - - vector m_nodes_to_ghost; - SetOfEntities m_pseudo_entities; - bool m_useCustomGhosting; - bool m_useAddNodeSharing; - bool m_checkForGhostedNodes; - - void resizeMap(); - - public: - int m_gee_cnt; - int m_gen_cnt; - - bool m_debug; - - NodeRegistryState m_state; - - ElementFunctionPrototype curr_funct; - inline bool do_curr_funct(const stk::mesh::Entity element, NeededEntityType& needed_entity_rank, unsigned iSubDimOrd, bool needNodes, - const CellTopologyData * const bucket_topo_data) - { return (this->*curr_funct)(element,needed_entity_rank,iSubDimOrd,needNodes,bucket_topo_data); } - private: - const unsigned m_waste_tolerance; - }; - - struct doForAllSubEntitiesFunctor - { - NodeRegistry_KOKKOS * m_nr; - const stk::mesh::Entity m_element; - std::vector & m_needed_entity_ranks; - const CellTopologyData * const m_bucket_topo_data; - const unsigned m_numSubDimNeededEntities; - const unsigned m_ineed_ent; - - doForAllSubEntitiesFunctor(NodeRegistry_KOKKOS * nr, const stk::mesh::Entity element, - vector& needed_entity_ranks,const CellTopologyData * const bucket_topo_data,const unsigned numSubDimNeededEntities, const unsigned ineed_ent) : - m_nr(nr), m_element(element), m_needed_entity_ranks(needed_entity_ranks), m_bucket_topo_data(bucket_topo_data), m_numSubDimNeededEntities(numSubDimNeededEntities), m_ineed_ent(ineed_ent) - { - } - - void run() - { - Kokkos::parallel_for(Kokkos::RangePolicy(0,m_numSubDimNeededEntities), *this); - } - - KOKKOS_INLINE_FUNCTION - void operator()(const unsigned& iSubDimOrd) const - { - m_nr->do_curr_funct(m_element, m_needed_entity_ranks[m_ineed_ent], iSubDimOrd, true, m_bucket_topo_data); - } - - }; - } // namespace percept - -#endif diff --git a/packages/percept/src/adapt/Refiner.cpp b/packages/percept/src/adapt/Refiner.cpp index 7d22b3526cb0..c782ad1fc777 100644 --- a/packages/percept/src/adapt/Refiner.cpp +++ b/packages/percept/src/adapt/Refiner.cpp @@ -1561,7 +1561,7 @@ } } } - stk::mesh::Selector owned = stk::mesh::MetaData::get(*m_eMesh.get_bulk_data()).locally_owned_part(); + stk::mesh::Selector owned = m_eMesh.get_bulk_data()->mesh_meta_data().locally_owned_part(); boundarySelector = boundarySelector & owned; smoothGeometry(0, &boundarySelector, option, use_ref_mesh); } diff --git a/packages/percept/src/adapt/UniformRefinerPattern.cpp b/packages/percept/src/adapt/UniformRefinerPattern.cpp index 3801225c467a..c6c9f079d804 100644 --- a/packages/percept/src/adapt/UniformRefinerPattern.cpp +++ b/packages/percept/src/adapt/UniformRefinerPattern.cpp @@ -721,7 +721,7 @@ << " with topo= " << getToTopoPartName() << std::endl; stk::mesh::set_topology(*block_to, stk::mesh::get_topology(shards::CellTopology(getToTopology()), eMesh.get_fem_meta_data()->spatial_dimension())); - if (block_to->attribute() == 0) { + if (!stk::io::is_part_io_part(block_to)) { stk::io::put_io_part_attribute(*block_to); } @@ -890,7 +890,7 @@ std::string ioss_elem_topo = ""; convert_stk_topology_to_ioss_name(elem_topo, ioss_elem_topo); - bool is_io_part = from_superset->attribute() != NULL; + bool is_io_part = stk::io::is_part_io_part(*from_superset); // need to create independent of toPart topology if (is_io_part && from_superset->id()>0) diff --git a/packages/percept/src/adapt/main/MeshAdapt.cpp b/packages/percept/src/adapt/main/MeshAdapt.cpp index 12afac0f9699..d7eabdaa3d9b 100644 --- a/packages/percept/src/adapt/main/MeshAdapt.cpp +++ b/packages/percept/src/adapt/main/MeshAdapt.cpp @@ -2279,6 +2279,8 @@ void MeshAdapt::initialize_m2g_geometry(std::string input_geometry) bd->modification_begin(); + std::vector nodesToCheck; + std::vector procsSharedTo; for (unsigned i = 0; i < curveIDs.size(); i++) { //create beams and put them into corresponding curve parts std::vector add_parts_beams(1, static_cast(0)); @@ -2293,22 +2295,20 @@ void MeshAdapt::initialize_m2g_geometry(std::string input_geometry) bool toDeclare = true; std::vector beamNodeIDs; - std::vector procsSharedTo; - std::vector keysToCheck; int lowestRank = std::numeric_limits::max(); for (unsigned j = 0; j < edge_node_ids[ii].size(); j++) beamNodeIDs.push_back((stk::mesh::EntityId) edge_node_ids[ii][j]); + nodesToCheck.clear(); for (unsigned j = 0; j < edge_node_ids[ii].size(); j++) { stk::mesh::Entity cur_node = bd->get_entity(stk::topology::NODE_RANK, beamNodeIDs[j]); - stk::mesh::EntityKey key = bd->entity_key(cur_node); - keysToCheck.push_back(key); + nodesToCheck.push_back(cur_node); } - bd->shared_procs_intersection(keysToCheck, procsSharedTo); + bd->shared_procs_intersection(nodesToCheck, procsSharedTo); procsSharedTo.push_back(THIS_PROC_NUM); //find all processes that own or have this node shared to it for (size_t iii = 0; iii < procsSharedTo.size(); iii++) { if (procsSharedTo[iii] < lowestRank) @@ -2353,7 +2353,8 @@ void MeshAdapt::initialize_m2g_geometry(std::string input_geometry) int lowestRank = std::numeric_limits::max(); std::vector keysToCheck; - std::vector procsSharedTo; + + procsSharedTo.clear(); //std::vector procsSharedTo; stk::mesh::Entity cur_node; for (unsigned j = 0; j < face_node_ids[ii].size(); j++) { diff --git a/packages/percept/src/adapt/sierra_element/CellTopology.cpp b/packages/percept/src/adapt/sierra_element/CellTopology.cpp index 04d37d7eab47..8747247bc00b 100644 --- a/packages/percept/src/adapt/sierra_element/CellTopology.cpp +++ b/packages/percept/src/adapt/sierra_element/CellTopology.cpp @@ -103,98 +103,6 @@ namespace percept { return s_cellTopologyMap; } - - CellTopologyIdMap & - get_cell_topology_id_map() - { - static CellTopologyIdMap s_cellTopologyMap; - - if (s_cellTopologyMap.empty()) { - s_cellTopologyMap[NODE_0] = getCellTopology(); - s_cellTopologyMap[PARTICLE_1] = getCellTopology(); - s_cellTopologyMap[EDGE_2] = getCellTopology >(); - s_cellTopologyMap[EDGE_3] = getCellTopology >(); - s_cellTopologyMap[SHELL_LINE_2] = getCellTopology >(); - s_cellTopologyMap[SHELL_LINE_3] = getCellTopology >(); - s_cellTopologyMap[ROD_2] = getCellTopology >(); - s_cellTopologyMap[ROD_3] = getCellTopology >(); - s_cellTopologyMap[FACE_TRI_3] = getCellTopology >(); - s_cellTopologyMap[FACE_TRI_4] = getCellTopology >(); - s_cellTopologyMap[FACE_TRI_6] = getCellTopology >(); - s_cellTopologyMap[SHELL_TRI_3] = getCellTopology >(); - // s_cellTopologyMap[SHELL_TRI_4] = getCellTopology >(); - s_cellTopologyMap[SHELL_TRI_6] = getCellTopology >(); - s_cellTopologyMap[FACE_QUAD_4] = getCellTopology >(); - s_cellTopologyMap[FACE_QUAD_8] = getCellTopology >(); - s_cellTopologyMap[FACE_QUAD_9] = getCellTopology >(); - s_cellTopologyMap[SHELL_QUAD_4] = getCellTopology >(); - s_cellTopologyMap[SHELL_QUAD_8] = getCellTopology >(); - s_cellTopologyMap[SHELL_QUAD_9] = getCellTopology >(); - s_cellTopologyMap[SOLID_TET_4] = getCellTopology >(); - s_cellTopologyMap[SOLID_TET_8] = getCellTopology >(); - s_cellTopologyMap[SOLID_TET_10] = getCellTopology >(); - s_cellTopologyMap[SOLID_HEX_8] = getCellTopology >(); - s_cellTopologyMap[SOLID_HEX_20] = getCellTopology >(); - s_cellTopologyMap[SOLID_HEX_27] = getCellTopology >(); - s_cellTopologyMap[SOLID_PYRAMID_5] = getCellTopology >(); - s_cellTopologyMap[SOLID_PYRAMID_13] = getCellTopology >(); - s_cellTopologyMap[SOLID_PYRAMID_14] = getCellTopology >(); - s_cellTopologyMap[SOLID_WEDGE_6] = getCellTopology >(); - s_cellTopologyMap[SOLID_WEDGE_15] = getCellTopology >(); - s_cellTopologyMap[SOLID_WEDGE_18] = getCellTopology >(); - s_cellTopologyMap[FACE_PENT_5] = getCellTopology >(); - s_cellTopologyMap[FACE_HEX_6] = getCellTopology >(); - } - - return s_cellTopologyMap; - } - - - IdCellTopologyMap & - get_id_cell_topology_map() - { - static IdCellTopologyMap s_cellTopologyMap; - - if (s_cellTopologyMap.empty()) { - s_cellTopologyMap[getCellTopology()] = NODE_0; - s_cellTopologyMap[getCellTopology()] = PARTICLE_1; - s_cellTopologyMap[getCellTopology >()] = EDGE_2; - s_cellTopologyMap[getCellTopology >()] = EDGE_3; - s_cellTopologyMap[getCellTopology >()] = SHELL_LINE_2; - s_cellTopologyMap[getCellTopology >()] = SHELL_LINE_3; - s_cellTopologyMap[getCellTopology >()] = ROD_2; - s_cellTopologyMap[getCellTopology >()] = ROD_3; - s_cellTopologyMap[getCellTopology >()] = FACE_TRI_3; - s_cellTopologyMap[getCellTopology >()] = FACE_TRI_4; - s_cellTopologyMap[getCellTopology >()] = FACE_TRI_6; - s_cellTopologyMap[getCellTopology >()] = SHELL_TRI_3; - // s_cellTopologyMap[getCellTopology >()] = SHELL_TRI_4; - s_cellTopologyMap[getCellTopology >()] = SHELL_TRI_6; - s_cellTopologyMap[getCellTopology >()] = FACE_QUAD_4; - s_cellTopologyMap[getCellTopology >()] = FACE_QUAD_8; - s_cellTopologyMap[getCellTopology >()] = FACE_QUAD_9; - s_cellTopologyMap[getCellTopology >()] = SHELL_QUAD_4; - s_cellTopologyMap[getCellTopology >()] = SHELL_QUAD_8; - s_cellTopologyMap[getCellTopology >()] = SHELL_QUAD_9; - s_cellTopologyMap[getCellTopology >()] = SOLID_TET_4; - s_cellTopologyMap[getCellTopology >()] = SOLID_TET_8; - s_cellTopologyMap[getCellTopology >()] = SOLID_TET_10; - s_cellTopologyMap[getCellTopology >()] = SOLID_HEX_8; - s_cellTopologyMap[getCellTopology >()] = SOLID_HEX_20; - s_cellTopologyMap[getCellTopology >()] = SOLID_HEX_27; - s_cellTopologyMap[getCellTopology >()] = SOLID_PYRAMID_5; - s_cellTopologyMap[getCellTopology >()] = SOLID_PYRAMID_13; - s_cellTopologyMap[getCellTopology >()] = SOLID_PYRAMID_14; - s_cellTopologyMap[getCellTopology >()] = SOLID_WEDGE_6; - s_cellTopologyMap[getCellTopology >()] = SOLID_WEDGE_15; - s_cellTopologyMap[getCellTopology >()] = SOLID_WEDGE_18; - s_cellTopologyMap[getCellTopology >()] = FACE_PENT_5; - s_cellTopologyMap[getCellTopology >()] = FACE_HEX_6; - } - - return s_cellTopologyMap; - } - } // namespace @@ -218,53 +126,6 @@ namespace percept { } - CellTopology - getBasicCellTopology( - TopologyId id) - { - if (id == INVALID) - return CellTopology(); - - CellTopologyIdMap &cell_topology_map = get_cell_topology_id_map(); - - CellTopologyIdMap::const_iterator it = cell_topology_map.find(id); - - if (it == cell_topology_map.end()) - { - //throw RuntimeError() << "Cell topology " << id << " is not defined"; - std::ostringstream msg; - msg << "Cell topology " << id << " is not defined"; - throw std::runtime_error(msg.str()); - } - - return (*it).second; - } - - - TopologyId - getCellTopologyId( - const CellTopology & cell_topology) - { - if (cell_topology.getCellTopologyData() == 0) - return INVALID; - - IdCellTopologyMap &cell_topology_map = get_id_cell_topology_map(); - - IdCellTopologyMap::const_iterator it = cell_topology_map.find(cell_topology); - - if (it == cell_topology_map.end()) - { - //throw RuntimeError() << "Cell topology " << cell_topology.getName() << " is not defined"; - std::ostringstream msg; - msg << "Cell topology " << cell_topology.getName() << " is not defined"; - throw std::runtime_error(msg.str()); - - } - - return (*it).second; - } - - Elem::CellTopology edgeCellTopology( const Elem::CellTopology & cell_topology, @@ -298,112 +159,6 @@ namespace percept { return cell_topology.getDimension() == 3 ? faceCellTopology(cell_topology, ordinal) : edgeCellTopology(cell_topology, ordinal); } - - bool isElement(const Elem::CellTopology &cell_topology, unsigned spatial_dimension) { - bool is_element = isShellElement(cell_topology) - || isRodElement(cell_topology) || isParticleElement(cell_topology) - || cell_topology.getDimension() == spatial_dimension; // isSolidElement(cell_topology, spatial_dimension); - - return is_element; - } - - - bool isSolidElement(const Elem::CellTopology &cell_topology, unsigned spatial_dimension) { - bool is_solid = cell_topology.getDimension() == spatial_dimension - && !isShellElement(cell_topology) && !isRodElement(cell_topology) && !isParticleElement(cell_topology); - - return is_solid; - } - - - bool isShellElement(const Elem::CellTopology &cell_topology) { - return cell_topology.getSideCount() == 2; - } - - - bool isRodElement(const Elem::CellTopology &cell_topology) { - bool is_rod = cell_topology.getDimension() == 2 - && cell_topology.getSideCount() == 1; - - return is_rod; - } - - - bool isParticleElement(const Elem::CellTopology &cell_topology) { - bool is_particle = cell_topology.getDimension() == 1 - && cell_topology.getNodeCount() == 1; - - return is_particle; - } - - - const unsigned * - getNodesOfEdge( - const Elem::CellTopology & cell_topology, - unsigned edge) - { - if (edge >= cell_topology.getEdgeCount()) - return NULL; - - // Get the topology to test the bounds of subcell_dim and subcell_ord. - cell_topology.getCellTopologyData(1, edge); - return cell_topology.getCellTopologyData()->subcell[1][edge].node; - } - - - const unsigned * - getNodesOfFace( - const Elem::CellTopology & cell_topology, - unsigned face) - { - if (face >= cell_topology.getFaceCount()) - return NULL; - - // Get the topology to test the bounds of subcell_dim and subcell_ord. - cell_topology.getCellTopologyData(2, face); - return cell_topology.getCellTopologyData()->subcell[2][face].node; - } - - - const unsigned * - getNodesOfSide( - const Elem::CellTopology & cell_topology, - unsigned side) - { - return cell_topology.getDimension() == 3 ? getNodesOfFace(cell_topology, side) : getNodesOfEdge(cell_topology, side); - } - - - int - getEdgeNode( - const Elem::CellTopology & cell_topology, - unsigned edge, - unsigned node_of_edge) - { - return cell_topology.getNodeMap(1, edge, node_of_edge); - } - - - int - getFaceNode( - const Elem::CellTopology & cell_topology, - unsigned face, - unsigned node_of_face) - { - return cell_topology.getNodeMap(2, face, node_of_face); - } - - - int - getSideNode( - const Elem::CellTopology & cell_topology, - unsigned side, - unsigned node_of_side) - { - return cell_topology.getDimension() == 3 ? getFaceNode(cell_topology, side, node_of_side) : getEdgeNode(cell_topology, side, node_of_side); - } - - int getFaceEdge( const Elem::CellTopology & cell_topology, @@ -413,82 +168,6 @@ namespace percept { return ::mapCellFaceEdge(cell_topology.getBaseCellTopologyData(), face, edge_of_face); } - - Elem::CellTopology - nodeCellTopology( - const Elem::CellTopology & cell_topology, - UInt ordinal) - { - return (ordinal < cell_topology.getNodeCount()) ? cell_topology.getCellTopologyData(0, ordinal) : NULL; - } - - - bool - isCellTopologySubsetOf( - const Elem::CellTopology & cell_topology, - const Elem::CellTopology & richer) - { - bool result = false; - UInt i, j; - - if (0 < cell_topology.getVertexCount()) { - result = - cell_topology.getDimension() == richer.getDimension() && - cell_topology.getVertexCount() == richer.getVertexCount() && - cell_topology.getEdgeCount() == richer.getEdgeCount() && - cell_topology.getFaceCount() == richer.getFaceCount() && - cell_topology.getNodeCount() <= richer.getNodeCount(); - - for (i = 0; result && i < cell_topology.getEdgeCount(); ++i) { - result = isCellTopologySubsetOf(edgeCellTopology(cell_topology, i), edgeCellTopology(richer, i)); - - for (j = 0; result && j < edgeCellTopology(cell_topology, i).getNodeCount(); ++j) - result = getEdgeNode(cell_topology, i, j) == getEdgeNode(richer, i, j); - } - - for (i = 0; i < cell_topology.getFaceCount() && result; ++i) { - result = isCellTopologySubsetOf(faceCellTopology(cell_topology, i), faceCellTopology(richer, i)); - - for (j = 0; result && j < faceCellTopology(cell_topology, i).getNodeCount(); ++j) - result = getFaceNode(cell_topology, i, j) == getFaceNode(richer, i, j); - } - } - - return result; - } - - - unsigned - getParametricDimension( - const Elem::CellTopology & cell_topology) - { - unsigned parametric_dimension = cell_topology.getDimension(); - if (isShellElement(cell_topology) || isRodElement(cell_topology) || isParticleElement(cell_topology)) - --parametric_dimension; - - return parametric_dimension; - } - - int - findReversePermutation( - const CellTopologyData & top, - int permutation_ord) - { - const int nv = top.vertex_count ; - const int np = top.permutation_count ; - const unsigned * const perm_node = top.permutation[permutation_ord].node ; - int p = 0 ; - for ( ; p < np ; ++p ) { - const unsigned * const reverse_perm_node = top.permutation[p].node ; - int j = 0 ; - for ( ; j < nv && reverse_perm_node[j] == perm_node[nv - 1 - j] ; ++j ) - ; - if ( nv == j ) - break ; - } - return p ; - } - } // namespace Elem } // namespace percept diff --git a/packages/percept/src/adapt/sierra_element/CellTopology.hpp b/packages/percept/src/adapt/sierra_element/CellTopology.hpp index 37c040b10878..6b1bdc4501ab 100644 --- a/packages/percept/src/adapt/sierra_element/CellTopology.hpp +++ b/packages/percept/src/adapt/sierra_element/CellTopology.hpp @@ -75,53 +75,16 @@ TopologyId getCellTopologyId(const CellTopology &cell_topology); - CellTopology getBasicCellTopology(TopologyId id); CellTopology getBasicCellTopology(const char *name); - CellTopology nodeCellTopology(const CellTopology &cell_topology, UInt ordinal); CellTopology edgeCellTopology(const CellTopology &cell_topology, UInt ordinal); CellTopology faceCellTopology(const CellTopology &cell_topology, UInt ordinal); CellTopology sideCellTopology(const CellTopology &cell_topology, UInt ordinal); - bool isElement(const CellTopology &cell_topology, unsigned spatial_dimension); - bool isSolidElement(const CellTopology &cell_topology, unsigned spatial_dimension); - bool isShellElement(const CellTopology &cell_topology); - bool isRodElement(const CellTopology &cell_topology); - bool isParticleElement(const CellTopology &cell_topology); - - int getEdgeNode(const CellTopology &cell_topology, unsigned edge, unsigned node_of_edge); - int getFaceNode(const CellTopology &cell_topology, unsigned face, unsigned node_of_face); - int getSideNode(const CellTopology &cell_topology, unsigned side, unsigned node_of_side); int getFaceEdge(const CellTopology &cell_topology, unsigned face, unsigned edge_of_face); - const unsigned *getNodesOfEdge(const CellTopology &cell_topology, unsigned edge); - const unsigned *getNodesOfFace(const CellTopology &cell_topology, unsigned face); - const unsigned *getNodesOfSide(const CellTopology &cell_topology, unsigned side); - - bool isCellTopologySubsetOf(const CellTopology &cell_topology, const CellTopology &richer); - - unsigned getParametricDimension(const CellTopology &cell_topology) ; - - int findReversePermutation(const CellTopologyData &top, int permutation_ord); - } // namespace Elem } // namespace percept -// namespace shards { - -// inline bool operator<(const CellTopology &c1, const CellTopology &c2) { -// return c1.getTopology() < c2.getTopology(); -// } - -// inline bool operator==(const CellTopology &c1, const CellTopology &c2) { -// return c1.getTopology() == c2.getTopology(); -// } - -// inline bool operator!=(const CellTopology &c1, const CellTopology &c2) { -// return !(c1 == c2); -// } - -// } // namespace shards - #endif // adapt_sierra_element_CellTopology_hpp diff --git a/packages/percept/src/adapt/sierra_element/RefinementKey.cpp b/packages/percept/src/adapt/sierra_element/RefinementKey.cpp deleted file mode 100644 index 4f4874924992..000000000000 --- a/packages/percept/src/adapt/sierra_element/RefinementKey.cpp +++ /dev/null @@ -1,99 +0,0 @@ -// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering -// Solutions of Sandia, LLC (NTESS). Under the terms of Contract -// DE-NA0003525 with NTESS, the U.S. Government retains certain rights -// in this software. -// -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - - - -#include - -#include - -namespace percept { -namespace Elem { - -/* - Derek Gaston: I've pretty much disabled this class so it always returns 0x0. - This is because I've removed heterogeneous refinement from the framework. - Hopefully this class will get completely removed one day... -*/ - -RefinementKey::RefinementKey() -{ - //Default constructor will assign 0x0 to value_ - //This will correspond to normal homogeneous refinement patterns - value_=0x0; -} -RefinementKey::RefinementKey(UInt valueArgue) -{ - //constructor where value is known -// value_=valueArgue; - value_=0x0; - -} - -RefinementKey::~RefinementKey() -{ -//Destructor -//Do Nothing -} - -UInt Elem::RefinementKey::value() -{ - return value_; -} - -void Elem::RefinementKey::assign_value(UInt valueArgue) -{ - //Assigns value_ based on a given id # - //Usefull when changing value to default during unrefinement - //as well as when assigning the unrefinement template id # - -// value_ = valueArgue; - -} - -void Elem::RefinementKey::assign_value(std::vector & edge_order) -{ - //Assigns value_ based on a given edge order cut vector - //Used durring heterogeneous marker resolution - -/* value_ = 0x0; - - for (UInt position = 0; position < edge_order.size(); ++position) { - value_ |= (edge_order[position]+1)<<(position*4); - } -*/ -} - -std::vector Elem::RefinementKey::ordered_cut_edges( UInt numEdges ) const -{ - //Given the number of edges of a mesh object populates a vector - //of edges to be refined in the correct order based on value_ - - std::vector edge_order; - edge_order.clear(); - - for (UInt iedge = 0; iedge>(iedge*4) & 0xF; - if(edge != 0x0) edge_order.push_back(edge-1); - } - return edge_order; - -} - - -bool Elem::RefinementKey::full_refinement( UInt numEdges ) -{ - //if number of edges on refinement key equal numEdges than full refinement - std::vector edgeOrder = ordered_cut_edges(numEdges); - UInt edgeOrderSize = edgeOrder.size(); - return edgeOrderSize == numEdges; -} - -} // namespace Elem -} // namespace percept - diff --git a/packages/percept/src/adapt/sierra_element/RefinementKey.hpp b/packages/percept/src/adapt/sierra_element/RefinementKey.hpp deleted file mode 100644 index 7dd457e28111..000000000000 --- a/packages/percept/src/adapt/sierra_element/RefinementKey.hpp +++ /dev/null @@ -1,79 +0,0 @@ -// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering -// Solutions of Sandia, LLC (NTESS). Under the terms of Contract -// DE-NA0003525 with NTESS, the U.S. Government retains certain rights -// in this software. -// -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -/** - * @file - * @specifications for heterogeneous key id to give refinement topologies - * @author Warren B. Tauber - * - * @par Description: - * - * In order to support removal of hanging nodes elements can be refined - * according to topologies that do not refine every edge and or face. To - * account for the different possibilities a key id is given to each mesh - * object. This key id will contain the edges to be cut information . - * Given the key id number member functions of this class will be able to - * give a vector of ordered edges to be cut. Also the reverse will be - * available as well. These functions will be different for each base - * topology class. - * - * value_ will be such that for 2d solid triangular elements - * a value of 0x021 would be edge 0 cut first then edge 1 - * a value 0f 0x321 would be edge 0 then edge 1 then edge 2 - * value of digits in key_id corresponds to (edge# +1 ) in left to right - * ordering - */ - -/*--------------------------------------------------------------------*/ - -#ifndef adapt_sierra_element_Elem_RefinementKey_hpp -#define adapt_sierra_element_Elem_RefinementKey_hpp - -#include -#include - -namespace percept { - namespace Elem { - - class RefinementKey { - private: - /** Unsigned integer that will store the correct template **/ - UInt value_; - - public: - - /** Default Constructer that assigns the value_ to 0x0 **/ - RefinementKey(); - - /** Constructer where value is known **/ - explicit RefinementKey(UInt valueArg); - - /** Destructor **/ - ~RefinementKey(); - - /** Returns the value_ **/ - UInt value(); - - /** Assigns value_**/ - void assign_value(UInt valueArg) ; - - /** function that assigns the value_ based on edges to be cut **/ - void assign_value( std::vector & edge_order ) ; - - /** function that returns ordered edges to be cut **/ - /** for this method need the objTopologyes number of edges **/ - std::vector ordered_cut_edges( UInt numEdges ) const ; - - /** Function returns boolean for full refinement **/ - bool full_refinement( UInt numEdges ); - }; - - } // namespace Elem -} // namespace percept - -#endif // adapt_sierra_element_Elem_RefinementKey_hpp diff --git a/packages/percept/src/adapt/sierra_element/RefinementTopology.cpp b/packages/percept/src/adapt/sierra_element/RefinementTopology.cpp index bec38dab40dc..d2cb31691239 100644 --- a/packages/percept/src/adapt/sierra_element/RefinementTopology.cpp +++ b/packages/percept/src/adapt/sierra_element/RefinementTopology.cpp @@ -139,705 +139,6 @@ using namespace percept; } - /*--------------------------------------------------------------------*/ - /*----MeshObjRefinementTopology----------------------------------*/ - /*--------------------------------------------------------------------*/ - MeshObjRefinementTopology::MeshObjRefinementTopology() - : m_numChild(0), - m_numChildNodes(0), - m_childCellTopology(0), - m_childNode(0), - m_homogeneous(true), - m_fullRefinement(false) - {} - - - MeshObjRefinementTopology::~MeshObjRefinementTopology() - { - //Destructor - //if homogeneous do not delete data as they refer to static variables - //if not homogeneous release refinement pattern data - - if(!m_homogeneous){ - - //Delete 2 dimensional array of childNode - for(UInt i = 0; igetCellTopology() : CellTopology(); - return (child < m_numChild) ? m_childCellTopology[child] : CellTopology(); - } - - - const UInt * - MeshObjRefinementTopology::child_node( - UInt child) const - { - return (child < m_numChild) ? m_childNode[child] : NULL; - } - - bool - MeshObjRefinementTopology::homogeneous_refinement() const - { - return m_homogeneous; - } - - - bool - MeshObjRefinementTopology::full_refinement() const - { - return m_fullRefinement; - } - - - std::pair - MeshObjRefinementTopology::child_face( - const UInt face_ordinal, - const UInt face_child_ordinal, - const CellTopology & objTop, - const RefinementKey & objDesiredKey) const - { - std::pair result(0,0); - - bool flag = (face_ordinal < objTop.getFaceCount()) && - (face_child_ordinal < getRefinementTopology(faceCellTopology(objTop, face_ordinal))->num_child()); - - if (flag) { - - // Match nodes of - // face_topology(face_ordinal)->child_nodes(face_child_node)[*] - // to - // child_topology(face_child.first)->face_nodes(face_child.second)[*] - - - - const CellTopology faceTop = Elem::faceCellTopology(objTop, face_ordinal); - - const UInt * const faceN = getRefinementTopology(objTop)->face_node(face_ordinal); - - //const UInt * const faceChildN = faceTop.child_node(face_child_ordinal); - - //faceChildN needs to be replaced with faceDesiredTop.child_node - //Need face desired key - - std::vector sideEdgePattern; - sideEdgePattern.clear(); - std::vector objEdgePattern = objDesiredKey.ordered_cut_edges(objTop.getEdgeCount()); - for(UInt jEdge = 0; jEdgequery_refinement_topology(faceDesiredKey,faceDesiredTop)) { - //throw RuntimeError() << "Failed query of face refinement topology." - // << std::endl << StackTrace; - THROW( "Failed query of face refinement topology."); - } - - const UInt * const faceChildN = faceDesiredTop.child_node(face_child_ordinal); - - const CellTopology faceChildTop = getRefinementTopology(faceTop)->child_cell_topology(face_child_ordinal); - - flag = false; - std::pair cf; - - for (cf.first = 0; cf.first < num_child(); ++cf.first) { - const UInt * const childN = child_node( cf.first); - const CellTopology childTop = child_cell_topology(cf.first); - - for (cf.second = 0; cf.second < childTop.getFaceCount(); ++cf.second) { - if (Elem::faceCellTopology(childTop, cf.second) == faceChildTop) { - - const UInt * const childFaceN = getRefinementTopology(childTop)->face_node(cf.second); - - UInt i = 0; - - for (; i < faceChildTop.getNodeCount() && - faceN[ faceChildN[i] ] == childN[ childFaceN[i] ]; ++i); - - if (i == faceChildTop.getNodeCount()) { - VERIFY_TRUE_ON(! flag); - flag = true; - result.first = cf.first; - result.second = cf.second; - } - } - } - } - } - - if (! flag) - { - // throw RuntimeError() << "MeshObjRefinementTopology[" << objTop.getName() << "]::child_face(" - // << face_ordinal << "," << face_child_ordinal - // << ") could not find a match." << std::endl << StackTrace; - THROW("MeshObjRefinementTopology[" << objTop.getName() << "]::child_face(" - << face_ordinal << "," << face_child_ordinal - << ") could not find a match." ); - - } - - return result; - } - - - std::pair - MeshObjRefinementTopology::child_edge( - const UInt edge_ordinal, - const UInt edge_child_ordinal, - const CellTopology & objTop) const - { - std::pair result(0,0); - - bool flag = (edge_ordinal < objTop.getEdgeCount()) && - (edge_child_ordinal < getRefinementTopology(edgeCellTopology(objTop, edge_ordinal))->num_child()); - - if (flag) { - - // Match nodes of - // edge_topology(edge_ordinal)->child_nodes(edge_child_node)[*] - // to - // child_topology(edge_child.first)->edge_nodes(edge_child.second)[*] - - const CellTopology edgeTop = edgeCellTopology(objTop, edge_ordinal); - - const UInt * const edgeN = getRefinementTopology(objTop)->edge_node(edge_ordinal); - const UInt * const edgeChildN = getRefinementTopology(edgeTop)->child_node(edge_child_ordinal); - - const CellTopology edgeChildTop = getRefinementTopology(edgeTop)->child_cell_topology(edge_child_ordinal); - - flag = false; - std::pair ce; - for (ce.first = 0; ce.first < num_child(); ++ce.first) { - - const UInt * const childN = child_node( ce.first); - const CellTopology childTop = child_cell_topology(ce.first); - - for (ce.second = 0; ce.second < childTop.getEdgeCount(); ++ce.second) { - - if (Elem::edgeCellTopology(childTop, ce.second) == edgeChildTop) { - - const UInt * const childEdgeN = getRefinementTopology(childTop)->edge_node(ce.second); - UInt i = 0; - - for (; i < edgeChildTop.getNodeCount() && - edgeN[ edgeChildN[i] ] == childN[ childEdgeN[i] ]; ++i); - - - if (i == edgeChildTop.getNodeCount()) { - flag = true; - result.first = ce.first; - result.second = ce.second; - - } - } - } - } - } - - if (! flag) - { - // throw RuntimeError() << "MeshObjRefinementTopology[" << objTop.getName() << "]::child_edge(" - // << edge_ordinal << "," << edge_child_ordinal - // << ") could not find a match." << std::endl << StackTrace; - THROW("MeshObjRefinementTopology[" << objTop.getName() << "]::child_edge(" - << edge_ordinal << "," << edge_child_ordinal - << ") could not find a match." ); - } - - return result; - } - - - std::pair - RefinementTopology::child_face( - const UInt face_ordinal, - const UInt face_child_ordinal) const - { - std::pair result(0,0); - - bool flag = (face_ordinal < m_cellTopology.getFaceCount()) && - (face_child_ordinal < getRefinementTopology(faceCellTopology(m_cellTopology, face_ordinal))->num_child()); - - if (flag) { - - // Match nodes of - // face_topology(face_ordinal)->child_nodes(face_child_node)[*] - // to - // child_topology(face_child.first)->face_nodes(face_child.second)[*] - - const CellTopology & faceTop = faceCellTopology(m_cellTopology, face_ordinal); - - const UInt * const faceN = face_node(face_ordinal); - - const UInt * const faceChildN = Elem::getRefinementTopology(faceTop)->child_node(face_child_ordinal); - - const CellTopology faceChildTop = getRefinementTopology(faceTop)->child_cell_topology(face_child_ordinal); - - flag = false; - std::pair cf; - - for (cf.first = 0; cf.first < num_child(); ++cf.first) { - - const UInt * const childN = child_node( cf.first); - const CellTopology childTop = child_cell_topology(cf.first); - - for (cf.second = 0; cf.second < childTop.getFaceCount(); ++cf.second) { - - if (faceCellTopology(childTop, cf.second) == faceChildTop) { - - const UInt * const childFaceN = Elem::getRefinementTopology(childTop)->face_node(cf.second); - - UInt i = 0; - - for (; i < faceChildTop.getNodeCount() && - faceN[ faceChildN[i] ] == childN[ childFaceN[i] ]; ++i); - - if (i == faceChildTop.getNodeCount()) { - VERIFY_TRUE(! flag); - flag = true; - result.first = cf.first; - result.second = cf.second; - } - } - } - } - } - - if (! flag) { - // throw RuntimeError() << "CellTopology[" << m_cellTopology.getName() << "]::child_face(" - // << face_ordinal << "," << face_child_ordinal - // << ") could not find a match." << std::endl << StackTrace; - THROW( "CellTopology[" << m_cellTopology.getName() << "]::child_face(" - << face_ordinal << "," << face_child_ordinal - << ") could not find a match." ); - } - - return result; - } - - - std::pair - RefinementTopology::child_edge( - const UInt edge_ordinal, - const UInt edge_child_ordinal) const - { - std::pair result(0,0); - - bool flag = (edge_ordinal < m_cellTopology.getEdgeCount()) && - (edge_child_ordinal < getRefinementTopology(edgeCellTopology(m_cellTopology, edge_ordinal))->num_child()); - - if (flag) { - - // Match nodes of - // edge_topology(edge_ordinal)->child_nodes(edge_child_node)[*] - // to - // child_topology(edge_child.first)->edge_nodes(edge_child.second)[*] - - const CellTopology & edgeTop = edgeCellTopology(m_cellTopology, edge_ordinal); - - const UInt * const edgeN = edge_node(edge_ordinal); - const UInt * const edgeChildN = getRefinementTopology(edgeTop)->child_node(edge_child_ordinal); - - const CellTopology edgeChildTop = getRefinementTopology(edgeTop)->child_cell_topology(edge_child_ordinal); - - flag = false; - std::pair ce; - - for (ce.first = 0; ce.first < num_child(); ++ce.first) { - - const UInt * const childN = child_node( ce.first); - const CellTopology childTop = child_cell_topology(ce.first); - - for (ce.second = 0; ce.second < childTop.getEdgeCount(); ++ce.second) { - - if (edgeCellTopology(childTop, ce.second) == edgeChildTop) { - - const UInt * const childEdgeN = getRefinementTopology(childTop)->edge_node(ce.second); - - UInt i = 0; - - for (; i < edgeChildTop.getNodeCount() && - edgeN[ edgeChildN[i] ] == childN[ childEdgeN[i] ]; ++i); - - if (i == edgeChildTop.getNodeCount()) { - VERIFY_TRUE(! flag); - flag = true; - result.first = ce.first; - result.second = ce.second; - } - } - } - } - } - - if (! flag) { - // throw RuntimeError() << "CellTopology[" << m_cellTopology.getName() << "]::child_edge(" - // << edge_ordinal << "," << edge_child_ordinal - // << ") could not find a match." << std::endl << StackTrace; - THROW("CellTopology[" << m_cellTopology.getName() << "]::child_edge(" - << edge_ordinal << "," << edge_child_ordinal - << ") could not find a match." ); - } - - return result; - } - - - bool - RefinementTopology::query_refinement_topology( - RefinementKey & objKey, - MeshObjRefinementTopology & objRefTop) const - { - // check if partial refinement key defined - if (0x0 != objKey.value()) { - // check if tri or tet element - if (this == getRefinementTopology(getCellTopology >())) { - return refine_rivara_tri(objKey, objRefTop); - } - else if (this == getRefinementTopology(getCellTopology >())) { - return refine_rivara_tet(objKey, objRefTop); - } - } - else { - //Return normal homogeneous patterns due to default value of 0x0 - //or the fact that obj topology name is currently not supported for partial refinement - //Need to return differnet pattens for all normal homogeneous objects - //these refinement patterns will not be dynamically created but will - //point to or be equal to the current static objects; - objRefTop.m_homogeneous = true; - objRefTop.m_fullRefinement = true; - objRefTop.m_numChild = num_child(); - objRefTop.m_numChildNodes= num_child_nodes(); - objRefTop.m_childCellTopology = child_cell_topology(); - objRefTop.m_childNode = (UInt * *) child_nodes(); - return true; - } - return false; - } - - - // Performs logic for bisection of triangle based on refinement key - // Bisections are performed in order of largest to smallest edge to be refined - bool - RefinementTopology::refine_rivara_tri( - RefinementKey & objKey, - MeshObjRefinementTopology & objRefTop) const - { - //partial refinement key defined so create pattern - - //refTop is not based on standard homogeneous patterns - objRefTop.m_homogeneous = false; - - UInt objNumNodes = m_cellTopology.getNodeCount(); - - //Vector of edges and order in which to be refined - std::vector edgePattern = objKey.ordered_cut_edges(m_cellTopology.getEdgeCount()); - - UInt num_edges_cut = edgePattern.size(); - - //number of child equals num_edges_cut +1 - objRefTop.m_numChild =num_edges_cut+1; - - //number of child nodes = m_cellTopology.getNodeCount() + num_edges_cut for 2D triangle - objRefTop.m_numChildNodes = objNumNodes+num_edges_cut; - - //Determine childNode for different templates based on largest edge bisection - UInt longEdge = edgePattern[0]; - - //For first bisection pNode is new child node need to find the node - //that is the opposite vertex of pNode - UInt node0 = edge_node(longEdge)[0]; - UInt node1 = edge_node(longEdge)[1]; - UInt pNode = edge_node(longEdge)[2]; - UInt oppVertex = objNumNodes; - for (UInt iNode = 0; iNode < objNumNodes; ++iNode) { - if(iNode != node0 && iNode != node1) oppVertex = iNode; - } - VERIFY_TRUE(oppVertex != m_cellTopology.getNodeCount()); - - //Define a vector of of vectors that hold child node maps - //Store all possible child node maps to be genereated - //First Bisection gives 2 child elements - //Then each of those elements gives possibly 2 more child - - std::vector > childNodeVec(6, std::vector(3)); - - //An array that tells wether a given possible child node map is active - bool childActive[6]={false,false,false,false,false,false}; - - //First Bisection - - childNodeVec[0][node0]=node0; - childNodeVec[0][node1]=pNode; - childNodeVec[0][oppVertex]=oppVertex; - - childNodeVec[1][node0]=pNode; - childNodeVec[1][node1]=node1; - childNodeVec[1][oppVertex]=oppVertex; - - childActive[0]=true; - childActive[1]=true; - - UInt iChildFilled = 2; - - //Now determine further bisection if needed - //If one of the elements formed in first bisection is refined than - //They will no longer be active - - for (UInt iEdge=1;iEdge edgePattern = objKey.ordered_cut_edges(objNumEdges); - - objRefTop.m_numChild = 1; - objRefTop.m_numChildNodes = objNumNodes; - UInt longEdge = edgePattern[0]; - UInt node0 = edge_node(longEdge)[0]; - UInt node1 = edge_node(longEdge)[1]; - UInt pNode = edge_node(longEdge)[2]; - - UInt oppEdge = objNumEdges; - for(UInt i=0;i > childNodeVec(14, std::vector(4)); - - bool childActive[14]={false,false,false,false,false,false, - false,false,false,false,false,false,false,false}; - - //First Bisection - objRefTop.m_numChild = objRefTop.m_numChild+1; - objRefTop.m_numChildNodes = objRefTop.m_numChildNodes +1; - - childNodeVec[0][node0]=node0; - childNodeVec[0][node1]=pNode; - childNodeVec[0][oppEdgeNode0]=oppEdgeNode0; - childNodeVec[0][oppEdgeNode1]=oppEdgeNode1; - - childNodeVec[1][node0]=pNode; - childNodeVec[1][node1]=node1; - childNodeVec[1][oppEdgeNode0]=oppEdgeNode0; - childNodeVec[1][oppEdgeNode1]=oppEdgeNode1; - - childActive[0]=true; - childActive[1]=true; - - UInt iChildFilled = 2; - - //After first bisection iterate over edges to be refined - for(UInt iEdge=1;iEdge > RefinementTopology::get_children_on_ordinal(const UInt face_ordinal) const { diff --git a/packages/percept/src/adapt/sierra_element/RefinementTopology.hpp b/packages/percept/src/adapt/sierra_element/RefinementTopology.hpp index 7d8278ad9dcb..80ad5ea230d8 100644 --- a/packages/percept/src/adapt/sierra_element/RefinementTopology.hpp +++ b/packages/percept/src/adapt/sierra_element/RefinementTopology.hpp @@ -14,15 +14,12 @@ #include -#include #include namespace percept { namespace Elem { class MeshObjTopology; - class RefinementKey; - class MeshObjRefinementTopology; class RefinementTopology { @@ -111,28 +108,6 @@ namespace percept { */ const UInt *edge_permutation(UInt orientation) const; - bool query_refinement_topology(RefinementKey &object_key, MeshObjRefinementTopology & refTop) const; - - bool refine_rivara_tri(RefinementKey &, MeshObjRefinementTopology & refTop) const; - bool refine_rivara_tet(RefinementKey &, MeshObjRefinementTopology & refTop) const; - - /*------------------------------------------------------------------*/ - /** Mapping of parent->face->child to parent->child->face - * @pre face_ordinal < num_faces() - * @pre face_child_ordinal < face_topology(face_ordinal)->num_child() - * @return (child_ordinal, child_face_ordinal) - */ - std::pair child_face(const UInt face_ordinal , - const UInt face_child_ordinal) const; - - /** Mapping of parent->edge->child to parent->child->edge - * @pre edge_ordinal < getEdgeCount() - * @pre edge_child_ordinal < edge_topology(edge_ordinal)->num_child() - * @return (child_ordinal, child_edge_ordinal) - */ - std::pair child_edge(const UInt edge_ordinal , - const UInt edge_child_ordinal) const; - std::vector< std::pair > get_children_on_ordinal(const UInt face_ordinal) const; @@ -158,63 +133,6 @@ namespace percept { }; - /** - * Class to hold refinement information for partial and heterogeneous - */ - class MeshObjRefinementTopology - { - public: - MeshObjRefinementTopology(); - ~MeshObjRefinementTopology(); - - UInt num_child() const; - - UInt num_child_nodes() const ; - - CellTopology child_cell_topology(UInt child) const ; - - const UInt * child_node(UInt child) const ; - - bool homogeneous_refinement() const; - - bool full_refinement() const; - /*------------------------------------------------------------------*/ - /** Mapping of parent->face->child to parent->child->face - * @pre face_ordinal < num_faces() - * @pre face_child_ordinal < face_topology(face_ordinal)->num_child() - * @arg objTop need to have objects regular topology as well - * @return (child_ordinal, child_face_ordinal) - */ - std::pair child_face(const UInt face_ordinal , - const UInt face_child_ordinal, - const Elem::CellTopology & objTop , - const RefinementKey &objDesiredKey) const ; - - /** Mapping of parent->edge->child to parent->child->edge - * @pre edge_ordinal < getEdgeCount() - * @pre edge_child_ordinal < edge_topology(edge_ordinal)->num_child() - * @arg objTop need to have objects regular topology as well - * @return (child_ordinal, child_edge_ordinal) - */ - std::pair child_edge(const UInt edge_ordinal , - const UInt edge_child_ordinal, - const Elem::CellTopology & objTop) const; - - - public: // private: - UInt m_numChild; - UInt m_numChildNodes; - const CellTopology * m_childCellTopology; - UInt * * m_childNode; - bool m_homogeneous; /*whether refinement is self simmilar*/ - bool m_fullRefinement; /*whether all edges are to be refined or not*/ - - private: - MeshObjRefinementTopology(const MeshObjRefinementTopology&); - MeshObjRefinementTopology&operator=(const MeshObjRefinementTopology&); - }; - - const RefinementTopology *getRefinementTopology(const Elem::CellTopology &cell_topology); const UInt *getRefinementEdgeNode(const Elem::CellTopology &cell_topology, UInt edge); const UInt *getRefinementFaceNode(const Elem::CellTopology &cell_topology, UInt face); diff --git a/packages/percept/src/adapt/sierra_element/StdMeshObjTopologies.cpp b/packages/percept/src/adapt/sierra_element/StdMeshObjTopologies.cpp index d50c8bc2106c..5b3b043075e1 100644 --- a/packages/percept/src/adapt/sierra_element/StdMeshObjTopologies.cpp +++ b/packages/percept/src/adapt/sierra_element/StdMeshObjTopologies.cpp @@ -64,60 +64,6 @@ static const_top_ptr wedge( UInt nnode); static const_top_ptr pyramid( UInt nnode); - // Node topologies - const MeshObjTopology * Node_0(); - - // Edge topologies - const MeshObjTopology * Edge_2(); // Linear 3D edge - const MeshObjTopology * Edge_3(); // Quadratic 3D edge - - // Triangular face topology, has three edges - const MeshObjTopology * Face_Tri_3(); // Triangle face 3N - const MeshObjTopology * Face_Tri_4(); // Triangle face 4N - const MeshObjTopology * Face_Tri_6(); // Triangle face 6N - - // Quadrilateral face topology, has four edges - const MeshObjTopology * Face_Quad_4(); // Quad face 4N - const MeshObjTopology * Face_Quad_8(); // Quad face 8N - const MeshObjTopology * Face_Quad_9(); // Quad face 9N - - // 3D Particle Element topologies - const MeshObjTopology * Particle_1(); - - // 3D Rod Element topology, has ZERO faces and ONE edge - const MeshObjTopology * Rod_2(); - const MeshObjTopology * Rod_3(); - - // 2D Shell Element topology, has TWO edges - const MeshObjTopology * Shell_Line_2(); - const MeshObjTopology * Shell_Line_3(); - - // 3D Triangular Shell Element topology, has TWO faces and THREE edges - const MeshObjTopology * Shell_Tri_3(); - const MeshObjTopology * Shell_Tri_6(); - - // 3D Quadrilateral Shell Element topology, has TWO faces and FOUR edges - const MeshObjTopology * Shell_Quad_4(); - const MeshObjTopology * Shell_Quad_9(); - - // 3D Tetrahedral Solid Element topology, has FOUR faces, 6 EDGES - const MeshObjTopology * Solid_Tet_4(); // 4-Nodes - const MeshObjTopology * Solid_Tet_8(); // 8-Nodes - const MeshObjTopology * Solid_Tet_10(); // 10-Nodes - - // 3D Wedge Solid Element topology, has FIVE faces and NINE edges - const MeshObjTopology * Solid_Wedge_6(); // 6-Nodes - const MeshObjTopology * Solid_Wedge_15(); // 15-Nodes - - // 3D Hexahedral Solid Element topology, has SIX faces and TWELVE edges - const MeshObjTopology * Solid_Hex_8(); // 8-Nodes - const MeshObjTopology * Solid_Hex_20(); // 20-Nodes - const MeshObjTopology * Solid_Hex_27(); // 27-Nodes - - // 3D Pyramid Solid Element topology, has FIVE faces and EIGHT edges - const MeshObjTopology * Solid_Pyramid_5(); // 5-Nodes - const MeshObjTopology * Solid_Pyramid_13(); // 13-Nodes - const MeshObjTopology * node() { @@ -2936,158 +2882,6 @@ return top ; } - const MeshObjTopology * Node_0() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Node_0()"); /* %TRACE% */ - return node(); - } - - const MeshObjTopology * Particle_1() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Particle_1()"); /* %TRACE% */ - return point(); - } - - /*--------------------------------------------------------------------*/ - /* Line topologies (edges, rods, 2D shell) */ - const MeshObjTopology * Edge_2() { - // // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Edge_2()"); /* %TRACE% */ - return line(NOT_ELEMENT, 2); - } - - const MeshObjTopology * Edge_3() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Edge_3()"); /* %TRACE% */ - return line(NOT_ELEMENT, 3); - } - - const MeshObjTopology * Rod_2() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Rod_2()"); /* %TRACE% */ - return line(ROD, 2); - } - - const MeshObjTopology * Rod_3() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Rod_3()"); /* %TRACE% */ - return line(ROD, 3); - } - - const MeshObjTopology * Shell_Line_2() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Line_2()"); /* %TRACE% */ - return line(SHELL, 2); - } - - const MeshObjTopology * Shell_Line_3() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Line_3()"); /* %TRACE% */ - return line(SHELL, 3); - } - - const MeshObjTopology * Face_Tri_3() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Tri_3()"); /* %TRACE% */ - return tri(NOT_ELEMENT, 3); - } - - const MeshObjTopology * Face_Tri_6() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Tri_6()"); /* %TRACE% */ - return tri(NOT_ELEMENT, 6); - } - - const MeshObjTopology * Face_Tri_4() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Tri_4()"); /* %TRACE% */ - return tri4(NOT_ELEMENT); - } - - const MeshObjTopology * Shell_Tri_3() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Tri_3()"); /* %TRACE% */ - return tri(SHELL, 3); - } - - const MeshObjTopology * Shell_Tri_6() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Tri_6()"); /* %TRACE% */ - return tri(SHELL, 6); - } - - - /*--------------------------------------------------------------------*/ - /* Quadrilaterals with 4, 8, or 9 nodes */ - const MeshObjTopology * Face_Quad_4() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Quad_4()"); /* %TRACE% */ - return quad(NOT_ELEMENT, 4); - } - - const MeshObjTopology * Face_Quad_8() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Quad_8()"); /* %TRACE% */ - return quad(NOT_ELEMENT, 8); - } - - const MeshObjTopology * Face_Quad_9() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Face_Quad_9()"); /* %TRACE% */ - return quad(NOT_ELEMENT, 9); - } - - const MeshObjTopology * Shell_Quad_4() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Quad_4()"); /* %TRACE% */ - return quad(SHELL, 4); - } - - const MeshObjTopology * Shell_Quad_9() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Shell_Quad_9()"); /* %TRACE% */ - return quad(SHELL, 9); - } - - /*--------------------------------------------------------------------*/ - /* Hexahedrons with 8, 20, or 27 nodes */ - const MeshObjTopology * Solid_Hex_8() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Hex_8()"); /* %TRACE% */ - return hex(8); - } - - const MeshObjTopology * Solid_Hex_20() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Hex_20()"); /* %TRACE% */ - return hex(20); - } - - const MeshObjTopology * Solid_Hex_27() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Hex_27()"); /* %TRACE% */ - return hex(27); - } - - /*--------------------------------------------------------------------*/ - /* Tetrahedrons with 4, 8, or 10 nodes */ - const MeshObjTopology * Solid_Tet_4() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Tet_4()"); /* %TRACE% */ - return tet(4); - } - - const MeshObjTopology * Solid_Tet_8() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Tet_8()"); /* %TRACE% */ - return tet(8); - } - - const MeshObjTopology * Solid_Tet_10() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Tet_10()"); /* %TRACE% */ - return tet(10); - } - - /*--------------------------------------------------------------------*/ - /* Pentahedrons (a.k.a. "wedges") with 6 or 15 nodes */ - const MeshObjTopology * Solid_Wedge_6() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Wedge_6()"); /* %TRACE% */ - return wedge(6); - } - - const MeshObjTopology * Solid_Wedge_15() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Wedge_15()"); /* %TRACE% */ - return wedge(15); - } - - /*--------------------------------------------------------------------*/ - /* Pyramids with 5 or 13 nodes */ - const MeshObjTopology * Solid_Pyramid_5() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Pyramid_5()"); /* %TRACE% */ - return pyramid(5); - } - - const MeshObjTopology * Solid_Pyramid_13() { - // /* %TRACE[SPEC]% */ Tracespec trace__("sierra::Solid_Pyramid_13()"); /* %TRACE% */ - return pyramid(13); - } void bootstrap() diff --git a/packages/percept/src/percept/FieldBLAS.hpp b/packages/percept/src/percept/FieldBLAS.hpp index d4eaa3ea6d2f..34ad4bb082cd 100644 --- a/packages/percept/src/percept/FieldBLAS.hpp +++ b/packages/percept/src/percept/FieldBLAS.hpp @@ -24,7 +24,7 @@ void field_axpby( const double beta , stk::mesh::Field *Y) { - const stk::mesh::MetaData & meta = stk::mesh::MetaData::get(bulkdata); + const stk::mesh::MetaData & meta = bulkdata.mesh_meta_data(); stk::mesh::Selector select_used = meta.locally_owned_part() | meta.globally_shared_part() ; @@ -55,7 +55,7 @@ void field_copy_component( stk::mesh::Field *Y, const int index) { - const stk::mesh::MetaData & meta = stk::mesh::MetaData::get(bulkdata); + const stk::mesh::MetaData & meta = bulkdata.mesh_meta_data(); const unsigned nDim = meta.spatial_dimension(); const unsigned field_size_y = Y->max_size(stk::topology::NODE_RANK); diff --git a/packages/percept/src/percept/GeometryVerifier.cpp b/packages/percept/src/percept/GeometryVerifier.cpp index b931b140401d..bbd1c1933c74 100644 --- a/packages/percept/src/percept/GeometryVerifier.cpp +++ b/packages/percept/src/percept/GeometryVerifier.cpp @@ -249,7 +249,7 @@ using namespace Intrepid; */ bool GeometryVerifier::isGeometryBad(stk::mesh::BulkData& bulk, bool printTable, std::vector *volume_histogram) //, stk::mesh::Part& mesh_part ) { - const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk); + const stk::mesh::MetaData& meta = bulk.mesh_meta_data(); const unsigned p_rank = bulk.parallel_rank(); bool checkLocalJacobians = m_checkLocalJacobians; PerceptMesh eMesh(&meta, &bulk, true); @@ -299,7 +299,7 @@ using namespace Intrepid; if ( stk::mesh::is_auto_declared_part(*part) ) continue; - const stk::topology part_cell_topo_data = stk::mesh::MetaData::get(bulk).get_topology(*part); + const stk::topology part_cell_topo_data = bulk.mesh_meta_data().get_topology(*part); //std::cout << "P[" << p_rank << "] part = " << part->name() << " part_cell_topo_data= " << part_cell_topo_data << " topo-name= " // << (part_cell_topo_data ? part_cell_topo_data->name : "null") << std::endl; @@ -348,8 +348,6 @@ using namespace Intrepid; unsigned topoDim = spaceDim; if (isShell) topoDim = 2; - //unsigned spatialDimMeta = stk::mesh::MetaData::get(bulk).spatial_dimension(); - // Rank-3 array with dimensions (C,N,D) for the node coordinates of 3 traingle cells FieldContainer cellNodes(numCells, numNodes, spaceDim); PerceptMesh::fillCellNodes(bulk, bucket, coord_field, cellNodes, spaceDim); diff --git a/packages/percept/src/percept/PerceptMesh.cpp b/packages/percept/src/percept/PerceptMesh.cpp index 5aa503804b5e..8aae58c44844 100644 --- a/packages/percept/src/percept/PerceptMesh.cpp +++ b/packages/percept/src/percept/PerceptMesh.cpp @@ -881,7 +881,7 @@ ost << sep << pv[ii]->name(); if (extra_info) { - bool is_io_part = pv[ii]->attribute() != NULL; + bool is_io_part = stk::io::is_part_io_part(pv[ii]); ost << " " << pv[ii]->primary_entity_rank() << " " << pv[ii]->topology() << " IO: " << is_io_part; } @@ -1866,9 +1866,6 @@ m_comm = comm; m_ownData = true; -#if HAVE_CGNS - m_iocgns_init.reset(new Iocgns::Initializer()); -#endif set_io_broker( new stk::io::StkMeshIoBroker(comm) ); m_isCommitted = false; @@ -2445,34 +2442,11 @@ //checkOmit(in_region. get_commsets(), omit_part ) ; /*const CommSetContainer& */ } - void PerceptMesh::read_metaDataNoCommitCGNSStructured(const std::string& in_filename) - { - Ioss::PropertyManager properties; - Ioss::DatabaseIO * dbi = Ioss::IOFactory::create("cgns", in_filename, Ioss::READ_MODEL, - (MPI_Comm)m_comm, properties); - if (dbi == nullptr || !dbi->ok(true)) { - std::exit(EXIT_FAILURE); - } - - // NOTE: 'region' owns 'db' pointer at this time... - m_cgns_structured_region.reset( new Ioss::Region(dbi, "region_1") ); - - if (m_cgns_structured_region->get_element_blocks().size()) - VERIFY_MSG("specified read of structured grid but this mesh has unstructured blocks"); - - } - // ======================================================================== void PerceptMesh::read_metaDataNoCommit( const std::string& in_filename, const std::string& type) { EXCEPTWATCH; - if (type == "cgns_structured") - { - read_metaDataNoCommitCGNSStructured(in_filename); - return; - } - stk::io::StkMeshIoBroker* mesh_data = m_iossMeshData.get(); mesh_data->add_mesh_database(in_filename, type, stk::io::READ_MESH); @@ -2528,7 +2502,6 @@ this->get_ioss_aliases(part.name(), aliases); for (auto alias : aliases) { - //std::cout << " alias= " << alias << std::endl; if (alias == bname) { return true; @@ -2548,14 +2521,7 @@ { return; } -#if !STK_PERCEPT_LITE - if (getProperty("file_type") == "cgns_structured") - { - m_block_structured_grid.reset(new BlockStructuredGrid(this->parallel(), m_cgns_structured_region.get())); - m_block_structured_grid->read_cgns(); - return; - } -#endif + int step = 1; if (get_database_time_step_count() == 0) @@ -2705,8 +2671,7 @@ if (name.find(PerceptMesh::s_omit_part) != std::string::npos) { //if (!get_rank()) std::cout << "tmp srk checkForPartsToAvoidWriting found omitted part= " << name << std::endl; - const Ioss::GroupingEntity *entity = part.attribute(); - if (entity) + if (stk::io::is_part_io_part(part)) { stk::io::remove_io_part_attribute(part); } @@ -2722,8 +2687,7 @@ std::string name = part.name(); //std::cout << "tmp srk checkForPartsToAvoidWriting found part from get_io_omitted_parts() = " << name << " s_omit_part= " << s_omit_part << std::endl; { - const Ioss::GroupingEntity *entity = part.attribute(); - if (entity) + if (stk::io::is_part_io_part(part)) { //std::cout << "tmp srk checkForPartsToAvoidWriting found part from get_io_omitted_parts() omitted part= " << name << std::endl; stk::io::remove_io_part_attribute(part); @@ -4359,29 +4323,6 @@ void PerceptMesh::add_coordinate_state_fields(const bool output_fields) { -#if !STK_PERCEPT_LITE - if (m_block_structured_grid) - { - int scalarDimension = get_spatial_dim(); // a scalar - - m_block_structured_grid->register_field("coordinates_N", scalarDimension); - m_block_structured_grid->register_field("coordinates_NM1", scalarDimension); - m_block_structured_grid->register_field("coordinates_lagged", scalarDimension); -// m_block_structured_grid->register_field("coordinates_0", scalarDimension); //madbrew : grepping for this coord field reveals that it only gets used for the stkmesh case - - m_block_structured_grid->register_field("cg_g", scalarDimension); - m_block_structured_grid->register_field("cg_r", scalarDimension); - m_block_structured_grid->register_field("cg_d", scalarDimension); - m_block_structured_grid->register_field("cg_s", scalarDimension); - - // edge length and adjacency - m_block_structured_grid->register_field("cg_edge_length", 1); - m_block_structured_grid->register_field("num_adj_elems", 1); //number of elements adjacent to a node, across block boundaries - - return; - } -#endif - m_num_coordinate_field_states = 3; int scalarDimension = get_spatial_dim(); // a scalar @@ -5453,7 +5394,7 @@ void PerceptMesh::add_part(const std::string& part_name, bool make_part_io_part) { stk::mesh::Part& part = get_fem_meta_data()->declare_part(part_name, stk::topology::NODE_RANK); - if (make_part_io_part && part.attribute() == NULL) { + if (make_part_io_part && !stk::io::is_part_io_part(part)) { stk::io::put_io_part_attribute(part); } } @@ -5539,7 +5480,7 @@ stk::mesh::EntitySideVector boundary; // select owned - stk::mesh::Selector owned = stk::mesh::MetaData::get(*get_bulk_data()).locally_owned_part(); + stk::mesh::Selector owned = get_bulk_data()->mesh_meta_data().locally_owned_part(); const stk::mesh::PartVector parts = get_fem_meta_data()->get_parts(); for (unsigned ip=0; ip < parts.size(); ip++) diff --git a/packages/percept/src/percept/PerceptMesh.hpp b/packages/percept/src/percept/PerceptMesh.hpp index c1040a531a9e..6bb852d0f88f 100644 --- a/packages/percept/src/percept/PerceptMesh.hpp +++ b/packages/percept/src/percept/PerceptMesh.hpp @@ -58,14 +58,6 @@ #include #include -#if HAVE_CGNS -# if defined(STK_BUILT_IN_SIERRA) -# include -# else -# include -# endif -#endif - #include #include #include @@ -1057,7 +1049,6 @@ /// read with no commit void read_metaDataNoCommit( const std::string& in_filename, const std::string &type = "exodus" ); - void read_metaDataNoCommitCGNSStructured( const std::string& in_filename ); /// create with no commit void create_metaDataNoCommit( const std::string& gmesh_spec); @@ -1108,7 +1099,6 @@ Teuchos::RCP m_iossMeshData; Teuchos::RCP m_iossMeshDataOut; - std::shared_ptr m_cgns_structured_region; #if !STK_PERCEPT_LITE std::shared_ptr m_block_structured_grid; #endif @@ -1198,9 +1188,6 @@ bool m_markNone; private: -#if HAVE_CGNS - std::shared_ptr m_iocgns_init; -#endif void checkStateSpec(const std::string& function, bool cond1=true, bool cond2=true, bool cond3=true); void checkState(const std::string& function) { diff --git a/packages/percept/src/percept/TopologyVerifier.cpp b/packages/percept/src/percept/TopologyVerifier.cpp index e6e47ede9d7d..ee83114f406b 100644 --- a/packages/percept/src/percept/TopologyVerifier.cpp +++ b/packages/percept/src/percept/TopologyVerifier.cpp @@ -102,7 +102,7 @@ */ bool TopologyVerifier::isTopologyBad(stk::mesh::BulkData& bulk) //, stk::mesh::Part& mesh_part ) { - const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk); + const stk::mesh::MetaData& meta = bulk.mesh_meta_data(); stk::mesh::Field *coord_field = meta.get_field >(stk::topology::NODE_RANK, "coordinates"); diff --git a/packages/percept/src/percept/eigen_verify/EigenVerify.cpp b/packages/percept/src/percept/eigen_verify/EigenVerify.cpp index 2d56797e6d65..471c0acc3f93 100644 --- a/packages/percept/src/percept/eigen_verify/EigenVerify.cpp +++ b/packages/percept/src/percept/eigen_verify/EigenVerify.cpp @@ -141,7 +141,7 @@ void compute_field_error( stk::mesh::Field * xferField, stk::mesh::Field * errorField) { - const stk::mesh::MetaData & meta = stk::mesh::MetaData::get(bulkdata); + const stk::mesh::MetaData & meta = bulkdata.mesh_meta_data(); const unsigned field_size = field->max_size(stk::topology::NODE_RANK); const unsigned nDim = meta.spatial_dimension(); diff --git a/packages/percept/src/percept/function/FieldFunction.hpp b/packages/percept/src/percept/function/FieldFunction.hpp index a49e948637ff..d80ced0f6a67 100644 --- a/packages/percept/src/percept/function/FieldFunction.hpp +++ b/packages/percept/src/percept/function/FieldFunction.hpp @@ -109,7 +109,6 @@ using namespace Intrepid; m_deriv_spec = deriv_spec; Dimensions domain_dimensions = getDomainDimensions(); Dimensions codomain_dimensions = getCodomainDimensions(); - //int meta_dimension = stk::mesh::MetaData::get_meta_data(stk::mesh::MetaData::get(*m_bulkData)).get_spatial_dimension(); int num_grad = deriv_spec.dimension(0); //domain_dimensions.push_back(num_grad); codomain_dimensions.back() = num_grad; @@ -128,7 +127,7 @@ using namespace Intrepid; virtual Teuchos::RCP gradient(int spatialDim=3) { - int meta_dimension = stk::mesh::MetaData::get(*m_bulkData).spatial_dimension(); + int meta_dimension = m_bulkData->mesh_meta_data().spatial_dimension(); VERIFY_OP_ON(meta_dimension, ==, spatialDim, "gradient: mismatch in spatial dimensions"); std::string xyz[] = {"x", "y", "z"}; MDArrayString mda(meta_dimension); @@ -201,7 +200,7 @@ using namespace Intrepid; VERIFY_OP_ON(output_field_values.rank(), <=, 3, "FieldFunction::operator() output_field_values bad rank"); int numInterpPoints = parametric_coordinates.dimension(0); - int spatialDim = stk::mesh::MetaData::get(*m_bulkData).spatial_dimension(); + int spatialDim = m_bulkData->mesh_meta_data().spatial_dimension(); int num_grad = getCodomainDimensions().back(); int numCells = PerceptMesh::size1(bucket_or_element); diff --git a/packages/percept/src/percept/function/FunctionOperator.hpp b/packages/percept/src/percept/function/FunctionOperator.hpp index 1a74dc92c87f..15c77cac26f6 100644 --- a/packages/percept/src/percept/function/FunctionOperator.hpp +++ b/packages/percept/src/percept/function/FunctionOperator.hpp @@ -55,7 +55,7 @@ using namespace Intrepid; } if (!part) { - m_selector = new stk::mesh::Selector(stk::mesh::MetaData::get(m_bulkData).universal_part()); + m_selector = new stk::mesh::Selector(m_bulkData.mesh_meta_data().universal_part()); } else { @@ -73,7 +73,7 @@ using namespace Intrepid; VERIFY_OP_ON(m_selector, !=, 0, "FunctionOperator::init"); delete m_selector; } - m_selector = new stk::mesh::Selector(stk::mesh::MetaData::get(m_bulkData).universal_part()); + m_selector = new stk::mesh::Selector(m_bulkData.mesh_meta_data().universal_part()); m_own_selector = true; } } diff --git a/packages/percept/src/percept/function/internal/IntegratedOp.hpp b/packages/percept/src/percept/function/internal/IntegratedOp.hpp index ee24dd31e2dc..4ca6eb51ef8d 100644 --- a/packages/percept/src/percept/function/internal/IntegratedOp.hpp +++ b/packages/percept/src/percept/function/internal/IntegratedOp.hpp @@ -120,7 +120,7 @@ return false; int cell_dimension = cell_topo.getDimension(); - int meta_dimension = stk::mesh::MetaData::get(bulkData).spatial_dimension(); + int meta_dimension = bulkData.mesh_meta_data().spatial_dimension(); if (cell_dimension == meta_dimension - 1) { @@ -129,7 +129,8 @@ VERIFY_OP_ON(cell_dimension, ==, meta_dimension, "Dimensions don't match"); - CoordinatesFieldType& coord_field = *(stk::mesh::MetaData::get(bulkData)).get_field(stk::topology::NODE_RANK, "coordinates"); + const stk::mesh::MetaData& meta = bulkData.mesh_meta_data(); + CoordinatesFieldType& coord_field = *meta.get_field(stk::topology::NODE_RANK, "coordinates"); // FIXME for fields not on a Node unsigned nDOF = m_nDOFs; @@ -318,12 +319,13 @@ const CellTopologyData * const child_cell_topo_data = stk::mesh::get_cell_topology(mybucket(bulkData, child_element).topology()).getCellTopologyData(); CellTopology child_cell_topo(child_cell_topo_data); int child_cell_dimension = child_cell_topo.getDimension(); - int meta_dimension = stk::mesh::MetaData::get(bulkData).spatial_dimension(); + int meta_dimension = bulkData.mesh_meta_data().spatial_dimension(); // for now, only allow face (or edge) VERIFY_OP_ON(child_cell_dimension, ==, meta_dimension - 1, "Dimensions don't match"); - CoordinatesFieldType& coord_field = *(stk::mesh::MetaData::get(bulkData)).get_field(stk::topology::NODE_RANK, "coordinates"); + const stk::mesh::MetaData& meta = bulkData.mesh_meta_data(); + CoordinatesFieldType& coord_field = *meta.get_field(stk::topology::NODE_RANK, "coordinates"); // FIXME for fields not on a Node unsigned nDOF = m_nDOFs; diff --git a/packages/percept/src/percept/mesh/geometry/recovery/GeometryRecoverySplineFit.cpp b/packages/percept/src/percept/mesh/geometry/recovery/GeometryRecoverySplineFit.cpp index 4190793340d7..4f1737b465a7 100644 --- a/packages/percept/src/percept/mesh/geometry/recovery/GeometryRecoverySplineFit.cpp +++ b/packages/percept/src/percept/mesh/geometry/recovery/GeometryRecoverySplineFit.cpp @@ -28,7 +28,7 @@ stk::mesh::Part& GeometryRecoverySplineFit::create_clone_part(const std::string& clone_name, const stk::mesh::EntityRank rank_to_clone, bool make_part_io_part ) { stk::mesh::Part& clone = m_eMesh.get_fem_meta_data()->declare_part(clone_name, rank_to_clone); - if (make_part_io_part && clone.attribute() == NULL) { + if (make_part_io_part && !stk::io::is_part_io_part(clone)) { stk::io::put_io_part_attribute(clone); } return clone; diff --git a/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.cpp b/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.cpp index e416060c1a64..ab3b941aef55 100644 --- a/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.cpp +++ b/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.cpp @@ -41,7 +41,6 @@ A(2,0) = 0; // (x[1][2] - x[0][2]); A(2,1) = 0; // (x[2][2] - x[0][2]); A(2,2) = 1.0; - //if (m_scale_to_unit) scale_to_unit(A); detJ = det(A); return detJ < 0.0; diff --git a/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.hpp b/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.hpp index 0638e1ab0a11..ddda7147813a 100644 --- a/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.hpp +++ b/packages/percept/src/percept/mesh/geometry/volume/VolumeUtil.hpp @@ -32,12 +32,10 @@ DenseMatrix<3,3> m_dMetric_dA[NNODES_MAX]; double m_grad[NNODES_MAX][NNODES_MAX][3]; int m_num_nodes; - bool m_scale_to_unit; bool m_use_approximate_quadratic_volume; VolumeUtil(bool use_approximate_quadratic_volume=true) : m_num_nodes(0), - m_scale_to_unit(false), m_use_approximate_quadratic_volume(use_approximate_quadratic_volume) { } @@ -78,7 +76,6 @@ A(2,0) = (x1[2] - x0[2]); A(2,1) = (x2[2] - x0[2]); A(2,2) = (x3[2] - x0[2]); - //if (m_scale_to_unit) scale_to_unit(A); detJ = det(A); return detJ < 0.0; @@ -212,69 +209,6 @@ return detJ < 0.0; } - void grad_util_pyramid_3d_new(int ibasis, const DenseMatrix<3,3>& dMdA, double grad[NNODES_MAX][3], const int nnode, const int spd, const int *indices, const int nind) - { - double rst[5][3] = {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{.5,.5,.49999}}; //note: avoid singularity at top vertex - - double r = rst[ibasis][0]; - double s = rst[ibasis][1]; - double t = rst[ibasis][2]; - - for (int i=0; i < nnode; i++) - for (int j=0; j < spd; j++) - grad[i][j]=0.0; - -#define Rule(x,y) x += y -#define GRAD(i,j) grad[i-1][j-1] -#define DMDA(i,j) dMdA(i-1,j-1) -#define Sqrt(x) std::sqrt(x) - Rule(GRAD(1,1),-(((-1 + s + t)*DMDA(1,1))/(-1 + 2*t))); - Rule(GRAD(1,1),-(((-1 + r + t)*DMDA(1,2))/(-1 + 2*t))); - Rule(GRAD(1,1),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(1,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(1,2),-(((-1 + s + t)*DMDA(2,1))/(-1 + 2*t))); - Rule(GRAD(1,2),-(((-1 + r + t)*DMDA(2,2))/(-1 + 2*t))); - Rule(GRAD(1,2),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(1,3),-(((-1 + s + t)*DMDA(3,1))/(-1 + 2*t))); - Rule(GRAD(1,3),-(((-1 + r + t)*DMDA(3,2))/(-1 + 2*t))); - Rule(GRAD(1,3),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(2,1),((-1 + s + t)*DMDA(1,1))/(-1 + 2*t)); - Rule(GRAD(2,1),((r - t)*DMDA(1,2))/(-1 + 2*t)); - Rule(GRAD(2,1),-(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(1,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(2,2),((-1 + s + t)*DMDA(2,1))/(-1 + 2*t)); - Rule(GRAD(2,2),((r - t)*DMDA(2,2))/(-1 + 2*t)); - Rule(GRAD(2,2), - -(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2)))); Rule(GRAD(2,3),((-1 + s + t)*DMDA(3,1))/(-1 + 2*t)); - Rule(GRAD(2,3),((r - t)*DMDA(3,2))/(-1 + 2*t)); - Rule(GRAD(2,3), - -(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2)))); Rule(GRAD(3,1),((-s + t)*DMDA(1,1))/(-1 + 2*t)); - Rule(GRAD(3,1),((-r + t)*DMDA(1,2))/(-1 + 2*t)); - Rule(GRAD(3,1), - -(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))*DMDA(1,3))/ - (Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(3,2),((-s + t)*DMDA(2,1))/(-1 + 2*t)); Rule(GRAD(3,2),((-r + t)*DMDA(2,2))/(-1 + 2*t)); - Rule(GRAD(3,2),-(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))* - DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(3,3),((-s + t)*DMDA(3,1))/(-1 + 2*t)); - Rule(GRAD(3,3),((-r + t)*DMDA(3,2))/(-1 + 2*t)); - Rule(GRAD(3,3),-(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))* - DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - - - -#undef List -#undef Power -#undef xi -#undef GRAD -#undef DMDA -#undef Rule - } - inline bool volume_matrix_tri_2D(double &detJ, DenseMatrix<3,3>& A, const double *x[3]) { return volume_matrix_2D(detJ, A, x); diff --git a/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.cpp b/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.cpp index b9502abb33b6..1ec3cb833f4e 100644 --- a/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.cpp +++ b/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.cpp @@ -8,15 +8,6 @@ #include -#if defined(STK_PERCEPT_LITE) -#if STK_PERCEPT_LITE -//d1 STK_PERCEPT_LITE -#endif -#else -//ud STK_PERCEPT_LITE -#endif - - #if !defined(NO_GEOM_SUPPORT) #include "JacobianUtil.hpp" @@ -34,38 +25,6 @@ {4, 7, 5, 0}, {5, 4, 6, 1}, {6, 5, 7, 2}, {7, 6, 4, 3}}; -#if 0 - static inline void set2d(JacobianUtil::Vec3D v, const double *x) - { - v[0] = x[0]; - v[1] = x[1]; - v[2] = 0.0; - } - static inline void set3d(JacobianUtil::Vec3D v, const double *x) - { - v[0] = x[0]; - v[1] = x[1]; - v[2] = x[2]; - } -#endif - - void scale_to_unit(DenseMatrix<3,3>& A) - { - for (int jvert=0; jvert < 3; jvert++) - { - double sum=0.0; - for (int ixyz=0; ixyz < 3; ixyz++) - { - sum += A(ixyz, jvert)*A(ixyz, jvert); - } - sum = std::max(1.e-10, std::sqrt(sum)); - for (int ixyz=0; ixyz < 3; ixyz++) - { - A(ixyz, jvert) /= sum; - } - } - } - template<> bool JacobianUtilImpl::jacobian_matrix_1D_in_2D(double &detJ, DenseMatrix<3,3>& A, const double *x[2]) { A.zero(); @@ -101,7 +60,6 @@ A(2,0) = 0; // (x[1][2] - x[0][2]); A(2,1) = 0; // (x[2][2] - x[0][2]); A(2,2) = 1.0; - //if (m_scale_to_unit) scale_to_unit(A); detJ = det(A); return detJ < 0.0; @@ -400,7 +358,6 @@ check_unhandled_topo(eMesh, topology_data); case shards::Pyramid<5>::key: { - bool err=false; for (i = 0; i < 5; ++i) { metric_valid = jacobian_matrix_pyramid_3D_new(i, m_detJ[i], m_J[i], @@ -409,21 +366,9 @@ VERTEX(v_i[2]), VERTEX(v_i[3]), VERTEX(v_i[4])); - if (m_detJ[i] < 1.e-12) err=true; } - // FIXME m = average_metrics(m_detJ, 5); - if (m < 1.e-12 || err) - { - std::cout << "pyramid detJ= " << m << std::endl; - for (i = 0; i < 5; ++i) { - std::cout << " detJ[" << i << "]= " << m_detJ[i] << std::endl; - } - for (i = 0; i < 5; ++i) { - std::cout << " J[" << i << "]= " << m_J[i] << std::endl; - } - } } break; diff --git a/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.hpp b/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.hpp index 9c12a53ef53c..68286da49763 100644 --- a/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.hpp +++ b/packages/percept/src/percept/mesh/mod/smoother/JacobianUtil.hpp @@ -36,12 +36,10 @@ DenseMatrix<3,3> m_dMetric_dA[NNODES_MAX]; double m_grad[NNODES_MAX][NNODES_MAX][3]; int m_num_nodes; - bool m_scale_to_unit; bool m_use_approximate_quadratic_jacobian; JacobianUtilBase(bool use_approximate_quadratic_jacobian=true) : m_num_nodes(0), - m_scale_to_unit(false), m_use_approximate_quadratic_jacobian(use_approximate_quadratic_jacobian) { } @@ -71,7 +69,6 @@ using Base::m_dMetric_dA; using Base::m_grad; using Base::m_num_nodes; - using Base::m_scale_to_unit; using Base::m_use_approximate_quadratic_jacobian; @@ -143,7 +140,6 @@ A(2,0) = (x1[2] - x0[2]); A(2,1) = (x2[2] - x0[2]); A(2,2) = (x3[2] - x0[2]); - //if (m_scale_to_unit) scale_to_unit(A); detJ = det(A); return detJ < 0.0; @@ -373,69 +369,6 @@ return detJ < 0.0; } - void grad_util_pyramid_3d_new(int ibasis, const DenseMatrix<3,3>& dMdA, double grad[NNODES_MAX][3], const int nnode, const int spd, const int *indices, const int nind) - { - double rst[5][3] = {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{.5,.5,.49999}}; //note: avoid singularity at top vertex - - double r = rst[ibasis][0]; - double s = rst[ibasis][1]; - double t = rst[ibasis][2]; - - for (int i=0; i < nnode; i++) - for (int j=0; j < spd; j++) - grad[i][j]=0.0; - -#define Rule(x,y) x += y -#define GRAD(i,j) grad[i-1][j-1] -#define DMDA(i,j) dMdA(i-1,j-1) -#define Sqrt(x) std::sqrt(x) - Rule(GRAD(1,1),-(((-1 + s + t)*DMDA(1,1))/(-1 + 2*t))); - Rule(GRAD(1,1),-(((-1 + r + t)*DMDA(1,2))/(-1 + 2*t))); - Rule(GRAD(1,1),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(1,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(1,2),-(((-1 + s + t)*DMDA(2,1))/(-1 + 2*t))); - Rule(GRAD(1,2),-(((-1 + r + t)*DMDA(2,2))/(-1 + 2*t))); - Rule(GRAD(1,2),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(1,3),-(((-1 + s + t)*DMDA(3,1))/(-1 + 2*t))); - Rule(GRAD(1,3),-(((-1 + r + t)*DMDA(3,2))/(-1 + 2*t))); - Rule(GRAD(1,3),((-2*(1 + Sqrt(2) - 2*t)*(-1 + t) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))* - DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2))); - Rule(GRAD(2,1),((-1 + s + t)*DMDA(1,1))/(-1 + 2*t)); - Rule(GRAD(2,1),((r - t)*DMDA(1,2))/(-1 + 2*t)); - Rule(GRAD(2,1),-(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(1,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(2,2),((-1 + s + t)*DMDA(2,1))/(-1 + 2*t)); - Rule(GRAD(2,2),((r - t)*DMDA(2,2))/(-1 + 2*t)); - Rule(GRAD(2,2), - -(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2)))); Rule(GRAD(2,3),((-1 + s + t)*DMDA(3,1))/(-1 + 2*t)); - Rule(GRAD(2,3),((r - t)*DMDA(3,2))/(-1 + 2*t)); - Rule(GRAD(2,3), - -(((1 + 3*Sqrt(2) - 2*t - 6*Sqrt(2)*t + 4*Sqrt(2)*Power(t,2) + s*(-1 - 3*Sqrt(2) + 2*(1 + Sqrt(2))*t) + - r*(-1 - 3*Sqrt(2) + 4*Sqrt(2)*s + 2*(1 + Sqrt(2))*t))*DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2)))); Rule(GRAD(3,1),((-s + t)*DMDA(1,1))/(-1 + 2*t)); - Rule(GRAD(3,1),((-r + t)*DMDA(1,2))/(-1 + 2*t)); - Rule(GRAD(3,1), - -(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))*DMDA(1,3))/ - (Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(3,2),((-s + t)*DMDA(2,1))/(-1 + 2*t)); Rule(GRAD(3,2),((-r + t)*DMDA(2,2))/(-1 + 2*t)); - Rule(GRAD(3,2),-(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))* - DMDA(2,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - Rule(GRAD(3,3),((-s + t)*DMDA(3,1))/(-1 + 2*t)); - Rule(GRAD(3,3),((-r + t)*DMDA(3,2))/(-1 + 2*t)); - Rule(GRAD(3,3),-(((s*(1 + 3*Sqrt(2) - 2*(1 + Sqrt(2))*t) + r*(1 + 3*Sqrt(2) - 4*Sqrt(2)*s - 2*(1 + Sqrt(2))*t) + 2*t*(-1 - 3*Sqrt(2) + (2 + 4*Sqrt(2))*t))* - DMDA(3,3))/(Sqrt(2)*Power(1 - 2*t,2)))); - - - -#undef List -#undef Power -#undef xi -#undef GRAD -#undef DMDA -#undef Rule - } - inline bool jacobian_matrix_tri_2D(double &detJ, DenseMatrix<3,3>& A, const double *x[3]) { diff --git a/packages/percept/src/percept/mesh/mod/smoother/SGridJacobianUtil.hpp b/packages/percept/src/percept/mesh/mod/smoother/SGridJacobianUtil.hpp index e9d555b22674..11779e8f650b 100644 --- a/packages/percept/src/percept/mesh/mod/smoother/SGridJacobianUtil.hpp +++ b/packages/percept/src/percept/mesh/mod/smoother/SGridJacobianUtil.hpp @@ -58,14 +58,11 @@ inline bool jacobian_matrix_3D(double &detJ, using Base::m_dMetric_dA; using Base::m_grad; using Base::m_num_nodes; - using Base::m_scale_to_unit; using Base::m_use_approximate_quadratic_jacobian; using MeshType = StructuredGrid; enum { NELEM_TYPES = 1 }; - //using Base::NELEM_TYPES; - JacobianUtilImpl(bool use_approximate_quadratic_jacobian=true) : Base(use_approximate_quadratic_jacobian) { diff --git a/packages/percept/src/percept/mesh/mod/smoother/SmootherMetric.hpp b/packages/percept/src/percept/mesh/mod/smoother/SmootherMetric.hpp index 795eff7408d1..f2ab73c17f85 100644 --- a/packages/percept/src/percept/mesh/mod/smoother/SmootherMetric.hpp +++ b/packages/percept/src/percept/mesh/mod/smoother/SmootherMetric.hpp @@ -1062,7 +1062,6 @@ { valid = true; JacobianUtil jacA, jacW; - //jacA.m_scale_to_unit = true; double A_ = 0.0, W_ = 0.0; // current and reference detJ jacA(A_, *m_eMesh, element, m_coord_field_current, m_topology_data); @@ -1110,7 +1109,6 @@ VERIFY_OP_ON(m_node, !=, stk::mesh::Entity(), "must set a node"); valid = true; JacobianUtil jacA, jacSA, jacW; - jacSA.m_scale_to_unit = true; double SA_ = 0.0, A_ = 0.0, W_ = 0.0; // current and reference detJ jacA(A_, *m_eMesh, element, m_coord_field_current, m_topology_data); @@ -1162,7 +1160,6 @@ { valid = true; JacobianUtil jacA, jacSA, jacW; - jacSA.m_scale_to_unit = true; double SA_ = 0.0, A_ = 0.0, W_ = 0.0; // current and reference detJ jacA(A_, *m_eMesh, element, m_coord_field_current, m_topology_data); @@ -2003,100 +2000,6 @@ }; -///////////////////////////////////// -//mesh specific implementations -///////////////////////////////////// -// template<> -// SmootherMetricUntangleImpl:: -// SmootherMetricUntangleImpl(PerceptMesh *eMesh) : SmootherMetricImpl(eMesh) { -// m_beta_mult = 0.05; -// } -// -// template<> -// SmootherMetricUntangleImpl:: -// SmootherMetricUntangleImpl(PerceptMesh *eMesh) : SmootherMetricImpl(eMesh) { -// std::shared_ptr bsg = eMesh->get_block_structured_grid(); -// m_beta_mult = 0.05; -// StructuredGrid::MTField *m_coord_field_current = bsg->m_fields["coordinates"].get(); -// StructuredGrid::MTField *m_coord_field_original = bsg->m_fields["coordinates_NM1"].get(); -// m_coords_current = *m_coord_field_current->m_block_fields[0]; -// m_coords_original = *m_coord_field_original->m_block_fields[0]; -// } -// -// // ETI -// template -// SmootherMetricUntangleImpl:: -// SmootherMetricUntangleImpl(PerceptMesh *eMesh); -// -// template -// SmootherMetricUntangleImpl:: -// SmootherMetricUntangleImpl(PerceptMesh *eMesh); -// -// -// template<> -// double SmootherMetricUntangleImpl:: -// metric(typename STKMesh::MTElement element, bool& valid) -// { -// valid = true; -// -// JacobianUtilImpl jacA, jacW; -// -// double A_ = 0.0, W_ = 0.0; // current and reference detJ -// jacA(A_, *Base::m_eMesh, element, Base::m_coord_field_current, Base::m_topology_data); -// jacW(W_, *Base::m_eMesh, element, Base::m_coord_field_original, Base::m_topology_data); -// double val_untangle=0.0; -// -// for (int i=0; i < jacA.m_num_nodes; i++) -// { -// double detAi = jacA.m_detJ[i]; -// double detWi = jacW.m_detJ[i]; -// -// if (detAi <= 0.) valid = false; -// -// double vv = m_beta_mult*detWi - detAi; -// vv = std::max(vv, 0.0); -// val_untangle += vv*vv; -// } -// return val_untangle; -// } -// -// template<> -// double SmootherMetricUntangleImpl:: -// metric(typename StructuredGrid::MTElement element, bool& valid) -// { -// valid = true; -// -// //SGridJacobianUtilImpl jacA, jacW; -// -// double A_ = 0.0, W_ = 0.0; // current and reference detJ -// double nodal_A[8], nodal_W[8]; -// SGridJacobianUtil(A_, nodal_A, m_coords_current, element); -// SGridJacobianUtil(W_, nodal_W, m_coords_original, element); -// double val_untangle=0.0; -// -// for (int i=0; i < 8; i++) -// { -// double detAi = nodal_A[i]; -// double detWi = nodal_W[i]; -// -// if (detAi <= 0.) valid = false; -// -// double vv = m_beta_mult*detWi - detAi; -// vv = std::max(vv, 0.0); -// val_untangle += vv*vv; -// } -// return val_untangle; -// } -// -// // ETI -// template -// double SmootherMetricUntangleImpl:: -// metric(typename STKMesh::MTElement element, bool& valid); -// -// template -// double SmootherMetricUntangleImpl:: -// metric(typename StructuredGrid::MTElement element, bool& valid); - }//percept #endif diff --git a/packages/percept/src/percept/mesh_transfer/RotationTranslation.hpp b/packages/percept/src/percept/mesh_transfer/RotationTranslation.hpp index 8b9a8263b9c5..fe3e2bed19a3 100644 --- a/packages/percept/src/percept/mesh_transfer/RotationTranslation.hpp +++ b/packages/percept/src/percept/mesh_transfer/RotationTranslation.hpp @@ -14,7 +14,7 @@ void applyRotation(stk::mesh::Field *vectorField, double&xrot, double&yrot, double&zrot) // angles in degrees { stk::mesh::BulkData& bulkdata = vectorField->get_mesh(); - const stk::mesh::MetaData & meta = stk::mesh::MetaData::get(bulkdata); + const stk::mesh::MetaData & meta = bulkdata.mesh_meta_data(); stk::mesh::EntityRank rank = vectorField->entity_rank(); // TODO support 2D? @@ -87,7 +87,7 @@ void applyTranslation(stk::mesh::Field *vectorFiel double&xtrans, double&ytrans, double&ztrans) { stk::mesh::BulkData& bulkdata = vectorField->get_mesh(); - const stk::mesh::MetaData & meta = stk::mesh::MetaData::get(bulkdata); + const stk::mesh::MetaData & meta = bulkdata.mesh_meta_data(); stk::mesh::EntityRank rank = vectorField->entity_rank(); stk::mesh::Selector select_used = diff --git a/packages/percept/src/percept/norm/H1Norm.hpp b/packages/percept/src/percept/norm/H1Norm.hpp index c382f7ca008c..19a2d58291bb 100644 --- a/packages/percept/src/percept/norm/H1Norm.hpp +++ b/packages/percept/src/percept/norm/H1Norm.hpp @@ -155,14 +155,15 @@ { // FIXME - make all percept code const-correct - PerceptMesh eMesh(&stk::mesh::MetaData::get(m_bulkData), &m_bulkData); + const stk::mesh::MetaData& meta = m_bulkData.mesh_meta_data(); + PerceptMesh eMesh(&meta, &m_bulkData); int spatialDim = eMesh.get_spatial_dim(); H1_NormOp H1_op(integrand, spatialDim); //CompositeFunction H1_of_integrand("H1_of_integrand", integrand, H1_op); IntegratedOp integrated_H1_op(H1_op, m_turboOpt); integrated_H1_op.setCubDegree(m_cubDegree); - const stk::mesh::Part& locally_owned_part = stk::mesh::MetaData::get(m_bulkData).locally_owned_part(); + const stk::mesh::Part& locally_owned_part = meta.locally_owned_part(); stk::mesh::Selector selector(*m_selector & locally_owned_part); //eMesh.print_info("Norm"); if (m_turboOpt == TURBO_NONE || m_turboOpt == TURBO_ELEMENT) diff --git a/packages/percept/src/percept/norm/IntrepidManager.cpp b/packages/percept/src/percept/norm/IntrepidManager.cpp index fd1e17f6106a..ad3514e592ea 100644 --- a/packages/percept/src/percept/norm/IntrepidManager.cpp +++ b/packages/percept/src/percept/norm/IntrepidManager.cpp @@ -561,7 +561,7 @@ found_it = 0; // FIXME consider caching the coords_field in FieldFunction - const stk::mesh::MetaData& metaData = stk::mesh::MetaData::get(bulkData); + const stk::mesh::MetaData& metaData = bulkData.mesh_meta_data(); CoordinatesFieldType *coords_field = metaData.get_field(stk::topology::NODE_RANK, "coordinates"); const stk::mesh::Bucket & bucket = bulkData.bucket(element); diff --git a/packages/percept/src/percept/norm/Norm.hpp b/packages/percept/src/percept/norm/Norm.hpp index 7e969e46b375..91fcae475a6e 100644 --- a/packages/percept/src/percept/norm/Norm.hpp +++ b/packages/percept/src/percept/norm/Norm.hpp @@ -129,7 +129,7 @@ Norm(stk::mesh::BulkData& bulkData, std::string partName, TurboOption turboOpt=TURBO_BUCKET, bool is_surface_norm=false) : FunctionOperator(bulkData, (stk::mesh::Part*)0), m_is_surface_norm(is_surface_norm), m_turboOpt(turboOpt), m_cubDegree(2), m_norm_field(0) { - stk::mesh::Part * part = stk::mesh::MetaData::get(bulkData).get_part(partName); + stk::mesh::Part * part = bulkData.mesh_meta_data().get_part(partName); if (!part) throw std::runtime_error(std::string("No part named ") +partName); init(part); error_check_is_surface_norm(); @@ -145,7 +145,7 @@ m_selector = new stk::mesh::Selector; for (int i = 0; i < partNames.dimension(0); i++) { - stk::mesh::Part * part = stk::mesh::MetaData::get(bulkData).get_part(partNames(i)); + stk::mesh::Part * part = bulkData.mesh_meta_data().get_part(partNames(i)); if (!part) throw std::runtime_error(std::string("No part named ") +partNames(i)); *m_selector = (*m_selector) | (*part); } @@ -161,7 +161,7 @@ stk::mesh::EntityRank rank = stk::topology::ELEMENT_RANK; if (m_is_surface_norm) { - rank = stk::mesh::MetaData::get(bulkData).side_rank(); + rank = bulkData.mesh_meta_data().side_rank(); } if (!m_selector) m_selector = new stk::mesh::Selector(); @@ -206,8 +206,8 @@ void error_check_is_surface_norm() { stk::mesh::EntityRank element_rank = stk::topology::ELEMENT_RANK; - stk::mesh::EntityRank side_rank = stk::mesh::MetaData::get(m_bulkData).side_rank(); - const stk::mesh::PartVector& parts = stk::mesh::MetaData::get(m_bulkData).get_parts(); + stk::mesh::EntityRank side_rank = m_bulkData.mesh_meta_data().side_rank(); + const stk::mesh::PartVector& parts = m_bulkData.mesh_meta_data().get_parts(); stk::mesh::EntityRank all_ranks = stk::topology::NODE_RANK; unsigned nparts = parts.size(); for (unsigned ipart=0; ipart < nparts; ipart++) @@ -278,7 +278,7 @@ integrated_LN_op.setAccumulationType(IntegratedOp::ACCUMULATE_MAX); } - const stk::mesh::Part& locally_owned_part = stk::mesh::MetaData::get(m_bulkData).locally_owned_part(); + const stk::mesh::Part& locally_owned_part = m_bulkData.mesh_meta_data().locally_owned_part(); stk::mesh::Selector selector(*m_selector & locally_owned_part); if (m_turboOpt == TURBO_NONE || m_turboOpt == TURBO_ELEMENT) { @@ -294,7 +294,7 @@ if (Power == -1) { - MaxOfNodeValues maxOfNodeValues(stk::mesh::MetaData::get(m_bulkData).spatial_dimension(), integrand); + MaxOfNodeValues maxOfNodeValues(m_bulkData.mesh_meta_data().spatial_dimension(), integrand); nodalOpLoop(m_bulkData, maxOfNodeValues, 0, &selector); for (unsigned iDim = 0; iDim < local.size(); iDim++) local[iDim] = std::max(local[iDim], maxOfNodeValues.maxVal[iDim]); diff --git a/packages/percept/src/percept/stk_rebalance_utils/RebalanceUtils.cpp b/packages/percept/src/percept/stk_rebalance_utils/RebalanceUtils.cpp index f86bc006b85c..aca208a44909 100644 --- a/packages/percept/src/percept/stk_rebalance_utils/RebalanceUtils.cpp +++ b/packages/percept/src/percept/stk_rebalance_utils/RebalanceUtils.cpp @@ -29,7 +29,7 @@ double stk::rebalance::check_balance(mesh::BulkData & bulk_data, const ParallelMachine &comm = bulk_data.parallel(); double my_load = 0.0; - const mesh::MetaData & meta_data = stk::mesh::MetaData::get(bulk_data); + const mesh::MetaData & meta_data = bulk_data.mesh_meta_data(); mesh::EntityVector local_elems; if (provided_elements)