diff --git a/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.cpp b/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.cpp index 1af33b41579e..2d61b7a8cc84 100644 --- a/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.cpp +++ b/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.cpp @@ -41,14 +41,14 @@ void filter_refinement_marker(const RefinementInterface & refinement, const stk: if (do_not_refine_or_unrefine_selector(*bucketPtr)) { for (int i = 0; i < size; ++i) - if (markers[i] == Refinement_Marker::REFINE || markers[i] == Refinement_Marker::COARSEN) - markers[i] = Refinement_Marker::NOTHING; + if (markers[i] == static_cast(Refinement_Marker::REFINE) || markers[i] == static_cast(Refinement_Marker::COARSEN)) + markers[i] = static_cast(Refinement_Marker::NOTHING); } else if (bucketPtr->member(parentPart)) { for (int i = 0; i < size; ++i) - if (markers[i] == Refinement_Marker::REFINE) - markers[i] = Refinement_Marker::NOTHING; + if (markers[i] == static_cast(Refinement_Marker::REFINE)) + markers[i] = static_cast(Refinement_Marker::NOTHING); } } } @@ -84,7 +84,7 @@ void perform_multilevel_adaptivity(RefinementInterface & refinement, int * markers = field_data(elem_marker, *b_ptr); for (size_t i = 0; i < b_ptr->size(); ++i) { - if (markers[i] == Refinement_Marker::REFINE) ++num_marked_refine; + if (markers[i] == static_cast(Refinement_Marker::REFINE)) ++num_marked_refine; } } diff --git a/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.hpp b/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.hpp index 74c4a91bb24c..4053e6e5d52f 100644 --- a/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.hpp +++ b/packages/krino/krino/krino_lib/Akri_AdaptivityHelpers.hpp @@ -21,12 +21,12 @@ class RefinementInterface; stk::mesh::Selector cdfem_do_not_refine_or_unrefine_selector(const CDFEM_Support & cdfem_support); -enum Refinement_Marker - { - COARSEN = -1, - NOTHING = 0, - REFINE = 1 - }; +enum class Refinement_Marker +{ + COARSEN = -1, + NOTHING = 0, + REFINE = 1 +}; void perform_multilevel_adaptivity(RefinementInterface & refinement, diff --git a/packages/krino/krino/krino_lib/Akri_AnalyticSurfaceInterfaceGeometry.cpp b/packages/krino/krino/krino_lib/Akri_AnalyticSurfaceInterfaceGeometry.cpp index dc4f768d8f90..86d1db640b31 100644 --- a/packages/krino/krino/krino_lib/Akri_AnalyticSurfaceInterfaceGeometry.cpp +++ b/packages/krino/krino/krino_lib/Akri_AnalyticSurfaceInterfaceGeometry.cpp @@ -165,8 +165,6 @@ static void append_surface_edge_intersection_points(const stk::mesh::BulkData & } } - - FieldRef AnalyticSurfaceInterfaceGeometry::get_coordinates_field(const stk::mesh::BulkData & mesh) const { FieldRef coordsField = myCdfemSupport.get_coords_field(); diff --git a/packages/krino/krino/krino_lib/Akri_CDMesh.cpp b/packages/krino/krino/krino_lib/Akri_CDMesh.cpp index 0a8576ad01a9..557a4d4e1d12 100644 --- a/packages/krino/krino/krino_lib/Akri_CDMesh.cpp +++ b/packages/krino/krino/krino_lib/Akri_CDMesh.cpp @@ -1033,46 +1033,6 @@ static void fill_nodes_of_elements_with_subelements_or_changed_phase(const stk:: } } -static -void pack_shared_nodes_for_sharing_procs(const stk::mesh::BulkData & mesh, - const std::set & nodes, - stk::CommSparse &commSparse) -{ - std::vector nodeSharedProcs; - stk::pack_and_communicate(commSparse,[&]() - { - for (auto node : nodes) - { - if (mesh.bucket(node).shared()) - { - mesh.comm_shared_procs(node, nodeSharedProcs); - for (int procId : nodeSharedProcs) - commSparse.send_buffer(procId).pack(mesh.identifier(node)); - } - } - }); -} - -static -void unpack_shared_nodes(const stk::mesh::BulkData & mesh, - std::set & nodes, - stk::CommSparse &commSparse) -{ - stk::unpack_communications(commSparse, [&](int procId) - { - stk::CommBuffer & buffer = commSparse.recv_buffer(procId); - - while ( buffer.remaining() ) - { - stk::mesh::EntityId nodeId; - commSparse.recv_buffer(procId).unpack(nodeId); - stk::mesh::Entity node = mesh.get_entity(stk::topology::NODE_RANK, nodeId); - STK_ThrowRequire(mesh.is_valid(node)); - nodes.insert(node); - } - }); -} - static std::set get_nodes_of_elements_with_subelements_or_have_changed_phase(const stk::mesh::BulkData & mesh, const Phase_Support & phaseSupport, const std::vector> & meshElements) @@ -1080,9 +1040,7 @@ static std::set get_nodes_of_elements_with_subelements_or_hav std::set nodesOfElements; fill_nodes_of_elements_with_subelements_or_changed_phase(mesh, phaseSupport, meshElements, nodesOfElements); - stk::CommSparse commSparse(mesh.parallel()); - pack_shared_nodes_for_sharing_procs(mesh, nodesOfElements, commSparse); - unpack_shared_nodes(mesh, nodesOfElements, commSparse); + communicate_shared_nodes_to_sharing_procs(mesh, nodesOfElements); return nodesOfElements; } @@ -1426,7 +1384,7 @@ static void find_nearest_matching_prolong_facet(const stk::mesh::BulkData & mesh const std::vector & requiredFields, const SubElementNode & targetNode, const ProlongationFacet *& nearestProlongFacet, - FacetDistanceQuery & nearestFacetQuery, + FacetDistanceQuery & nearestFacetQuery, bool & haveMissingRemoteProlongFacets) { const stk::math::Vector3d & targetCoordinates = targetNode.coordinates(); @@ -1450,7 +1408,7 @@ static void find_nearest_matching_prolong_facet(const stk::mesh::BulkData & mesh for (auto && prolong_facet : nearest_prolong_facets) { - FacetDistanceQuery facet_query(*prolong_facet->get_facet(), targetCoordinates); + FacetDistanceQuery facet_query(*prolong_facet->get_facet(), targetCoordinates); if (nearestFacetQuery.empty() || facet_query.distance_squared() < nearestFacetQuery.distance_squared()) { nearestProlongFacet = prolong_facet; @@ -1526,7 +1484,7 @@ CDMesh::find_prolongation_node(const SubElementNode & targetNode) const STK_ThrowRequire(need_facets_for_prolongation()); const ProlongationFacet * nearestProlongFacet = nullptr; - FacetDistanceQuery nearestFacetQuery; + FacetDistanceQuery nearestFacetQuery; find_nearest_matching_prolong_facet(stk_bulk(), my_phase_prolong_tree_map, requiredFields, targetNode, nearestProlongFacet, nearestFacetQuery, my_missing_remote_prolong_facets); @@ -1923,7 +1881,8 @@ CDMesh::snap_nearby_intersections_to_nodes(const InterfaceGeometry & interfaceGe { /* %TRACE[ON]% */ Trace trace__("krino::Mesh::snap_nearby_intersections_to_nodes(void)"); /* %TRACE% */ stk::diag::TimeBlock timer__(my_timer_snap); - snap_to_node(stk_bulk(), interfaceGeometry, get_snapper(), domainsAtNodes); + const stk::mesh::Selector parentElementSelector = get_cdfem_parent_element_selector(get_active_part(), my_cdfem_support, my_phase_support); + snap_to_node(stk_bulk(), parentElementSelector, interfaceGeometry, get_snapper(), domainsAtNodes); for (auto && entry : domainsAtNodes) { const SubElementNode * node = get_mesh_node(stk_bulk().identifier(entry.first)); diff --git a/packages/krino/krino/krino_lib/Akri_CDMesh_Refinement.cpp b/packages/krino/krino/krino_lib/Akri_CDMesh_Refinement.cpp index 452a424f4fbb..b062b83e66cc 100644 --- a/packages/krino/krino/krino_lib/Akri_CDMesh_Refinement.cpp +++ b/packages/krino/krino/krino_lib/Akri_CDMesh_Refinement.cpp @@ -40,11 +40,11 @@ namespace { if (refine_level >= interface_max_refine_level ) { - marker = Refinement_Marker::NOTHING; + marker = static_cast(Refinement_Marker::NOTHING); } else { - marker = Refinement_Marker::REFINE; + marker = static_cast(Refinement_Marker::REFINE); } } } @@ -232,7 +232,7 @@ mark_nearest_node_on_cut_edges(const stk::mesh::BulkData& mesh, int determine_refinement_marker(const bool isElementIndicated, const int refinementIterCount, const int interfaceMinRefineLevel, const int elementRefineLevel) { - int marker = Refinement_Marker::NOTHING; + auto marker = Refinement_Marker::NOTHING; const int targetRefineLevel = isElementIndicated ? interfaceMinRefineLevel : 0; if (elementRefineLevel < targetRefineLevel) { @@ -246,7 +246,7 @@ int determine_refinement_marker(const bool isElementIndicated, const int refinem // with refinement marker = Refinement_Marker::COARSEN; } - return marker; + return static_cast(marker); } void diff --git a/packages/krino/krino/krino_lib/Akri_Compute_Surface_Distance.cpp b/packages/krino/krino/krino_lib/Akri_Compute_Surface_Distance.cpp index 1ccb86d5531a..1a726a02171d 100644 --- a/packages/krino/krino/krino_lib/Akri_Compute_Surface_Distance.cpp +++ b/packages/krino/krino/krino_lib/Akri_Compute_Surface_Distance.cpp @@ -33,7 +33,7 @@ Compute_Surface_Distance::calculate( } static void compute_distance_to_facets(const stk::mesh::BulkData & mesh, - const MeshSurface & facet_list, + const FacetedSurfaceBase & facets, const stk::mesh::Field& coordinates, const stk::mesh::Field& distance, const stk::mesh::Selector & volume_selector, @@ -53,18 +53,18 @@ static void compute_distance_to_facets(const stk::mesh::BulkData & mesh, STK_ThrowAssert(&(dist[n]) != NULL); const stk::math::Vector3d xvec(coord+n*spatial_dimension, spatial_dimension); - dist[n] = facet_list.point_unsigned_distance(xvec, narrowBandSize, farFieldValue); + dist[n] = facets.point_unsigned_distance(xvec, narrowBandSize, farFieldValue); } } } -void print_facet_info(const MeshSurface & facet_list, stk::ParallelMachine parallel) +void print_facet_info(const FacetedSurfaceBase & facets, stk::ParallelMachine parallel) { constexpr int vec_size = 3; std::array local_sizes, global_min, global_max; - local_sizes[0] = facet_list.storage_size(); - local_sizes[1] = facet_list.size(); - local_sizes[2] = facet_list.size()+facet_list.nonlocal_size(); + local_sizes[0] = facets.storage_size(); + local_sizes[1] = facets.size(); + local_sizes[2] = facets.size()+facets.nonlocal_size(); stk::all_reduce_min( parallel, local_sizes.data(), global_min.data(), vec_size ); stk::all_reduce_max( parallel, local_sizes.data(), global_max.data(), vec_size ); @@ -75,7 +75,6 @@ void print_facet_info(const MeshSurface & facet_list, stk::ParallelMachine paral krinolog << " Memory usage (mb): min=" << global_min[0]/(1024.0*1024.0) << ", max=" << global_max[0]/(1024.0*1024.0) << stk::diag::dendl; } - void Compute_Surface_Distance::calculate( const stk::mesh::BulkData & mesh, @@ -83,7 +82,7 @@ Compute_Surface_Distance::calculate( const stk::mesh::Field& coordinates, const stk::mesh::Field& distance, const stk::mesh::Selector & volume_selector, - const stk::mesh::Selector & surface_selector, + const stk::mesh::Selector & surfaceSelector, const double narrowBandSize, const double farFieldValue) { /* %TRACE[ON]% */ Trace trace__("krino::Compute_Surface_Distance::compute_surface_distance(void)"); /* %TRACE% */ @@ -91,14 +90,15 @@ Compute_Surface_Distance::calculate( stk::diag::Timer timer( "Compute Surface Distance", parent_timer ); stk::diag::TimeBlock timer_(timer); - MeshSurface facet_list(mesh.mesh_meta_data(), coordinates, surface_selector, +1); + std::unique_ptr facets = build_mesh_surface(mesh.mesh_meta_data(), coordinates, surfaceSelector, +1); + FieldRef coordsField(coordinates); const BoundingBox nodeBbox = compute_nodal_bbox(mesh, volume_selector & stk::mesh::selectField(distance), coordsField); - facet_list.prepare_to_compute(0.0, nodeBbox, narrowBandSize); // Setup including communication of facets that are within this processors narrow band + facets->prepare_to_compute(0.0, nodeBbox, narrowBandSize); // Setup including communication of facets that are within this processors narrow band - print_facet_info(facet_list, mesh.parallel()); + print_facet_info(*facets, mesh.parallel()); - compute_distance_to_facets(mesh, facet_list, coordinates, distance, volume_selector, narrowBandSize, farFieldValue); + compute_distance_to_facets(mesh, *facets, coordinates, distance, volume_selector, narrowBandSize, farFieldValue); } } // namespace krino diff --git a/packages/krino/krino/krino_lib/Akri_ContourElement.cpp b/packages/krino/krino/krino_lib/Akri_ContourElement.cpp index 3c8248cb926a..06e1ca84b5e3 100644 --- a/packages/krino/krino/krino_lib/Akri_ContourElement.cpp +++ b/packages/krino/krino/krino_lib/Akri_ContourElement.cpp @@ -527,7 +527,7 @@ ContourElement::dump_subelement_details() const } void -ContourElement::build_subelement_facets( Faceted_Surface & facets ) +ContourElement::build_subelement_facets( FacetedSurfaceBase & facets ) { if (my_sign == 0) { diff --git a/packages/krino/krino/krino_lib/Akri_ContourElement.hpp b/packages/krino/krino/krino_lib/Akri_ContourElement.hpp index 5ec1fd1efe7b..5fb03b8a38c1 100644 --- a/packages/krino/krino/krino_lib/Akri_ContourElement.hpp +++ b/packages/krino/krino/krino_lib/Akri_ContourElement.hpp @@ -74,7 +74,7 @@ class ContourElement { void dump_subelement_structure( void ) const; void dump_subelement_details( void ) const; - void build_subelement_facets( Faceted_Surface & facets ); + void build_subelement_facets( FacetedSurfaceBase & facets ); int gather_intg_pts( const int intg_pt_sign, sierra::ArrayContainer & intg_pt_locations, diff --git a/packages/krino/krino/krino_lib/Akri_ContourSubElement.cpp b/packages/krino/krino/krino_lib/Akri_ContourSubElement.cpp index 295abbe11a46..c68125639dd1 100644 --- a/packages/krino/krino/krino_lib/Akri_ContourSubElement.cpp +++ b/packages/krino/krino/krino_lib/Akri_ContourSubElement.cpp @@ -43,7 +43,7 @@ ContourSubElement::ContourSubElement( const ContourElement *in_owner, } int -ContourSubElement::build_facets( Faceted_Surface & facets ) +ContourSubElement::build_facets( FacetedSurfaceBase & facets ) { int start_size = facets.size(); @@ -118,7 +118,7 @@ ContourSubElement::dump_details() const } int -ContourSubElement::side_facets( Faceted_Surface & facets, int side ) const +ContourSubElement::side_facets( FacetedSurfaceBase & facets, int side ) const { const std::string & owner_type = my_owner->dist_topology().name(); const std::string & sub_type = topology().name(); @@ -1094,7 +1094,7 @@ ContourSubElement_Tri_3::is_degenerate( const std::array & edge_node_ids, } int -ContourSubElement_Tri_3::side_facets( Faceted_Surface & facets, +ContourSubElement_Tri_3::side_facets( FacetedSurfaceBase & facets, int side ) const { STK_ThrowAssert( get_side_ids()[side] == -2 ); @@ -1105,15 +1105,9 @@ ContourSubElement_Tri_3::side_facets( Faceted_Surface & facets, const unsigned * const lnn = get_side_node_ordinals(topology(), side); if ( LevelSet::sign_change(0.0, (double) my_sign) ) - { - std::unique_ptr facet = std::make_unique( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[1]]) ); - facets.add( std::move(facet) ); - } + facets.emplace_back_2d( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[1]]) ); else - { - std::unique_ptr facet = std::make_unique( my_owner->coordinates(myCoords[lnn[1]]), my_owner->coordinates(myCoords[lnn[0]]) ); - facets.add( std::move(facet) ); - } + facets.emplace_back_2d( my_owner->coordinates(myCoords[lnn[1]]), my_owner->coordinates(myCoords[lnn[0]]) ); return( num_facets ); } @@ -2086,7 +2080,7 @@ ContourSubElement_Tet_4::is_degenerate( const std::array & edge_node_ids } int -ContourSubElement_Tet_4::side_facets( Faceted_Surface & facets, +ContourSubElement_Tet_4::side_facets( FacetedSurfaceBase & facets, int side ) const { STK_ThrowAssert( mySideIds[side] == -2 ); @@ -2097,15 +2091,9 @@ ContourSubElement_Tet_4::side_facets( Faceted_Surface & facets, const unsigned * const lnn = get_side_node_ordinals(topology(), side); if ( LevelSet::sign_change(0.0, (double) my_sign) ) - { - std::unique_ptr facet = std::make_unique( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[1]]), my_owner->coordinates(myCoords[lnn[2]]) ); - facets.add( std::move(facet) ); - } + facets.emplace_back_3d( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[1]]), my_owner->coordinates(myCoords[lnn[2]]) ); else - { - std::unique_ptr facet = std::make_unique( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[2]]), my_owner->coordinates(myCoords[lnn[1]]) ); - facets.add( std::move(facet) ); - } + facets.emplace_back_3d( my_owner->coordinates(myCoords[lnn[0]]), my_owner->coordinates(myCoords[lnn[2]]), my_owner->coordinates(myCoords[lnn[1]]) ); return( num_facets ); } diff --git a/packages/krino/krino/krino_lib/Akri_ContourSubElement.hpp b/packages/krino/krino/krino_lib/Akri_ContourSubElement.hpp index d6c00e17afed..18a43f6cbc4e 100644 --- a/packages/krino/krino/krino_lib/Akri_ContourSubElement.hpp +++ b/packages/krino/krino/krino_lib/Akri_ContourSubElement.hpp @@ -64,10 +64,10 @@ class ContourSubElement { sierra::ArrayContainer & determinants, int index = 0 ); - int build_facets( Faceted_Surface & facets ); + int build_facets( FacetedSurfaceBase & facets ); // default implementation - virtual int side_facets( Faceted_Surface & facets, int side ) const; + virtual int side_facets( FacetedSurfaceBase & facets, int side ) const; virtual double side_area( int side ) const; stk::topology topology() const { return get_master_element().get_topology(); } @@ -228,7 +228,7 @@ class ContourSubElement_Tri_3 : public ContourSubElementWithTopology create_interface_geometry(const stk::mesh::Me if (manager.get_bounding_surfaces().size() > 0) return create_bounding_surface_geometry(manager, activePart, cdfemSupport, phaseSupport); - return create_levelset_geometry(activePart, cdfemSupport, phaseSupport, Phase_Support::get_levelset_fields(meta)); + return create_levelset_geometry(meta.spatial_dimension(), activePart, cdfemSupport, phaseSupport, Phase_Support::get_levelset_fields(meta)); } std::unique_ptr create_bounding_surface_geometry(Surface_Manager & manager, @@ -43,7 +43,8 @@ std::unique_ptr create_bounding_surface_geometry(Surface_Mana return geom; } -std::unique_ptr create_levelset_geometry(const stk::mesh::Part & activePart, +std::unique_ptr create_levelset_geometry(const int dim, + const stk::mesh::Part & activePart, const CDFEM_Support & cdfemSupport, const Phase_Support & phaseSupport, const std::vector & LSFields) @@ -53,7 +54,7 @@ std::unique_ptr create_levelset_geometry(const stk::mesh::Par (RESNAP_USING_INTERFACE_ON_PREVIOUS_SNAPPED_MESH == cdfemSupport.get_resnap_method() || RESNAP_USING_INTERPOLATION == cdfemSupport.get_resnap_method())); if (facetsRequested && !phaseSupport.has_one_levelset_per_phase()) - return std::make_unique(activePart, cdfemSupport, phaseSupport, LSFields); + return std::make_unique(dim, activePart, cdfemSupport, phaseSupport, LSFields); return std::make_unique(activePart, cdfemSupport, phaseSupport, LSFields); } diff --git a/packages/krino/krino/krino_lib/Akri_CreateInterfaceGeometry.hpp b/packages/krino/krino/krino_lib/Akri_CreateInterfaceGeometry.hpp index 5a02c5ed248c..1691b7420469 100644 --- a/packages/krino/krino/krino_lib/Akri_CreateInterfaceGeometry.hpp +++ b/packages/krino/krino/krino_lib/Akri_CreateInterfaceGeometry.hpp @@ -26,7 +26,8 @@ std::unique_ptr create_bounding_surface_geometry(Surface_Mana const CDFEM_Support & cdfemSupport, const Phase_Support & phaseSupport); -std::unique_ptr create_levelset_geometry(const stk::mesh::Part & activePart, +std::unique_ptr create_levelset_geometry(const int dim, + const stk::mesh::Part & activePart, const CDFEM_Support & cdfemSupport, const Phase_Support & phaseSupport, const std::vector & LSFields); diff --git a/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.cpp b/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.cpp index cfdb7ecb178d..857e6a152b8a 100644 --- a/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.cpp +++ b/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.cpp @@ -109,11 +109,13 @@ std::array get_oriented_nodes_tetrahedron(const int nodeToUpdate) return lnn; } -static double eikonal_solve_triangle_2d_or_3d(const std::array & x, const std::array & d, const int dim, const double far, const double speed) +static void eikonal_solve_triangle_2d_or_3d_for_distance_and_edge_position(const std::array & x, const std::array & d, const int dim, const double far, const double speed, double & dist, double & edgePos) { if (d[0] == far || d[1] == far) { - return far; + dist = far; + edgePos = 0.; + return; } const double sqrSpeed = speed*speed; @@ -148,38 +150,66 @@ static double eikonal_solve_triangle_2d_or_3d(const std::array= 0.0 && loc <= 1.0) - return d[0] + d20; + if (edgePos >= 0.0 && edgePos <= 1.0) + { + dist = d[0] + d20; + return; + } } const double h20 = std::sqrt(sqr20); const double h21 = (x[2] - x[1]).length(); - return std::min(d[0]+h20/speed,d[1]+h21/speed); + const double distSide0 = d[0]+h20/speed; + const double distSide1 = d[1]+h21/speed; + if (distSide0 < distSide1) + { + dist = distSide0; + edgePos = 0.; + } + else + { + dist = distSide1; + edgePos = 1.; + } } double eikonal_solve_triangle(const std::array & x, const std::array & dSigned, const int sign, const double far, const double speed) { const std::array d{{ std::abs(dSigned[0]), std::abs(dSigned[1]) }}; static constexpr int dim = 2; - return sign*eikonal_solve_triangle_2d_or_3d(x, d, dim, far, speed); + double dist = far; + double edgePos = 0.; + eikonal_solve_triangle_2d_or_3d_for_distance_and_edge_position(x, d, dim, far, speed, dist, edgePos); + return sign*dist; +} + +std::pair eikonal_solve_triangle_for_distance_and_extension_speed(const std::array & x, const std::array & d, const std::array & extSpeed, const double far) +{ + static constexpr int dim = 2; + double nodeDist = far; + double edgePos = 0.; + eikonal_solve_triangle_2d_or_3d_for_distance_and_edge_position(x, d, dim, far, 1.0, nodeDist, edgePos); + const double nodeSpeed = (1.-edgePos)*extSpeed[0] + edgePos*extSpeed[1]; + return {nodeDist, nodeSpeed}; } -double eikonal_solve_triangle3d(const std::array & x, const std::array & dSigned, const double far, const double speed) +std::pair eikonal_solve_triangle_for_distance_and_extension_speed(const std::array & x, const std::array & dSigned, const std::array & extSpeed, const int sign, const double far) { const std::array d{{ std::abs(dSigned[0]), std::abs(dSigned[1]) }}; - static constexpr int dim = 3; - return eikonal_solve_triangle_2d_or_3d(x, d, dim, far, speed); + const auto [unsignedNodeDist, nodeSpeed] = eikonal_solve_triangle_for_distance_and_extension_speed(x, d, extSpeed, far); + return {sign*unsignedNodeDist, nodeSpeed}; } -double eikonal_solve_tetrahedron_side(const std::array & x, const std::array & d, const double far, const double speed, const int iSide) +void eikonal_solve_tetrahedron_side_for_distance_and_barycentric_coords(const std::array & x, const std::array & d, const double far, const double speed, const int iSide, double & dist, double & edgePos) { static const std::array,3> sides{{ {{0,1,3}}, {{1,2,3}}, {{2,0,3}} }}; const std::array & sideLNN = sides[iSide]; const std::array xSide{x[sideLNN[0]],x[sideLNN[1]],x[sideLNN[2]]}; const std::array dSide{d[sideLNN[0]],d[sideLNN[1]]}; - return eikonal_solve_triangle3d(xSide,dSide,far,speed); + static constexpr int dim = 3; + eikonal_solve_triangle_2d_or_3d_for_distance_and_edge_position(xSide, dSide, dim, far, speed, dist, edgePos); } double eikonal_solve_tetrahedron(const std::array & x, const std::array & dSigned, const int sign, const double far, const double speed) @@ -188,11 +218,37 @@ double eikonal_solve_tetrahedron(const std::array & x, co return sign*eikonal_solve_tetrahedron(x, d, far, speed); } -double eikonal_solve_tetrahedron(const std::array & x, const std::array & d, const double far, const double speed) +void eikonal_solve_tetrahedron_sides_for_distance_and_barycentric_coords(const std::array & x, const std::array & d, const double far, const double speed, double & dist, std::array & baseBarycentricCoords) +{ + std::array sideDist; + std::array sideEdgePos; + eikonal_solve_tetrahedron_side_for_distance_and_barycentric_coords(x, d, far, speed, 0, sideDist[0], sideEdgePos[0]); + eikonal_solve_tetrahedron_side_for_distance_and_barycentric_coords(x, d, far, speed, 1, sideDist[1], sideEdgePos[1]); + eikonal_solve_tetrahedron_side_for_distance_and_barycentric_coords(x, d, far, speed, 2, sideDist[2], sideEdgePos[2]); + const int minSide = (sideDist[0] < sideDist[1]) ? (sideDist[0] < sideDist[2] ? 0 : 2) : (sideDist[1] < sideDist[2] ? 1 : 2); + dist = sideDist[minSide]; + if (0 == minSide) + { + baseBarycentricCoords[0] = sideEdgePos[0]; baseBarycentricCoords[1] = 0.; + } + else if (1 == minSide) + { + baseBarycentricCoords[0] = 1.-sideEdgePos[1]; baseBarycentricCoords[1] = sideEdgePos[1]; + } + else + { + baseBarycentricCoords[0] = 0.; baseBarycentricCoords[1] = 1.-sideEdgePos[2]; + } + baseBarycentricCoords[2] = 1.0 - baseBarycentricCoords[0] - baseBarycentricCoords[1]; +} + +void eikonal_solve_tetrahedron_for_distance_and_barycentric_coords(const std::array & x, const std::array & d, const double far, const double speed, double & dist, std::array & baseBarycentricCoords) { if (d[0] == far || d[1] == far || d[2] == far) { - return far; + dist = far; + baseBarycentricCoords = {0.,0.,0.}; + return; } const double sqrSpeed = speed*speed; @@ -235,17 +291,42 @@ double eikonal_solve_tetrahedron(const std::array & x, co const double B2 = x20.length_squared() - sqrSpeed*d20*d20; const double C2 = Dot(x20,x30) - sqrSpeed*d20*d30; const double denom = A2*B1 - A1*B2; - const double loc_r = (C2*B1 - C1*B2)/denom; - const double loc_s = (A2*C1 - A1*C2)/denom; - const double loc_t = 1.0 - loc_r - loc_s; - - if (loc_r >= 0.0 && loc_s >= 0.0 && loc_t >= 0.0) - return d[0] + d30; + baseBarycentricCoords[0] = (C2*B1 - C1*B2)/denom; + baseBarycentricCoords[1] = (A2*C1 - A1*C2)/denom; + baseBarycentricCoords[2] = 1.0 - baseBarycentricCoords[0] - baseBarycentricCoords[1]; + + if (baseBarycentricCoords[0] >= 0.0 && baseBarycentricCoords[1] >= 0.0 && baseBarycentricCoords[2] >= 0.0) + { + dist = d[0] + d30; + return; + } } - return std::min({ - eikonal_solve_tetrahedron_side(x,d,far,speed,0), - eikonal_solve_tetrahedron_side(x,d,far,speed,1), - eikonal_solve_tetrahedron_side(x,d,far,speed,2)}); + eikonal_solve_tetrahedron_sides_for_distance_and_barycentric_coords(x, d, far, speed, dist, baseBarycentricCoords); +} + +double eikonal_solve_tetrahedron(const std::array & x, const std::array & d, const double far, const double speed) +{ + double dist = far; + std::array baseBarycentricCoords; + eikonal_solve_tetrahedron_for_distance_and_barycentric_coords(x, d, far, speed, dist, baseBarycentricCoords); + return dist; } + +std::pair eikonal_solve_tetrahedron_for_distance_and_extension_speed(const std::array & x, const std::array & d, const std::array & extSpeed, const double far) +{ + double nodeDist = far; + std::array baseBarycentricCoords; + eikonal_solve_tetrahedron_for_distance_and_barycentric_coords(x, d, far, 1.0, nodeDist, baseBarycentricCoords); + const double nodeSpeed = baseBarycentricCoords[2]*extSpeed[0] + baseBarycentricCoords[0]*extSpeed[1] + + baseBarycentricCoords[1]*extSpeed[2]; + return {nodeDist, nodeSpeed}; +} + +std::pair eikonal_solve_tetrahedron_for_distance_and_extension_speed(const std::array & x, const std::array & dSigned, const std::array & extSpeed, const int sign, const double far) +{ + const std::array d{{ std::abs(dSigned[0]), std::abs(dSigned[1]), std::abs(dSigned[2]) }}; + const auto [unsignedNodeDist, nodeSpeed] = eikonal_solve_tetrahedron_for_distance_and_extension_speed(x, d, extSpeed, far); + return {sign*unsignedNodeDist, nodeSpeed}; +} + } diff --git a/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.hpp b/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.hpp index eabd27de6fc7..7ed52044098c 100644 --- a/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.hpp +++ b/packages/krino/krino/krino_lib/Akri_Eikonal_Calc.hpp @@ -40,7 +40,7 @@ double eikonal_solve_triangle(const std::array & x, const double speed); double eikonal_solve_tetrahedron(const std::array & x, - const std::array & d, + const std::array & dSigned, const int sign, const double far, const double speed); @@ -50,6 +50,28 @@ double eikonal_solve_tetrahedron(const std::array & x, const double far, const double speed); +std::pair eikonal_solve_triangle_for_distance_and_extension_speed(const std::array & x, + const std::array & d, + const std::array & extSpeed, + const double far); + +std::pair eikonal_solve_triangle_for_distance_and_extension_speed(const std::array & x, + const std::array & dSigned, + const std::array & extSpeed, + const int sign, + const double far); + +std::pair eikonal_solve_tetrahedron_for_distance_and_extension_speed(const std::array & x, + const std::array & d, + const std::array & extSpeed, + const double far); + +std::pair eikonal_solve_tetrahedron_for_distance_and_extension_speed(const std::array & x, + const std::array & dSigned, + const std::array & extSpeed, + const int sign, + const double far); + std::array get_oriented_nodes_triangle(const int nodeToUpdate); std::array get_oriented_nodes_tetrahedron(const int nodeToUpdate); diff --git a/packages/krino/krino/krino_lib/Akri_Fast_Marching.cpp b/packages/krino/krino/krino_lib/Akri_Fast_Marching.cpp index 65720a59b477..f54501057bc5 100644 --- a/packages/krino/krino/krino_lib/Akri_Fast_Marching.cpp +++ b/packages/krino/krino/krino_lib/Akri_Fast_Marching.cpp @@ -12,25 +12,33 @@ #include #include -#include #include #include #include +#include #include #include "Akri_Eikonal_Calc.hpp" namespace krino{ -Fast_Marching::Fast_Marching(LevelSet & ls, const stk::mesh::Selector & selector, stk::diag::Timer & parent_timer) - : my_ls(ls), - my_selector(selector), - my_timer("Fast Marching", parent_timer), - my_fm_node_less(ls.mesh()), +Fast_Marching::Fast_Marching(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & activeElementSelector, + const FieldRef& coordinates, + const FieldRef& distance, + const std::function & get_interface_speed, + stk::diag::Timer & parentTimer) + : myMesh(mesh), + mySelector(activeElementSelector), + myCoordinates(coordinates), + myDistance(distance), + my_get_interface_speed(get_interface_speed), + myTimer("Fast Marching", parentTimer), + my_fm_node_less(mesh), trial_nodes(my_fm_node_less) { - ParallelErrorMessage err(mesh().parallel()); - stk::mesh::Selector fieldSelector(my_ls.get_distance_field()); - const stk::mesh::BucketVector& buckets = selector.get_buckets(stk::topology::ELEMENT_RANK); + ParallelErrorMessage err(mesh.parallel()); + stk::mesh::Selector fieldSelector(myDistance); + const stk::mesh::BucketVector& buckets = activeElementSelector.get_buckets(stk::topology::ELEMENT_RANK); for (auto && bucket : buckets) { if (fieldSelector(*bucket) && @@ -43,6 +51,21 @@ Fast_Marching::Fast_Marching(LevelSet & ls, const stk::mesh::Selector & selector check_error(err, "Checking topology"); } +double & Fast_Marching::get_node_distance(const stk::mesh::Entity node) +{ + return get_scalar_field(myMesh, myDistance, node); +} + +double Fast_Marching::get_node_distance(const stk::mesh::Entity node) const +{ + return get_scalar_field(myMesh, myDistance, node); +} + +double Fast_Marching::get_element_interface_speed(ParallelErrorMessage& err, const stk::mesh::Entity elem) const +{ + return my_get_interface_speed ? my_get_interface_speed(err, elem) : 1.0; +} + void Fast_Marching::check_error(const ParallelErrorMessage& err, const std::string & context) const { @@ -55,21 +78,104 @@ Fast_Marching::check_error(const ParallelErrorMessage& err, const std::string & STK_ThrowRequireMsg(!globalError.first, "Error in " << context << "."); } +void Fast_Marching::redistance(const std::vector & initialNodes) +{ + stk::diag::TimeBlock timer__(myTimer); + + initialize_algorithm(); + initialize_given_nodes(initialNodes); + propagate_distance_to_convergence(); +} + +void Fast_Marching::initialize_given_nodes(const std::vector & initialNodes) +{ + for (auto initialNode : initialNodes) + { + Fast_Marching_Node * fm_node = get_fm_node(initialNode); + STK_ThrowAssert(nullptr != fm_node && fm_node->status() != STATUS_UNUSED); + const double nodeDist = get_node_distance(initialNode); + fm_node->set_signed_dist(nodeDist); + fm_node->set_status(STATUS_INITIAL); + } +} + void Fast_Marching::redistance() { - stk::diag::TimeBlock timer__(my_timer); - const FieldRef& dRef = my_ls.get_distance_field(); - const FieldRef& olddRef = my_ls.get_old_distance_field(); + stk::diag::TimeBlock timer__(myTimer); + + initialize_algorithm(); + initialize_nodes_of_crossed_elements(); + propagate_distance_to_convergence(); +} +void Fast_Marching::initialize_algorithm() +{ // make sure field is parallel consistent to start out { - std::vector parallel_fields(1, &dRef.field()); - stk::mesh::copy_owned_to_shared(mesh(), parallel_fields); + std::vector parallel_fields(1, &myDistance.field()); + stk::mesh::copy_owned_to_shared(myMesh, parallel_fields); } - ParallelErrorMessage err(mesh().parallel()); + fm_nodes.clear(); + fm_nodes.resize(stk::mesh::count_selected_entities(myMesh.mesh_meta_data().universal_part(), myMesh.buckets(stk::topology::NODE_RANK))); + + std::vector field_nodes; + stk::mesh::get_selected_entities( selected_with_field_not_ghost_selector(), myMesh.buckets( stk::topology::NODE_RANK ), field_nodes ); + + for ( auto&& node : field_nodes ) + initialize_node(node); +} + +void Fast_Marching::initialize_nodes_of_crossed_elements() +{ + ParallelErrorMessage err(myMesh.parallel()); + + std::vector field_elems; + stk::mesh::get_selected_entities( selected_with_field_not_ghost_selector(), myMesh.buckets( stk::topology::ELEMENT_RANK ), field_elems ); + for (auto&& elem : field_elems) + if (have_crossing(elem)) + initialize_element(elem, get_element_interface_speed(err, elem)); + + if (myMesh.parallel_size() > 1) + { + stk::mesh::Selector globally_shared_selector = selected_with_field_selector() & myMesh.mesh_meta_data().globally_shared_part(); + std::vector< stk::mesh::Entity> shared_nodes; + stk::mesh::get_selected_entities( globally_shared_selector, myMesh.buckets( stk::topology::NODE_RANK ), shared_nodes ); + + for ( auto && shared_node : shared_nodes ) + { + Fast_Marching_Node * fm_node = get_fm_node(shared_node); + if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) + { + double & node_dist = get_node_distance(fm_node->node()); + node_dist = fm_node->signed_dist()*fm_node->sign(); + } + } + stk::mesh::parallel_min(myMesh, {&myDistance.field()}); + + for ( auto && shared_node : shared_nodes ) + { + Fast_Marching_Node * fm_node = get_fm_node(shared_node); + if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) + { + stk::mesh::Entity node = fm_node->node(); + const double min_node_unsigned_dist = get_node_distance(node); + const double fm_node_unsigned_dist = fm_node->signed_dist()*fm_node->sign(); + fm_node->set_signed_dist(min_node_unsigned_dist*fm_node->sign()); + if (min_node_unsigned_dist < fm_node_unsigned_dist) + { + fm_node->set_status(STATUS_INITIAL); + } + } + } + } + + check_error(err, "Fast Marching Initialization"); +} - initialize(err); +void Fast_Marching::propagate_distance_to_convergence() +{ + ParallelErrorMessage err(myMesh.parallel()); // neighbors of initial nodes are trial nodes for ( auto && fm_node : fm_nodes ) @@ -80,17 +186,15 @@ void Fast_Marching::redistance() } } - check_error(err, "Fast Marching Initialization"); - - stk::mesh::Selector globally_shared_selector = my_selector & mesh().mesh_meta_data().globally_shared_part(); + stk::mesh::Selector globally_shared_selector = selected_with_field_selector() & myMesh.mesh_meta_data().globally_shared_part(); std::vector< stk::mesh::Entity> shared_nodes; - stk::mesh::get_selected_entities( globally_shared_selector, mesh().buckets( stk::topology::NODE_RANK ), shared_nodes ); + stk::mesh::get_selected_entities( globally_shared_selector, myMesh.buckets( stk::topology::NODE_RANK ), shared_nodes ); bool done = false; const size_t local_num_nodes = fm_nodes.size(); size_t max_num_nodes = 0; - stk::all_reduce_max(mesh().parallel(), &local_num_nodes, &max_num_nodes, 1); + stk::all_reduce_max(myMesh.parallel(), &local_num_nodes, &max_num_nodes, 1); const unsigned max_outer_steps = max_num_nodes; // should be overkill unsigned num_outer_steps = 0; @@ -123,18 +227,18 @@ void Fast_Marching::redistance() check_error(err, "Fast Marching Iteration"); unsigned num_locally_updated = 0; - if (mesh().parallel_size() > 1) + if (myMesh.parallel_size() > 1) { for ( auto && shared_node : shared_nodes ) { Fast_Marching_Node * fm_node = get_fm_node(shared_node); if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) { - double & node_dist = *field_data(dRef, fm_node->node()); + double & node_dist = get_node_distance(fm_node->node()); node_dist = fm_node->signed_dist()*fm_node->sign(); } } - stk::mesh::parallel_min(mesh(), {&dRef.field()}); + stk::mesh::parallel_min(myMesh, {&myDistance.field()}); for ( auto && shared_node : shared_nodes ) { @@ -142,7 +246,7 @@ void Fast_Marching::redistance() if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) { stk::mesh::Entity node = fm_node->node(); - const double min_node_unsigned_dist = *field_data(dRef, node); + const double min_node_unsigned_dist = get_node_distance(node); const double fm_node_unsigned_dist = fm_node->signed_dist()*fm_node->sign(); if(min_node_unsigned_dist < fm_node_unsigned_dist) { @@ -156,7 +260,7 @@ void Fast_Marching::redistance() } unsigned num_globally_updated = 0; - stk::all_reduce_sum(mesh().parallel(), &num_locally_updated, &num_globally_updated, 1); + stk::all_reduce_sum(myMesh.parallel(), &num_locally_updated, &num_globally_updated, 1); done = (num_globally_updated == 0); } @@ -165,24 +269,23 @@ void Fast_Marching::redistance() { if (fm_node.status() == STATUS_TRIAL || fm_node.status() == STATUS_FAR) { - err << "Node " << mesh().identifier(fm_node.node()) << " with status " << fm_node.status() << " with distance " << fm_node.signed_dist() << " did not get updated!\n"; + err << "Node " << myMesh.identifier(fm_node.node()) << " with status " << fm_node.status() << " with distance " << fm_node.signed_dist() << " did not get updated!\n"; } if (fm_node.status() != STATUS_UNUSED) { - double & node_dist = *field_data(dRef, fm_node.node()); + double & node_dist = get_node_distance(fm_node.node()); node_dist = fm_node.signed_dist(); } } check_error(err, "Fast Marching Update"); - stk::mesh::field_copy(dRef, olddRef); } Fast_Marching_Node * Fast_Marching::get_fm_node(stk::mesh::Entity node) { - if (mesh().is_valid(node) && mesh().local_id(node) < fm_nodes.size()) + if (myMesh.is_valid(node) && myMesh.local_id(node) < fm_nodes.size()) { - return &fm_nodes[mesh().local_id(node)]; + return &fm_nodes[myMesh.local_id(node)]; } else { @@ -190,102 +293,40 @@ Fast_Marching_Node * Fast_Marching::get_fm_node(stk::mesh::Entity node) } } -void Fast_Marching::initialize(ParallelErrorMessage& err) +stk::mesh::Selector Fast_Marching::selected_with_field_not_ghost_selector() const { - stk::mesh::BulkData& stk_mesh = mesh(); - const FieldRef xRef = my_ls.get_coordinates_field(); - const FieldRef dRef = my_ls.get_distance_field(); - - const int dim = mesh().mesh_meta_data().spatial_dimension(); - - fm_nodes.clear(); - fm_nodes.resize(stk::mesh::count_selected_entities(stk_mesh.mesh_meta_data().universal_part(), stk_mesh.buckets(stk::topology::NODE_RANK))); - - std::vector field_nodes; - stk::mesh::Selector field_not_ghost = aux_meta().active_not_ghost_selector() & my_selector & stk::mesh::selectField(dRef); - stk::mesh::get_selected_entities( field_not_ghost, stk_mesh.buckets( stk::topology::NODE_RANK ), field_nodes ); - - for ( auto&& node : field_nodes ) - { - const double * curr_node_dist = field_data(dRef, node); - STK_ThrowAssert(nullptr != curr_node_dist); - - stk::math::Vector3d coords = stk::math::Vector3d::ZERO; - const double * xptr = field_data(xRef, node); - STK_ThrowAssert(nullptr != xptr); - for ( int d = 0; d < dim; d++ ) - { - coords[d] = xptr[d]; - } - - Fast_Marching_Node * fm_node = get_fm_node(node); - STK_ThrowAssert(nullptr != fm_node); - *fm_node = Fast_Marching_Node(node,STATUS_FAR,LevelSet::sign(*curr_node_dist) * std::numeric_limits::max(),LevelSet::sign(*curr_node_dist),coords); - } - - { - std::vector field_elems; - stk::mesh::get_selected_entities( field_not_ghost, stk_mesh.buckets( stk::topology::ELEMENT_RANK ), field_elems ); - for (auto&& elem : field_elems) - { - if (have_crossing(elem)) - { - const double speed = my_ls.get_time_of_arrival_speed(elem, err); - initialize_element(elem, speed); - } - } + stk::mesh::Selector selector = mySelector & stk::mesh::selectField(myDistance) & (myMesh.mesh_meta_data().locally_owned_part() | myMesh.mesh_meta_data().globally_shared_part()); + return selector; +} - if (mesh().parallel_size() > 1) - { - stk::mesh::Selector globally_shared_selector = my_selector & mesh().mesh_meta_data().globally_shared_part(); - std::vector< stk::mesh::Entity> shared_nodes; - stk::mesh::get_selected_entities( globally_shared_selector, stk_mesh.buckets( stk::topology::NODE_RANK ), shared_nodes ); +stk::mesh::Selector Fast_Marching::selected_with_field_selector() const +{ + stk::mesh::Selector selector = mySelector & stk::mesh::selectField(myDistance); + return selector; +} - for ( auto && shared_node : shared_nodes ) - { - Fast_Marching_Node * fm_node = get_fm_node(shared_node); - if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) - { - double & node_dist = *field_data(dRef, fm_node->node()); - node_dist = fm_node->signed_dist()*fm_node->sign(); - } - } - stk::mesh::parallel_min(mesh(), {&dRef.field()}); +void Fast_Marching::initialize_node(const stk::mesh::Entity node) +{ + const int currentSign = sign(get_node_distance(node)); + const stk::math::Vector3d coords = get_vector_field(myMesh, myCoordinates, node, myMesh.mesh_meta_data().spatial_dimension()); - for ( auto && shared_node : shared_nodes ) - { - Fast_Marching_Node * fm_node = get_fm_node(shared_node); - if (nullptr != fm_node && fm_node->status() != STATUS_UNUSED) - { - stk::mesh::Entity node = fm_node->node(); - const double min_node_unsigned_dist = *field_data(dRef, node); - const double fm_node_unsigned_dist = fm_node->signed_dist()*fm_node->sign(); - fm_node->set_signed_dist(min_node_unsigned_dist*fm_node->sign()); - if (min_node_unsigned_dist < fm_node_unsigned_dist) - { - fm_node->set_status(STATUS_INITIAL); - } - } - } - } - } + Fast_Marching_Node * fm_node = get_fm_node(node); + STK_ThrowAssert(nullptr != fm_node); + *fm_node = Fast_Marching_Node(node,STATUS_FAR, currentSign*std::numeric_limits::max(), currentSign, coords); } bool Fast_Marching::have_crossing(const stk::mesh::Entity & elem) const { - const FieldRef dRef = my_ls.get_distance_field(); - const unsigned npe = mesh().bucket(elem).topology().num_nodes(); + const unsigned npe = myMesh.bucket(elem).topology().num_nodes(); STK_ThrowAssert(npe > 0); - const stk::mesh::Entity * elem_nodes = mesh().begin(elem, stk::topology::NODE_RANK); - const double * dist0 = field_data(dRef, elem_nodes[0]); - STK_ThrowAssert(nullptr != dist0); + const stk::mesh::Entity * elem_nodes = myMesh.begin(elem, stk::topology::NODE_RANK); + const double dist0 = get_node_distance(elem_nodes[0]); for (unsigned n=1; n(dRef, elem_nodes[n]); - STK_ThrowAssert(nullptr != dist); - if (LevelSet::sign_change(*dist0, *dist)) + const double dist = get_node_distance(elem_nodes[n]); + if (sign_change(dist0, dist)) { return true; } @@ -326,19 +367,18 @@ Fast_Marching::initialize_element(const stk::mesh::Entity & elem, const double s // Initialize using method #5 (element rescaling) - const FieldRef dRef = my_ls.get_distance_field(); - const stk::mesh::Entity * elem_nodes = mesh().begin(elem, stk::topology::NODE_RANK); - const int npe = mesh().bucket(elem).topology().num_nodes(); + const stk::mesh::Entity * elem_nodes = myMesh.begin(elem, stk::topology::NODE_RANK); + const int npe = myMesh.bucket(elem).topology().num_nodes(); auto get_coordinates = build_get_fm_node_coordinates(this); - const double mag_grad = calculate_gradient_magnitude(npe, elem_nodes, dRef, get_coordinates); + const double mag_grad = calculate_gradient_magnitude(npe, elem_nodes, myDistance, get_coordinates); for (int inode=0; inodestatus() != STATUS_UNUSED); - const double elem_node_dist = *field_data(dRef, elem_nodes[inode]) / (mag_grad * speed); + const double elem_node_dist = get_node_distance(elem_nodes[inode]) / (mag_grad * speed); const int sign = fm_node->sign(); fm_node->set_signed_dist(sign * std::min(std::abs(fm_node->signed_dist()), std::abs(elem_node_dist))); fm_node->set_status(STATUS_INITIAL); @@ -348,34 +388,33 @@ Fast_Marching::initialize_element(const stk::mesh::Entity & elem, const double s void Fast_Marching::update_neighbors(Fast_Marching_Node & accepted_node, ParallelErrorMessage& err) { - const FieldRef dRef = my_ls.get_distance_field(); - const stk::mesh::Selector active_field_not_aura = my_selector & aux_meta().active_not_ghost_selector() & stk::mesh::selectField(dRef); + const stk::mesh::Selector elemSelector = selected_with_field_not_ghost_selector(); stk::mesh::Entity node = accepted_node.node(); STK_ThrowAssertMsg(STATUS_ACCEPTED == accepted_node.status() || STATUS_INITIAL == accepted_node.status(), "Expected ACCEPTED OR INITIAL status"); - const int dim = mesh().mesh_meta_data().spatial_dimension(); + const int dim = myMesh.mesh_meta_data().spatial_dimension(); STK_ThrowAssert(2 == dim || 3 == dim); const int npe_dist = (2==dim) ? 3 : 4; std::vector elem_nodes(npe_dist); - const unsigned num_node_elems = mesh().num_elements(node); - const stk::mesh::Entity* node_elems = mesh().begin_elements(node); + const unsigned num_node_elems = myMesh.num_elements(node); + const stk::mesh::Entity* node_elems = myMesh.begin_elements(node); for (unsigned node_elem_index=0; node_elem_index #include #include +#include namespace krino { @@ -65,9 +65,15 @@ class Fast_Marching_Node_Less class Fast_Marching { public: - Fast_Marching(LevelSet & ls, const stk::mesh::Selector & selector, stk::diag::Timer & parent_timer); + Fast_Marching(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & activeElementSelector, + const FieldRef& coordinates, + const FieldRef& distance, + const std::function & get_interface_speed, + stk::diag::Timer & parentTimer); void redistance(); + void redistance(const std::vector & initialNodes); void initialize(ParallelErrorMessage& err); void update_neighbors(Fast_Marching_Node & accepted_node, ParallelErrorMessage & err); void update_node(std::vector & elem_nodes, int node_to_update, const double speed); @@ -82,14 +88,26 @@ class Fast_Marching { void update_trial_node(Fast_Marching_Node & add_trial_node, const double dist); Fast_Marching_Node * get_fm_node(stk::mesh::Entity node); - const AuxMetaData& aux_meta() const { return my_ls.aux_meta(); } - AuxMetaData& aux_meta() { return my_ls.aux_meta(); } - const stk::mesh::BulkData& mesh() const { return my_ls.mesh(); } - stk::mesh::BulkData& mesh() { return my_ls.mesh(); } private: - LevelSet & my_ls; - stk::mesh::Selector my_selector; - stk::diag::Timer my_timer; + void initialize_algorithm(); + void initialize_nodes_of_crossed_elements(); + void initialize_given_nodes(const std::vector & initialNodes); + void propagate_distance_to_convergence(); + + void initialize_node(const stk::mesh::Entity node); + double & get_node_distance(const stk::mesh::Entity node); + double get_node_distance(const stk::mesh::Entity node) const; + double get_element_interface_speed(ParallelErrorMessage& err, const stk::mesh::Entity elem) const; + stk::mesh::Selector selected_with_field_not_ghost_selector() const; + stk::mesh::Selector selected_with_field_selector() const; + + const stk::mesh::BulkData & myMesh; + stk::mesh::Selector mySelector; + const FieldRef myCoordinates; + const FieldRef myDistance; + const std::function my_get_interface_speed; + stk::diag::Timer myTimer; + std::vector fm_nodes; Fast_Marching_Node_Less my_fm_node_less; std::set trial_nodes; //set sorted by distance, then by id to break "ties" diff --git a/packages/krino/krino/krino_lib/Akri_Intersection_Points.cpp b/packages/krino/krino/krino_lib/Akri_Intersection_Points.cpp index 23b6b9673ff8..5381c3c4b485 100644 --- a/packages/krino/krino/krino_lib/Akri_Intersection_Points.cpp +++ b/packages/krino/krino/krino_lib/Akri_Intersection_Points.cpp @@ -176,12 +176,13 @@ static void communicate_all_intersection_points(const stk::mesh::BulkData & mesh } static std::vector build_intersection_points(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const NodeToCapturedDomainsMap & nodesToCapturedDomains, const IntersectionPointFilter & intersectionPointFilter) { std::vector elementsToIntersect; - stk::mesh::get_entities( mesh, stk::topology::ELEMENT_RANK, mesh.mesh_meta_data().locally_owned_part(), elementsToIntersect, false); + stk::mesh::get_entities( mesh, stk::topology::ELEMENT_RANK, mesh.mesh_meta_data().locally_owned_part() & elementSelector, elementsToIntersect, false); std::vector intersectionPoints; geometry.append_element_intersection_points(mesh, nodesToCapturedDomains, elementsToIntersect, intersectionPointFilter, intersectionPoints); @@ -191,11 +192,12 @@ static std::vector build_intersection_points(const stk::mesh: } std::vector build_all_intersection_points(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const NodeToCapturedDomainsMap & nodesToCapturedDomains) { const IntersectionPointFilter intersectionPointFilter = keep_all_intersection_points_filter(); - return build_intersection_points(mesh, geometry, nodesToCapturedDomains, intersectionPointFilter); + return build_intersection_points(mesh, elementSelector, geometry, nodesToCapturedDomains, intersectionPointFilter); } static IntersectionPointFilter @@ -210,12 +212,13 @@ filter_intersection_points_to_those_not_already_handled(const NodeToCapturedDoma } std::vector build_uncaptured_intersection_points(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const NodeToCapturedDomainsMap & nodesToCapturedDomains) { IntersectionPointFilter intersectionPointFilter = filter_intersection_points_to_those_not_already_handled(nodesToCapturedDomains); - return build_intersection_points(mesh, geometry, nodesToCapturedDomains, intersectionPointFilter); + return build_intersection_points(mesh, elementSelector, geometry, nodesToCapturedDomains, intersectionPointFilter); } static IntersectionPointFilter @@ -231,17 +234,26 @@ filter_intersection_points_to_those_using_previous_iteration_snap_nodes_but_not_ } std::vector get_owned_elements_using_nodes_knowing_that_nodes_dont_have_common_elements(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const std::vector & nodes) { std::vector nodeElements; for (auto node : nodes) + { for (auto element : StkMeshEntities{mesh.begin_elements(node), mesh.end_elements(node)}) - if (mesh.bucket(element).owned()) + { + const auto & bucket = mesh.bucket(element); + if (bucket.owned() && elementSelector(bucket)) + { nodeElements.push_back(element); + } + } + } return nodeElements; } std::vector update_intersection_points_after_snap_iteration(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const std::vector & iterationSortedSnapNodes, const NodeToCapturedDomainsMap & nodesToCapturedDomains, @@ -269,7 +281,7 @@ std::vector update_intersection_points_after_snap_iteration(const stk::m if (geometry.snapped_elements_may_have_new_intersections()) { - const std::vector updateElements = get_owned_elements_using_nodes_knowing_that_nodes_dont_have_common_elements(mesh, iterationSortedSnapNodes); + const std::vector updateElements = get_owned_elements_using_nodes_knowing_that_nodes_dont_have_common_elements(mesh, elementSelector, iterationSortedSnapNodes); geometry.append_element_intersection_points(mesh, nodesToCapturedDomains, updateElements, intersectionPointFilter, intersectionPoints); communicate_intersection_points(mesh, newSize, intersectionPoints); } diff --git a/packages/krino/krino/krino_lib/Akri_Intersection_Points.hpp b/packages/krino/krino/krino_lib/Akri_Intersection_Points.hpp index 108d138cbef0..ac241e0b2f60 100644 --- a/packages/krino/krino/krino_lib/Akri_Intersection_Points.hpp +++ b/packages/krino/krino/krino_lib/Akri_Intersection_Points.hpp @@ -76,14 +76,17 @@ std::string debug_output(const stk::mesh::BulkData & mesh, const IntersectionPoi IntersectionPointFilter keep_all_intersection_points_filter(); std::vector build_all_intersection_points(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const NodeToCapturedDomainsMap & nodesToCapturedDomains); std::vector build_uncaptured_intersection_points(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const NodeToCapturedDomainsMap & nodesToCapturedDomains); std::vector update_intersection_points_after_snap_iteration(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & geometry, const std::vector & iterationSortedSnapNodes, const NodeToCapturedDomainsMap & nodesToCapturedDomains, diff --git a/packages/krino/krino/krino_lib/Akri_LevelSet.cpp b/packages/krino/krino/krino_lib/Akri_LevelSet.cpp index 221ea8ac1ba5..de07769574f2 100644 --- a/packages/krino/krino/krino_lib/Akri_LevelSet.cpp +++ b/packages/krino/krino/krino_lib/Akri_LevelSet.cpp @@ -8,6 +8,7 @@ #include +#include #include #include #include @@ -38,6 +39,7 @@ #include #include #include +#include namespace krino { @@ -139,8 +141,8 @@ void LevelSet::setup(void) { register_fields(); - facets.reset(new Faceted_Surface("current_facets")); - facets_old.reset(new Faceted_Surface("old_facets")); + facets = FacetedSurfaceBase::build(my_meta.spatial_dimension()); + facets_old = FacetedSurfaceBase::build(my_meta.spatial_dimension()); // initializes CDFEM_Support if (!meta().is_commit()) @@ -280,9 +282,11 @@ void LevelSet::write_facets(void) { /* %TRACE[ON]% */ Trace trace__("krino::LevelSet::facets_exoii(void)"); /* %TRACE% */ - Faceted_Surface & f = *facets; const std::string fileBaseName = "facets_" + name(); - krino::write_facets(spatial_dimension, f, fileBaseName, my_facetFileIndex++, mesh().parallel()); + if (2 == mesh().mesh_meta_data().spatial_dimension()) + krino::write_facets(facets->get_facets_2d(), fileBaseName, my_facetFileIndex++, mesh().parallel()); + else + krino::write_facets(facets->get_facets_3d(), fileBaseName, my_facetFileIndex++, mesh().parallel()); } //----------------------------------------------------------------------------------- @@ -600,10 +604,8 @@ void LevelSet::locally_conserved_redistance() } } // end bucket loop - const double localSumCorrectionSquared = sumCorrectionSquared; - stk::all_reduce_sum(mesh().parallel(), &localSumCorrectionSquared, &sumCorrectionSquared, 1); - const size_t localCountCorrection = countCorrection; - stk::all_reduce_sum(mesh().parallel(), &localCountCorrection, &countCorrection, 1); + all_reduce_sum(mesh().parallel(), sumCorrectionSquared); + all_reduce_sum(mesh().parallel(), countCorrection); const double correctionNorm = std::sqrt(sumCorrectionSquared/countCorrection/avgElemSize); if(krinolog.shouldPrint(LOG_DEBUG)) @@ -763,6 +765,25 @@ void LevelSet::redistance_using_existing_facets(const stk::mesh::Selector & volu facets->swap( *facets_old ); } +void LevelSet::redistance_nodes_using_existing_facets(const std::vector & nodesToRedistance) +{ + const FieldRef coordsField = get_coordinates_field(); + const FieldRef distField = get_distance_field(); + + const BoundingBox nodeBbox = krino::compute_nodal_bbox( mesh(), coordsField, nodesToRedistance ); + + const double narrowBandSize = 0.; + facets->prepare_to_compute(nodeBbox, narrowBandSize); + + for ( auto && node : nodesToRedistance ) + { + double & d = get_scalar_field(mesh(), distField, node); + const stk::math::Vector3d coords = get_vector_field(mesh(), coordsField, node, spatial_dimension); + + d = facets->point_signed_distance(coords); + } +} + void LevelSet::redistance() { redistance(my_meta.universal_part()); } void @@ -792,19 +813,62 @@ LevelSet::redistance(const stk::mesh::Selector & volumeSelector) redistance_using_existing_facets(volumeSelector); } +static std::vector get_owned_and_shared_interface_and_child_element_nodes(const stk::mesh::BulkData & mesh, + const AuxMetaData & auxMeta, + const CDFEM_Support & cdfemSupport, + const stk::mesh::Selector interfaceSelector) +{ + std::vector initialNodes; + const stk::mesh::Selector ownedOrSharedInterfaceSelector = interfaceSelector & auxMeta.active_part() & (mesh.mesh_meta_data().locally_owned_part() | mesh.mesh_meta_data().globally_shared_part()); + stk::mesh::get_selected_entities( ownedOrSharedInterfaceSelector, mesh.buckets(stk::topology::NODE_RANK), initialNodes, false ); + + stk::mesh::Selector ownedChildSelector = cdfemSupport.get_child_part() & mesh.mesh_meta_data().locally_owned_part(); + for (auto * bucketPtr : mesh.get_buckets(stk::topology::ELEMENT_RANK, ownedChildSelector)) + for (auto elem : *bucketPtr) + for (auto node : StkMeshEntities{mesh.begin_nodes(elem), mesh.end_nodes(elem)}) + initialNodes.push_back(node); + stk::util::sort_and_unique(initialNodes, stk::mesh::EntityLess(mesh)); + return initialNodes; +} + void LevelSet::interface_conforming_redistance() { - krinolog << "Redistancing the level set field..." << stk::diag::dendl; + if (FAST_MARCHING == my_redistance_method) + krinolog << "Redistancing the level set field using CDFEM fast marching method..." << stk::diag::dendl; + else + krinolog << "Redistancing the level set field using CDFEM method..." << stk::diag::dendl; + + build_interface_conforming_facets(); sync_all_fields_to_host(); + if (FAST_MARCHING == my_redistance_method) + fast_marching_interface_conforming_redistance_using_existing_facets(); + else + redistance_using_existing_facets(my_meta.universal_part()); +} + +void +LevelSet::fast_marching_interface_conforming_redistance_using_existing_facets() +{ const auto & phaseSupport = Phase_Support::get(meta()); const stk::mesh::Selector interfaceSelector = phaseSupport.get_negative_levelset_interface_selector(my_identifier); - const stk::mesh::Selector negativeBlockSelector = phaseSupport.get_negative_levelset_block_selector(my_identifier); - build_interface_conforming_facets(interfaceSelector, negativeBlockSelector); + const CDFEM_Support & cdfemSupport = CDFEM_Support::get(meta()); + const std::vector initialNodes = get_owned_and_shared_interface_and_child_element_nodes(mesh(), aux_meta(), cdfemSupport, interfaceSelector); + + redistance_nodes_using_existing_facets(initialNodes); - redistance_using_existing_facets(my_meta.universal_part()); + const stk::mesh::Selector activeVolumeSelector = aux_meta().active_part(); + + std::function get_interface_speed; + Fast_Marching fm(mesh(), + activeVolumeSelector, + get_coordinates_field(), + get_distance_field(), + get_interface_speed, + my_redistance_timer); + fm.redistance(initialNodes); } static bool determine_polarity_for_negative_side_of_interface(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side) @@ -832,13 +896,12 @@ static std::array get_oriented_triangle_side_nodes(const st return {{sideNodes[0], sideNodes[2], sideNodes[1]}}; } -static void append_facets_from_triangle_side(const stk::mesh::BulkData & mesh, const FieldRef coords, const stk::mesh::Selector & interfaceSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side, Faceted_Surface & facets) +static void append_facets_from_triangle_side(const stk::mesh::BulkData & mesh, const FieldRef coords, const stk::mesh::Selector & interfaceSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side, FacetedSurfaceBase & facets) { const std::array orientedSideNodes = get_oriented_triangle_side_nodes(mesh, negativeSideElementSelector, side); const std::array sideNodeCoords{{stk::math::Vector3d(field_data(coords, orientedSideNodes[0]), 3), stk::math::Vector3d(field_data(coords, orientedSideNodes[1]), 3), stk::math::Vector3d(field_data(coords, orientedSideNodes[2]), 3)}}; - std::unique_ptr facet = std::make_unique( sideNodeCoords[0], sideNodeCoords[1], sideNodeCoords[2] ); - facets.add( std::move(facet) ); + facets.emplace_back_3d( sideNodeCoords[0], sideNodeCoords[1], sideNodeCoords[2] ); } static std::array get_oriented_line_side_nodes(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side) @@ -851,13 +914,12 @@ static std::array get_oriented_line_side_nodes(const stk::m return {{sideNodes[1], sideNodes[0]}}; } -static void append_facets_from_line_side(const stk::mesh::BulkData & mesh, const FieldRef coords, const stk::mesh::Selector & sideSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side, Faceted_Surface & facets) +static void append_facets_from_line_side(const stk::mesh::BulkData & mesh, const FieldRef coords, const stk::mesh::Selector & sideSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side, FacetedSurfaceBase & facets) { const std::array orientedSideNodes = get_oriented_line_side_nodes(mesh, negativeSideElementSelector, side); const std::array sideNodeCoords{{stk::math::Vector3d(field_data(coords, orientedSideNodes[0]), 2), stk::math::Vector3d(field_data(coords, orientedSideNodes[1]), 2)}}; - std::unique_ptr facet = std::make_unique(sideNodeCoords[0], sideNodeCoords[1]); - facets.add( std::move(facet) ); + facets.emplace_back_2d(sideNodeCoords[0], sideNodeCoords[1]); } void LevelSet::append_facets_from_side(const stk::mesh::Selector & sideSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side) @@ -869,8 +931,12 @@ void LevelSet::append_facets_from_side(const stk::mesh::Selector & sideSelector, } void -LevelSet::build_interface_conforming_facets(const stk::mesh::Selector & interfaceSelector, const stk::mesh::Selector & negativeSideBlockSelector) +LevelSet::build_interface_conforming_facets() { + const auto & phaseSupport = Phase_Support::get(meta()); + const stk::mesh::Selector interfaceSelector = phaseSupport.get_negative_levelset_interface_selector(my_identifier); + const stk::mesh::Selector negativeSideBlockSelector = phaseSupport.get_negative_levelset_block_selector(my_identifier); + const stk::mesh::Selector sideSelector = interfaceSelector & aux_meta().active_part(); const stk::mesh::Selector ownedSideSelector = sideSelector & meta().locally_owned_part(); @@ -913,7 +979,7 @@ LevelSet::get_time_of_arrival_speed(stk::mesh::Entity elem, ParallelErrorMessage } void -LevelSet::fast_methods_redistance(const stk::mesh::Selector & selector, const bool compute_time_of_arrival) +LevelSet::fast_methods_redistance(const stk::mesh::Selector & volumeSelector, const bool compute_time_of_arrival) { /* %TRACE[ON]% */ Trace trace__("krino::LevelSet::fast_marching_redistance(const stk::mesh::Selector & selector)"); /* %TRACE% */ // Unlike redistance() this method provides an approximate (not exact) distance to the isosurface. @@ -921,6 +987,12 @@ LevelSet::fast_methods_redistance(const stk::mesh::Selector & selector, const bo // different than a pure distance function because it provides the distance through domain. It can't see // through walls like the redistance() method does. + std::function get_interface_speed; + if (compute_time_of_arrival) + get_interface_speed = [&](ParallelErrorMessage& err, stk::mesh::Entity elem) { return get_time_of_arrival_speed(elem, err); }; + + const stk::mesh::Selector activeVolumeSelector = aux_meta().active_part() & volumeSelector; + if (my_redistance_method == FAST_MARCHING) { if (compute_time_of_arrival) @@ -928,25 +1000,25 @@ LevelSet::fast_methods_redistance(const stk::mesh::Selector & selector, const bo else krinolog << "Redistancing the level set field using a fast marching method..." << stk::diag::dendl; - Fast_Marching fm(*this, selector, my_redistance_timer); + Fast_Marching fm(mesh(), + activeVolumeSelector, + get_coordinates_field(), + get_distance_field(), + get_interface_speed, + my_redistance_timer); fm.redistance(); + stk::mesh::field_copy(get_distance_field(), get_old_distance_field()); } else { STK_ThrowRequire(my_redistance_method == FAST_ITERATIVE); - std::function get_interface_speed; if (compute_time_of_arrival) - { krinolog << "Initializing the level set field to be the time-of-arrival using a fast iterative method..." << stk::diag::dendl; - get_interface_speed = [&](ParallelErrorMessage& err, stk::mesh::Entity elem) { return get_time_of_arrival_speed(elem, err); }; - } else - { krinolog << "Redistancing the level set field using a fast iterative method..." << stk::diag::dendl; - } FastIterativeMethod fim(mesh(), - selector, + activeVolumeSelector, get_coordinates_field(), get_distance_field(), get_interface_speed, @@ -1485,10 +1557,7 @@ LevelSet::remove_wall_features() const } } - int global_small_feature_removed = false; - stk::all_reduce_sum(mesh().parallel(), &small_feature_removed, &global_small_feature_removed, 1); - - if (global_small_feature_removed) + if (stk::is_true_on_any_proc(mesh().parallel(), small_feature_removed)) { stk::mesh::communicate_field_data(mesh(), {&dField.field()}); return true; @@ -1584,7 +1653,6 @@ LevelSet::elem_on_interface(stk::mesh::Entity e) const double *d = field_data(isoField, nodes[0]); double first_value = *d - my_threshold; - bool inside_narrow_band = fabs(first_value) < 0.5*my_narrow_band_size; for (unsigned i = 1; i < nnodes; ++i ) { @@ -1593,7 +1661,6 @@ LevelSet::elem_on_interface(stk::mesh::Entity e) const double value = *d - my_threshold; have_crossing |= sign_change(value, first_value); - inside_narrow_band |= (fabs(value) < 0.5*my_narrow_band_size); } // It is the user's job to make sure that the narrow band is sufficiently @@ -1625,20 +1692,7 @@ LevelSet::compute_levelset_sizes( double & area, double & negVol, double & posVo for (auto && elem : objs) accumulate_levelset_integrals_on_element(mesh(), elem, area, negVol, posVol, h_avg, xField, isovar, isoval); - // communicate global sums - const int vec_length = 3; - std::vector local_sum( vec_length ); - std::vector global_sum( vec_length ); - local_sum[0] = area; - local_sum[1] = negVol; - local_sum[2] = posVol; - - stk::all_reduce_sum(mesh().parallel(), &local_sum[0], &global_sum[0], vec_length); - - area = global_sum[0]; - negVol = global_sum[1]; - posVol = global_sum[2]; - + all_reduce_sum(mesh().parallel(), area, negVol, posVol); } //-------------------------------------------------------------------------------- void @@ -1647,28 +1701,6 @@ LevelSet::compute_sizes( double & area, double & negVol, double & posVol, const compute_levelset_sizes(area, negVol, posVol, get_isovar_field(), isoval); } -template -void all_reduce_sum(stk::ParallelMachine comm, Args&&... args) -{ - typedef typename std::common_type::type T; - const std::array local = {{ args... }}; - std::array global; - stk::all_reduce_sum(comm, local.data(), global.data(), local.size()); - const T * data = global.data(); - ((std::forward(args) = *(data++)), ...); -} - -template -void all_reduce_max(stk::ParallelMachine comm, Args&&... args) -{ - typedef typename std::common_type::type T; - const std::array local = {{ args... }}; - std::array global; - stk::all_reduce_max(comm, local.data(), global.data(), local.size()); - const T * data = global.data(); - ((std::forward(args) = *(data++)), ...); -} - static double get_gradient_magnitude_at_ip(const sierra::ArrayContainer & gradDist, const int ip) { double mag2GradPhi = 0.; @@ -1834,7 +1866,7 @@ LevelSet::build_facets_locally(const stk::mesh::Selector & selector) } void -LevelSet::build_facets_for_elements(const stk::mesh::BulkData & mesh, const FieldRef xField, const FieldRef isoField, const std::vector & elementsToIntersect, const double avgEdgeLength, Faceted_Surface & facets) +LevelSet::build_facets_for_elements(const stk::mesh::BulkData & mesh, const FieldRef xField, const FieldRef isoField, const std::vector & elementsToIntersect, const double avgEdgeLength, FacetedSurfaceBase & facets) { facets.clear(); @@ -1931,49 +1963,12 @@ LevelSet::gather_nodal_field( //-------------------------------------------------------------------------------- std::string print_sizes(const LevelSet & ls) -{ /* %TRACE[ON]% */ Trace trace__("krino::print_sizes(const LevelSet & levelSet)"); /* %TRACE% */ - - // - // find sizes of facets stored in vector of facet descriptions - // - - const auto & ls_surfaces = ls.get_facets().get_facets(); - unsigned local_facet_num = ls_surfaces.size(); - unsigned global_facet_num = 0; - stk::all_reduce_sum(ls.mesh().parallel(), &local_facet_num , &global_facet_num, 1); - - // iterate local facets and calc area - - double local_facets_totalArea = 0.0; // total surface area of interface from facets - double local_facets_maxArea = -1.0; // area of largest facet on interface - double local_facets_minArea = std::numeric_limits::max(); // area of smallest facet on interface - - // loop over facets - for ( auto&& surface : ls_surfaces ) - { - double area = surface->facet_area(); - local_facets_totalArea += area; - - local_facets_maxArea = std::min(area, local_facets_maxArea); - local_facets_minArea = std::max(area, local_facets_maxArea); - } - - double global_facets_totalArea = 0.0; - double global_facets_maxArea = 0.0; - double global_facets_minArea = 0.0; - - stk::all_reduce_min(ls.mesh().parallel(), &local_facets_minArea, &global_facets_minArea, 1); - stk::all_reduce_max(ls.mesh().parallel(), &local_facets_maxArea, &global_facets_maxArea, 1); - stk::all_reduce_sum(ls.mesh().parallel(), &local_facets_totalArea, &global_facets_totalArea, 1); - +{ std::ostringstream out; - // facet info - out << "P" << ls.mesh().parallel_rank() << ": facets for level set '" << ls.name() << "' " << std::endl ; - out << "\t " << "Global sizes: { " << global_facet_num << " facets on }" << std::endl; - out << "\t " << "Local sizes: { " << local_facet_num << " facets on }" << std::endl; - out << "\t Global areas: { Min = " << global_facets_minArea << ", Max = " << global_facets_maxArea - << ", Total = " << global_facets_totalArea << " }" << std::endl; + out << "P" << ls.mesh().parallel_rank() << ": facets for level set '" << ls.name() << "' " << std::endl; + out << ls.get_facets().print_sizes() << std::endl; + return out.str(); } //-------------------------------------------------------------------------------- diff --git a/packages/krino/krino/krino_lib/Akri_LevelSet.hpp b/packages/krino/krino/krino_lib/Akri_LevelSet.hpp index dc35e823200b..e31d8c558569 100644 --- a/packages/krino/krino/krino_lib/Akri_LevelSet.hpp +++ b/packages/krino/krino/krino_lib/Akri_LevelSet.hpp @@ -81,7 +81,7 @@ friend class LevelSet_Size; const FieldRef & field_ref, double * field); - static void build_facets_for_elements(const stk::mesh::BulkData & mesh, const FieldRef xField, const FieldRef isoField, const std::vector & elementsToIntersect, const double avgEdgeLength, Faceted_Surface & facets); + static void build_facets_for_elements(const stk::mesh::BulkData & mesh, const FieldRef xField, const FieldRef isoField, const std::vector & elementsToIntersect, const double avgEdgeLength, FacetedSurfaceBase & facets); double compute_average_edge_length() const; void build_facets_locally(const stk::mesh::Selector & selector); @@ -152,8 +152,8 @@ friend class LevelSet_Size; void set_redistance_method( const Redistance_Method type ) { my_redistance_method = type; } void set_time_of_arrival_element_speed_field_name( const std::string & time_of_arrival_speed_field_name) { my_time_of_arrival_element_speed_field_name = time_of_arrival_speed_field_name; } void set_time_of_arrival_block_speed(const std::string & blockName, const double blockSpeed); - Faceted_Surface & get_facets() { return *facets; } - const Faceted_Surface & get_facets() const { return *facets; } + FacetedSurfaceBase & get_facets() { return *facets; } + const FacetedSurfaceBase & get_facets() const { return *facets; } void narrow_band_multiplier( double multiplier ) { my_narrow_band_multiplier = multiplier; } const double & narrow_band_size() const { return my_narrow_band_size; } @@ -202,6 +202,7 @@ friend class LevelSet_Size; void redistance(const stk::mesh::Selector & selector); void fast_methods_redistance(const stk::mesh::Selector & selector, const bool compute_time_of_arrival = false); void interface_conforming_redistance(); + void fast_marching_interface_conforming_redistance_using_existing_facets(); std::pair get_conserved_negative_volume_and_time() const; void set_conserved_negative_volume_and_time(const double vol, const double time); @@ -233,8 +234,9 @@ friend class LevelSet_Size; LevelSet(stk::mesh::MetaData & in_meta, const std::string & in_name, const stk::diag::Timer & parent_timer); void sync_all_fields_to_host(); void redistance_using_existing_facets(const stk::mesh::Selector & volumeSelector); + void redistance_nodes_using_existing_facets(const std::vector & nodesToRedistance); void append_facets_from_side(const stk::mesh::Selector & interfaceSelector, const stk::mesh::Selector & negativeSideElementSelector, const stk::mesh::Entity side); - void build_interface_conforming_facets(const stk::mesh::Selector & interfaceSelector, const stk::mesh::Selector & negativeSideBlockSelector); + void build_interface_conforming_facets(); stk::mesh::MetaData & my_meta; AuxMetaData & my_aux_meta; @@ -283,10 +285,10 @@ friend class LevelSet_Size; std::unique_ptr my_IC_alg; // vector of previous facets - std::unique_ptr facets_old; + std::unique_ptr facets_old; // vector of current facets - std::unique_ptr facets; + std::unique_ptr facets; stk::math::Vector3d my_extension_velocity; const double epsilon; diff --git a/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.cpp b/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.cpp index eaa7c597d340..c693c52f8fb6 100644 --- a/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.cpp +++ b/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.cpp @@ -12,7 +12,8 @@ namespace krino { -LevelSetSurfaceInterfaceGeometry::LevelSetSurfaceInterfaceGeometry(const stk::mesh::Part & activePart, +LevelSetSurfaceInterfaceGeometry::LevelSetSurfaceInterfaceGeometry(const int dim, + const stk::mesh::Part & activePart, const CDFEM_Support & cdfemSupport, const Phase_Support & phaseSupport, const std::vector & LSFields) @@ -22,11 +23,11 @@ LevelSetSurfaceInterfaceGeometry::LevelSetSurfaceInterfaceGeometry(const stk::me for (auto && lsField : myLSFields) { mySurfaceIdentifiers.push_back(lsField.identifier); - myLSSurfaces.emplace_back("My Facets"); + myLSSurfaces.emplace_back(std::move(FacetedSurfaceBase::build(dim))); } for (size_t i=0; i LevelSetSurfaceInterfaceGeometry::get_levelset_element_selectors() const @@ -101,7 +102,7 @@ void LevelSetSurfaceInterfaceGeometry::build_levelset_facets(const stk::mesh::Bu const stk::mesh::Selector elemFieldSelector = elementSelector & levelsetElementSelectors[i]; stk::mesh::get_selected_entities( elemFieldSelector, mesh.buckets(stk::topology::ELEMENT_RANK), elementsToIntersect, false ); - LevelSet::build_facets_for_elements(mesh, coordsField, myLSFields[i].isovar, elementsToIntersect, avgEdgeLength, myLSSurfaces[i]); + LevelSet::build_facets_for_elements(mesh, coordsField, myLSFields[i].isovar, elementsToIntersect, avgEdgeLength, *myLSSurfaces[i]); } } @@ -111,11 +112,5 @@ std::vector LevelSetSurfaceInterfaceGeometry::get_possibly_cu return LevelSetInterfaceGeometry::get_active_elements_that_may_be_cut_by_levelsets(mesh, get_active_part(), myLSFields); } -std::vector LevelSetSurfaceInterfaceGeometry::get_elements_that_intersect_interval(const stk::mesh::BulkData & mesh, const std::array loAndHi) const -{ - // NOTE: Uses levelset field directly, not facetted interface, because it needs the analog, not discrete version - return LevelSetInterfaceGeometry::get_active_elements_that_intersect_levelset_interval(mesh, get_active_part(), myLSFields, loAndHi); -} - } diff --git a/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.hpp b/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.hpp index 1b72a45ee9aa..9087c431754b 100644 --- a/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.hpp +++ b/packages/krino/krino/krino_lib/Akri_LevelSetSurfaceInterfaceGeometry.hpp @@ -14,7 +14,8 @@ class Phase_Support; class LevelSetSurfaceInterfaceGeometry : public AnalyticSurfaceInterfaceGeometry { public: - LevelSetSurfaceInterfaceGeometry(const stk::mesh::Part & activePart, + LevelSetSurfaceInterfaceGeometry(const int dim, + const stk::mesh::Part & activePart, const CDFEM_Support & cdfemSupport, const Phase_Support & phaseSupport, const std::vector & LSFields); @@ -31,7 +32,6 @@ class LevelSetSurfaceInterfaceGeometry : public AnalyticSurfaceInterfaceGeometry // Methods that use levelset directly and not facetted levelset surface, because they need the analog, not discrete version virtual std::vector get_possibly_cut_elements(const stk::mesh::BulkData & mesh) const override; - virtual std::vector get_elements_that_intersect_interval(const stk::mesh::BulkData & mesh, const std::array loAndHi) const override; private: std::vector get_levelset_element_selectors() const; @@ -40,7 +40,7 @@ class LevelSetSurfaceInterfaceGeometry : public AnalyticSurfaceInterfaceGeometry std::vector myLSFields; std::vector mySurfaceIdentifiers; - mutable std::vector myLSSurfaces; + mutable std::vector> myLSSurfaces; mutable size_t myLastMeshSyncCount{0}; mutable bool myDoUpdateFacetsWhenMeshChanges{true}; }; diff --git a/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.cpp b/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.cpp index a5c2097f9483..5e8c034805b5 100644 --- a/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.cpp +++ b/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.cpp @@ -30,4 +30,19 @@ compute_nodal_bbox( const stk::mesh::BulkData & mesh, const stk::mesh::Selector return bbox; } +BoundingBox +compute_nodal_bbox( const stk::mesh::BulkData & mesh, const FieldRef coordsField, const std::vector & nodes ) +{ + const int ndim = mesh.mesh_meta_data().spatial_dimension(); + BoundingBox bbox; + + for ( auto && node : nodes ) + { + stk::math::Vector3d x(field_data(coordsField, node), ndim); + bbox.accommodate( x ); + } + + return bbox; +} + } diff --git a/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.hpp b/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.hpp index df898b128057..c3fca1b9734c 100644 --- a/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.hpp +++ b/packages/krino/krino/krino_lib/Akri_NodalBoundingBox.hpp @@ -5,12 +5,14 @@ namespace stk { namespace mesh { class BulkData; } } namespace stk { namespace mesh { class Selector; } } +namespace stk { namespace mesh { class Entity; } } namespace krino { class FieldRef; BoundingBox compute_nodal_bbox( const stk::mesh::BulkData & mesh, const stk::mesh::Selector & selector, const FieldRef coordsField ); +BoundingBox compute_nodal_bbox( const stk::mesh::BulkData & mesh, const FieldRef coordsField, const std::vector & nodes ); } diff --git a/packages/krino/krino/krino_lib/Akri_OutputUtils.cpp b/packages/krino/krino/krino_lib/Akri_OutputUtils.cpp index 4306e17b2e50..57acf9b21827 100644 --- a/packages/krino/krino/krino_lib/Akri_OutputUtils.cpp +++ b/packages/krino/krino/krino_lib/Akri_OutputUtils.cpp @@ -86,30 +86,32 @@ std::string create_file_name(const std::string & fileBaseName, const int fileInd return filename; } -void -write_facets( const int dim, const Faceted_Surface & facetedSurface, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm) +template +void write_facets( const std::vector & facets, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm) { - int nfaces = facetedSurface.size(); - const int nodes_per_elem = dim; - const int nnodes = nfaces * nodes_per_elem; - const int nelems = nfaces * 1; //writing out faces as elements + int numFacets = facets.size(); + const int nodesPerFacet = FACET::DIM; + const int numNodes = numFacets * nodesPerFacet; + const unsigned numProcs = stk::parallel_machine_size(comm); + const int procId = stk::parallel_machine_rank(comm); // Type (ExodusII) is hard-wired at this time. - Ioss::DatabaseIO *db = Ioss::IOFactory::create("exodusII", create_file_name(fileBaseName, fileIndex), Ioss::WRITE_RESULTS); + Ioss::PropertyManager properties; + properties.add(Ioss::Property("COMPOSE_RESULTS", 1)); + const std::string fileName = create_file_name(fileBaseName, fileIndex); + Ioss::DatabaseIO *db = Ioss::IOFactory::create("exodusII", create_file_name(fileBaseName, fileIndex), Ioss::WRITE_RESULTS, comm, properties); + STK_ThrowRequireMsg(db != nullptr && db->ok(), "ERROR: Could not open output database '" << fileName << "' of type 'exodus'\n"); Ioss::Region io(db, "FacetRegion"); // find offsets for consistent global numbering int elem_start = 0; - if ( stk::parallel_machine_size(comm) > 1 ) { - std::vector< int > num_faces; - num_faces.resize( stk::parallel_machine_size(comm) ); - - // Gather everyone's facet list sizes - // I don't think there is a framework call for this... - MPI_Allgather( &nfaces, 1, MPI_INT, - &(num_faces[0]), 1, MPI_INT, comm ); - for (int i = 0; i< stk::parallel_machine_rank(comm); ++i) { - elem_start += num_faces[i]; + if ( numProcs > 1 ) { + std::vector< int > numFaces(numProcs); + // Gather everyone's facet list sizes. I don't think there is a stk call for this... + MPI_Allgather( &numFacets, 1, MPI_INT, + &(numFaces[0]), 1, MPI_INT, comm ); + for (int i = 0; i 0) +// { +// io.end_mode(Ioss::STATE_DEFINE_MODEL); +// return; +// } - Ioss::NodeBlock *nb = new Ioss::NodeBlock(db, "nodeblock_1", nnodes, dim); + Ioss::NodeBlock *nb = new Ioss::NodeBlock(db, "nodeblock_1", numNodes, nodesPerFacet); io.add(nb); + nb->property_add(Ioss::Property("locally_owned_count", numNodes)); std::string el_type; - if (dim == 3) + if constexpr (FACET::DIM == 3) el_type = "trishell3"; else el_type = "shellline2d2"; - Ioss::ElementBlock *eb = new Ioss::ElementBlock(db, "block_1", el_type, nelems); + Ioss::ElementBlock *eb = new Ioss::ElementBlock(db, "block_1", el_type, numFacets); io.add(eb); io.end_mode(Ioss::STATE_DEFINE_MODEL); io.begin_mode(Ioss::STATE_MODEL); - const FacetOwningVec & facets = facetedSurface.get_facets(); - // loop over elements and node to create maps { std::vector< int > nmap, emap; - nmap.reserve(nnodes); - emap.reserve(nelems); + nmap.reserve(numNodes); + emap.reserve(numFacets); for ( unsigned n=0, e=0; eput_field_data("ids", emap); } + if (nb->get_database()->needs_shared_node_information()) + { + std::vector owningProcessor(numNodes, procId); + nb->put_field_data("owning_processor", owningProcessor); + } + // generate coordinates { std::vector< double > xyz; - xyz.reserve(dim*nnodes); + xyz.reserve(FACET::DIM*numNodes); for ( auto&& facet : facets ) { - for ( int j = 0; j < nodes_per_elem; ++j ) { - const stk::math::Vector3d & vert = facet->facet_vertex(j); + for ( int j = 0; j < nodesPerFacet; ++j ) { + const stk::math::Vector3d & vert = facet.facet_vertex(j); xyz.push_back(vert[0]); xyz.push_back(vert[1]); - if (3 == dim) xyz.push_back(vert[2]); + if (3 == FACET::DIM) xyz.push_back(vert[2]); } } nb->put_field_data("mesh_model_coordinates", xyz); @@ -176,9 +184,9 @@ write_facets( const int dim, const Faceted_Surface & facetedSurface, const std:: // generate connectivity { std::vector< int > conn; - conn.reserve(nnodes); - for ( int n = 0; n < nnodes; ++n ) { - conn.push_back(nodes_per_elem*elem_start + n+1); + conn.reserve(numNodes); + for ( int n = 0; n < numNodes; ++n ) { + conn.push_back(nodesPerFacet*elem_start + n+1); } eb->put_field_data("connectivity", conn); } @@ -207,6 +215,11 @@ stk::mesh::PartVector turn_off_output_for_empty_io_parts(const stk::mesh::BulkDa return emptyParts; } +// Explicit template instantiation + +template void write_facets( const std::vector & facets, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm); +template void write_facets( const std::vector & facets, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm); + } diff --git a/packages/krino/krino/krino_lib/Akri_OutputUtils.hpp b/packages/krino/krino/krino_lib/Akri_OutputUtils.hpp index fd04670fe4cb..2b8f347c2642 100644 --- a/packages/krino/krino/krino_lib/Akri_OutputUtils.hpp +++ b/packages/krino/krino/krino_lib/Akri_OutputUtils.hpp @@ -14,16 +14,19 @@ #include #include -namespace krino { class Faceted_Surface; } - namespace krino { +template +class Faceted_Surface; + void output_mesh_with_fields(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & outputSelector, const std::string & fileName, int step, double time, stk::io::DatabasePurpose purpose = stk::io::WRITE_RESULTS); void output_composed_mesh_with_fields(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & outputSelector, const std::string & fileName, int step, double time, stk::io::DatabasePurpose purpose = stk::io::WRITE_RESULTS); std::string create_file_name(const std::string & fileBaseName, const int fileIndex); -void write_facets( const int dim, const Faceted_Surface & facetedSurface, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm); stk::mesh::PartVector turn_off_output_for_empty_io_parts(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & outputSelector); +template +void write_facets( const std::vector & facetedSurface, const std::string & fileBaseName, const int fileIndex, const stk::ParallelMachine comm); + } diff --git a/packages/krino/krino/krino_lib/Akri_ProlongationData.cpp b/packages/krino/krino/krino_lib/Akri_ProlongationData.cpp index e59f774228c0..732399ac0449 100644 --- a/packages/krino/krino/krino_lib/Akri_ProlongationData.cpp +++ b/packages/krino/krino/krino_lib/Akri_ProlongationData.cpp @@ -575,14 +575,14 @@ static bool does_any_node_have_field(const FieldRef field, const std::vector & facetDistanceQuery, const std::vector & facetNodes) { interpolate_to_point(mesh.stk_meta(), facetDistanceQuery, facetNodes); } void ProlongationFacetPointData::interpolate_to_point(const stk::mesh::MetaData & meta, - const FacetDistanceQuery & facetDistanceQuery, + const FacetDistanceQuery & facetDistanceQuery, const std::vector & facetNodes) { const stk::math::Vector3d node_wts = facetDistanceQuery.closest_point_weights(); @@ -871,7 +871,7 @@ void ProlongationFacet::build_and_append_edge_facets(ProlongFacetVec & facetVec) } } -std::unique_ptr ProlongationFacet::get_prolongation_point_data(const FacetDistanceQuery & dist_query) const +std::unique_ptr ProlongationFacet::get_prolongation_point_data(const FacetDistanceQuery & dist_query) const { STK_ThrowAssert(&dist_query.facet() == my_facet.get()); return std::make_unique(my_mesh, dist_query, my_prolong_nodes); diff --git a/packages/krino/krino/krino_lib/Akri_ProlongationData.hpp b/packages/krino/krino/krino_lib/Akri_ProlongationData.hpp index 59d0176e2bac..061f9b787d1e 100644 --- a/packages/krino/krino/krino_lib/Akri_ProlongationData.hpp +++ b/packages/krino/krino/krino_lib/Akri_ProlongationData.hpp @@ -113,9 +113,9 @@ class ProlongationPointData : public ProlongationData { class ProlongationFacetPointData : public ProlongationPointData { public: - ProlongationFacetPointData(const CDMesh & mesh, const FacetDistanceQuery & facetDistanceQuery, const std::vector & facetNodes); + ProlongationFacetPointData(const CDMesh & mesh, const FacetDistanceQuery & facetDistanceQuery, const std::vector & facetNodes); private: - void interpolate_to_point(const stk::mesh::MetaData & meta, const FacetDistanceQuery & facetDistanceQuery, const std::vector & facetNodes); + void interpolate_to_point(const stk::mesh::MetaData & meta, const FacetDistanceQuery & facetDistanceQuery, const std::vector & facetNodes); std::vector myFieldStorageIndices; }; @@ -203,7 +203,7 @@ class ProlongationFacet { Facet * get_facet() const { return my_facet.get(); } std::vector compute_common_fields() const; - std::unique_ptr get_prolongation_point_data(const FacetDistanceQuery & dist_query) const; + std::unique_ptr get_prolongation_point_data(const FacetDistanceQuery & dist_query) const; const std::vector & get_prolongation_nodes() const { return my_prolong_nodes; } bool communicate_me(const BoundingBox & proc_target_bbox) const; static void insert_into_bounding_box(const ProlongationFacet * prolongFacet, BoundingBox & bbox) { return prolongFacet->get_facet()->insert_into(bbox); } @@ -225,7 +225,7 @@ class ProlongationFacet { class ProlongationQuery { public: ProlongationQuery() = default; - ProlongationQuery(const ProlongationFacet & facet, const FacetDistanceQuery & distQuery) { myProlongFacetPointData = facet.get_prolongation_point_data(distQuery); } + ProlongationQuery(const ProlongationFacet & facet, const FacetDistanceQuery & distQuery) { myProlongFacetPointData = facet.get_prolongation_point_data(distQuery); } ProlongationQuery(const ProlongationNodeData * prolongNodeData) : myProlongNodeData(prolongNodeData) {} const ProlongationPointData * get_prolongation_point_data() const { return myProlongFacetPointData ? get_facet_point_data() : myProlongNodeData; } private: diff --git a/packages/krino/krino/krino_lib/Akri_RefineNearLevelSets.cpp b/packages/krino/krino/krino_lib/Akri_RefineNearLevelSets.cpp index 9349c02ac380..f6b5db2c0037 100644 --- a/packages/krino/krino/krino_lib/Akri_RefineNearLevelSets.cpp +++ b/packages/krino/krino/krino_lib/Akri_RefineNearLevelSets.cpp @@ -40,7 +40,7 @@ void refine_elements_that_intersect_distance_interval_from_levelsets(stk::mesh:: const CDFEM_Support & cdfemSupport = CDFEM_Support::get(mesh.mesh_meta_data()); const Phase_Support & phaseSupport = Phase_Support::get(mesh.mesh_meta_data()); const auto lsFields = build_levelset_fields(stkLevelSetFields); - const LevelSetSurfaceInterfaceGeometry interfaceGeom(activePart, cdfemSupport, phaseSupport, lsFields); + const LevelSetSurfaceInterfaceGeometry interfaceGeom(mesh.mesh_meta_data().spatial_dimension(), activePart, cdfemSupport, phaseSupport, lsFields); const int numRefinementSteps = 2*refinementSupport.get_interface_maximum_refinement_level(); // Make sure refinement completes so that elements touching interval are fully refined diff --git a/packages/krino/krino/krino_lib/Akri_RefinementInterface.cpp b/packages/krino/krino/krino_lib/Akri_RefinementInterface.cpp index df112ae5875f..2c7d991f3766 100644 --- a/packages/krino/krino/krino_lib/Akri_RefinementInterface.cpp +++ b/packages/krino/krino/krino_lib/Akri_RefinementInterface.cpp @@ -24,14 +24,14 @@ namespace krino { void clear_refinement_marker(const RefinementInterface & refinement) { FieldRef markerField = refinement.get_marker_field(); - stk::mesh::field_fill(static_cast(Refinement::NOTHING), markerField, stk::mesh::selectField(markerField)); + stk::mesh::field_fill(static_cast(Refinement::RefinementMarker::NOTHING), markerField, stk::mesh::selectField(markerField)); } void mark_selected_elements_for_refinement(const RefinementInterface & refinement, const stk::mesh::Selector & selector) { FieldRef markerField = refinement.get_marker_field(); clear_refinement_marker(refinement); - stk::mesh::field_fill(static_cast(Refinement::REFINE), markerField, selector); + stk::mesh::field_fill(static_cast(Refinement::RefinementMarker::REFINE), markerField, selector); } void mark_selected_elements_for_refinement(const RefinementInterface & refinement, const int current_refinement_level, const int max_refinement_levels, const stk::mesh::Selector & selector) @@ -50,7 +50,7 @@ void mark_elements_for_refinement(const RefinementInterface & refinement, const for (auto && elem : elemsToRefine) { int * elemMarker = field_data(markerField, elem); - *elemMarker = Refinement::REFINE; + *elemMarker = static_cast(Refinement::RefinementMarker::REFINE); } } @@ -131,7 +131,9 @@ void mark_based_on_indicator_field(const stk::mesh::BulkData & mesh, const int size = bucket.size(); for(int i=0; i < size; ++i) { - markerData[i] = indicatorData[i] > threshold_val ? Refinement::REFINE : Refinement::NOTHING; + markerData[i] = indicatorData[i] > threshold_val + ? static_cast(Refinement::RefinementMarker::REFINE) + : static_cast(Refinement::RefinementMarker::NOTHING); } } } @@ -189,19 +191,29 @@ KrinoRefinement::get(const stk::mesh::MetaData & meta) } KrinoRefinement & -KrinoRefinement::create(stk::mesh::MetaData & meta) +KrinoRefinement::create(stk::mesh::MetaData & meta, stk::diag::Timer & timer) { KrinoRefinement * refinement = const_cast(meta.get_attribute()); STK_ThrowRequireMsg(nullptr == refinement, "KrinoRefinement::create should be called only once per MetaData."); if (nullptr == refinement) { AuxMetaData & auxMeta = AuxMetaData::get(meta); - refinement = new KrinoRefinement(meta, &auxMeta.active_part(), false/*auxMeta.get_force_64bit_flag()*/, auxMeta.get_assert_32bit_flag()); + refinement = new KrinoRefinement(meta, + &auxMeta.active_part(), + false /*auxMeta.get_force_64bit_flag()*/, + auxMeta.get_assert_32bit_flag(), + timer); meta.declare_attribute_with_delete(refinement); } return *refinement; } +KrinoRefinement & +KrinoRefinement::create(stk::mesh::MetaData & meta) +{ + return KrinoRefinement::create(meta, sierra::Diag::sierraTimer()); +} + bool KrinoRefinement::is_created(const stk::mesh::MetaData & meta) { @@ -224,9 +236,14 @@ KrinoRefinement::register_parts_and_fields_via_aux_meta_for_fmwk(stk::mesh::Meta } } -KrinoRefinement::KrinoRefinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart, const bool force64Bit, const bool assert32Bit) : - myMeta(meta), - myRefinement(meta, activePart, force64Bit, assert32Bit) +KrinoRefinement::KrinoRefinement(stk::mesh::MetaData & meta, + stk::mesh::Part * activePart, + const bool force64Bit, + const bool assert32Bit, + stk::diag::Timer & parent_timer) + : myMeta(meta), + myRefinementTimer("Refinement", parent_timer), + myRefinement(meta, activePart, force64Bit, assert32Bit, myRefinementTimer) { myElementMarkerField = meta.declare_field(stk::topology::ELEMENT_RANK, "REFINEMENT_ELEMENT_MARKER", 1); stk::mesh::put_field_on_mesh(myElementMarkerField.field(), meta.universal_part(), 1, 1, nullptr); @@ -378,9 +395,9 @@ std::pair KrinoRefinement::get_marked_element_counts() const for (int i = 0; i < length; ++i) { - if (marker[i] == Refinement::COARSEN) + if (marker[i] == static_cast(Refinement::RefinementMarker::COARSEN)) ++numUnrefine; - else if (marker[i] == Refinement::REFINE) + else if (marker[i] == static_cast(Refinement::RefinementMarker::REFINE)) ++numRefine; } } @@ -408,8 +425,11 @@ bool KrinoRefinement::do_refinement(const int debugLevel) krinolog << "Adapt: before refine, mesh has " << counts[0] << " nodes, " << counts[1] << " edges, " << counts[2] << " faces, " << counts[3] << " elements" << stk::diag::dendl; - const bool didMakeAnyChanges = myRefinement.do_refinement(get_marker()); - + bool didMakeAnyChanges = false; + { + stk::diag::TimeBlock timer_(myRefinementTimer); + didMakeAnyChanges = myRefinement.do_refinement(get_marker()); + } stk::mesh::comm_mesh_counts(mesh, counts); krinolog << "Adapt: after refine, mesh has " << counts[0] << " nodes, " << counts[1] diff --git a/packages/krino/krino/krino_lib/Akri_RefinementInterface.hpp b/packages/krino/krino/krino_lib/Akri_RefinementInterface.hpp index cc7438a95065..8bc04748726e 100644 --- a/packages/krino/krino/krino_lib/Akri_RefinementInterface.hpp +++ b/packages/krino/krino/krino_lib/Akri_RefinementInterface.hpp @@ -11,10 +11,12 @@ #include #include "Akri_Refinement.hpp" +#include "stk_util/diag/Timer.hpp" namespace stk { namespace mesh { class Entity; } } namespace stk { namespace mesh { class MetaData; } } namespace stk { namespace mesh { class Part; } } +namespace stk { namespace diag { class Timer; } } namespace krino { @@ -74,6 +76,7 @@ class KrinoRefinement : public RefinementInterface public: static KrinoRefinement & get(const stk::mesh::MetaData & meta); static KrinoRefinement & create(stk::mesh::MetaData & meta); + static KrinoRefinement & create(stk::mesh::MetaData & meta, stk::diag::Timer & timer); static bool is_created(const stk::mesh::MetaData & meta); static void register_parts_and_fields_via_aux_meta_for_fmwk(stk::mesh::MetaData & meta); @@ -108,14 +111,15 @@ class KrinoRefinement : public RefinementInterface void set_marker_field(const std::string & markerFieldName); private: - KrinoRefinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart, const bool force64Bit, const bool assert32Bit); + KrinoRefinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart, const bool force64Bit, const bool assert32Bit, stk::diag::Timer & parent_timer); std::pair get_marked_element_counts() const; TransitionElementEdgeMarker & setup_marker() const; TransitionElementEdgeMarker & get_marker() const; stk::mesh::MetaData & myMeta; - Refinement myRefinement; FieldRef myElementMarkerField; mutable std::unique_ptr myMarker; + mutable stk::diag::Timer myRefinementTimer; + Refinement myRefinement; }; void check_leaf_children_have_parents_on_same_proc(const stk::ParallelMachine comm, const RefinementInterface & refinement); diff --git a/packages/krino/krino/krino_lib/Akri_Snap.cpp b/packages/krino/krino/krino_lib/Akri_Snap.cpp index b240240944c1..addc5104df26 100644 --- a/packages/krino/krino/krino_lib/Akri_Snap.cpp +++ b/packages/krino/krino/krino_lib/Akri_Snap.cpp @@ -822,7 +822,6 @@ static void fill_interpolation_nodes_and_weights_at_node_location_in_previous_co static std::vector build_interpolation_points_for_snapping(const stk::mesh::BulkData & mesh, const stk::mesh::Part & activePart, const FieldRef coordsField, const FieldRef cdfemSnapField, const std::vector & snapNodes) { stk::mesh::Entity containingElem; - stk::math::Vector3d containingElementParametricCoords; std::vector interpNodes; std::vector interpWeights; @@ -1057,7 +1056,7 @@ void update_intersection_points_and_snap_infos_after_snap_iteration(const stk::m std::vector & intersectionPoints, std::vector & snapInfos) { - const std::vector oldToNewIntPts = update_intersection_points_after_snap_iteration(mesh, geometry, iterationSortedSnapNodes, nodesToCapturedDomains, intersectionPoints); + const std::vector oldToNewIntPts = update_intersection_points_after_snap_iteration(mesh, elementSelector, geometry, iterationSortedSnapNodes, nodesToCapturedDomains, intersectionPoints); const std::vector sortedIdsOfNodesThatNeedNewSnapInfos = get_sorted_ids_of_owned_nodes_of_elements_of_nodes(mesh, elementSelector, iterationSortedSnapNodes); @@ -1092,7 +1091,7 @@ NodeToCapturedDomainsMap snap_as_much_as_possible_while_maintaining_quality(cons std::vector intersectionPoints; geometry.store_phase_for_uncut_elements(mesh); - intersectionPoints = build_all_intersection_points(mesh, geometry, nodesToCapturedDomains); + intersectionPoints = build_all_intersection_points(mesh, elementSelector, geometry, nodesToCapturedDomains); std::vector snapInfos = build_snap_infos_from_intersection_points(mesh, sharpFeatureInfo.get(), elementSelector, nodesToCapturedDomains, intersectionPoints, qualityMetric, minIntPtWeightForEstimatingCutQuality, maxSnapForEdges, globalIDsAreParallelConsistent); while (true) diff --git a/packages/krino/krino/krino_lib/Akri_SnapToNode.cpp b/packages/krino/krino/krino_lib/Akri_SnapToNode.cpp index 3ed6212987b4..17ba7fa2ac13 100644 --- a/packages/krino/krino/krino_lib/Akri_SnapToNode.cpp +++ b/packages/krino/krino/krino_lib/Akri_SnapToNode.cpp @@ -44,11 +44,12 @@ void determine_node_snapping_from_intersection_points(const stk::mesh::BulkData } void snap_to_node(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & interfaceGeometry, const CDFEM_Snapper & snapper, NodeToCapturedDomainsMap & nodesToCapturedDomains) { - const std::vector intersectionPoints = build_uncaptured_intersection_points(mesh, interfaceGeometry, nodesToCapturedDomains); + const std::vector intersectionPoints = build_uncaptured_intersection_points(mesh, elementSelector, interfaceGeometry, nodesToCapturedDomains); determine_node_snapping_from_intersection_points(mesh, intersectionPoints, snapper, nodesToCapturedDomains); communicate_node_captured_domains_for_all_nodes(mesh, nodesToCapturedDomains); diff --git a/packages/krino/krino/krino_lib/Akri_SnapToNode.hpp b/packages/krino/krino/krino_lib/Akri_SnapToNode.hpp index dcd106704e8d..7ea1fcfd2616 100644 --- a/packages/krino/krino/krino_lib/Akri_SnapToNode.hpp +++ b/packages/krino/krino/krino_lib/Akri_SnapToNode.hpp @@ -19,6 +19,7 @@ class InterfaceGeometry; class Phase_Support; void snap_to_node(const stk::mesh::BulkData & mesh, + const stk::mesh::Selector & elementSelector, const InterfaceGeometry & interfaceGeometry, const CDFEM_Snapper & snapper, NodeToCapturedDomainsMap & nodesToCapturedDomains); diff --git a/packages/krino/krino/krino_lib/CMakeLists.txt b/packages/krino/krino/krino_lib/CMakeLists.txt index 2b26e6d79d54..8824d52a1da8 100644 --- a/packages/krino/krino/krino_lib/CMakeLists.txt +++ b/packages/krino/krino/krino_lib/CMakeLists.txt @@ -5,8 +5,11 @@ TRIBITS_INCLUDE_DIRECTORIES(${${PACKAGE_NAME}_SOURCE_DIR}) TRIBITS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}) TRIBITS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) -set(KRINO_HAVE_YAML "${${PACKAGE_NAME}_ENABLE_yamlcpp}") - +TRIBITS_ADD_OPTION_AND_DEFINE(TPL_ENABLE_yamlcpp + KRINO_HAVE_YAML + "Enable support for YAML parser in ${PACKAGE_NAME}." + OFF ) + TRIBITS_ADD_OPTION_AND_DEFINE(TPL_ENABLE_ParMETIS KRINO_HAVE_PARMETIS "Enable support for ParMETIS parser in ${PACKAGE_NAME}." @@ -25,3 +28,4 @@ TRIBITS_ADD_LIBRARY( INSTALL(FILES ${HEADERS} DESTINATION ${${PROJECT_NAME}_INSTALL_INCLUDE_DIR}/krino_lib) + diff --git a/packages/krino/krino/mesh_surface/Akri_MeshSurface.cpp b/packages/krino/krino/mesh_surface/Akri_MeshSurface.cpp index b585b4586535..61b9b0fe5f0d 100644 --- a/packages/krino/krino/mesh_surface/Akri_MeshSurface.cpp +++ b/packages/krino/krino/mesh_surface/Akri_MeshSurface.cpp @@ -68,7 +68,7 @@ void Parallel_Facet_File_Reader::get_batch_size(const int local_num_facets, int } Faceted_Surface_From_File::Faceted_Surface_From_File(const std::string & surface_name, const stk::diag::Timer &parent_timer) -: Faceted_Surface(surface_name), +: Faceted_Surface(), my_timer("Facet File Reader", parent_timer), my_built_local_facets(false) { @@ -202,8 +202,7 @@ void FACSurface::read_facets(const int batch_size, const int num_facets, const s p1 = p2; p2 = tmp; } - std::unique_ptr facet = std::make_unique( points[p1], points[p2], points[p3] ); - add( std::move(facet) ); + emplace_back_3d( points[p1], points[p2], points[p3] ); } } @@ -346,8 +345,7 @@ void PLYSurface::read_facets(const int batch_size, const std::vector facet = std::make_unique( points[p1], points[p2], points[p3] ); - add( std::move(facet) ); + emplace_back_3d( points[p1], points[p2], points[p3] ); } } @@ -367,20 +365,14 @@ STLSurface::STLSurface(const std::string & surface_name, my_reader.close_file(); } -static std::unique_ptr build_facet(const std::array & pts, const int sign) +static void flip_facet_with_negative_sign(std::array & pts, const int sign) { - unsigned p0 = 0; - unsigned p1 = 1; - unsigned p2 = 2; if (sign == -1) { - // permute to flip normal direction - p1 = 2; - p2 = 1; + std::swap(pts[1][0], pts[2][0]); + std::swap(pts[1][1], pts[2][1]); + std::swap(pts[1][2], pts[2][2]); } - - std::unique_ptr facet = std::make_unique( pts[p0], pts[p1], pts[p2] ); - return facet; } static bool is_finite_normal(const stk::math::Vector3d normal) @@ -397,7 +389,7 @@ static bool is_finite_normal(const stk::math::Vector3d normal) return finiteFlag; } -static std::unique_ptr read_ascii_facet(std::ifstream & input, const stk::math::Vector3d & scale, const int sign) +void read_ascii_facet(std::ifstream & input, const stk::math::Vector3d & scale, const int sign, std::array & facetPts, bool & isFacetValid) { std::string symbol, nx, ny, nz; double X, Y, Z; @@ -415,13 +407,12 @@ static std::unique_ptr read_ascii_facet(std::ifstream & input, const stk: STK_ThrowRequire(symbol.compare("loop") == 0); // Read in the vertices - std::array pts; for (int i=0; i<3; i++) { input >> symbol; STK_ThrowRequire(symbol.compare("vertex") == 0); input >> X >> Y >> Z; - pts[i] = stk::math::Vector3d(X*scale[0], Y*scale[1], Z*scale[2]); + facetPts[i] = stk::math::Vector3d(X*scale[0], Y*scale[1], Z*scale[2]); } // Read in the "endloop" and "endfacet" lines @@ -430,15 +421,8 @@ static std::unique_ptr read_ascii_facet(std::ifstream & input, const stk: input >> symbol; STK_ThrowRequire(symbol.compare("endfacet") == 0); - std::unique_ptr facet; - if (is_finite_normal(normal)) - { - facet = build_facet(pts, sign); - if (facet->degenerate()) - facet.reset(); - } - - return facet; + isFacetValid = is_finite_normal(normal) && !Facet3d::is_degenerate(facetPts); + flip_facet_with_negative_sign(facetPts, sign); } static bool read_start_of_next_ascii_facet(std::ifstream & input) @@ -464,12 +448,11 @@ static stk::math::Vector3d read_binary_vector(std::ifstream& input) return stk::math::Vector3d(x, y, z); } -static std::unique_ptr read_binary_facet(std::ifstream & input, const stk::math::Vector3d & scale, const int sign) +static void read_binary_facet(std::ifstream & input, const stk::math::Vector3d & scale, const int sign, std::array & facetPts, bool & isFacetValid) { const stk::math::Vector3d normal = read_binary_vector(input); - std::array pts; - for (auto && pt : pts) + for (auto && pt : facetPts) { const stk::math::Vector3d vec = read_binary_vector(input); pt = stk::math::Vector3d(vec[0]*scale[0], vec[1]*scale[1], vec[2]*scale[2]); @@ -478,15 +461,8 @@ static std::unique_ptr read_binary_facet(std::ifstream & input, const stk char dummy[2]; input.read(dummy, 2); - std::unique_ptr facet; - if (is_finite_normal(normal)) - { - facet = build_facet(pts, sign); - if (facet->degenerate()) - facet.reset(); - } - - return facet; + isFacetValid = is_finite_normal(normal) && !Facet3d::is_degenerate(facetPts); + flip_facet_with_negative_sign(facetPts, sign); } BoundingBox @@ -603,14 +579,17 @@ unsigned STLSurface::read_ascii_facets(const unsigned max_batch_size) { std::ifstream & input = my_reader.input(); + bool isFacetValid = false; + std::array facetPts; + unsigned count = 0; bool done = false; while (!done) { - std::unique_ptr facet = read_ascii_facet(input, my_scale, my_dist_sign); - if (facet) + read_ascii_facet(input, my_scale, my_dist_sign, facetPts, isFacetValid); + if (isFacetValid) { - add(std::move(facet)); + emplace_back_3d(facetPts[0], facetPts[1], facetPts[2]); ++count; } done = (!read_start_of_next_ascii_facet(input) || count >= max_batch_size); @@ -622,14 +601,17 @@ unsigned STLSurface::read_binary_facets(const unsigned maxBatchSize, unsigned & { std::ifstream & input = my_reader.input(); + bool isFacetValid = false; + std::array facetPts; + unsigned count = 0; bool done = numRemainingInFile == 0; while (!done) { - std::unique_ptr facet = read_binary_facet(input, my_scale, my_dist_sign); - if (facet) + read_binary_facet(input, my_scale, my_dist_sign, facetPts, isFacetValid); + if (isFacetValid) { - add(std::move(facet)); + emplace_back_3d(facetPts[0], facetPts[1], facetPts[2]); ++count; } done = (--numRemainingInFile == 0 || count >= maxBatchSize); @@ -642,14 +624,20 @@ BoundingBox STLSurface::get_ascii_facet_bounding_box() { std::ifstream & input = my_reader.input(); + bool isFacetValid = false; + std::array facetPts; + BoundingBox bbox; bool done = false; while (!done) { - std::unique_ptr facet = read_ascii_facet(input, my_scale, my_dist_sign); - if (facet) - facet->insert_into(bbox); + read_ascii_facet(input, my_scale, my_dist_sign, facetPts, isFacetValid); + if (isFacetValid) + { + const Facet3d facet(facetPts[0], facetPts[1], facetPts[2]); + facet.insert_into(bbox); + } done = !read_start_of_next_ascii_facet(input); } return bbox; @@ -659,13 +647,19 @@ BoundingBox STLSurface::get_binary_facet_bounding_box(const unsigned numFacets) { std::ifstream & input = my_reader.input(); + bool isFacetValid = false; + std::array facetPts; + BoundingBox bbox; for (unsigned count=0; count facet = read_binary_facet(input, my_scale, my_dist_sign); - if (facet) - facet->insert_into(bbox); + read_binary_facet(input, my_scale, my_dist_sign, facetPts, isFacetValid); + if (isFacetValid) + { + const Facet3d facet(facetPts[0], facetPts[1], facetPts[2]); + facet.insert_into(bbox); + } } return bbox; } @@ -827,8 +821,7 @@ EXOSurface::read_file(const std::vector & proc_bboxes) p2 = 1; } - std::unique_ptr facet = std::make_unique( pt[p0], pt[p1], pt[p2] ); - add( std::move(facet) ); + emplace_back_3d( pt[p0], pt[p1], pt[p2] ); } } parallel_distribute_facets(current_batch_size, proc_bboxes); @@ -836,11 +829,22 @@ EXOSurface::read_file(const std::vector & proc_bboxes) } } -MeshSurface::MeshSurface(const stk::mesh::MetaData & meta, +std::unique_ptr build_mesh_surface(const stk::mesh::MetaData & meta, + const stk::mesh::Field& coordsField, + const stk::mesh::Selector & surfaceSelector, + const int sign) +{ + if (meta.spatial_dimension() == 2) + return std::make_unique>(meta, coordsField, surfaceSelector, sign); + return std::make_unique>(meta, coordsField, surfaceSelector, sign); +} + +template +MeshSurface::MeshSurface(const stk::mesh::MetaData & meta, const stk::mesh::Field& coord_ref, const stk::mesh::Selector & surface_selector, const int sign) - : Faceted_Surface("Mesh surface"), + : Faceted_Surface(), my_sign(sign), my_mesh_meta(meta), my_coord_ref(coord_ref), @@ -848,8 +852,8 @@ MeshSurface::MeshSurface(const stk::mesh::MetaData & meta, { } -void -MeshSurface::build_local_facets(const BoundingBox & proc_bbox) +template +void MeshSurface::build_local_facets(const BoundingBox & proc_bbox) { /* %TRACE[ON]% */ Trace trace__("krino::MeshSurface::build_local_facets()"); /* %TRACE% */ const stk::mesh::BulkData & mesh = my_mesh_meta.mesh_bulk_data(); @@ -862,7 +866,7 @@ MeshSurface::build_local_facets(const BoundingBox & proc_bbox) stk::mesh::BucketVector const& buckets = mesh.get_buckets( my_mesh_meta.side_rank(), active_locally_owned_part_selector); - clear(); + Faceted_Surface::clear(); stk::mesh::BucketVector::const_iterator ib = buckets.begin(); stk::mesh::BucketVector::const_iterator ib_end = buckets.end(); @@ -909,8 +913,8 @@ MeshSurface::build_local_facets(const BoundingBox & proc_bbox) } } -void -MeshSurface::add_facet2d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1) +template +void MeshSurface::add_facet2d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1) { STK_ThrowAssert(2 == my_mesh_meta.spatial_dimension()); stk::math::Vector3d pt[2]; @@ -938,8 +942,8 @@ MeshSurface::add_facet2d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side add_facet2d(pt[0], pt[1]); } -void -MeshSurface::add_facet3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1, unsigned p2) +template +void MeshSurface::add_facet3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1, unsigned p2) { STK_ThrowAssert(3 == my_mesh_meta.spatial_dimension()); stk::math::Vector3d pt[3]; @@ -961,8 +965,8 @@ MeshSurface::add_facet3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side add_facet3d(pt[0], pt[1], pt[2]); } -void -MeshSurface::add_quad3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1, unsigned p2, unsigned p3) +template +void MeshSurface::add_quad3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, unsigned p0, unsigned p1, unsigned p2, unsigned p3) { STK_ThrowAssert(3 == my_mesh_meta.spatial_dimension()); stk::math::Vector3d pt[5]; @@ -1005,22 +1009,22 @@ MeshSurface::add_quad3d(const stk::mesh::BulkData& mesh, stk::mesh::Entity side, } } -void -MeshSurface::add_facet2d(stk::math::Vector3d & p0, stk::math::Vector3d & p1) +template +void MeshSurface::add_facet2d(stk::math::Vector3d & p0, stk::math::Vector3d & p1) { STK_ThrowAssert(2 == my_mesh_meta.spatial_dimension()); - - std::unique_ptr facet = std::make_unique( p0, p1 ); - add( std::move(facet) ); + Faceted_Surface::emplace_back_2d( p0, p1 ); } -void -MeshSurface::add_facet3d(stk::math::Vector3d & p0, stk::math::Vector3d & p1, stk::math::Vector3d & p2) +template +void MeshSurface::add_facet3d(stk::math::Vector3d & p0, stk::math::Vector3d & p1, stk::math::Vector3d & p2) { STK_ThrowAssert(3 == my_mesh_meta.spatial_dimension()); - - std::unique_ptr facet = std::make_unique( p0, p1, p2 ); - add( std::move(facet) ); + Faceted_Surface::emplace_back_3d( p0, p1, p2 ); } +// Explicit template instantiation +template class MeshSurface; +template class MeshSurface; + } // namespace krino diff --git a/packages/krino/krino/mesh_surface/Akri_MeshSurface.hpp b/packages/krino/krino/mesh_surface/Akri_MeshSurface.hpp index 31b3a65be096..bae09d112aee 100644 --- a/packages/krino/krino/mesh_surface/Akri_MeshSurface.hpp +++ b/packages/krino/krino/mesh_surface/Akri_MeshSurface.hpp @@ -24,7 +24,8 @@ namespace stk { namespace mesh { class Selector; } } namespace krino { -class MeshSurface : public Faceted_Surface { +template +class MeshSurface : public Faceted_Surface { public: MeshSurface(const stk::mesh::MetaData & meta, const stk::mesh::Field& coord_ref, @@ -48,6 +49,11 @@ class MeshSurface : public Faceted_Surface { const stk::mesh::Selector my_surface_selector; }; +std::unique_ptr build_mesh_surface(const stk::mesh::MetaData & meta, + const stk::mesh::Field& coordsField, + const stk::mesh::Selector & surfaceSelector, + const int sign); + class Parallel_Facet_File_Reader { public: Parallel_Facet_File_Reader(const std::string & in_filename) : my_filename(in_filename) {} @@ -66,7 +72,7 @@ class Parallel_Facet_File_Reader { std::ifstream my_input; }; -class Faceted_Surface_From_File : public Faceted_Surface { +class Faceted_Surface_From_File : public Faceted_Surface { public: Faceted_Surface_From_File(const std::string & surface_name, const stk::diag::Timer &parent_timer); virtual ~Faceted_Surface_From_File() {} diff --git a/packages/krino/krino/mesh_utils/Akri_AllReduce.hpp b/packages/krino/krino/mesh_utils/Akri_AllReduce.hpp new file mode 100644 index 000000000000..75f31057285c --- /dev/null +++ b/packages/krino/krino/mesh_utils/Akri_AllReduce.hpp @@ -0,0 +1,66 @@ +#ifndef KRINO_KRINO_MESH_UTILS_AKRI_ALLREDUCE_HPP_ +#define KRINO_KRINO_MESH_UTILS_AKRI_ALLREDUCE_HPP_ + +#include +#include + +namespace krino { + +template std::array::type, sizeof...(Args)> make_array(Args&&... args) +{ + return {{ args... }}; +} + +template +std::array all_reduce_sum_array(stk::ParallelMachine comm, const std::array & local) +{ + std::array global; + stk::all_reduce_sum(comm, local.data(), global.data(), local.size()); + return global; +} + +template +std::array all_reduce_min_array(stk::ParallelMachine comm, const std::array & local) +{ + std::array global; + stk::all_reduce_min(comm, local.data(), global.data(), local.size()); + return global; +} + +template +std::array all_reduce_max_array(stk::ParallelMachine comm, const std::array & local) +{ + std::array global; + stk::all_reduce_max(comm, local.data(), global.data(), local.size()); + return global; +} + +template +void all_reduce_sum(stk::ParallelMachine comm, Args&&... args) +{ + const auto global = all_reduce_sum_array(comm, make_array(args...)); + const auto * data = global.data(); + ((std::forward(args) = *(data++)), ...); +} + +template +void all_reduce_min(stk::ParallelMachine comm, Args&&... args) +{ + const auto global = all_reduce_min_array(comm, make_array(args...)); + const auto * data = global.data(); + ((std::forward(args) = *(data++)), ...); +} + +template +void all_reduce_max(stk::ParallelMachine comm, Args&&... args) +{ + const auto global = all_reduce_max_array(comm, make_array(args...)); + const auto * data = global.data(); + ((std::forward(args) = *(data++)), ...); +} + +} + + + +#endif /* KRINO_KRINO_MESH_UTILS_AKRI_ALLREDUCE_HPP_ */ diff --git a/packages/krino/krino/mesh_utils/Akri_Edge.cpp b/packages/krino/krino/mesh_utils/Akri_Edge.cpp index 93ef49375456..0ba1b2210fee 100644 --- a/packages/krino/krino/mesh_utils/Akri_Edge.cpp +++ b/packages/krino/krino/mesh_utils/Akri_Edge.cpp @@ -105,4 +105,26 @@ std::vector get_edges_of_selected_elements(const stk::mesh::BulkData & mes return edges; } +std::vector get_edges_of_elements(const stk::mesh::BulkData & mesh, const std::vector & elements) +{ + std::vector edges; + + if (!elements.empty()) + { + const stk::topology elemTopology = mesh.bucket(elements[0]).topology(); // Assume all elements have same topology and check later in debug + const size_t edgeCount = elements.size() * elemTopology.num_edges(); + edges.reserve(edgeCount); + + for (const auto elem : elements) + { + STK_ThrowAssert(elemTopology == mesh.bucket(elem).topology()); + append_entity_edges(mesh, elemTopology, elem, edges); + } + + stk::util::sort_and_unique(edges); + } + + return edges; +} + } diff --git a/packages/krino/krino/mesh_utils/Akri_Edge.hpp b/packages/krino/krino/mesh_utils/Akri_Edge.hpp index 7363b6941760..0dff58f15240 100644 --- a/packages/krino/krino/mesh_utils/Akri_Edge.hpp +++ b/packages/krino/krino/mesh_utils/Akri_Edge.hpp @@ -60,6 +60,7 @@ int get_edge_parallel_owner_rank(const stk::mesh::BulkData & mesh, const Edge ed std::string debug_edge(const stk::mesh::BulkData & mesh, const Edge edge); std::vector get_edges_of_selected_elements(const stk::mesh::BulkData & mesh, const stk::mesh::Selector & elementSelector); +std::vector get_edges_of_elements(const stk::mesh::BulkData & mesh, const std::vector & elements); } diff --git a/packages/krino/krino/mesh_utils/Akri_MeshHelpers.cpp b/packages/krino/krino/mesh_utils/Akri_MeshHelpers.cpp index 305c7289045d..a62f2c134d65 100644 --- a/packages/krino/krino/mesh_utils/Akri_MeshHelpers.cpp +++ b/packages/krino/krino/mesh_utils/Akri_MeshHelpers.cpp @@ -47,7 +47,7 @@ stk::mesh::PartVector get_all_block_parts(const stk::mesh::MetaData & meta) return blockParts; } -static double * get_field_data(const stk::mesh::BulkData& mesh, const FieldRef field, const stk::mesh::Entity entity) +double * get_field_data(const stk::mesh::BulkData& mesh, const FieldRef field, const stk::mesh::Entity entity) { STK_ThrowRequireMsg(field.valid(), "Invalid field: " << field.name()); double * fieldData = field_data(field, entity); @@ -55,6 +55,11 @@ static double * get_field_data(const stk::mesh::BulkData& mesh, const FieldRef f return fieldData; } +double & get_scalar_field(const stk::mesh::BulkData& mesh, const FieldRef field, const stk::mesh::Entity entity) +{ + return *get_field_data(mesh, field, entity); +} + stk::math::Vector3d get_vector_field(const stk::mesh::BulkData& mesh, const FieldRef vecField, const stk::mesh::Entity entity) { return stk::math::Vector3d(get_field_data(mesh, vecField, entity)); @@ -500,10 +505,10 @@ void pack_entities_for_sharing_procs(const stk::mesh::BulkData & mesh, }); } -void unpack_shared_entities(const stk::mesh::BulkData & mesh, - std::vector & sharedEntities, +std::vector unpack_entities_from_other_procs(const stk::mesh::BulkData & mesh, stk::CommSparse &commSparse) { + std::vector entities; stk::unpack_communications(commSparse, [&](int procId) { stk::CommBuffer & buffer = commSparse.recv_buffer(procId); @@ -514,9 +519,10 @@ void unpack_shared_entities(const stk::mesh::BulkData & mesh, commSparse.recv_buffer(procId).unpack(entityKey); stk::mesh::Entity entity = mesh.get_entity(entityKey); STK_ThrowAssert(mesh.is_valid(entity)); - sharedEntities.push_back(entity); + entities.push_back(entity); } }); + return entities; } static @@ -526,8 +532,7 @@ void append_shared_entities_to_owned_ones(const stk::mesh::BulkData & mesh, stk::CommSparse commSparse(mesh.parallel()); pack_entities_for_sharing_procs(mesh, entities, commSparse); - std::vector sharedEntities; - unpack_shared_entities(mesh, sharedEntities, commSparse); + const std::vector sharedEntities = unpack_entities_from_other_procs(mesh, commSparse); entities.insert(entities.end(), sharedEntities.begin(), sharedEntities.end()); } @@ -2637,4 +2642,80 @@ void communicate_owned_entities_to_ghosting_procs(const stk::mesh::BulkData & me unpack_ghosted_entities_from_owners(mesh, entities, commSparse); } +template +void pack_shared_nodes_for_sharing_procs(const stk::mesh::BulkData & mesh, + const NODE_CONTAINER & nodes, + stk::CommSparse &commSparse) +{ + std::vector nodeSharedProcs; + stk::pack_and_communicate(commSparse,[&]() + { + for (auto node : nodes) + { + if (mesh.bucket(node).shared()) + { + mesh.comm_shared_procs(node, nodeSharedProcs); + for (int procId : nodeSharedProcs) + commSparse.send_buffer(procId).pack(mesh.identifier(node)); + } + } + }); +} + +static +void unpack_shared_nodes(const stk::mesh::BulkData & mesh, + std::set & nodes, + stk::CommSparse &commSparse) +{ + stk::unpack_communications(commSparse, [&](int procId) + { + stk::CommBuffer & buffer = commSparse.recv_buffer(procId); + + while ( buffer.remaining() ) + { + stk::mesh::EntityId nodeId; + commSparse.recv_buffer(procId).unpack(nodeId); + stk::mesh::Entity node = mesh.get_entity(stk::topology::NODE_RANK, nodeId); + STK_ThrowRequire(mesh.is_valid(node)); + nodes.insert(node); + } + }); +} + +static +void unpack_shared_nodes(const stk::mesh::BulkData & mesh, + std::vector & nodes, + stk::CommSparse &commSparse) +{ + stk::unpack_communications(commSparse, [&](int procId) + { + stk::CommBuffer & buffer = commSparse.recv_buffer(procId); + + while ( buffer.remaining() ) + { + stk::mesh::EntityId nodeId; + commSparse.recv_buffer(procId).unpack(nodeId); + stk::mesh::Entity node = mesh.get_entity(stk::topology::NODE_RANK, nodeId); + STK_ThrowRequire(mesh.is_valid(node)); + nodes.push_back(node); + } + }); +} + +void communicate_shared_nodes_to_sharing_procs(const stk::mesh::BulkData & mesh, std::set & nodes) +{ + stk::CommSparse commSparse(mesh.parallel()); + pack_shared_nodes_for_sharing_procs(mesh, nodes, commSparse); + unpack_shared_nodes(mesh, nodes, commSparse); +} + +void communicate_shared_nodes_to_sharing_procs_and_sort_and_unique(const stk::mesh::BulkData & mesh, std::vector & nodes) +{ + stk::CommSparse commSparse(mesh.parallel()); + pack_shared_nodes_for_sharing_procs(mesh, nodes, commSparse); + unpack_shared_nodes(mesh, nodes, commSparse); + + stk::util::sort_and_unique(nodes, stk::mesh::EntityLess(mesh)); +} + } // namespace krino diff --git a/packages/krino/krino/mesh_utils/Akri_MeshHelpers.hpp b/packages/krino/krino/mesh_utils/Akri_MeshHelpers.hpp index 5514ec00e814..3311f1d80a15 100644 --- a/packages/krino/krino/mesh_utils/Akri_MeshHelpers.hpp +++ b/packages/krino/krino/mesh_utils/Akri_MeshHelpers.hpp @@ -46,6 +46,8 @@ struct StkMeshEntities void fill_node_ids_for_nodes(const stk::mesh::BulkData & mesh, const std::vector & parentNodes, std::vector & parentNodeIds); stk::mesh::PartVector get_all_block_parts(const stk::mesh::MetaData & meta); +double * get_field_data(const stk::mesh::BulkData& mesh, const FieldRef field, const stk::mesh::Entity entity); +double & get_scalar_field(const stk::mesh::BulkData& mesh, const FieldRef field, const stk::mesh::Entity entity); stk::math::Vector3d get_vector_field(const stk::mesh::BulkData& mesh, const FieldRef vecField, const stk::mesh::Entity entity); stk::math::Vector3d get_vector_field(const stk::mesh::BulkData& mesh, const FieldRef vecField, const stk::mesh::Entity entity, const unsigned vecLen); bool is_less_than_in_x_then_y_then_z(const stk::math::Vector3d& A, const stk::math::Vector3d &B); @@ -74,7 +76,7 @@ void attach_entity_to_element(stk::mesh::BulkData & mesh, const stk::mesh::Entit void attach_entity_to_elements(stk::mesh::BulkData & mesh, stk::mesh::Entity entity); void unpack_entities_from_other_procs(const stk::mesh::BulkData & mesh, std::set & entities, stk::CommSparse &commSparse); void pack_entities_for_sharing_procs(const stk::mesh::BulkData & mesh, const std::vector & entities, stk::CommSparse &commSparse); -void unpack_shared_entities(const stk::mesh::BulkData & mesh, std::vector & sharedEntities, stk::CommSparse &commSparse); +std::vector unpack_entities_from_other_procs(const stk::mesh::BulkData & mesh, stk::CommSparse &commSparse); void activate_selected_entities_touching_active_elements(stk::mesh::BulkData & mesh, const stk::mesh::EntityRank entityRank, const stk::mesh::Selector & entitySelector, stk::mesh::Part & activePart); void activate_all_entities(stk::mesh::BulkData & mesh, stk::mesh::Part & active_part); void destroy_custom_ghostings(stk::mesh::BulkData & mesh); @@ -163,6 +165,8 @@ void batch_create_sides(stk::mesh::BulkData & mesh, const std::vector< SideDescr void make_side_ids_consistent_with_stk_convention(stk::mesh::BulkData & mesh); void communicate_owned_entities_to_ghosting_procs(const stk::mesh::BulkData & mesh, std::vector & entities); +void communicate_shared_nodes_to_sharing_procs(const stk::mesh::BulkData & mesh, std::set & nodes); +void communicate_shared_nodes_to_sharing_procs_and_sort_and_unique(const stk::mesh::BulkData & mesh, std::vector & nodes); double compute_element_volume_to_edge_ratio(stk::mesh::BulkData & mesh, stk::mesh::Entity element, const stk::mesh::Field * const coords_field); diff --git a/packages/krino/krino/parser/Akri_IC_Parser.cpp b/packages/krino/krino/parser/Akri_IC_Parser.cpp index 6829f7bed5ef..a32d8ed34ef5 100644 --- a/packages/krino/krino/parser/Akri_IC_Parser.cpp +++ b/packages/krino/krino/parser/Akri_IC_Parser.cpp @@ -13,7 +13,6 @@ #include #include #include -#include #include #include diff --git a/packages/krino/krino/parser/Akri_Surface_Parser.cpp b/packages/krino/krino/parser/Akri_Surface_Parser.cpp index 9e5f288cedd3..655466e3b70a 100644 --- a/packages/krino/krino/parser/Akri_Surface_Parser.cpp +++ b/packages/krino/krino/parser/Akri_Surface_Parser.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include @@ -257,7 +258,7 @@ parse_facets(const Parser::Node & ic_node, const stk::diag::Timer &parent_timer) return surface; } -MeshSurface * +Surface * parse_mesh_surface(const Parser::Node & ic_node, const stk::mesh::MetaData & meta) { std::string surface_name; @@ -279,8 +280,11 @@ parse_mesh_surface(const Parser::Node & ic_node, const stk::mesh::MetaData & met const stk::mesh::Field* coords = reinterpret_cast*>(&auxMeta.get_current_coordinates().field()); STK_ThrowRequire(nullptr != coords); - const stk::mesh::Selector surface_selector = stk::mesh::Selector(io_part); - return new MeshSurface(meta, *coords, surface_selector, sign); + const stk::mesh::Selector surfaceSelector = stk::mesh::Selector(io_part); + + if (meta.spatial_dimension() == 2) + return new MeshSurface(meta, *coords, surfaceSelector, sign); + return new MeshSurface(meta, *coords, surfaceSelector, sign); } LevelSet_String_Function * diff --git a/packages/krino/krino/refinement/Akri_Refinement.cpp b/packages/krino/krino/refinement/Akri_Refinement.cpp index c00485e901ac..daa873dd8a48 100644 --- a/packages/krino/krino/refinement/Akri_Refinement.cpp +++ b/packages/krino/krino/refinement/Akri_Refinement.cpp @@ -58,26 +58,34 @@ void Refinement::declare_refinement_fields() stk::mesh::put_field_on_mesh(*myOriginatingProcForParentElementField, myMeta.universal_part(), 1, nullptr); // needed everywhere for restart, otherwise could be parent_part } -Refinement::Refinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart, const bool force64Bit, const bool assert32Bit) - : myMeta(meta), - myForce64Bit(force64Bit), - myAssert32Bit(assert32Bit), - myNodeRefiner(force64Bit, assert32Bit), - myEntityIdPool(meta), - myActivePart(activePart) +Refinement::Refinement(stk::mesh::MetaData & meta, + stk::mesh::Part * activePart, + const bool force64Bit, + const bool assert32Bit, + stk::diag::Timer & parentTimer) + : myMeta(meta), + myForce64Bit(force64Bit), + myAssert32Bit(assert32Bit), + myNodeRefiner(force64Bit, assert32Bit), + myEntityIdPool(meta), + myActivePart(activePart), + refineTimer(parentTimer), + unrefineTimer(parentTimer), + myFixPartsandOwnersTimer("Fix Parts and Owners", parentTimer) { myCoordsField = static_cast*>(myMeta.coordinate_field()); declare_refinement_parts(); declare_refinement_fields(); } -Refinement::Refinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart) - : Refinement(meta, activePart, false, false) +Refinement::Refinement( + stk::mesh::MetaData & meta, stk::mesh::Part * activePart, stk::diag::Timer & parentTimer) + : Refinement(meta, activePart, false, false, parentTimer) { } -Refinement::Refinement(stk::mesh::MetaData & meta) - : Refinement(meta, nullptr, false, false) +Refinement::Refinement(stk::mesh::MetaData & meta, stk::diag::Timer & parentTimer) + : Refinement(meta, nullptr, false, false, parentTimer) { } @@ -906,28 +914,41 @@ void Refinement::update_element_rebalance_weights_incorporating_parallel_owner_c void Refinement::create_refined_nodes_elements_and_sides(const EdgeMarkerInterface & edgeMarker) { stk::mesh::BulkData & mesh = myMeta.mesh_bulk_data(); + std::vector sideRequests; - mesh.modification_begin(); + { + mesh.modification_begin(); - destroy_custom_ghostings(); + // Mesh modification cycle + destroy_custom_ghostings(); - myNodeRefiner.create_refined_edge_nodes(mesh, get_parts_for_new_refined_edge_nodes(), myRefinedEdgeNodeParentIdsField); + myNodeRefiner.create_refined_edge_nodes( + mesh, get_parts_for_new_refined_edge_nodes(), myRefinedEdgeNodeParentIdsField); - const std::vector bucketsData = get_buckets_data_for_candidate_elements_to_refine(edgeMarker); - const size_t numNewElements = count_new_child_elements(edgeMarker, bucketsData); - myEntityIdPool.reserve(stk::topology::ELEMENT_RANK, numNewElements, myAssert32Bit, myForce64Bit); + const std::vector bucketsData = + get_buckets_data_for_candidate_elements_to_refine(edgeMarker); + const size_t numNewElements = count_new_child_elements(edgeMarker, bucketsData); + myEntityIdPool.reserve( + stk::topology::ELEMENT_RANK, numNewElements, myAssert32Bit, myForce64Bit); - std::vector elementsToDelete; - std::vector sideRequests; - refine_elements_with_refined_edges_and_store_sides_to_create(edgeMarker, bucketsData, sideRequests, elementsToDelete); - stk::mesh::destroy_elements_no_mod_cycle(mesh, elementsToDelete, mesh.mesh_meta_data().universal_part()); + std::vector elementsToDelete; + refine_elements_with_refined_edges_and_store_sides_to_create( + edgeMarker, bucketsData, sideRequests, elementsToDelete); + stk::mesh::destroy_elements_no_mod_cycle( + mesh, elementsToDelete, mesh.mesh_meta_data().universal_part()); - mesh.modification_end(); + mesh.modification_end(); + } + //timer batch create sides if(stk::is_true_on_any_proc(mesh.parallel(), !sideRequests.empty())) + { batch_create_sides(mesh, sideRequests); + } - myNodeRefiner.prolong_refined_edge_nodes(mesh); + { + myNodeRefiner.prolong_refined_edge_nodes(mesh); + } } void Refinement::create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(const EdgeMarkerInterface & edgeMarker) @@ -1029,81 +1050,145 @@ void Refinement::destroy_custom_ghostings() krino::destroy_custom_ghostings(myMeta.mesh_bulk_data()); } -bool Refinement::do_unrefinement(const EdgeMarkerInterface & edgeMarker) +void Refinement::finalize() { + stk::diag::TimeBlock timer_(myFixPartsandOwnersTimer); stk::mesh::BulkData & mesh = myMeta.mesh_bulk_data(); + activate_selected_entities_touching_active_elements( + mesh, myMeta.side_rank(), myMeta.universal_part(), *myActivePart); + fix_node_owners_to_assure_active_owned_element_for_node(mesh, *myActivePart); +} - check_leaf_children_have_parents_on_same_proc(); +bool Refinement::refine_elements(const EdgeMarkerInterface & edgeMarker) +{ + stk::diag::TimeBlock timer_(refineTimer.rootTimer); + + stk::mesh::BulkData & mesh = myMeta.mesh_bulk_data(); + + { + stk::diag::TimeBlock timer_2(refineTimer.checkLeafChildren); + check_leaf_children_have_parents_on_same_proc(); + } + + { + stk::diag::TimeBlock timer_2(refineTimer.findEdgesToRefine); + find_edges_to_refine(edgeMarker); + } + bool didMakeAnyChanges = false; + + const bool haveEdgesToRefineLocally = get_num_edges_to_refine() > 0; + if (stk::is_true_on_any_proc(mesh.parallel(), haveEdgesToRefineLocally)) + { + didMakeAnyChanges = true; + { + stk::diag::TimeBlock timer_2(refineTimer.createRefinedEdges); + create_refined_nodes_elements_and_sides(edgeMarker); + } + + { + stk::diag::TimeBlock timer_2(refineTimer.elminateHangingNodes); + create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(edgeMarker); + } + + { + stk::diag::TimeBlock timer_2(refineTimer.prolongNodes); + myNodeRefiner.prolong_refined_edge_nodes(mesh); + } + } + + return didMakeAnyChanges; +} + +bool Refinement::unrefine_elements(const EdgeMarkerInterface & edgeMarker) +{ + stk::diag::TimeBlock timer_(unrefineTimer.rootTimer); + + stk::mesh::BulkData & mesh = myMeta.mesh_bulk_data(); + + { + stk::diag::TimeBlock timer_2(unrefineTimer.checkLeafChildren); + check_leaf_children_have_parents_on_same_proc(); + } bool didMakeAnyChanges = false; if(stk::is_true_on_any_proc(mesh.parallel(), edgeMarker.locally_have_elements_to_unrefine())) { std::vector childElementsToDeleteForUnrefinement; std::vector ownedParentElementsModifiedByUnrefinement; - edgeMarker.fill_elements_modified_by_unrefinement(ownedParentElementsModifiedByUnrefinement, childElementsToDeleteForUnrefinement); + { + stk::diag::TimeBlock timer_2(unrefineTimer.fillElements); + edgeMarker.fill_elements_modified_by_unrefinement( + ownedParentElementsModifiedByUnrefinement, childElementsToDeleteForUnrefinement); + } if(stk::is_true_on_any_proc(mesh.parallel(), !childElementsToDeleteForUnrefinement.empty())) { - for (auto parentElement : ownedParentElementsModifiedByUnrefinement) { + stk::diag::TimeBlock timer_2(unrefineTimer.restrictElementFields); + auto childElements = get_children(parentElement); restrict_element_fields(mesh, parentElement, childElements); } + std::vector originatingProcForParentsBeingModified; didMakeAnyChanges = true; +{ + stk::diag::TimeBlock timer_2(unrefineTimer.fillElements); - const std::vector originatingProcForParentsBeingModified = get_originating_procs_for_elements(ownedParentElementsModifiedByUnrefinement); + originatingProcForParentsBeingModified = + get_originating_procs_for_elements(ownedParentElementsModifiedByUnrefinement); +} + { + stk::diag::TimeBlock timer_2(unrefineTimer.meshMod); + + mesh.modification_begin(); + destroy_custom_ghostings(); + stk::mesh::destroy_elements_no_mod_cycle( + mesh, childElementsToDeleteForUnrefinement, mesh.mesh_meta_data().universal_part()); + remove_parent_parts(ownedParentElementsModifiedByUnrefinement); + mesh.modification_end(); + } - mesh.modification_begin(); - destroy_custom_ghostings(); - stk::mesh::destroy_elements_no_mod_cycle(mesh, childElementsToDeleteForUnrefinement, mesh.mesh_meta_data().universal_part()); - remove_parent_parts(ownedParentElementsModifiedByUnrefinement); - mesh.modification_end(); + { + stk::diag::TimeBlock timer_2(unrefineTimer.fixFaceEdgeOwnership); + fix_face_and_edge_ownership(mesh); + mark_already_refined_edges(); + } - fix_face_and_edge_ownership(mesh); - mark_already_refined_edges(); + { + stk::diag::TimeBlock timer_2(unrefineTimer.elminateHangingNodes); + create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(edgeMarker); + } - create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(edgeMarker); + { + stk::diag::TimeBlock timer_2(unrefineTimer.fixFaceEdgeOwnership); - respect_originating_proc_for_parents_modified_by_unrefinement(ownedParentElementsModifiedByUnrefinement, originatingProcForParentsBeingModified); + respect_originating_proc_for_parents_modified_by_unrefinement( + ownedParentElementsModifiedByUnrefinement, originatingProcForParentsBeingModified); + } } } - check_leaf_children_have_parents_on_same_proc(); + { + stk::diag::TimeBlock timer_2(unrefineTimer.checkLeafChildren); + check_leaf_children_have_parents_on_same_proc(); + } return didMakeAnyChanges; } bool Refinement::do_refinement(const EdgeMarkerInterface & edgeMarker) { - stk::mesh::BulkData & mesh = myMeta.mesh_bulk_data(); - - check_leaf_children_have_parents_on_same_proc(); - - find_edges_to_refine(edgeMarker); - bool didMakeAnyChanges = false; - const bool haveEdgesToRefineLocally = get_num_edges_to_refine() > 0; - if(stk::is_true_on_any_proc(mesh.parallel(), haveEdgesToRefineLocally)) - { - didMakeAnyChanges = true; - - create_refined_nodes_elements_and_sides(edgeMarker); - - create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(edgeMarker); - - myNodeRefiner.prolong_refined_edge_nodes(mesh); - } - - didMakeAnyChanges |= do_unrefinement(edgeMarker); + didMakeAnyChanges |= refine_elements(edgeMarker); + didMakeAnyChanges |= unrefine_elements(edgeMarker); if (didMakeAnyChanges && myActivePart) { - activate_selected_entities_touching_active_elements(mesh, myMeta.side_rank(), myMeta.universal_part(), *myActivePart); - fix_node_owners_to_assure_active_owned_element_for_node(mesh, *myActivePart); + finalize(); } STK_ThrowAssertMsg(!have_any_hanging_refined_nodes(), "Mesh has hanging refined node."); diff --git a/packages/krino/krino/refinement/Akri_Refinement.hpp b/packages/krino/krino/refinement/Akri_Refinement.hpp index eca86472fd80..377205704bad 100644 --- a/packages/krino/krino/refinement/Akri_Refinement.hpp +++ b/packages/krino/krino/refinement/Akri_Refinement.hpp @@ -8,6 +8,7 @@ #include #include #include "Akri_NodeRefiner.hpp" +#include "stk_util/diag/Timer.hpp" namespace stk { namespace mesh { class MetaData; } } namespace stk { class topology; } @@ -20,15 +21,60 @@ class EdgeMarkerInterface; class Refinement { public: - enum RefinementMarker + enum class RefinementMarker { COARSEN = -1, NOTHING = 0, REFINE = 1 }; - Refinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart, const bool force64Bit, const bool assert32Bit); - Refinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart); - Refinement(stk::mesh::MetaData & meta); + + struct RefineElementsTimers + { + RefineElementsTimers(stk::diag::Timer & parent_timer) + : rootTimer("Refine Elements", parent_timer), + checkLeafChildren("Check Leaf Children", rootTimer), + findEdgesToRefine("Find Edges to Refine", rootTimer), + createRefinedEdges("Create Refined Edges", rootTimer), + elminateHangingNodes("Eliminate Hanging Nodes", rootTimer), + prolongNodes("Prolong Edge Nodes", rootTimer) + {}; + + mutable stk::diag::Timer rootTimer; + mutable stk::diag::Timer checkLeafChildren; + mutable stk::diag::Timer findEdgesToRefine; + mutable stk::diag::Timer createRefinedEdges; + mutable stk::diag::Timer elminateHangingNodes; + mutable stk::diag::Timer prolongNodes; + }; + + struct UnrefineElementsTimers + { + UnrefineElementsTimers(stk::diag::Timer & parent_timer) + : rootTimer("Unrefine Elements", parent_timer), + checkLeafChildren("Check Leaf Children", rootTimer), + restrictElementFields("Restrict Element Fields", rootTimer), + fillElements("Fill elements", rootTimer), + meshMod("Mesh Modification", rootTimer), + elminateHangingNodes("Eliminate Hanging Nodes", rootTimer), + fixFaceEdgeOwnership("Fix ownership and mark", rootTimer) + {}; + + mutable stk::diag::Timer rootTimer; + mutable stk::diag::Timer checkLeafChildren; + mutable stk::diag::Timer restrictElementFields; + mutable stk::diag::Timer fillElements; + mutable stk::diag::Timer meshMod; + mutable stk::diag::Timer elminateHangingNodes; + mutable stk::diag::Timer fixFaceEdgeOwnership; + }; + + Refinement(stk::mesh::MetaData & meta, + stk::mesh::Part * activePart, + const bool force64Bit, + const bool assert32Bit, + stk::diag::Timer & parentTimer); + Refinement(stk::mesh::MetaData & meta, stk::mesh::Part * activePart,stk::diag::Timer & parentTimer); + Refinement(stk::mesh::MetaData & meta, stk::diag::Timer & parentTimer); Refinement ( const Refinement & ) = delete; Refinement & operator= ( const Refinement & ) = delete; @@ -100,7 +146,9 @@ class Refinement bool locally_have_any_hanging_refined_nodes() const; void create_refined_nodes_elements_and_sides(const EdgeMarkerInterface & edgeMarker); void create_another_layer_of_refined_elements_and_sides_to_eliminate_hanging_nodes(const EdgeMarkerInterface & edgeMarker); - bool do_unrefinement(const EdgeMarkerInterface & edgeMarker); + bool unrefine_elements(const EdgeMarkerInterface & edgeMarker); + bool refine_elements(const EdgeMarkerInterface & edgeMarker); + void finalize(); void mark_already_refined_edges(); void destroy_custom_ghostings(); @@ -144,6 +192,10 @@ class Refinement stk::mesh::Field * myChildElementIds8Field{nullptr}; stk::mesh::Field * myRefinedEdgeNodeParentIdsField{nullptr}; stk::mesh::Field * myOriginatingProcForParentElementField{nullptr}; + + mutable RefineElementsTimers refineTimer; + mutable UnrefineElementsTimers unrefineTimer; + mutable stk::diag::Timer myFixPartsandOwnersTimer; }; } // namespace krino diff --git a/packages/krino/krino/refinement/Akri_TransitionElementEdgeMarker.cpp b/packages/krino/krino/refinement/Akri_TransitionElementEdgeMarker.cpp index 247dcb58912a..cf18ecd9329c 100644 --- a/packages/krino/krino/refinement/Akri_TransitionElementEdgeMarker.cpp +++ b/packages/krino/krino/refinement/Akri_TransitionElementEdgeMarker.cpp @@ -237,7 +237,7 @@ static bool is_any_element_marked(FieldRef elementMarkerField, const std::vector for(const auto & elem : elems) { auto * elemMarker = field_data(elementMarkerField, elem); - if (elemMarker && *elemMarker == Refinement::REFINE) + if (elemMarker && *elemMarker == static_cast(Refinement::RefinementMarker::REFINE)) return true; } return false; @@ -246,7 +246,7 @@ static bool is_any_element_marked(FieldRef elementMarkerField, const std::vector static bool is_element_marked_for_unrefinement(FieldRef elementMarkerField, const stk::mesh::Entity elem) { auto * elemMarker = field_data(elementMarkerField, elem); - if (elemMarker && *elemMarker == Refinement::COARSEN) + if (elemMarker && *elemMarker == static_cast(Refinement::RefinementMarker::COARSEN)) return true; return false; } @@ -361,7 +361,7 @@ void TransitionElementEdgeMarker::locally_mark_edges_of_marked_non_transition_el { const stk::mesh::Entity elem = bucket[i]; - if (elemMarker[i*markerFieldLength] == Refinement::REFINE && !is_transition(elem)) + if (elemMarker[i*markerFieldLength] == static_cast(Refinement::RefinementMarker::REFINE) && !is_transition(elem)) { fill_entity_edges(myMesh, elem, edgesToRefineForElement); for (auto && edgeToRefineForElement : edgesToRefineForElement) @@ -454,8 +454,7 @@ static void communicate_to_get_sorted_edges_nodes_that_will_be_removed_by_unrefi stk::CommSparse commSparse(mesh.parallel()); pack_entities_for_sharing_procs(mesh, sharedEdgeNodesThatMustBeKept, commSparse); - std::vector sharedEdgeNodesThatMustBeKeptFromOtherProcs; - unpack_shared_entities(mesh, sharedEdgeNodesThatMustBeKeptFromOtherProcs, commSparse); + std::vector sharedEdgeNodesThatMustBeKeptFromOtherProcs = unpack_entities_from_other_procs(mesh, commSparse); stk::util::remove_intersection_from_first(ownedOrSharedEdgeNodesThatMayBeUnrefined, sharedEdgeNodesThatMustBeKeptFromOtherProcs); } diff --git a/packages/krino/krino/region/Akri_Region.cpp b/packages/krino/krino/region/Akri_Region.cpp index 5f1af46d7d9a..12b8929194f7 100644 --- a/packages/krino/krino/region/Akri_Region.cpp +++ b/packages/krino/krino/region/Akri_Region.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include diff --git a/packages/krino/krino/region/Akri_Startup.cpp b/packages/krino/krino/region/Akri_Startup.cpp index 614b2df5ed45..01bafb6c0a47 100644 --- a/packages/krino/krino/region/Akri_Startup.cpp +++ b/packages/krino/krino/region/Akri_Startup.cpp @@ -266,7 +266,8 @@ void Startup::handle_exception(const char * what, const bool is_parsing) std::cout << std::flush; std::cerr << std::flush; - std::abort(); + const int errorCode{-1}; + MPI_Abort(stk::EnvData::parallel_comm(), errorCode); } void Startup::report_handler(const char * message, int type) diff --git a/packages/krino/krino/surface/Akri_AnalyticSurf.cpp b/packages/krino/krino/surface/Akri_AnalyticSurf.cpp index fcef0f8a821f..127498eea624 100644 --- a/packages/krino/krino/surface/Akri_AnalyticSurf.cpp +++ b/packages/krino/krino/surface/Akri_AnalyticSurf.cpp @@ -255,7 +255,7 @@ Ellipsoid::point_signed_distance(const stk::math::Vector3d &x) const return mySign*(mySemiAxesNorm*std::sqrt(dx*dx+dy*dy+dz*dz)-mySemiAxesNorm); } -Plane::Plane(const double normal[3], +Plane::Plane(const stk::math::Vector3d & normal, const double offset, const double multiplier) : SurfaceThatDoesntTakeAdvantageOfNarrowBandAndThereforeHasCorrectSign(), diff --git a/packages/krino/krino/surface/Akri_AnalyticSurf.hpp b/packages/krino/krino/surface/Akri_AnalyticSurf.hpp index b35fd64fc0ea..decf25db4547 100644 --- a/packages/krino/krino/surface/Akri_AnalyticSurf.hpp +++ b/packages/krino/krino/surface/Akri_AnalyticSurf.hpp @@ -108,7 +108,9 @@ class Ellipsoid : public SurfaceThatDoesntTakeAdvantageOfNarrowBandAndThereforeH class Plane: public SurfaceThatDoesntTakeAdvantageOfNarrowBandAndThereforeHasCorrectSign { public: - Plane(const double normal[3], + Plane(const stk::math::Vector3d & normal, + const double offset) : Plane(normal, offset, 1.) {} + Plane(const stk::math::Vector3d & normal, const double offset, const double multiplier); diff --git a/packages/krino/krino/surface/Akri_Composite_Surface.cpp b/packages/krino/krino/surface/Akri_Composite_Surface.cpp index ed3a1df8f061..3f18e519a141 100644 --- a/packages/krino/krino/surface/Akri_Composite_Surface.cpp +++ b/packages/krino/krino/surface/Akri_Composite_Surface.cpp @@ -39,7 +39,8 @@ Composite_Surface::prepare_to_compute(const double time, const BoundingBox & poi for ( auto&& surface : my_subsurfaces ) { // keep nearby or moving surfaces - if (FACETED_SURFACE == surface->type() || + if (FACETED_SURFACE_2D == surface->type() || + FACETED_SURFACE_3D == surface->type() || nullptr != surface->get_transformation() || surface->does_intersect(point_bbox)) { diff --git a/packages/krino/krino/surface/Akri_Facet.cpp b/packages/krino/krino/surface/Akri_Facet.cpp index 92ab4ef36acf..67277c5beeb7 100644 --- a/packages/krino/krino/surface/Akri_Facet.cpp +++ b/packages/krino/krino/surface/Akri_Facet.cpp @@ -17,23 +17,18 @@ namespace krino{ // //-------------------------------------------------------------------------------- -std::unique_ptr -Facet::unpack_from_buffer( stk::CommBuffer & b ) +static void unpack_vector3d_from_buffer( stk::CommBuffer & b, stk::math::Vector3d & coords ) { - std::unique_ptr facet; - - int dim = 0; - b.unpack(dim); - - STK_ThrowRequire(dim == 2 || dim == 3); - - if (dim == 3) - facet = Facet3d::unpack_from_buffer(b); - else - facet = Facet2d::unpack_from_buffer(b); + b.unpack(coords[0]); + b.unpack(coords[1]); + b.unpack(coords[2]); +} - STK_ThrowRequire(facet); - return facet; +static void unpack_vector2d_from_buffer( stk::CommBuffer & b, stk::math::Vector3d & coords ) +{ + b.unpack(coords[0]); + b.unpack(coords[1]); + coords[2] = 0.; } //-------------------------------------------------------------------------------- @@ -41,68 +36,37 @@ Facet::unpack_from_buffer( stk::CommBuffer & b ) void Facet2d::pack_into_buffer(stk::CommBuffer & b) const { - const int dim = 2; - b.pack(dim); for (int n = 0; n < 2; ++n ) { - const stk::math::Vector3d & pt = facet_vertex(n); - b.pack(pt[0]); - b.pack(pt[1]); + b.pack(myCoords[n][0]); + b.pack(myCoords[n][1]); } } -std::unique_ptr -Facet2d::unpack_from_buffer( stk::CommBuffer & b ) +void +Facet2d::unpack_facet_data_from_buffer( stk::CommBuffer & b, std::array & facetCoords) { - std::unique_ptr facet; - - double vx, vy; - b.unpack(vx); - b.unpack(vy); - stk::math::Vector3d x0( vx, vy, 0.0 ); - b.unpack(vx); - b.unpack(vy); - stk::math::Vector3d x1( vx, vy, 0.0 ); - - facet = std::make_unique( x0, x1 ); - return facet; + unpack_vector2d_from_buffer(b, facetCoords[0]); + unpack_vector2d_from_buffer(b, facetCoords[1]); } void Facet3d::pack_into_buffer(stk::CommBuffer & b) const { - const int dim = 3; - b.pack(dim); for (int n = 0; n < 3; ++n ) { - const stk::math::Vector3d & pt = facet_vertex(n); - b.pack(pt[0]); - b.pack(pt[1]); - b.pack(pt[2]); + b.pack(myCoords[n][0]); + b.pack(myCoords[n][1]); + b.pack(myCoords[n][2]); } } -std::unique_ptr -Facet3d::unpack_from_buffer( stk::CommBuffer & b ) +void +Facet3d::unpack_facet_data_from_buffer( stk::CommBuffer & b, std::array & facetCoords) { - std::unique_ptr facet; - - double vx, vy, vz; - b.unpack(vx); - b.unpack(vy); - b.unpack(vz); - stk::math::Vector3d x0( vx, vy, vz ); - b.unpack(vx); - b.unpack(vy); - b.unpack(vz); - stk::math::Vector3d x1( vx, vy, vz ); - b.unpack(vx); - b.unpack(vy); - b.unpack(vz); - stk::math::Vector3d x2( vx, vy, vz ); - - facet = std::make_unique( x0, x1, x2 ); - return facet; + unpack_vector3d_from_buffer(b, facetCoords[0]); + unpack_vector3d_from_buffer(b, facetCoords[1]); + unpack_vector3d_from_buffer(b, facetCoords[2]); } std::ostream & Facet3d::put( std::ostream & os ) const @@ -159,25 +123,12 @@ Facet3d::Facet3d( const stk::math::Vector3d & pt0, { } -Facet3d::Facet3d( const double * pt0, - const double * pt1, - const double * pt2) - : myCoords{stk::math::Vector3d(pt0), stk::math::Vector3d(pt1), stk::math::Vector3d(pt2)} -{ -} - Facet2d::Facet2d( const stk::math::Vector3d & pt0, const stk::math::Vector3d & pt1 ) : myCoords{pt0, pt1} { } -Facet2d::Facet2d( const double * pt0, - const double * pt1 ) - : myCoords{stk::math::Vector3d(pt0,2), stk::math::Vector3d(pt1,2)} -{ -} - template bool are_all_components_lo(const VecType & bboxMin, const std::array & points, const unsigned comp) { diff --git a/packages/krino/krino/surface/Akri_Facet.hpp b/packages/krino/krino/surface/Akri_Facet.hpp index 489862c0b0a8..fcd70559870e 100644 --- a/packages/krino/krino/surface/Akri_Facet.hpp +++ b/packages/krino/krino/surface/Akri_Facet.hpp @@ -25,9 +25,6 @@ namespace krino { class Facet; class Transformation; -typedef std::vector< Facet * > FacetVec; -typedef std::vector< std::unique_ptr > FacetOwningVec; - class Facet { public: Facet() {} @@ -39,10 +36,6 @@ class Facet { virtual std::ostream & put( std::ostream& os ) const = 0; friend std::ostream & operator << ( std::ostream &os , const Facet &f ) { return f.put(os); } - // methods for off-processor communication - static std::unique_ptr unpack_from_buffer( stk::CommBuffer & b ); // static method that builds facet from data in buffer for off-processor communication - virtual void pack_into_buffer(stk::CommBuffer & b) const = 0; // pack into buffer for off-processor communication - static stk::math::Vector3d get_centroid(const Facet * facet) { return facet->centroid(); } static void insert_into_bounding_box(const Facet * facet, BoundingBox & bbox) { return facet->insert_into(bbox); } @@ -76,19 +69,23 @@ class Facet3d : public Facet { Facet3d() = delete; virtual ~Facet3d() {} + static constexpr int DIM=3; typedef CalcTriangle3 Calc; + static bool is_degenerate(const std::array & coords) { return Calc::normal_dir(coords).zero_length(); } + static stk::math::Vector3d get_centroid(const Facet3d * facet) { return facet->centroid(); } + static void insert_into_bounding_box(const Facet3d * facet, BoundingBox & bbox) { return facet->insert_into(bbox); } virtual size_t storage_size() const override { return sizeof(Facet3d); } virtual std::ostream & put( std::ostream& os ) const override; - static std::unique_ptr unpack_from_buffer( stk::CommBuffer & b ); // static method that builds facet from data in buffer for off-processor communication - virtual void pack_into_buffer(stk::CommBuffer & b) const override; + static void unpack_facet_data_from_buffer( stk::CommBuffer & b, std::array & facetCoords); + virtual void pack_into_buffer(stk::CommBuffer & b) const; virtual const stk::math::Vector3d & facet_vertex(const int i) const override { return myCoords[i]; } virtual stk::math::Vector3d & facet_vertex(const int i) override { return myCoords[i]; } virtual double facet_area() const override { return Calc::area(myCoords); } virtual stk::math::Vector3d facet_normal() const override { return Calc::normal(myCoords); } - virtual bool degenerate() const override { return Calc::normal_dir(myCoords).zero_length(); } + virtual bool degenerate() const override { return is_degenerate(myCoords); } virtual double point_distance_squared( const stk::math::Vector3d & x ) const override { return Calc::distance_squared( myCoords, x ); } virtual void closest_point( const stk::math::Vector3d & queryPt, stk::math::Vector3d & closestPt, stk::math::Vector2d & paramAtClosestPt ) const override @@ -116,13 +113,16 @@ class Facet2d : public Facet { Facet2d() = delete; virtual ~Facet2d() {} + static constexpr int DIM=2; typedef CalcSegment3 Calc; + static stk::math::Vector3d get_centroid(const Facet2d * facet) { return facet->centroid(); } + static void insert_into_bounding_box(const Facet2d * facet, BoundingBox & bbox) { return facet->insert_into(bbox); } virtual size_t storage_size() const override { return sizeof(Facet2d); } virtual std::ostream & put( std::ostream& os ) const override; - static std::unique_ptr unpack_from_buffer( stk::CommBuffer & b ); // static method that builds facet from data in buffer for off-processor communication - virtual void pack_into_buffer(stk::CommBuffer & b) const override; + static void unpack_facet_data_from_buffer( stk::CommBuffer & b, std::array & facetCoords); + virtual void pack_into_buffer(stk::CommBuffer & b) const; virtual const stk::math::Vector3d & facet_vertex(const int i) const override { return myCoords[i]; } virtual stk::math::Vector3d & facet_vertex(const int i) override { return myCoords[i]; } @@ -148,24 +148,25 @@ class Facet2d : public Facet { std::array myCoords; }; +template class FacetDistanceQuery { public: FacetDistanceQuery() : myFacet(nullptr), mySqrDistance(0.0) {} - FacetDistanceQuery(const Facet & facet, const stk::math::Vector3d & queryPt) : myFacet(&facet) + FacetDistanceQuery(const FACET & facet, const stk::math::Vector3d & queryPt) : myFacet(&facet) { myFacet->closest_point(queryPt, myClosestPt, myParamCoords); mySqrDistance = (queryPt-myClosestPt).length_squared(); } bool empty() const { return myFacet == nullptr; } - const Facet & facet() const { return *myFacet; } + const FACET & facet() const { return *myFacet; } double distance_squared() const { return mySqrDistance; } const stk::math::Vector3d & closest_point() const { return myClosestPt; } double signed_distance(const stk::math::Vector3d & queryPt) const { return std::sqrt(mySqrDistance)*myFacet->point_distance_sign(queryPt); } stk::math::Vector3d closest_point_weights() const { return myFacet->weights(myParamCoords); } private: - const Facet* myFacet; + const FACET* myFacet; stk::math::Vector3d myClosestPt; double mySqrDistance; stk::math::Vector2d myParamCoords; diff --git a/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.cpp b/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.cpp index 5540a7435842..d8a14e4e0c65 100644 --- a/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.cpp +++ b/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.cpp @@ -53,7 +53,8 @@ std::vector fill_processor_bounding_boxes(const BoundingBox & local return procPaddedQueryBboxes; } -static double compute_point_distance_squared(const stk::math::Vector3d &x, const FacetVec & nearestFacets) +template +static double compute_point_distance_squared(const stk::math::Vector3d &x, const std::vector & nearestFacets) { double minSqrDist = std::numeric_limits::max(); for ( auto&& facet : nearestFacets ) @@ -65,8 +66,8 @@ static double compute_point_distance_squared(const stk::math::Vector3d &x, const return minSqrDist; } -double -point_distance_given_nearest_facets(const stk::math::Vector3d &x, const FacetVec & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) +template +double point_distance_given_nearest_facets(const stk::math::Vector3d &x, const std::vector & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) { if (nearestFacets.empty()) { @@ -98,42 +99,93 @@ point_distance_given_nearest_facets(const stk::math::Vector3d &x, const FacetVec return dist; } -double -compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const FacetVec & facets) +template +std::vector> build_distance_queries(const stk::math::Vector3d &x, const std::vector & facets) { - - // If the closest_point weights are all larger than this value, then the closest point - // is considered to be on the face of the closest facet rather than on the edges of the facet, and - // therefore only the closest facet is considered in the distance calculation. Otherwise, all of the - // facets are considered to compute an average normal in order to compute the distance. - const double edge_tol = 1.e-6; - - std::vector facet_queries; - facet_queries.reserve(facets.size()); + std::vector> facetDistQueries; + facetDistQueries.reserve(facets.size()); for ( auto&& facet : facets ) { if (facet->degenerate()) continue; // Skip zero-sized facets - facet_queries.emplace_back(*facet, x); + facetDistQueries.emplace_back(*facet, x); } + return facetDistQueries; +} - STK_ThrowRequireMsg(!facet_queries.empty(), "All facets are degenerate in compute_point_to_facets_distance_by_average_normal."); +template +unsigned find_index_of_closest_facet(const std::vector> & facetDistQueries) +{ + unsigned closest = 0; + for ( unsigned index=0; index +stk::math::Vector3d compute_pseudo_normal(const std::vector> & facet_queries, const unsigned nearest) +{ + const double tol = 1.e-6; + stk::math::Vector3d pseudo_normal = stk::math::Vector3d::ZERO; + stk::math::Vector3d average_normal = stk::math::Vector3d::ZERO; + + const stk::math::Vector3d nearest_closest_point = facet_queries[nearest].closest_point(); + const double nearest_size2 = facet_queries[nearest].facet().mean_squared_edge_length(); + unsigned close_count = 0; + for ( auto&& query : facet_queries ) { - if ( facet_queries[index].distance_squared() < facet_queries[nearest].distance_squared() ) + const stk::math::Vector3d closest_point = query.closest_point(); + const double dist2_from_nearest = (closest_point-nearest_closest_point).length_squared(); + if (dist2_from_nearest < tol*tol*nearest_size2) { - nearest = index; + ++close_count; + const FACET & facet = query.facet(); + + average_normal += facet.facet_normal(); + + if constexpr (3 == FACET::DIM) + { + const stk::math::Vector3d closest_pt_wts = query.closest_point_weights(); + const int closest_node = (closest_pt_wts[0] > closest_pt_wts[1]) ? ((closest_pt_wts[0] > closest_pt_wts[2]) ? 0 : 2) : ((closest_pt_wts[1] > closest_pt_wts[2]) ? 1 : 2); + + const int n0 = closest_node; + const int n1 = (n0<2) ? (n0+1) : 0; + const int n2 = (n1<2) ? (n1+1) : 0; + const stk::math::Vector3d edge0 = facet.facet_vertex(n1) - facet.facet_vertex(n0); + const stk::math::Vector3d edge1 = facet.facet_vertex(n2) - facet.facet_vertex(n0); + const double facet_angle = std::acos(Dot(edge0, edge1)/(edge0.length()*edge1.length())); + + pseudo_normal += facet.facet_normal()*facet_angle; + } } } + STK_ThrowRequireMsg(close_count>0,"Issue with tolerance in compute_pseudo_normal. No facet found within tolerance of closest point."); + + return (3 == FACET::DIM && close_count > 2) ? pseudo_normal : average_normal; +} + +template +double compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const std::vector & facets) +{ + + // If the closest_point weights are all larger than this value, then the closest point + // is considered to be on the face of the closest facet rather than on the edges of the facet, and + // therefore only the closest facet is considered in the distance calculation. Otherwise, all of the + // facets are considered to compute an average normal in order to compute the distance. + const double edge_tol = 1.e-6; - if ( facet_queries[nearest].distance_squared() == 0. ) + const std::vector> facetDistQueries = build_distance_queries(x, facets); + STK_ThrowRequireMsg(!facetDistQueries.empty(), "All facets are degenerate in compute_point_to_facets_distance_by_average_normal."); + + const unsigned nearest = find_index_of_closest_facet(facetDistQueries); + + if ( facetDistQueries[nearest].distance_squared() == 0. ) { return 0.0; } - const stk::math::Vector3d closest_pt_wts = facet_queries[nearest].closest_point_weights(); - const int dim = dynamic_cast(facets[nearest]) ? 3 : 2; + const stk::math::Vector3d closest_pt_wts = facetDistQueries[nearest].closest_point_weights(); + const int dim = FACET::DIM; bool closest_point_on_edge = false; for (int d=0; d 0) + if (Dot(pseudo_normal, x-facetDistQueries[nearest].closest_point()) > 0) { return std::sqrt(min_sqr_dist); } @@ -172,48 +224,6 @@ compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, } } -stk::math::Vector3d -compute_pseudo_normal(const unsigned dim, const std::vector & facet_queries, const unsigned nearest) -{ - const double tol = 1.e-6; - stk::math::Vector3d pseudo_normal = stk::math::Vector3d::ZERO; - stk::math::Vector3d average_normal = stk::math::Vector3d::ZERO; - - const stk::math::Vector3d nearest_closest_point = facet_queries[nearest].closest_point(); - const double nearest_size2 = facet_queries[nearest].facet().mean_squared_edge_length(); - unsigned close_count = 0; - for ( auto&& query : facet_queries ) - { - const stk::math::Vector3d closest_point = query.closest_point(); - const double dist2_from_nearest = (closest_point-nearest_closest_point).length_squared(); - if (dist2_from_nearest < tol*tol*nearest_size2) - { - ++close_count; - const Facet & facet = query.facet(); - - average_normal += facet.facet_normal(); - - if (3 == dim) - { - const stk::math::Vector3d closest_pt_wts = query.closest_point_weights(); - const int closest_node = (closest_pt_wts[0] > closest_pt_wts[1]) ? ((closest_pt_wts[0] > closest_pt_wts[2]) ? 0 : 2) : ((closest_pt_wts[1] > closest_pt_wts[2]) ? 1 : 2); - - const int n0 = closest_node; - const int n1 = (n0<2) ? (n0+1) : 0; - const int n2 = (n1<2) ? (n1+1) : 0; - const stk::math::Vector3d edge0 = facet.facet_vertex(n1) - facet.facet_vertex(n0); - const stk::math::Vector3d edge1 = facet.facet_vertex(n2) - facet.facet_vertex(n0); - const double facet_angle = std::acos(Dot(edge0, edge1)/(edge0.length()*edge1.length())); - - pseudo_normal += facet.facet_normal()*facet_angle; - } - } - } - STK_ThrowRequireMsg(close_count>0,"Issue with tolerance in compute_pseudo_normal. No facet found within tolerance of closest point."); - - return (3 == dim && close_count > 2) ? pseudo_normal : average_normal; -} - bool is_projection_of_point_inside_enlarged_triangle(const stk::math::Vector3d & triPt0, const stk::math::Vector3d & triPt1, const stk::math::Vector3d & triPt2, const stk::math::Vector3d& p) { constexpr double expand {1e-10}; @@ -233,16 +243,18 @@ bool is_projection_of_point_inside_enlarged_segment(const stk::math::Vector3d & return Facet2d::Calc::is_projection_of_point_inside_segment(p0, p1, p); } -bool is_projection_of_point_inside_enlarged_facet(const Facet & facet, const stk::math::Vector3d& p) +bool is_projection_of_point_inside_enlarged_facet(const Facet3d & facet, const stk::math::Vector3d& p) { - const int dim = (dynamic_cast(&facet)) ? 3 : 2; + return is_projection_of_point_inside_enlarged_triangle(facet.facet_vertex(0), facet.facet_vertex(1), facet.facet_vertex(2), p); +} - if (3 == dim) - return is_projection_of_point_inside_enlarged_triangle(facet.facet_vertex(0), facet.facet_vertex(1), facet.facet_vertex(2), p); +bool is_projection_of_point_inside_enlarged_facet(const Facet2d & facet, const stk::math::Vector3d& p) +{ return is_projection_of_point_inside_enlarged_segment(facet.facet_vertex(0), facet.facet_vertex(1), p); } -std::pair compute_facet_edge_intersection(const Facet & facet, +template +std::pair compute_facet_edge_intersection(const FACET & facet, const stk::math::Vector3d& edgePt0, const stk::math::Vector3d& edgePt1) { @@ -259,26 +271,16 @@ std::pair compute_facet_edge_intersection(const Facet & facet, return {0, -1.}; } -std::pair compute_intersection_between_and_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1) +template +double compute_intersection_between_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1) { if (candidates.empty()) - return {0, -1.}; - - const double dist0 = compute_point_to_facets_distance_by_average_normal(edgePt0, candidates); - const double dist1 = compute_point_to_facets_distance_by_average_normal(edgePt1, candidates); - - if (!sign_change(dist0, dist1)) - return {0, -1.}; - - if (0. == dist0) - return {-1, 0.}; - if (0. == dist1) - return {1, 1.}; + return -1.; bool haveCrossing = false; double intersectionLoc = -1.; - for (const Facet * surfFacet : candidates) + for (const FACET * surfFacet : candidates) { const auto [facetCrossingSign, facetIntersectionLoc] = compute_facet_edge_intersection(*surfFacet, edgePt0, edgePt1); if (facetCrossingSign != 0) @@ -289,23 +291,31 @@ std::pair compute_intersection_between_and_surface_facets_and_edge( } } - if (haveCrossing) - return {sign(dist1), intersectionLoc}; - - // Sign change, but no crossing. This could be because there is a small gap in the facets, or there could be - // an unterminated surface far away. Decide based on magnitude of the distance at the linear crossing location. - - const double linearCrossingLoc = dist0 / (dist0 - dist1); + return intersectionLoc; +} - const stk::math::Vector3d linearCrossingPt = (1.-linearCrossingLoc) * edgePt0 + linearCrossingLoc * edgePt1; - const double minSqrDist = compute_point_distance_squared(linearCrossingPt, candidates); - const double edgeSqrLen = (edgePt1 - edgePt0).length_squared(); +template +stk::math::Vector3d compute_pseudo_normal(const stk::math::Vector3d &x, const std::vector & nearestFacets) +{ + const std::vector> facetDistQueries = build_distance_queries(x, nearestFacets); + STK_ThrowRequireMsg(!facetDistQueries.empty(), "All facets are degenerate in compute_pseudo_normal."); - constexpr double sqrTol = 1.e-4; // Large tol here is fine, right? - if (minSqrDist < sqrTol*edgeSqrLen) - return {sign(dist1), linearCrossingLoc}; + const unsigned nearest = find_index_of_closest_facet(facetDistQueries); - return {0, -1.}; + return compute_pseudo_normal(facetDistQueries, nearest); } +// Explicit template instantiation + +template double point_distance_given_nearest_facets(const stk::math::Vector3d &x, const std::vector & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance); +template double point_distance_given_nearest_facets(const stk::math::Vector3d &x, const std::vector & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance); +template double compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const std::vector & facets); +template double compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const std::vector & facets); +template stk::math::Vector3d compute_pseudo_normal(const stk::math::Vector3d &x, const std::vector & nearestFacets); +template stk::math::Vector3d compute_pseudo_normal(const stk::math::Vector3d &x, const std::vector & nearestFacets); +template std::pair compute_facet_edge_intersection(const Facet2d & facet, const stk::math::Vector3d& edgePt0, const stk::math::Vector3d& edgePt1); +template std::pair compute_facet_edge_intersection(const Facet3d & facet, const stk::math::Vector3d& edgePt0, const stk::math::Vector3d& edgePt1); +template double compute_intersection_between_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1); +template double compute_intersection_between_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1); + } diff --git a/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.hpp b/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.hpp index b6ce00b86591..365dd493c09e 100644 --- a/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.hpp +++ b/packages/krino/krino/surface/Akri_FacetedSurfaceCalcs.hpp @@ -8,15 +8,25 @@ #ifndef KRINO_KRINO_SURFACE_AKRI_FACETEDSURFACECALCS_HPP_ #define KRINO_KRINO_SURFACE_AKRI_FACETEDSURFACECALCS_HPP_ #include -#include +#include namespace krino { std::vector fill_processor_bounding_boxes(const BoundingBox & localFacetBbox, const BoundingBox & localQueryBbox, const double truncationLength); - double point_distance_given_nearest_facets(const stk::math::Vector3d &x, const FacetVec & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance); - double compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const FacetVec & facets); - stk::math::Vector3d compute_pseudo_normal(const unsigned dim, const std::vector & facet_queries, const unsigned nearest); - std::pair compute_facet_edge_intersection(const Facet & facet, const stk::math::Vector3d& edgePt0, const stk::math::Vector3d& edgePt1); - std::pair compute_intersection_between_and_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1); + + template + double point_distance_given_nearest_facets(const stk::math::Vector3d &x, const std::vector & nearestFacets, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance); + + template + double compute_point_to_facets_distance_by_average_normal(const stk::math::Vector3d &x, const std::vector & facets); + + template + stk::math::Vector3d compute_pseudo_normal(const stk::math::Vector3d &x, const std::vector & nearestFacets); + + template + std::pair compute_facet_edge_intersection(const FACET & facet, const stk::math::Vector3d& edgePt0, const stk::math::Vector3d& edgePt1); + + template + double compute_intersection_between_surface_facets_and_edge(const std::vector & candidates, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1); } diff --git a/packages/krino/krino/surface/Akri_Faceted_Surface.cpp b/packages/krino/krino/surface/Akri_Faceted_Surface.cpp index 9b28f9ea6de4..7878e1e666dc 100644 --- a/packages/krino/krino/surface/Akri_Faceted_Surface.cpp +++ b/packages/krino/krino/surface/Akri_Faceted_Surface.cpp @@ -12,30 +12,33 @@ #include #include #include +#include +#include #include #include #include #include #include +#include namespace krino{ -Faceted_Surface::Faceted_Surface(const std::string & sn) - : SurfaceThatTakesAdvantageOfNarrowBandAndThereforeMightHaveWrongSign(), - my_name(sn) {} - -static int find_destination_proc_for_facet(const std::vector & procBboxes, const Facet * facet) +template +static int find_destination_proc_for_facet(const std::vector & procBboxes, const FACET & facet) { int numProcs = procBboxes.size(); for ( int destProc = 1; destProc < numProcs; ++destProc ) - if ( facet->does_intersect(procBboxes[destProc]) ) + if ( facet.does_intersect(procBboxes[destProc]) ) return destProc; return 0; } -static void unpack_and_append_facets_from_proc(stk::CommSparse & commSparse, const int recvProc, FacetOwningVec & facetVec) +template +static void unpack_and_append_facets_from_proc(stk::CommSparse & commSparse, const int recvProc, std::vector & facetVec) { + std::array facetCoords; + stk::CommBuffer & b = commSparse.recv_buffer(recvProc); if (b.remaining()) { @@ -47,21 +50,37 @@ static void unpack_and_append_facets_from_proc(stk::CommSparse & commSparse, con for ( size_t n = 0; n < numProcFacets; ++n ) { - std::unique_ptr facet = Facet::unpack_from_buffer( b ); - facetVec.emplace_back( std::move(facet) ); + FACET::unpack_facet_data_from_buffer(b, facetCoords); + if constexpr (FACET::DIM == 3) + facetVec.emplace_back(facetCoords[0], facetCoords[1], facetCoords[2]); + else + facetVec.emplace_back(facetCoords[0], facetCoords[1]); } STK_ThrowAssert( 0 == b.remaining() ); } } -void -Faceted_Surface::parallel_distribute_facets(const size_t batch_size, const std::vector & procBboxes) +const std::vector & FacetedSurfaceBase::get_facets_2d() const { return as_derived_type().get_facets(); } +const std::vector & FacetedSurfaceBase::get_facets_3d() const { return as_derived_type().get_facets(); } +std::vector & FacetedSurfaceBase::get_facets_2d() { return as_derived_type().get_facets(); } +std::vector & FacetedSurfaceBase::get_facets_3d() { return as_derived_type().get_facets(); } + +template +void Faceted_Surface::swap(FacetedSurfaceBase & rhs) +{ + Faceted_Surface * other = dynamic_cast*>(&rhs); + STK_ThrowRequire(nullptr != other); + myLocalFacets.swap(other->myLocalFacets); +} + +template +void Faceted_Surface::parallel_distribute_facets(const size_t batch_size, const std::vector & procBboxes) { const int num_procs = stk::EnvData::parallel_size(); if ( num_procs == 1 ) return; const int me = stk::EnvData::parallel_rank(); - STK_ThrowRequire(me != 0 || batch_size <= my_local_facets.size()); + STK_ThrowRequire(me != 0 || batch_size <= myLocalFacets.size()); stk::CommSparse comm_sparse(stk::EnvData::parallel_comm()); @@ -72,10 +91,10 @@ Faceted_Surface::parallel_distribute_facets(const size_t batch_size, const std:: { dest_procs.resize(batch_size, 0); proc_facet_counts.resize(num_procs, 0); - start = my_local_facets.size() - batch_size; + start = myLocalFacets.size() - batch_size; for ( size_t index = 0; index < batch_size; ++index ) { - dest_procs[index] = find_destination_proc_for_facet(procBboxes, my_local_facets[start+index].get()); + dest_procs[index] = find_destination_proc_for_facet(procBboxes, myLocalFacets[start+index]); ++(proc_facet_counts[dest_procs[index]]); } } @@ -100,8 +119,8 @@ Faceted_Surface::parallel_distribute_facets(const size_t batch_size, const std:: { if (dest_procs[index] == dest_proc) { - const Facet * facet = my_local_facets[start+index].get(); - facet->pack_into_buffer(b); + const FACET & facet = myLocalFacets[start+index]; + facet.pack_into_buffer(b); } } } @@ -117,75 +136,36 @@ Faceted_Surface::parallel_distribute_facets(const size_t batch_size, const std:: if (me == 0) { // delete facets that have been sent away (even if they are headed back to 0) - my_local_facets.erase(my_local_facets.begin()+start, my_local_facets.end()); - } - - unpack_and_append_facets_from_proc(comm_sparse, 0, my_local_facets); -} - -void -Faceted_Surface::prepare_to_compute(const double time, const BoundingBox & point_bbox, const double truncation_length) -{ - - build_local_facets(point_bbox); - - if (my_transformation != nullptr) - { - for ( auto&& facet : my_local_facets ) - { - facet->apply_transformation(*my_transformation); - } - } - - my_bounding_box.clear(); - for (auto && facet : my_local_facets) - { - facet->insert_into(my_bounding_box); + myLocalFacets.erase(myLocalFacets.begin()+start, myLocalFacets.end()); } - my_all_facets.clear(); - for ( auto&& facet : my_local_facets ) - { - my_all_facets.push_back(facet.get()); - } - - // Get all remote facets that might be closest to this processors query points - gather_nonlocal_facets(point_bbox, truncation_length); - - if(krinolog.shouldPrint(LOG_DEBUG)) - krinolog << "P" << stk::EnvData::parallel_rank() << ":" << " Building facet tree for " << my_all_facets.size() << " facets." << stk::diag::dendl; - - my_facet_tree = std::make_unique>( my_all_facets, Facet::get_centroid, Facet::insert_into_bounding_box ); - - if ( krinolog.shouldPrint(LOG_DEBUG) ) - krinolog << "P" << stk::EnvData::parallel_rank() << ": After building search tree, storage size is " << storage_size()/(1024.*1024.) << " Mb." << stk::diag::dendl; -} - -void Faceted_Surface::prepare_to_compute(const BoundingBox & point_bbox, const double truncation_length) -{ - STK_ThrowAssert(nullptr == my_transformation); - prepare_to_compute(0.0, point_bbox, truncation_length); + unpack_and_append_facets_from_proc(comm_sparse, 0, myLocalFacets); } -static void append_intersecting_facets(const BoundingBox & sourceBbox, const std::vector & sourceFacets, const BoundingBox & targetBbox, std::vector & targetFacets) +template +static void append_intersecting_facets(const BoundingBox & sourceBbox, const std::vector & sourceFacets, const BoundingBox & targetBbox, std::vector & targetFacets) { - targetFacets.reserve(targetFacets.size() + sourceFacets.size()); if (targetBbox.intersects(sourceBbox)) { for ( auto&& facet : sourceFacets ) - if ( facet->does_intersect(targetBbox) ) - targetFacets.push_back(facet); + if ( facet.does_intersect(targetBbox) ) + targetFacets.push_back(&facet); } } -static void retain_intersecting_facets(const BoundingBox & bbox, std::vector & searchFacets) +template +static void store_pointers_to_local_intersecting_facets_and_nonlocal_facets(const BoundingBox & bbox, const std::vector & localFacets, const std::vector & nonLocalFacets, std::vector & allFacetPtrs) { - std::vector targetFacets; - append_intersecting_facets(bbox, searchFacets, bbox, targetFacets); - searchFacets.swap(targetFacets); + allFacetPtrs.clear(); + allFacetPtrs.reserve(localFacets.size()+nonLocalFacets.size()); + append_intersecting_facets(bbox, localFacets, bbox, allFacetPtrs); + + for ( auto&& facet : nonLocalFacets ) + allFacetPtrs.push_back(&facet); } -static size_t get_global_num_facets(const stk::ParallelMachine comm, const std::vector & localFacets) +template +static size_t get_global_num_facets(const stk::ParallelMachine comm, const std::vector & localFacets) { const size_t localNumFacets = localFacets.size(); size_t globalNumFacets = 0; @@ -193,26 +173,23 @@ static size_t get_global_num_facets(const stk::ParallelMachine comm, const std:: return globalNumFacets; } -void -Faceted_Surface::gather_nonlocal_facets(const BoundingBox & point_bbox, const double truncation_length) +template +static void fill_intersecting_nonlocal_facets(const BoundingBox & localFacetBBox, const std::vector & procPaddedQueryBboxes, const std::vector & localFacets, std::vector & nonLocalFacetsThatIntersectProcPaddedQueryBbox) { + nonLocalFacetsThatIntersectProcPaddedQueryBbox.clear(); - // If truncation length is specified, get all of the facets that are within - // our padded processor's bounding box. - // To do this, see if any local facets lie in the padded nodal bounding box - // of another proc. if so, send them a copy of those facets. - - my_nonlocal_facets.clear(); - const stk::ParallelMachine comm = stk::EnvData::parallel_comm(); const int numProcs = stk::parallel_machine_size(comm); if ( numProcs == 1) return; - const std::vector procPaddedQueryBboxes = fill_processor_bounding_boxes(my_bounding_box, point_bbox, truncation_length); + // If truncation length is specified, get all of the facets that are within + // our padded processor's bounding box. + // To do this, see if any local facets lie in the padded nodal bounding box + // of another proc. if so, send them a copy of those facets. const int me = stk::parallel_machine_rank(comm); - const size_t numFacetsPerBatch = 1 + get_global_num_facets(comm, my_all_facets)/numProcs; // min size of 1 - std::vector intersectingFacets; + const size_t numFacetsPerBatch = 1 + get_global_num_facets(comm, localFacets)/numProcs; // min size of 1 + std::vector intersectingFacets; std::vector> procAndNumFacets; // Perform communication in batches, communicating with as many neighbors as possible until the batch size is reached. @@ -230,7 +207,7 @@ Faceted_Surface::gather_nonlocal_facets(const BoundingBox & point_bbox, const do const int destProc = (me+commPartner) % numProcs; const size_t numPreviousIntersectingFacets = intersectingFacets.size(); - append_intersecting_facets(my_bounding_box, my_all_facets, procPaddedQueryBboxes[destProc], intersectingFacets); + append_intersecting_facets(localFacetBBox, localFacets, procPaddedQueryBboxes[destProc], intersectingFacets); const size_t procNumIntersectingFacets = intersectingFacets.size() - numPreviousIntersectingFacets; procAndNumFacets.emplace_back(destProc, procNumIntersectingFacets); @@ -263,7 +240,7 @@ Faceted_Surface::gather_nonlocal_facets(const BoundingBox & point_bbox, const do for(int recvProc=0; recvProc +void Faceted_Surface::prepare_to_compute(const double time, const BoundingBox & point_bbox, const double truncation_length) +{ + build_local_facets(point_bbox); + + if (my_transformation != nullptr) + { + for ( auto&& facet : myLocalFacets ) + { + facet.apply_transformation(*my_transformation); + } + } + + my_bounding_box.clear(); + for (auto && facet : myLocalFacets) + facet.insert_into(my_bounding_box); + + const std::vector procPaddedQueryBboxes = fill_processor_bounding_boxes(my_bounding_box, point_bbox, truncation_length); + + // Get all remote facets that might be closest to this processors query points + fill_intersecting_nonlocal_facets(my_bounding_box, procPaddedQueryBboxes, myLocalFacets, myNonLocalFacets); + + const int me = stk::EnvData::parallel_rank(); + store_pointers_to_local_intersecting_facets_and_nonlocal_facets(procPaddedQueryBboxes[me], myLocalFacets, myNonLocalFacets, myAllFacetPtrs); + + if(krinolog.shouldPrint(LOG_DEBUG)) + krinolog << "P" << stk::EnvData::parallel_rank() << ":" << " Building facet tree for " << myAllFacetPtrs.size() << " facets." << stk::diag::dendl; - retain_intersecting_facets(procPaddedQueryBboxes[me], my_all_facets); + my_facet_tree = std::make_unique>( myAllFacetPtrs, FACET::get_centroid, FACET::insert_into_bounding_box ); - // copy pointers to nonlocal facets into my_all_facets - my_all_facets.reserve(my_all_facets.size()+ my_nonlocal_facets.size()); - for (auto && nonlocal_descendant : my_nonlocal_facets) { my_all_facets.push_back(nonlocal_descendant.get()); } + if ( krinolog.shouldPrint(LOG_DEBUG) ) + krinolog << "P" << stk::EnvData::parallel_rank() << ": After building search tree, storage size is " << storage_size()/(1024.*1024.) << " Mb." << stk::diag::dendl; } -double -Faceted_Surface::point_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) const +template +double Faceted_Surface::point_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) const { STK_ThrowAssertMsg(my_facet_tree, "ERROR: Empty facet tree"); // get all facets we need to check - FacetVec nearestFacets; + std::vector nearestFacets; my_facet_tree->find_closest_entities( x, nearestFacets, narrow_band_size ); STK_ThrowRequire( !nearestFacets.empty() || 0.0 != narrow_band_size || my_facet_tree->empty() ); return point_distance_given_nearest_facets(x, nearestFacets, narrow_band_size, far_field_value, compute_signed_distance); } -size_t -Faceted_Surface::storage_size() const +template +stk::math::Vector3d Faceted_Surface::pseudo_normal_at_closest_point(const stk::math::Vector3d &x) const { - size_t store_size = sizeof(Faceted_Surface); + if (my_facet_tree->empty()) + return stk::math::Vector3d::ZERO; + + std::vector nearestFacets; + my_facet_tree->find_closest_entities( x, nearestFacets, 0. ); + STK_ThrowRequire( !nearestFacets.empty() ); + + return compute_pseudo_normal(x, nearestFacets); +} + +template +size_t Faceted_Surface::storage_size() const +{ + size_t store_size = sizeof(Faceted_Surface); if (my_facet_tree) { store_size += my_facet_tree->storage_size(); } - store_size += krino::storage_size(my_local_facets); - store_size += krino::storage_size(my_nonlocal_facets); - store_size += krino::storage_size(my_all_facets); + store_size += krino::storage_size(myLocalFacets); + store_size += krino::storage_size(myNonLocalFacets); + store_size += krino::storage_size(myAllFacetPtrs); return store_size; } -static std::vector find_candidate_surface_facets_for_intersection_with_segment(SearchTree & facetTree, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1) +template +static std::vector find_candidate_surface_facets_for_intersection_with_segment(SearchTree & facetTree, const stk::math::Vector3d & edgePt0, const stk::math::Vector3d & edgePt1) { BoundingBox facetBbox; facetBbox.accommodate(edgePt0); facetBbox.accommodate(edgePt1); facetBbox.pad_epsilon(); - std::vector results; + std::vector results; facetTree.get_intersecting_entities(facetBbox, results); return results; } -std::pair Faceted_Surface::compute_intersection_with_segment(const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double edgeCrossingTol) const +template +std::pair Faceted_Surface::compute_intersection_with_segment(const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double edgeCrossingTol) const +{ + const double dist0 = point_signed_distance(pt0); + const double dist1 = point_signed_distance(pt1); + if (sign_change(dist0, dist1)) + { + constexpr double tightEnoughTolToCalcDistanceAtSignedDistCrossing = 0.001; + constexpr double looseEnoughTolToHandleSmallGaps = 0.05; + const std::vector candidates = find_candidate_surface_facets_for_intersection_with_segment(*my_facet_tree, pt0, pt1); + const double intersectionLoc = compute_intersection_between_surface_facets_and_edge(candidates, pt0, pt1); + if (intersectionLoc >= 0.) + return {sign(dist1), intersectionLoc}; + auto [crossingSign, linearCrossingLoc] = compute_surface_intersection_with_crossed_segment_from_signed_distance(*this, pt0, pt1, dist0, dist1, tightEnoughTolToCalcDistanceAtSignedDistCrossing); + const stk::math::Vector3d linearCrossingCoords = (1.-linearCrossingLoc)*pt0 + linearCrossingLoc*pt1; + const double crossingDist = point_signed_distance(linearCrossingCoords); + if (std::abs(crossingDist) < looseEnoughTolToHandleSmallGaps*(pt1-pt0).length()) + return {sign(dist1), linearCrossingLoc}; + +// krinolog << "Kind of surprised " << crossingDist << " " << looseEnoughTolToHandleSmallGaps*(pt1-pt0).length() << " " << (pt1-pt0).length() << " " << dist0 << " " << dist1 << " " << linearCrossingLoc << stk::diag::dendl; +// krinolog << pt0 << " " << pt1 << stk::diag::dendl; + } + return {0, -1.}; +} + +template +std::string Faceted_Surface::print_sizes() const { - const std::vector candidates = find_candidate_surface_facets_for_intersection_with_segment(*my_facet_tree, pt0, pt1); - return compute_intersection_between_and_surface_facets_and_edge(candidates, pt0, pt1); + const stk::ParallelMachine comm = stk::EnvData::parallel_comm(); + unsigned numFacets = size(); + all_reduce_sum(comm, numFacets); + + // iterate local facets and calc area + + double facets_totalArea = 0.0; // total surface area of interface from facets + double facets_maxArea = -1.0; // area of largest facet on interface + double facets_minArea = std::numeric_limits::max(); // area of smallest facet on interface + + // loop over facets + for ( auto&& facet : myLocalFacets ) + { + double area = facet.facet_area(); + facets_totalArea += area; + + facets_minArea = std::min(area, facets_minArea); + facets_maxArea = std::max(area, facets_maxArea); + } + + all_reduce_min(comm, facets_minArea); + all_reduce_max(comm, facets_maxArea); + all_reduce_sum(comm, facets_totalArea); + + std::ostringstream out; + out << "\t " << "Global sizes: { " << numFacets << " facets on }" << std::endl; + out << "\t " << "Local sizes: { " << size() << " facets on }" << std::endl; + out << "\t Global areas: { Min = " << facets_minArea << ", Max = " << facets_maxArea + << ", Total = " << facets_totalArea << " }" << std::endl; + return out.str(); } +// Explicit template instantiation +template class Faceted_Surface; +template class Faceted_Surface; + } // namespace krino diff --git a/packages/krino/krino/surface/Akri_Faceted_Surface.hpp b/packages/krino/krino/surface/Akri_Faceted_Surface.hpp index d2aa7b9e57d0..b7674ad14d2a 100644 --- a/packages/krino/krino/surface/Akri_Faceted_Surface.hpp +++ b/packages/krino/krino/surface/Akri_Faceted_Surface.hpp @@ -12,15 +12,78 @@ #include #include #include +#include +#include namespace krino { -class Faceted_Surface : public SurfaceThatTakesAdvantageOfNarrowBandAndThereforeMightHaveWrongSign { +template +class Faceted_Surface; +class FacetedSurfaceBase : public SurfaceThatTakesAdvantageOfNarrowBandAndThereforeMightHaveWrongSign +{ public: - Faceted_Surface(const std::string & sn); + FacetedSurfaceBase(const int dim) + : SurfaceThatTakesAdvantageOfNarrowBandAndThereforeMightHaveWrongSign(), myDim(dim) + { STK_ThrowRequire(myDim==2 || myDim==3); } + + static std::unique_ptr build(const int dim); + + virtual Surface_Type type() const override { return (myDim == 3) ? FACETED_SURFACE_3D : FACETED_SURFACE_2D; } + + virtual void swap(FacetedSurfaceBase & other) = 0; + virtual void clear() = 0; + virtual size_t size() const = 0; + virtual size_t nonlocal_size() const = 0; + virtual double point_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) const = 0; + virtual void prepare_to_compute(const double time, const BoundingBox & point_bbox, const double truncation_length) = 0; + virtual std::string print_sizes() const = 0; + + void prepare_to_compute(const BoundingBox & point_bbox, const double truncation_length) + { + STK_ThrowAssert(nullptr == my_transformation); + prepare_to_compute(0.0, point_bbox, truncation_length); + } + + double point_unsigned_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value) const + { + return point_distance(x, narrow_band_size, far_field_value, false); + } + + template + const Faceted_Surface & as_derived_type() const + { + const auto * derived = dynamic_cast *>(this); + STK_ThrowRequireMsg(derived, "Can't access FacetedSurfaceBase as Faceted_Surface<" << typeid(this).name() << ">"); + return *derived; + } + + template + Faceted_Surface & as_derived_type() + { + auto * derived = dynamic_cast *>(this); + STK_ThrowRequireMsg(derived, "Can't access FacetedSurfaceBase as Faceted_Surface<" << typeid(this).name() << ">"); + return *derived; + } + + const std::vector & get_facets_2d() const; + const std::vector & get_facets_3d() const; + std::vector & get_facets_2d(); + std::vector & get_facets_3d(); + + void emplace_back_2d(const stk::math::Vector3d & x0, const stk::math::Vector3d & x1) { get_facets_2d().emplace_back(x0,x1); } + void emplace_back_3d(const stk::math::Vector3d & x0, const stk::math::Vector3d & x1, const stk::math::Vector3d & x2) { get_facets_3d().emplace_back(x0,x1,x2); } + +private: + int myDim; +}; + +template +class Faceted_Surface : public FacetedSurfaceBase +{ +public: + Faceted_Surface() : FacetedSurfaceBase(FACET::DIM) {} - virtual Surface_Type type() const override { return FACETED_SURFACE; } virtual size_t storage_size() const override; virtual void prepare_to_compute(const double time, const BoundingBox & point_bbox, const double truncation_length) override; virtual double truncated_point_signed_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value) const override @@ -30,37 +93,40 @@ class Faceted_Surface : public SurfaceThatTakesAdvantageOfNarrowBandAndTherefore virtual std::pair compute_intersection_with_segment(const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double edgeCrossingTol) const override; // query/modify facets - void add( std::unique_ptr facet ) { my_local_facets.emplace_back(std::move(facet)); } - void reserve(unsigned size_) { my_local_facets.reserve(size_); } - unsigned size() const { return my_local_facets.size(); } - unsigned nonlocal_size() const { return my_nonlocal_facets.size(); } - Facet * operator()( const unsigned index ) const { return my_local_facets[index].get(); } - void clear() { my_local_facets.clear(); } - void swap(Faceted_Surface & other) { my_local_facets.swap(other.my_local_facets); } - const FacetOwningVec & get_facets() const { return my_local_facets; } + virtual void clear() override { myLocalFacets.clear(); } + virtual size_t size() const override { return myLocalFacets.size(); } + virtual void swap(FacetedSurfaceBase & other) override; + + virtual size_t nonlocal_size() const override { return myNonLocalFacets.size(); } + virtual std::string print_sizes() const; + void parallel_distribute_facets(const size_t batch_size, const std::vector & proc_bboxes); + double point_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) const override; + const std::vector & get_facets() const { return myLocalFacets; } + std::vector & get_facets() { return myLocalFacets; } + public: - void prepare_to_compute(const BoundingBox & point_bbox, const double truncation_length); - double point_unsigned_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value) const - { - return point_distance(x, narrow_band_size, far_field_value, false); - } + stk::math::Vector3d pseudo_normal_at_closest_point(const stk::math::Vector3d &x) const; private: virtual void build_local_facets(const BoundingBox & proc_bbox) {} - double point_distance(const stk::math::Vector3d &x, const double narrow_band_size, const double far_field_value, const bool compute_signed_distance) const; - void gather_nonlocal_facets(const BoundingBox & local_bbox, const double truncation_length); - std::string my_name; - FacetOwningVec my_local_facets; + std::vector myLocalFacets; - mutable std::unique_ptr> my_facet_tree; - mutable FacetOwningVec my_nonlocal_facets; - mutable FacetVec my_all_facets; + mutable std::unique_ptr> my_facet_tree; + mutable std::vector myNonLocalFacets; + mutable std::vector myAllFacetPtrs; BoundingBox my_bounding_box; }; +inline std::unique_ptr FacetedSurfaceBase::build(const int dim) +{ + if (dim == 2) + return std::make_unique>(); + return std::make_unique>(); +} + } // namespace krino #endif // Akri_Faceted_Surface_h diff --git a/packages/krino/krino/surface/Akri_Surface.hpp b/packages/krino/krino/surface/Akri_Surface.hpp index f70ff8d331b0..9ffda9d5c662 100644 --- a/packages/krino/krino/surface/Akri_Surface.hpp +++ b/packages/krino/krino/surface/Akri_Surface.hpp @@ -18,17 +18,7 @@ namespace krino { template inline size_t -storage_size(const std::vector & vec) {return (sizeof(std::vector) + vec.size()*sizeof(T*)); } - -// linear cost, requires T::storage_size -template -inline size_t -storage_size(const std::vector> & autovec) -{ - size_t store_size = sizeof(std::vector>) + autovec.size()*sizeof(std::unique_ptr); - for (auto && entry : autovec) store_size += entry->storage_size(); - return store_size; -} +storage_size(const std::vector & vec) {return (sizeof(std::vector) + vec.capacity()*sizeof(T)); } class Transformation; @@ -42,7 +32,8 @@ enum Surface_Type PLANE, RANDOM, STRING_FUNCTION, - FACETED_SURFACE, + FACETED_SURFACE_2D, + FACETED_SURFACE_3D, // Never, ever, ever add an entry after MAX_SURFACE_TYPE. Never. MAX_SURFACE_TYPE }; diff --git a/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.cpp b/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.cpp index d5431f04fab6..2fc2b07b525c 100644 --- a/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.cpp +++ b/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.cpp @@ -13,25 +13,30 @@ namespace krino { // Note: Seems pretty strange, but the boost::toms748_solve in find_root appears to do better on the interval -1->1 than from 0->1 -static std::function build_edge_distance_function(const Surface & surface, const std::array & edgeNodeCoords) +static std::function build_edge_distance_function(const Surface & surface, const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1) { std::function distanceFunction = - [&surface, &edgeNodeCoords](const double x) + [&surface, &pt0, &pt1](const double x) { - return surface.point_signed_distance(0.5*(1.-x)*edgeNodeCoords[0] + 0.5*(1.+x)*edgeNodeCoords[1]); + return surface.point_signed_distance(0.5*(1.-x)*pt0 + 0.5*(1.+x)*pt1); }; return distanceFunction; } -static double find_crossing_position(const Surface & surface, const std::array & edgeNodeCoords, const double edgeTol) +static double find_edge_crossing_position(const std::function & edgeDistanceFunc, const double dist0, const double dist1, const double edgeTol) { - const double phi0 = surface.point_signed_distance(edgeNodeCoords[0]); - const double phi1 = surface.point_signed_distance(edgeNodeCoords[1]); const int maxIters = 100; - const auto result = find_root(build_edge_distance_function(surface, edgeNodeCoords), -1., 1., phi0, phi1, maxIters, 2.*edgeTol); + const auto result = find_root(edgeDistanceFunc, -1., 1., dist0, dist1, maxIters, edgeTol); STK_ThrowRequire(result.first); - return 0.5*(1.+result.second); + return result.second; +} + +std::pair compute_surface_intersection_with_crossed_segment_from_signed_distance(const Surface & surface, const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double dist0, const double dist1, const double edgeCrossingTol) +{ + STK_ThrowAssert(sign_change(dist0, dist1)); + const double location = find_edge_crossing_position(build_edge_distance_function(surface, pt0, pt1), dist0, dist1, 2*edgeCrossingTol); // 2x edgeTol because of 2x range (-1->1) + return {sign(dist1), 0.5*(1.+location)}; } std::pair compute_surface_intersection_with_segment_from_signed_distance(const Surface & surface, const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double edgeCrossingTol) @@ -39,11 +44,7 @@ std::pair compute_surface_intersection_with_segment_from_signed_dis const double phi0 = surface.point_signed_distance(pt0); const double phi1 = surface.point_signed_distance(pt1); if (sign_change(phi0, phi1)) - { - const std::array edgeNodeCoords{pt0, pt1}; - const double location = find_crossing_position(surface, edgeNodeCoords, edgeCrossingTol); - return {sign(phi1), location}; - } + return compute_surface_intersection_with_crossed_segment_from_signed_distance(surface, pt0, pt1, phi0, phi1, edgeCrossingTol); return {0, -1.}; } diff --git a/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.hpp b/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.hpp index 442ea29346b8..563d19c86c62 100644 --- a/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.hpp +++ b/packages/krino/krino/surface/Akri_SurfaceIntersectionFromSignedDistance.hpp @@ -6,6 +6,7 @@ namespace krino { class Surface; + std::pair compute_surface_intersection_with_crossed_segment_from_signed_distance(const Surface & surface, const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double dist0, const double dist1, const double edgeCrossingTol); std::pair compute_surface_intersection_with_segment_from_signed_distance(const Surface & surface, const stk::math::Vector3d &pt0, const stk::math::Vector3d &pt1, const double edgeCrossingTol); } diff --git a/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.cpp b/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.cpp index 69cd02f60073..7f2058d268b6 100644 --- a/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.cpp +++ b/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.cpp @@ -21,6 +21,12 @@ StkMeshBuilder::StkMeshBuilder(stk::mesh::BulkData & mesh, const stk::Para mMesh.mesh_meta_data().use_simple_fields(); } +template +const stk::mesh::FieldBase & StkMeshBuilder::get_coordinates_field() const +{ + return *mMesh.mesh_meta_data().coordinate_field(); +} + template void StkMeshBuilder::declare_coordinates() { diff --git a/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.hpp b/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.hpp index f6fbd3cf60d6..19fd433e40e9 100644 --- a/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.hpp +++ b/packages/krino/krino/unit_tests/Akri_StkMeshBuilder.hpp @@ -66,6 +66,7 @@ class StkMeshBuilder stk::mesh::Entity get_side_with_nodes(const std::vector &nodesOfSide) const; void create_block_parts(const std::vector &elementBlockIDs); stk::math::Vector3d get_node_coordinates(const stk::mesh::Entity node) const; + const stk::mesh::FieldBase & get_coordinates_field() const; const stk::mesh::PartVector & get_block_parts() const { return mBlockParts; } void write_mesh(const std::string & fileName); diff --git a/packages/krino/krino/unit_tests/Akri_StkMeshFixture.hpp b/packages/krino/krino/unit_tests/Akri_StkMeshFixture.hpp index e4112f1bd5fc..c8b21a1eae33 100644 --- a/packages/krino/krino/unit_tests/Akri_StkMeshFixture.hpp +++ b/packages/krino/krino/unit_tests/Akri_StkMeshFixture.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -33,6 +34,7 @@ class StkMeshAndBuilder std::vector get_ids_of_nodes_with_given_indices(const std::vector & nodeIndices) const { return mBuilder.get_ids_of_nodes_with_given_indices(nodeIndices); } const std::vector & get_owned_elements() const { return mBuilder.get_owned_elements(); } stk::math::Vector3d get_node_coordinates(const stk::mesh::Entity node) const { return mBuilder.get_node_coordinates(node); } + const stk::mesh::FieldBase & get_coordinates_field() const { return mBuilder.get_coordinates_field(); } template void build_full_np1_mesh(const MeshSpecType &meshSpec) diff --git a/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.cpp b/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.cpp new file mode 100644 index 000000000000..86b8942ab65a --- /dev/null +++ b/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.cpp @@ -0,0 +1,37 @@ +#include + +#include +#include +#include +#include +#include +#include + +namespace krino { + +stk::mesh::Entity find_local_node_closest_to_location(const stk::mesh::BulkData& mesh, const stk::mesh::Selector& nodeSelector, const FieldRef coordsField, const stk::math::Vector3d& location) +{ + const unsigned dim = mesh.mesh_meta_data().spatial_dimension(); + double minSqrDistance = std::numeric_limits::max(); + stk::mesh::Entity closest; + + for(const auto & bucketPtr : mesh.get_buckets(stk::topology::NODE_RANK, nodeSelector)) + { + for(const auto node : *bucketPtr) + { + const stk::math::Vector3d nodeLocation = get_vector_field(mesh, coordsField, node, dim); + const double nodeSqrDist = (nodeLocation-location).length_squared(); + if(nodeSqrDist < minSqrDistance) + { + minSqrDistance = nodeSqrDist; + closest = node; + } + } + } + + return closest; +} + +} + + diff --git a/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.hpp b/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.hpp new file mode 100644 index 000000000000..b69ac3dc5a25 --- /dev/null +++ b/packages/krino/krino/unit_tests/Akri_UnitMeshUtils.hpp @@ -0,0 +1,15 @@ +#ifndef KRINO_KRINO_UNIT_TESTS_AKRI_UNITMESHUTILS_HPP_ +#define KRINO_KRINO_UNIT_TESTS_AKRI_UNITMESHUTILS_HPP_ + +#include +#include + +namespace krino { + +class FieldRef; + +stk::mesh::Entity find_local_node_closest_to_location(const stk::mesh::BulkData& mesh, const stk::mesh::Selector& nodeSelector, const FieldRef coordsField, const stk::math::Vector3d& location); + +} + +#endif /* KRINO_KRINO_UNIT_TESTS_AKRI_UNITMESHUTILS_HPP_ */ diff --git a/packages/krino/krino/unit_tests/Akri_UnitTestUtils.cpp b/packages/krino/krino/unit_tests/Akri_UnitTestUtils.cpp index 6f202d04c033..a1e5696bc0bf 100644 --- a/packages/krino/krino/unit_tests/Akri_UnitTestUtils.cpp +++ b/packages/krino/krino/unit_tests/Akri_UnitTestUtils.cpp @@ -8,20 +8,72 @@ #include #include +#include #include namespace krino { -void expect_eq(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double relativeTol) +static bool is_near_absolute(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double absoluteTol) +{ + for (int i=0; i<3; ++i) + if (std::abs(gold[i]-result[i]) > absoluteTol) + return false; + return true; +} + +static bool is_near_relative(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double relativeTol) { const double absoluteTol = relativeTol * (gold.length() + result.length()); - expect_eq_absolute(gold, result, absoluteTol); + return is_near_absolute(gold, result, absoluteTol); +} + +void expect_eq(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double relativeTol) +{ + EXPECT_TRUE(is_near_relative(gold, result, relativeTol)) << "Failed vector comparison: gold: " << gold << " actual:" << result << " relative tol:" << relativeTol; } void expect_eq_absolute(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double absoluteTol) { - for (int i=0; i<3; ++i) - EXPECT_NEAR(gold[i], result[i], absoluteTol) <<"gold: " << gold << " actual:" << result; + EXPECT_TRUE(is_near_absolute(gold, result, absoluteTol)) << "Failed vector comparison: gold: " << gold << " actual:" << result << " absolute tol:" << absoluteTol;; +} + +bool is_near_relative(const Facet2d & gold, const Facet2d & result, const double relativeTol) +{ + return is_near_relative(gold.facet_vertex(0), result.facet_vertex(0), relativeTol) && + is_near_relative(gold.facet_vertex(1), result.facet_vertex(1), relativeTol); +} + +void expect_eq(const Facet2d & gold, const Facet2d & result, const double relativeTol) +{ + EXPECT_TRUE(is_near_relative(gold, result, relativeTol)) << "Failed Facet2d comparison: gold: " << gold << " actual:" << result << " relative tol:" << relativeTol;; +} + +bool is_near_relative(const Facet3d & gold, const Facet3d & result, const double relativeTol) +{ + const double err0 = (gold.facet_vertex(0) - result.facet_vertex(0)).length_squared(); + const double err1 = (gold.facet_vertex(1) - result.facet_vertex(0)).length_squared(); + const double err2 = (gold.facet_vertex(2) - result.facet_vertex(0)).length_squared(); + if (err1 < err0 && err1 < err2) + { + return is_near_relative(gold.facet_vertex(1), result.facet_vertex(0), relativeTol) && + is_near_relative(gold.facet_vertex(2), result.facet_vertex(1), relativeTol) && + is_near_relative(gold.facet_vertex(0), result.facet_vertex(2), relativeTol); + } + else if (err2 < err0 && err2 < err1) + { + return is_near_relative(gold.facet_vertex(2), result.facet_vertex(0), relativeTol) && + is_near_relative(gold.facet_vertex(0), result.facet_vertex(1), relativeTol) && + is_near_relative(gold.facet_vertex(1), result.facet_vertex(2), relativeTol); + } + + return is_near_relative(gold.facet_vertex(0), result.facet_vertex(0), relativeTol) && + is_near_relative(gold.facet_vertex(1), result.facet_vertex(1), relativeTol) && + is_near_relative(gold.facet_vertex(2), result.facet_vertex(2), relativeTol); +} + +void expect_eq(const Facet3d & gold, const Facet3d & result, const double relativeTol) +{ + EXPECT_TRUE(is_near_relative(gold, result, relativeTol)) << "Failed Facet3d comparison: gold: " << gold << " actual:" << result << " relative tol:" << relativeTol;; } } diff --git a/packages/krino/krino/unit_tests/Akri_UnitTestUtils.hpp b/packages/krino/krino/unit_tests/Akri_UnitTestUtils.hpp index 30e07dd43e14..221dc4e1ab14 100644 --- a/packages/krino/krino/unit_tests/Akri_UnitTestUtils.hpp +++ b/packages/krino/krino/unit_tests/Akri_UnitTestUtils.hpp @@ -12,8 +12,15 @@ namespace krino { +class Facet2d; +class Facet3d; + void expect_eq(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double relativeTol=1.e-6); void expect_eq_absolute(const stk::math::Vector3d & gold, const stk::math::Vector3d & result, const double absoluteTol=1.e-6); +void expect_eq(const Facet2d & gold, const Facet2d & result, const double relativeTol=1.e-6); +void expect_eq(const Facet3d & gold, const Facet3d & result, const double relativeTol=1.e-6); +bool is_near_relative(const Facet2d & gold, const Facet2d & result, const double relativeTol=1.e-6); +bool is_near_relative(const Facet3d & gold, const Facet3d & result, const double relativeTol=1.e-6); } diff --git a/packages/krino/krino/unit_tests/Akri_Unit_Analytic_CDMesh.cpp b/packages/krino/krino/unit_tests/Akri_Unit_Analytic_CDMesh.cpp index 3d606b3a0608..b7f585a4dbc5 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_Analytic_CDMesh.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_Analytic_CDMesh.cpp @@ -243,10 +243,7 @@ class CubeDecompositionFixture : public AnalyticDecompositionFixture facet = std::make_unique( cubeVerts[facetVerts[0]], cubeVerts[facetVerts[1]], cubeVerts[facetVerts[2]] ); - myCube.add( std::move(facet) ); - } + myCube.emplace_back_3d( cubeVerts[facetVerts[0]], cubeVerts[facetVerts[1]], cubeVerts[facetVerts[2]] ); myCubeGeometry.reset(new AnalyticSurfaceInterfaceGeometry({this->surfaceIdentifier}, {&myCube}, AuxMetaData::get(this->fixture.meta_data()).active_part(), this->cdfemSupport, Phase_Support::get(this->fixture.meta_data()))); } @@ -259,7 +256,7 @@ class CubeDecompositionFixture : public AnalyticDecompositionFixture myCube; std::unique_ptr myCubeGeometry; }; diff --git a/packages/krino/krino/unit_tests/Akri_Unit_CDMesh.cpp b/packages/krino/krino/unit_tests/Akri_Unit_CDMesh.cpp index d2ee508662cf..7a6c72cc693a 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_CDMesh.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_CDMesh.cpp @@ -319,7 +319,7 @@ class CompleteDecompositionFixture : public ::testing::Test void decompose_mesh() { NodeToCapturedDomainsMap nodesToSnappedDomains; - std::unique_ptr interfaceGeometry = create_levelset_geometry(krino_mesh->get_active_part(), cdfemSupport, Phase_Support::get(fixture.meta_data()), levelset_fields()); + std::unique_ptr interfaceGeometry = create_levelset_geometry(fixture.meta_data().spatial_dimension(), krino_mesh->get_active_part(), cdfemSupport, Phase_Support::get(fixture.meta_data()), levelset_fields()); if (cdfemSupport.get_cdfem_edge_degeneracy_handling() == SNAP_TO_INTERFACE_WHEN_QUALITY_ALLOWS_THEN_SNAP_TO_NODE) { const double minIntPtWeightForEstimatingCutQuality = cdfemSupport.get_snapper().get_edge_tolerance(); diff --git a/packages/krino/krino/unit_tests/Akri_Unit_Constrained_Redistance.cpp b/packages/krino/krino/unit_tests/Akri_Unit_Constrained_Redistance.cpp index d6b0ccbefda7..58f9240d23bb 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_Constrained_Redistance.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_Constrained_Redistance.cpp @@ -22,7 +22,7 @@ class ConstrainedRedistance : public StkMeshTriFixture ConstrainedRedistance() : logfile(), my_ls(LevelSet::build(mMesh.mesh_meta_data(), "LS", sierra::Diag::sierraTimer())), - myRefinement(mMesh.mesh_meta_data(), &this->get_aux_meta().active_part()) + myRefinement(mMesh.mesh_meta_data(), &this->get_aux_meta().active_part(), sierra::Diag::sierraTimer()) { my_ls.setup(); // registers field and sets field refs on the object @@ -40,11 +40,11 @@ class ConstrainedRedistance : public StkMeshTriFixture for (auto && bucket : mMesh.get_buckets(stk::topology::ELEMENT_RANK, mMesh.mesh_meta_data().locally_owned_part())) { - const int markerValue = - myRefinement.is_parent(*bucket) ? Refinement::NOTHING : Refinement::REFINE; + const auto markerValue = + myRefinement.is_parent(*bucket) ? Refinement::RefinementMarker::NOTHING : Refinement::RefinementMarker::REFINE; auto * elemMarker = field_data(myElementMarkerField, *bucket); for (size_t i = 0; i < bucket->size(); ++i) - elemMarker[i] = markerValue; + elemMarker[i] = static_cast(markerValue); } } diff --git a/packages/krino/krino/unit_tests/Akri_Unit_ContourElement.cpp b/packages/krino/krino/unit_tests/Akri_Unit_ContourElement.cpp index f4b49357bf7a..37d8f871942f 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_ContourElement.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_ContourElement.cpp @@ -61,13 +61,13 @@ TEST(SingleElementFixture, LS_Element_Tet4) const double length_scale = 1.0; // Used for snapping facets to vertices of element when distance is small compared to length_scale ls_elem.compute_subelement_decomposition(length_scale); - Faceted_Surface faceted_surface("tmp"); + Faceted_Surface faceted_surface; ls_elem.build_subelement_facets(faceted_surface); - const auto & facets = faceted_surface.get_facets(); + const auto & facets = faceted_surface.get_facets_3d(); ASSERT_EQ(1u, facets.size()); - Facet & facet = *facets[0]; + const Facet3d & facet = facets[0]; EXPECT_EQ(2.0, facet.facet_vertex(0)[2]); EXPECT_EQ(2.0, facet.facet_vertex(1)[2]); @@ -111,18 +111,18 @@ class ContourElementFixture : public StkMeshFixture *dist = nodalLevelset[n]; } } - void build_facets_in_element(const stk::mesh::Entity elem, Faceted_Surface & facetedSurface) + void build_facets_in_element(const stk::mesh::Entity elem, FacetedSurfaceBase & facetedSurface) { facetedSurface.clear(); ContourElement contourElement(mMesh, elem, myLs.get_coordinates_field(), myLs.get_distance_field(), 0.); contourElement.compute_subelement_decomposition(avgElemSize); contourElement.build_subelement_facets(facetedSurface); } - void expect_num_facets_in_element(const stk::mesh::Entity elem, const size_t goldNumFacets) + void expect_num_facets_in_element(const int dim, const stk::mesh::Entity elem, const size_t goldNumFacets) { - Faceted_Surface facetedSurface("tmp"); - build_facets_in_element(elem, facetedSurface); - EXPECT_EQ(goldNumFacets, facetedSurface.size()); + std::unique_ptr facetedSurface = FacetedSurfaceBase::build(dim); + build_facets_in_element(elem, *facetedSurface); + EXPECT_EQ(goldNumFacets, facetedSurface->size()); } protected: @@ -142,8 +142,8 @@ TEST_F(ContourElementTwoTri306090, twoElementsWhereSnappedSignPatternIsDifferent const std::vector & elems = get_owned_elements(); ASSERT_EQ(2u, elems.size()); - expect_num_facets_in_element(elems[0], 0); - expect_num_facets_in_element(elems[1], 1); + expect_num_facets_in_element(2, elems[0], 0); + expect_num_facets_in_element(2, elems[1], 1); } } @@ -158,8 +158,8 @@ TEST_F(ContourElementTwoQuads, twoQuadElementsWithSmallVariationAcrossElementAnd const std::vector & elems = get_owned_elements(); ASSERT_EQ(2u, elems.size()); - expect_num_facets_in_element(elems[0], 1); - expect_num_facets_in_element(elems[1], 0); + expect_num_facets_in_element(2, elems[0], 1); + expect_num_facets_in_element(2, elems[1], 0); } } diff --git a/packages/krino/krino/unit_tests/Akri_Unit_DecomposeWithSensitivities.cpp b/packages/krino/krino/unit_tests/Akri_Unit_DecomposeWithSensitivities.cpp index e566123323e1..081e9a0087f2 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_DecomposeWithSensitivities.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_DecomposeWithSensitivities.cpp @@ -36,7 +36,7 @@ static void decompose_mesh_to_conform_to_levelsets(stk::mesh::BulkData & mesh, c krino::CDFEM_Support & cdfemSupport = krino::CDFEM_Support::get(meta); krino::Phase_Support & phaseSupport = krino::Phase_Support::get(meta); - std::unique_ptr interfaceGeometry = krino::create_levelset_geometry(auxMeta.active_part(), cdfemSupport, phaseSupport, lsFields); + std::unique_ptr interfaceGeometry = krino::create_levelset_geometry(meta.spatial_dimension(), auxMeta.active_part(), cdfemSupport, phaseSupport, lsFields); auxMeta.clear_force_64bit_flag(); krino::CDMesh::decompose_mesh(mesh, *interfaceGeometry); } diff --git a/packages/krino/krino/unit_tests/Akri_Unit_DecompositionFixture.hpp b/packages/krino/krino/unit_tests/Akri_Unit_DecompositionFixture.hpp index 3e8688246e0c..bda6cee3a7b0 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_DecompositionFixture.hpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_DecompositionFixture.hpp @@ -151,7 +151,7 @@ void check_nonfatal_error(const std::string & baseName, const int iteration) void decompose_mesh() { NodeToCapturedDomainsMap nodesToSnappedDomains; - std::unique_ptr interfaceGeometry = create_levelset_geometry(cdmesh->get_active_part(), cdfem_support(), Phase_Support::get(mMesh.mesh_meta_data()), levelset_fields()); + std::unique_ptr interfaceGeometry = create_levelset_geometry(mMesh.mesh_meta_data().spatial_dimension(), cdmesh->get_active_part(), cdfem_support(), Phase_Support::get(mMesh.mesh_meta_data()), levelset_fields()); if (cdfem_support().get_cdfem_edge_degeneracy_handling() == SNAP_TO_INTERFACE_WHEN_QUALITY_ALLOWS_THEN_SNAP_TO_NODE) { const double minIntPtWeightForEstimatingCutQuality = cdfem_support().get_snapper().get_edge_tolerance(); diff --git a/packages/krino/krino/unit_tests/Akri_Unit_Eikonal.cpp b/packages/krino/krino/unit_tests/Akri_Unit_Eikonal.cpp index 80e941c5be9b..bdfc4dab6815 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_Eikonal.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_Eikonal.cpp @@ -45,6 +45,22 @@ void expect_tetrahedron_eikonal_solution(const std::array EXPECT_NEAR(gold, resultCase, 1.e-8); } +void expect_tetrahedron_distance_and_extension_speed(const std::array & x, const std::array & d, const std::array & extSpeed, const double goldDist, const double goldExtSpeed) +{ + const int sign = 1; + const auto [nodeDist, nodeSpeed] = eikonal_solve_tetrahedron_for_distance_and_extension_speed(x, d, extSpeed, sign, farDistance); + EXPECT_NEAR(goldDist, nodeDist, 1.e-8); + EXPECT_NEAR(goldExtSpeed, nodeSpeed, 1.e-8); +} + +void expect_triangle_distance_and_extension_speed(const std::array & x, const std::array & d, const std::array & extSpeed, const double goldDist, const double goldExtSpeed) +{ + const int sign = 1; + const auto [nodeDist, nodeSpeed] = eikonal_solve_triangle_for_distance_and_extension_speed(x, d, extSpeed, sign, farDistance); + EXPECT_NEAR(goldDist, nodeDist, 1.e-8); + EXPECT_NEAR(goldExtSpeed, nodeSpeed, 1.e-8); +} + void expect_tetrahedron_distance(const std::array & x, const std::array & d, const double gold) { const double speed = 1.0; @@ -140,6 +156,17 @@ TEST(Eikonal_Calc, regularTet_nbrsAllNearZero_resultIsTetHeight) expect_tetrahedron_distance(get_regular_tetrahedron_coordinates(), distance, gold); } +TEST(Eikonal_Calc, regularTet_nbrsAllZero_distanceIsHeightAndExtensionFromFaceCentroid) +{ + const std::array distance{{0.,0.,0.}}; + const std::array extSpeed{{0.,1.,2.}}; + + const double edgeLen = get_regular_tetrahedron_edge_length(); + const double goldDist = std::sqrt(2./3.)*edgeLen; + const double goldExtSpeed = (extSpeed[0]+extSpeed[1]+extSpeed[2])/3.; + expect_tetrahedron_distance_and_extension_speed(get_regular_tetrahedron_coordinates(), distance, extSpeed, goldDist, goldExtSpeed); +} + TEST(Eikonal_Calc, regularTet_gradienAlong0to3_resultIsEdgeLength) { const std::array coords = get_regular_tetrahedron_coordinates(); @@ -150,6 +177,18 @@ TEST(Eikonal_Calc, regularTet_gradienAlong0to3_resultIsEdgeLength) expect_tetrahedron_distance(coords, distance, gold); } +TEST(Eikonal_Calc, regularTet_gradienAlong0to3_distanceIsEdgeLengthAndExtensionFromNode0) +{ + const std::array coords = get_regular_tetrahedron_coordinates(); + const double edgeLen = get_regular_tetrahedron_edge_length(); + const std::array distance{{0.,Dot(coords[1]-coords[0],coords[3]-coords[0])/edgeLen,Dot(coords[2]-coords[0],coords[3]-coords[0])/edgeLen}}; + const std::array extSpeed{{1.,2.,3.}}; + + const double goldDist = edgeLen; + const double goldExtSpeed = extSpeed[0]; + expect_tetrahedron_distance_and_extension_speed(get_regular_tetrahedron_coordinates(), distance, extSpeed, goldDist, goldExtSpeed); +} + TEST(Eikonal_Calc, regularTri_nbrsAllZero_resultIsTriHeight) { const std::array distance{{0.,0.}}; @@ -158,6 +197,16 @@ TEST(Eikonal_Calc, regularTri_nbrsAllZero_resultIsTriHeight) expect_triangle_distance(get_regular_triangle_coordinates(), distance, gold); } +TEST(Eikonal_Calc, regularTri_nbrsAllZero_distanceIsHeightAndExtensionFromEdgeMidpt) +{ + const std::array distance{{0.,0.}}; + const std::array extSpeed{{0.,1.}}; + + const double goldDist = 0.5*std::sqrt(3.0)*get_regular_triangle_edge_length(); + const double goldExtSpeed = (extSpeed[0]+extSpeed[1])/2.; + expect_triangle_distance_and_extension_speed(get_regular_triangle_coordinates(), distance, extSpeed, goldDist, goldExtSpeed); +} + TEST(Eikonal_Calc, regularTri_gradienAlong0to2_resultIsEdgeLength) { const std::array coords = get_regular_triangle_coordinates(); @@ -168,6 +217,18 @@ TEST(Eikonal_Calc, regularTri_gradienAlong0to2_resultIsEdgeLength) expect_triangle_distance(coords, distance, gold); } +TEST(Eikonal_Calc, regularTri_gradienAlong0to2_distanceIsEdgeLengthAndExtensionFromNode0) +{ + const std::array coords = get_regular_triangle_coordinates(); + const double edgeLen = get_regular_triangle_edge_length(); + const std::array distance{{0.,Dot(coords[1]-coords[0],coords[2]-coords[0])/edgeLen}}; + const std::array extSpeed{{1.,2.}}; + + const double goldDist = edgeLen; + const double goldExtSpeed = extSpeed[0]; + expect_triangle_distance_and_extension_speed(coords, distance, extSpeed, goldDist, goldExtSpeed); +} + TEST(Eikonal_Calc, rightTriangleWithOneNeighborNodeAsClosestPointAndOtherNbrOnLongerEdge_resultIsDistanceFromClosestPoint) { const std::array coords{{ {0.,0.,0.}, {2.,0.,0.}, {0.,1.,0.} }}; @@ -199,6 +260,29 @@ TEST(Eikonal_Calc, regularTriangle_withFrontEmanatingFromOutsideBase_resultIsEdg } } +TEST(Eikonal_Calc, regularTriangle_withFrontEmanatingFromOutsideBase_distanceIsEdgeLengthAndExtensionFromNearestNode) +{ + const std::array coords = get_regular_triangle_coordinates(); + const double goldDist = 1.0; + const std::array extSpeed{{1.,2.}}; + + { + const std::array distance{{0., 0.9}}; + expect_triangle_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[0]); + } + + { + const std::array distance{{0.9, 0.}}; + expect_triangle_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[1]); + } + + // detJ = 0 case + { + const std::array distance{{0., 1.}}; + expect_triangle_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[0]); + } +} + TEST(Eikonal_Calc, rightTetrahedronWithOneNeighborNodeAsClosestPointAndAnotherNbrOnLongerEdge_resultIsDistanceFromClosestPoint) { const std::array coords{{ {0.,0.,0.}, {2.,0.,0.}, {0.,1.,0.}, {0.,0.,1.} }}; @@ -226,7 +310,6 @@ TEST(Eikonal_Calc, rightTetrahedron_nbrsAllZeroAndUnitSpeed_resultIsDistance) expect_tetrahedron_eikonal_solution(get_right_tetrahedron_coordinates(), distance, speed, gold); } - TEST(Eikonal_Calc, rightTetrahedron_nbrsAllZeroAndNonUnitSpeed_resultIsTimeOfArrival) { const std::array distance{{0.,0.,0.}}; @@ -354,6 +437,29 @@ TEST(Eikonal_Calc, regularTetrahedron_withFrontEmanatingFromOutsideBase_resultIs } } +TEST(Eikonal_Calc, regularTetrahedron_withFrontEmanatingFromOutsideBase_distanceIsEdgeLengthAndExtensionFromNearestNode) +{ + const std::array coords = get_regular_tetrahedron_coordinates(); + const double edgeLen = get_regular_tetrahedron_edge_length(); + const double goldDist = edgeLen; + const std::array extSpeed{{1.,2.,3.}}; + + { + const std::array distance{{0., 0.6*edgeLen, 0.6*edgeLen}}; + expect_tetrahedron_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[0]); + } + + { + const std::array distance{{0.6*edgeLen, 0., 0.6*edgeLen}}; + expect_tetrahedron_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[1]); + } + + { + const std::array distance{{0.6*edgeLen, 0.6*edgeLen, 0.}}; + expect_tetrahedron_distance_and_extension_speed(coords, distance, extSpeed, goldDist, extSpeed[2]); + } +} + TEST(Eikonal_Calc, rightTriangle_usingMultipleSpeeds_resultIsLengthOverSpeed) { const std::array distance{{0.,0.}}; diff --git a/packages/krino/krino/unit_tests/Akri_Unit_RebalanceUtils.cpp b/packages/krino/krino/unit_tests/Akri_Unit_RebalanceUtils.cpp index fc413ff69f6f..2d4eb56c4888 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_RebalanceUtils.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_RebalanceUtils.cpp @@ -257,7 +257,7 @@ class ParallelRebalanceForAdaptivityFixture3D : public ::testing::Test if(std::fabs(coords[meta->spatial_dimension()-1]-1.) <= 1e-6) { int * elemMarker = field_data(refinement.get_marker_field(), elem); - *elemMarker = Refinement::REFINE; + *elemMarker = static_cast(Refinement::RefinementMarker::REFINE); } } } @@ -519,7 +519,7 @@ TEST_F(ParallelRebalanceForAdaptivityFixture3D, ParentChildRebalanceRules) EXPECT_TRUE(parent != stk::mesh::Entity()); EXPECT_TRUE(refinement.is_parent(parent)); int * elemMarker = field_data(markerField, elem); - *elemMarker = Refinement::COARSEN; + *elemMarker = static_cast(Refinement::RefinementMarker::COARSEN); } } refinement.do_refinement(); diff --git a/packages/krino/krino/unit_tests/Akri_Unit_RefineInterval.cpp b/packages/krino/krino/unit_tests/Akri_Unit_RefineInterval.cpp index 0dcd8f13b2b5..64871468e252 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_RefineInterval.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_RefineInterval.cpp @@ -7,6 +7,7 @@ #include #include #include +#include stk::mesh::Field & create_levelset_field(stk::mesh::BulkData & mesh, const std::string & levelsetName) @@ -84,5 +85,12 @@ TEST(RefineWithinDistanceOfLevelSets, twoIntersectingLevelSets_correctElementsGe test_refinement_level_for_elements_attached_to_node(mesh, activePart, meshAndBuilder.get_assigned_node_for_index(2), 2); test_refinement_level_for_elements_attached_to_node(mesh, activePart, meshAndBuilder.get_assigned_node_for_index(4), numRefinementLevels); + + if (mesh.parallel_rank() == 0) // Whole mesh on 0 + { + stk::mesh::Entity nodeOffCenter = krino::find_local_node_closest_to_location(mesh, mesh.mesh_meta_data().universal_part(), meshAndBuilder.get_coordinates_field(), stk::math::Vector3d(0.05,-0.05,0)); + std::cout << "Checking element size for elements using node " << mesh.identifier(nodeOffCenter) << std::endl; + test_refinement_level_for_elements_attached_to_node(mesh, activePart, nodeOffCenter, numRefinementLevels); + } } diff --git a/packages/krino/krino/unit_tests/Akri_Unit_RefinementFixture.hpp b/packages/krino/krino/unit_tests/Akri_Unit_RefinementFixture.hpp index 335ea7f76125..a48c64ff42ba 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_RefinementFixture.hpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_RefinementFixture.hpp @@ -18,14 +18,14 @@ namespace krino { -inline void set_refinement_marker_field(FieldRef elementMarkerField, const int value) +inline void set_refinement_marker_field(FieldRef elementMarkerField, const Refinement::RefinementMarker value) { - stk::mesh::field_fill(value, elementMarkerField, stk::mesh::selectField(elementMarkerField)); + stk::mesh::field_fill(static_cast(value), elementMarkerField, stk::mesh::selectField(elementMarkerField)); } inline void clear_refinement_marker_field(FieldRef elementMarkerField) { - set_refinement_marker_field(elementMarkerField, Refinement::NOTHING); + set_refinement_marker_field(elementMarkerField, Refinement::RefinementMarker::NOTHING); } @@ -34,7 +34,7 @@ class RefinementFixture : public StkMeshFixture { public: RefinementFixture() - : myRefinement(mMesh.mesh_meta_data(), &this->get_aux_meta().active_part()) + : myRefinement(mMesh.mesh_meta_data(), &this->get_aux_meta().active_part(), sierra::Diag::sierraTimer()) { stk::mesh::MetaData & meta = mMesh.mesh_meta_data(); stk::mesh::FieldBase & elemMarkerField = meta.declare_field(stk::topology::ELEMENT_RANK, myElementMarkerFieldName, 1); @@ -55,7 +55,7 @@ class RefinementFixture : public StkMeshFixture for (auto && elem : elements) { int * elemMarker = field_data(myElementMarkerField, elem); - *elemMarker = Refinement::REFINE; + *elemMarker = static_cast(Refinement::RefinementMarker::REFINE); } } void mark_elements_for_unrefinement(const std::vector & elements) @@ -63,17 +63,17 @@ class RefinementFixture : public StkMeshFixture for (auto && elem : elements) { int * elemMarker = field_data(myElementMarkerField, elem); - *elemMarker = Refinement::COARSEN; + *elemMarker = static_cast(Refinement::RefinementMarker::COARSEN); } } void mark_nonparent_elements() { for ( auto && bucket : mMesh.get_buckets( stk::topology::ELEMENT_RANK, mMesh.mesh_meta_data().locally_owned_part() ) ) { - const int markerValue = myRefinement.is_parent(*bucket) ? Refinement::NOTHING : Refinement::REFINE; + const auto markerValue = myRefinement.is_parent(*bucket) ? Refinement::RefinementMarker::NOTHING : Refinement::RefinementMarker::REFINE; auto * elemMarker = field_data(myElementMarkerField, *bucket); for (size_t i=0; isize(); ++i) - elemMarker[i] = markerValue; + elemMarker[i] = static_cast(markerValue); } } @@ -167,7 +167,7 @@ class RefinementFixture : public StkMeshFixture auto doMarkAndFieldVal = tag_element_spans_z_equal_0_and_compute_element_field((*bucket)[iElem], flip); auto doMark = std::get<0>(doMarkAndFieldVal); auto centroidValue = std::get<1>(doMarkAndFieldVal); - if (doMark) elemMarker[iElem] = Refinement::REFINE; + if (doMark) elemMarker[iElem] = static_cast(Refinement::RefinementMarker::REFINE); if (centroidValue <= 0 ) elemField[iElem] = -1; else elemField[iElem] = 1; @@ -208,7 +208,7 @@ class RefinementFixture : public StkMeshFixture for ( size_t iElem=0; iElemsize(); ++iElem ) { if (element_spans_x_equal_0((*bucket)[iElem])) - elemMarker[iElem] = Refinement::REFINE; + elemMarker[iElem] = static_cast(Refinement::RefinementMarker::REFINE); } } } @@ -253,7 +253,7 @@ class RefinementFixture : public StkMeshFixture const double rand_val = rand_dist(rand_gen); if(rand_val < refine_prob) { - elemMarker[iElem] = Refinement::REFINE; + elemMarker[iElem] = static_cast(Refinement::RefinementMarker::REFINE); } } } @@ -266,10 +266,10 @@ class RefinementFixture : public StkMeshFixture const double rand_val = rand_dist(rand_gen); if(rand_val < unrefine_prob) { - elemMarker[iElem] = Refinement::COARSEN; + elemMarker[iElem] = static_cast(Refinement::RefinementMarker::COARSEN); for(auto && child : myRefinement.get_children((*bucket)[iElem])) { - *field_data(myElementMarkerField, child) = Refinement::COARSEN; + *field_data(myElementMarkerField, child) = static_cast(Refinement::RefinementMarker::COARSEN); } } } @@ -333,12 +333,12 @@ class RefinementFixture : public StkMeshFixture void mark_all_elements_for_unrefinement() { - set_refinement_marker_field(myElementMarkerField, Refinement::COARSEN); + set_refinement_marker_field(myElementMarkerField, Refinement::RefinementMarker::COARSEN); } void mark_all_elements_for_refinement() { - set_refinement_marker_field(myElementMarkerField, Refinement::REFINE); + set_refinement_marker_field(myElementMarkerField, Refinement::RefinementMarker::REFINE); } void test_refinement_of_transition_element_leads_to_refinement_of_parent(const int indexOfCenterElement) diff --git a/packages/krino/krino/unit_tests/Akri_Unit_Snap.cpp b/packages/krino/krino/unit_tests/Akri_Unit_Snap.cpp index fb530f4b341d..37577c9616fb 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_Snap.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_Snap.cpp @@ -442,7 +442,7 @@ class SnapTrisFixture : public StkMeshTriFixture IntersectionPointFromNodalLevelsetInterfaceGeometry geometry; geometry.set_nodal_levelset(mMesh, get_assigned_node_global_ids(), nodeLSValues); NodeToCapturedDomainsMap nodesToCapturedDomains; - return build_all_intersection_points(mMesh, geometry, nodesToCapturedDomains); + return build_all_intersection_points(mMesh, mMesh.mesh_meta_data().universal_part(), geometry, nodesToCapturedDomains); } std::vector get_allowed_snap_indices(const std::vector & intersectionPoints, const stk::mesh::Entity snapNode) diff --git a/packages/krino/krino/unit_tests/Akri_Unit_TriangleCalcs.cpp b/packages/krino/krino/unit_tests/Akri_Unit_TriangleCalcs.cpp index 893b6f09feb9..0138d1a834bb 100644 --- a/packages/krino/krino/unit_tests/Akri_Unit_TriangleCalcs.cpp +++ b/packages/krino/krino/unit_tests/Akri_Unit_TriangleCalcs.cpp @@ -10,7 +10,7 @@ #include #include #include - +#include namespace { diff --git a/packages/krino/krino_mesh_adapt/KrinoMeshAdaptMain.cpp b/packages/krino/krino_mesh_adapt/KrinoMeshAdaptMain.cpp new file mode 100644 index 000000000000..372af915e3e9 --- /dev/null +++ b/packages/krino/krino_mesh_adapt/KrinoMeshAdaptMain.cpp @@ -0,0 +1,122 @@ +/*--------------------------------------------------------------------*/ +/* Copyright 2002 - 2008, 2010, 2011 National Technology & */ +/* Engineering Solutions of Sandia, LLC (NTESS). Under the terms */ +/* of Contract DE-NA0003525 with NTESS, there is a */ +/* non-exclusive license for use of this work by or on behalf */ +/* of the U.S. Government. Export of this program may require */ +/* a license from the United States Government. */ +/*--------------------------------------------------------------------*/ + +// ####################### Start Clang Header Tool Managed Headers ######################## +// clang-format off +#include +#include +#include // for MPI_Comm_rank, MPI_Comm_size +#include +#include +#include // for operator<<, basic_ostream:... +#include // for unique_ptr +#include // for parallel_machine_finalize +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for ProductRegistry + +void krino_mesh_adapt(const krino::MeshAdaptInputData &inputData, const stk::ParallelMachine comm) +{ + std::shared_ptr bulk = stk::mesh::MeshBuilder(comm).create(); + stk::mesh::MetaData& meta = bulk->mesh_meta_data(); + meta.use_simple_fields(); + meta.enable_late_fields(); + + stk::io::fill_mesh_with_auto_decomp(inputData.meshIn, *bulk); + + krino::refine_mesh_with_params(*bulk, inputData.algorithmParams, bulk->parallel()); + + Ioss::PropertyManager properties; + if (inputData.algorithmParams.autoCompose) + properties.add(Ioss::Property("COMPOSE_RESULTS", 1)); + + stk::io::StkMeshIoBroker io_broker; + io_broker.set_bulk_data(bulk); + size_t resultFileIndex = io_broker.create_output_mesh(inputData.meshOut, stk::io::WRITE_RESULTS, properties); + io_broker.write_output_mesh(resultFileIndex); +} + +bool run_krino_mesh_adapt(int argc, char **argv, stk::ParallelMachine comm) +{ + const double startWallTime = stk::wall_time(); + const double startCpuTime = stk::cpu_time(); + + int proc = stk::parallel_machine_rank(comm); + if(proc != 0) + stk::EnvData::instance().m_outputP0 = &stk::EnvData::instance().m_outputNull; + + krino::MeshAdaptParser parser; + bool didParseOk = parser.read_command_line(argc, argv, comm); + + bool successfullyMeshed{false}; + double elapsedWallTime{0.0}; + double elapsedCpuTime{0.0}; + + if (didParseOk) + { + const krino::MeshAdaptInputData &inputData = parser.get_input_data(); + + krino_mesh_adapt(inputData, comm); + + stk::parallel_machine_barrier(comm); + + elapsedWallTime = stk::wall_time() - startWallTime; + elapsedCpuTime = stk::cpu_time() - startCpuTime; + + if (0 == stk::parallel_machine_rank(comm)) + std::cout << "Elapsed wall time: " + std::to_string(elapsedWallTime) + " seconds" << std::endl; + + if (0 == stk::parallel_machine_rank(comm)) + std::cout << "Elapsed cpu time: " + std::to_string(elapsedCpuTime) + " seconds" << std::endl; + + successfullyMeshed = true; + } + +// const double startTime{static_cast(::time(nullptr))}; +// common_utils::write_audit_file( +// "krino_mesh_adapt", +// "mesh refinement", +// startTime, +// elapsedWallTime, +// elapsedCpuTime, +// successfullyMeshed, +// [&](const auditdata &data) {OutputAuditLog(&data);}, +// comm); + return successfullyMeshed; +} + +int main(int argc, char **argv) +{ + stk::ParallelMachine comm{stk::parallel_machine_init(&argc, &argv)}; + + bool successfullyMeshed{false}; + + try + { + successfullyMeshed = run_krino_mesh_adapt(argc, argv, comm); + } + catch(std::exception &e) + { + std::cerr << "Proc " << stk::parallel_machine_rank(comm) << ": " << e.what() << std::endl; + const int errorCode{-1}; + MPI_Abort(comm, errorCode); + } + + stk::parallel_machine_finalize(); + + return successfullyMeshed ? 0 : 1; +} diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.cpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.cpp new file mode 100644 index 000000000000..5a4d1b632b12 --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.cpp @@ -0,0 +1,42 @@ +/*--------------------------------------------------------------------*/ +/* Copyright 2002 - 2008, 2010, 2011 National Technology & */ +/* Engineering Solutions of Sandia, LLC (NTESS). Under the terms */ +/* of Contract DE-NA0003525 with NTESS, there is a */ +/* non-exclusive license for use of this work by or on behalf */ +/* of the U.S. Government. Export of this program may require */ +/* a license from the United States Government. */ +/*--------------------------------------------------------------------*/ + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace krino { + +void refine_mesh_with_params(stk::mesh::BulkData & mesh, const MeshAdaptAlgorithmParameters &algParams, const stk::ParallelMachine comm) +{ + stk::diag::Timer refinementTimer("Refinement", sierra::Diag::sierraTimer()); + stk::mesh::Part * activePart = nullptr; + Refinement refinement(mesh.mesh_meta_data(), activePart, algParams.force64Bit, algParams.assert32Bit, refinementTimer); + + std::vector counts; + stk::mesh::comm_mesh_counts(mesh, counts); + + sierra::Env::outputP0() << "Uniform refinement: before refine, mesh has " << counts[0] << " nodes, " << counts[1] + << " edges, " << counts[2] << " faces, " << counts[3] << " elements" << std::endl; + + refinement.do_uniform_refinement(algParams.numUniformRefinementLevels); + + stk::mesh::comm_mesh_counts(mesh, counts); + + sierra::Env::outputP0() << "Uniform refinement: after refine, mesh has " << counts[0] << " nodes, " << counts[1] + << " edges, " << counts[2] << " faces, " << counts[3] << " elements" << std::endl; +} + +} + diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.hpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.hpp new file mode 100644 index 000000000000..c23346e2b711 --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdapt.hpp @@ -0,0 +1,18 @@ +#ifndef KrinoMeshAdapt_hpp +#define KrinoMeshAdapt_hpp + +#include +#include + +namespace krino { struct MeshAdaptAlgorithmParameters; } + +namespace krino +{ + +void refine_mesh_with_params(stk::mesh::BulkData & mesh, + const MeshAdaptAlgorithmParameters &algParams, + const stk::ParallelMachine comm); + +} + +#endif diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptAlgorithmParameters.hpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptAlgorithmParameters.hpp new file mode 100644 index 000000000000..40e499dacd92 --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptAlgorithmParameters.hpp @@ -0,0 +1,17 @@ +#ifndef KRINOMESHADAPTALGORITHMPARAMETERS_HPP_ +#define KRINOMESHADAPTALGORITHMPARAMETERS_HPP_ + +namespace krino +{ + +struct MeshAdaptAlgorithmParameters +{ + bool force64Bit{false}; + bool assert32Bit{false}; + bool autoCompose{true}; + unsigned numUniformRefinementLevels{1}; +}; + +} + +#endif /* KRINOMESHADAPTALGORITHMPARAMETERS_HPP_ */ diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptInputData.hpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptInputData.hpp new file mode 100644 index 000000000000..2b965e3a7fdc --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptInputData.hpp @@ -0,0 +1,27 @@ +#ifndef KrinoMeshAdaptInputData_hpp +#define KrinoMeshAdaptInputData_hpp + +#include +#include + +namespace krino +{ + +struct MeshAdaptInputData +{ + static constexpr const char * mDefaultString{"None"}; + std::string meshIn{mDefaultString}; + std::string meshOut{mDefaultString}; + MeshAdaptAlgorithmParameters algorithmParams; + + MeshAdaptInputData() {} + MeshAdaptInputData(const std::string& iMeshIn, + const std::string& iMeshOut) : + meshIn(iMeshIn), + meshOut(iMeshOut) + {} +}; + +} + +#endif diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.cpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.cpp new file mode 100644 index 000000000000..fa53229c8eae --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.cpp @@ -0,0 +1,70 @@ +// ####################### Start Clang Header Tool Managed Headers ######################## +// clang-format off +#include +#include +#include +#include +#include +#include // for exception +#include // for string +// clang-format on +// ####################### End Clang Header Tool Managed Headers ######################## + +namespace krino +{ + +static std::string make_flag(const stk::CommandLineOption &option) +{ + return option.name+","+option.abbreviation; +} + +bool MeshAdaptParser::read_command_line( int argc, char *argv[], stk::ParallelMachine comm) +{ + MeshAdaptInputData defaultValues; + + const stk::CommandLineOption inmeshOption{"input_mesh", "i", "Filename of input genesis mesh."}; + const stk::CommandLineOption outmeshOption{"output_mesh", "o", "Filename of output genesis refined mesh to write."}; + + const stk::CommandLineOption numUniformRefinementLevelsOption{"--number_refines", "n", "Number of uniform refinement passes. Default is 1."}; + const stk::CommandLineOption force64BitOption{"force_64_bit", "", "Force all new nodes and elements past 32bit limit."}; + const stk::CommandLineOption assert32BitOption{"assert_32_bit", "", "Assert that all nodes and elements are within the 32bit limit."}; + const stk::CommandLineOption noAutoComposeOption{"no_auto_compose", "", "Do not automatically recompose the parallel files."}; + + stk::CommandLineParserParallel commandLine(comm); + commandLine.disallow_unrecognized(); + + commandLine.add_required(inmeshOption); + commandLine.add_required(outmeshOption); + commandLine.add_optional(numUniformRefinementLevelsOption, mInputData.algorithmParams.numUniformRefinementLevels); + commandLine.add_flag(make_flag(force64BitOption), force64BitOption.description); + commandLine.add_flag(make_flag(assert32BitOption), assert32BitOption.description); + commandLine.add_flag(make_flag(noAutoComposeOption), noAutoComposeOption.description); + + stk::CommandLineParser::ParseState state = commandLine.parse(argc, const_cast(argv)); + if(state == stk::CommandLineParser::ParseComplete) + { + mInputData.meshIn = commandLine.get_option_value(inmeshOption.name); + mInputData.meshOut = commandLine.get_option_value(outmeshOption.name); + + mInputData.algorithmParams.force64Bit = commandLine.is_option_provided(force64BitOption.name); + mInputData.algorithmParams.assert32Bit = commandLine.is_option_provided(assert32BitOption.name); + mInputData.algorithmParams.autoCompose = !commandLine.is_option_provided(noAutoComposeOption.name); + + if(commandLine.is_option_parsed(numUniformRefinementLevelsOption.name)) + { + mInputData.algorithmParams.numUniformRefinementLevels = commandLine.get_option_value(numUniformRefinementLevelsOption.name); + if (mInputData.algorithmParams.numUniformRefinementLevels < 1) + { + sierra::Env::outputP0() << "ERROR: Number of uniform refinement passes specified via --"+numUniformRefinementLevelsOption.name+" must be >= 1." << std::endl; + state = stk::CommandLineParser::ParseError; + } + } + } + + if(state != stk::CommandLineParser::ParseComplete) + sierra::Env::outputP0() << commandLine.get_usage() << std::endl; + + return state == stk::CommandLineParser::ParseComplete; +} + +} diff --git a/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.hpp b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.hpp new file mode 100644 index 000000000000..8b56ba8c067c --- /dev/null +++ b/packages/krino/krino_mesh_adapt/mesh_adapt_lib/KrinoMeshAdaptParser.hpp @@ -0,0 +1,27 @@ +#ifndef KrinoMeshAdaptParser_hpp +#define KrinoMeshAdaptParser_hpp + +// ####################### Start Clang Header Tool Managed Headers ######################## +// clang-format off +#include +#include +// clang-format on +// ####################### End Clang Header Tool Managed Headers ######################## + +namespace krino +{ + +class MeshAdaptParser +{ +public: + const MeshAdaptInputData& get_input_data() const { return mInputData; } + bool read_command_line(int argc, char *argv[], stk::ParallelMachine comm); + +private: + MeshAdaptInputData mInputData; +}; + +} + + +#endif