diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index 65a8edb8..b380472a 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -114,7 +114,7 @@ void crackBranchingExample( const std::string filename ) // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, prenotch ); + inputs, particles, force_model ); // ==================================================== // Boundary conditions @@ -137,7 +137,7 @@ void crackBranchingExample( const std::string filename ) // ==================================================== // Simulation run // ==================================================== - cabana_pd->init(); + cabana_pd->init( bc, prenotch ); cabana_pd->run( bc ); } diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index 0c84a553..94dfaa28 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -85,7 +85,8 @@ void fragmentingCylinderExample( const std::string filename ) { double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + ( x[1] - y_center ) * ( x[1] - y_center ); - if ( rsq < Rin * Rin || rsq > Rout * Rout ) + if ( rsq < Rin * Rin || rsq > Rout * Rout || x[2] > 0.05 || + x[2] < -0.05 ) return false; return true; }; @@ -137,12 +138,29 @@ void fragmentingCylinderExample( const std::string filename ) particles->updateParticles( exec_space{}, init_functor ); // ==================================================== - // Simulation run + // Simulation run with contact physics // ==================================================== - auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model ); - cabana_pd->init(); - cabana_pd->run(); + if ( inputs["use_contact"] ) + { + double r_c = inputs["contact_horizon_factor"]; + r_c *= dx[0]; + CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model, contact_model ); + cabana_pd->init(); + cabana_pd->run(); + } + // ==================================================== + // Simulation run without contact + // ==================================================== + else + { + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model ); + cabana_pd->init(); + cabana_pd->run(); + } } // Initialize MPI+Kokkos. diff --git a/examples/mechanics/inputs/fragmenting_cylinder.json b/examples/mechanics/inputs/fragmenting_cylinder.json index a6b7041b..cc1deba3 100644 --- a/examples/mechanics/inputs/fragmenting_cylinder.json +++ b/examples/mechanics/inputs/fragmenting_cylinder.json @@ -1,12 +1,15 @@ { "num_cells" : {"value": [51, 51, 101]}, - "system_size" : {"value": [0.05, 0.05, 0.1], "unit": "m"}, + "dx": {"value": [0.001, 0.001, 0.001]}, + "system_size" : {"value": [0.1, 0.1, 0.15], "unit": "m"}, "grid_perturbation_factor" : {"value": 0.1}, "density" : {"value": 7800, "unit": "kg/m^3"}, "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, "critical_stretch" : {"value": 0.02}, "horizon" : {"value": 0.00417462, "unit": "m"}, + "use_contact" : {"value": true}, + "contact_horizon_factor" : {"value": 0.9}, "cylinder_outer_radius" : {"value": 0.025, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.02, "unit": "m"}, "max_radial_velocity" : {"value": 200, "unit": "m/s"}, diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index 78dab593..b6a3372b 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -116,7 +116,7 @@ void kalthoffWinklerExample( const std::string filename ) // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, prenotch ); + inputs, particles, force_model ); // ==================================================== // Boundary conditions @@ -132,7 +132,7 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== // Simulation run // ==================================================== - cabana_pd->init(); + cabana_pd->init( bc, prenotch ); cabana_pd->run( bc ); } diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 7da7ad22..f1014885 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -29,8 +29,10 @@ #include #include +#include #include #include +#include #include #include diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 6ab1508c..5ee56eb3 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -286,7 +286,9 @@ class Comm _init_timer.stop(); } - ~Comm() {} + + auto size() { return mpi_size; } + auto rank() { return mpi_rank; } // Determine which particles should be ghosted, reallocating and recounting // if needed. diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 6dab7f08..8299c78b 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -156,7 +156,20 @@ class Force { } - // Primary constructor: use positions and construct neighbors. + // General constructor (necessary for contact, but could be used by any + // force routine). + template + Force( const bool half_neigh, const double delta, + const PositionType& positions, const std::size_t num_local, + const double mesh_min[3], const double mesh_max[3], + const double tol = 1e-14 ) + : _half_neigh( half_neigh ) + , _neigh_list( neighbor_list_type( positions, 0, num_local, delta + tol, + 1.0, mesh_min, mesh_max ) ) + { + } + + // Constructor which stores existing neighbors. template Force( const bool half_neigh, const NeighborListType& neighbors ) : _half_neigh( half_neigh ) @@ -224,7 +237,7 @@ class Force ******************************************************************************/ template void computeForce( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag ) + const ParallelType& neigh_op_tag, const bool reset = true ) { auto n_local = particles.n_local; auto x = particles.sliceReferencePosition(); @@ -233,7 +246,8 @@ void computeForce( ForceType& force, ParticleType& particles, auto f_a = particles.sliceForceAtomic(); // Reset force. - Cabana::deep_copy( f, 0.0 ); + if ( reset ) + Cabana::deep_copy( f, 0.0 ); // if ( half_neigh ) // Forces must be atomic for half list @@ -278,7 +292,7 @@ double computeEnergy( ForceType& force, ParticleType& particles, template void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu, - const ParallelType& neigh_op_tag ) + const ParallelType& neigh_op_tag, const bool reset = true ) { auto n_local = particles.n_local; auto x = particles.sliceReferencePosition(); @@ -287,7 +301,8 @@ void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu, auto f_a = particles.sliceForceAtomic(); // Reset force. - Cabana::deep_copy( f, 0.0 ); + if ( reset ) + Cabana::deep_copy( f, 0.0 ); // if ( half_neigh ) // Forces must be atomic for half list diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index b615b15e..29a0f898 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -32,22 +32,7 @@ class Inputs inputs = parse( filename ); // Add additional derived inputs to json. System size. - auto size = inputs["system_size"]["value"]; - std::string size_unit = inputs["system_size"]["unit"]; - for ( std::size_t d = 0; d < size.size(); d++ ) - { - double s = size[d]; - double low = -0.5 * s; - double high = 0.5 * s; - inputs["low_corner"]["value"][d] = low; - inputs["high_corner"]["value"][d] = high; - - double nc = inputs["num_cells"]["value"][d]; - inputs["dx"]["value"][d] = ( high - low ) / nc; - } - inputs["low_corner"]["unit"] = size_unit; - inputs["high_corner"]["unit"] = size_unit; - inputs["dx"]["unit"] = size_unit; + setupSize(); // Number of steps. double tf = inputs["final_time"]["value"]; @@ -104,7 +89,73 @@ class Inputs // Not yet a user option. inputs["half_neigh"]["value"] = false; } - ~Inputs() {} + + void setupSize() + { + if ( inputs.contains( "system_size" ) ) + { + auto size = inputs["system_size"]["value"]; + std::string size_unit = inputs["system_size"]["unit"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double s = size[d]; + double low = -0.5 * s; + double high = 0.5 * s; + inputs["low_corner"]["value"][d] = low; + inputs["high_corner"]["value"][d] = high; + } + inputs["low_corner"]["unit"] = size_unit; + inputs["high_corner"]["unit"] = size_unit; + } + else if ( inputs.contains( "low_corner" ) && + ( inputs.contains( "high_corner" ) ) ) + { + auto low_corner = inputs["low_corner"]["value"]; + auto high_corner = inputs["high_corner"]["value"]; + std::string size_unit = inputs["low_corner"]["unit"]; + for ( std::size_t d = 0; d < low_corner.size(); d++ ) + { + double low = low_corner[d]; + double high = high_corner[d]; + inputs["system_size"]["value"][d] = high - low; + } + inputs["system_size"]["unit"] = size_unit; + } + else + { + throw std::runtime_error( "Must input either system_size or " + "both low_corner and high_corner." ); + } + + if ( inputs.contains( "dx" ) ) + { + auto size = inputs["system_size"]["value"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double system_size = inputs["system_size"]["value"][d]; + double dx = inputs["dx"]["value"][d]; + inputs["num_cells"]["value"][d] = + static_cast( system_size / dx ); + } + } + else if ( inputs.contains( "num_cells" ) ) + { + auto size = inputs["system_size"]["value"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double low = inputs["low_corner"]["value"][d]; + double high = inputs["low_corner"]["value"][d]; + double nc = inputs["num_cells"]["value"][d]; + inputs["dx"]["value"][d] = ( high - low ) / nc; + } + std::string size_unit = inputs["system_size"]["unit"]; + inputs["dx"]["unit"] = size_unit; + } + else + { + throw std::runtime_error( "Must input either num_cells or dx." ); + } + } void computeCriticalTimeStep() { diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index d7a99321..d8be3db4 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -92,7 +92,7 @@ class SolverBase }; template + class ForceModelType, class ContactModelType = NoContact> class SolverElastic { public: @@ -102,7 +102,7 @@ class SolverElastic // Core module types - required for all problems. using particle_type = ParticleType; using integrator_type = Integrator; - using force_model_type = ForceModel; + using force_model_type = ForceModelType; using force_type = Force; using comm_type = Comm; @@ -111,6 +111,8 @@ class SolverElastic // Optional module types. using heat_transfer_type = HeatTransfer; + using contact_type = Force; + using contact_model_type = ContactModelType; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -118,6 +120,27 @@ class SolverElastic : inputs( _inputs ) , particles( _particles ) , _init_time( 0.0 ) + { + setup( force_model ); + } + + SolverElastic( input_type _inputs, + std::shared_ptr _particles, + force_model_type force_model, + contact_model_type contact_model ) + : inputs( _inputs ) + , particles( _particles ) + , _init_time( 0.0 ) + { + setup( force_model ); + + _neighbor_timer.start(); + contact = std::make_shared( inputs["half_neigh"], + *particles, contact_model ); + _neighbor_timer.stop(); + } + + void setup( force_model_type force_model ) { num_steps = inputs["num_steps"]; output_frequency = inputs["output_frequency"]; @@ -130,12 +153,29 @@ class SolverElastic // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); + if constexpr ( is_contact::value ) + { + if ( comm->size() > 1 ) + throw std::runtime_error( + "Contact with MPI is currently disabled." ); + } + // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) force_model.update( particles->sliceTemperature() ); + // Create heat transfer if needed. + if constexpr ( is_heat_transfer< + typename force_model_type::thermal_type>::value ) + { + thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; + heat_transfer = std::make_shared( + inputs["half_neigh"], force->_neigh_list, force_model ); + } + _neighbor_timer.start(); + // This will either be PD or DEM forces. force = std::make_shared( inputs["half_neigh"], *particles, force_model ); _neighbor_timer.stop(); @@ -255,6 +295,9 @@ class SolverElastic // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeForce( *contact, *particles, neigh_iter_tag{}, false ); + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, step * dt ); @@ -287,6 +330,9 @@ class SolverElastic // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeForce( *contact, *particles, neigh_iter_tag{}, false ); + if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); @@ -422,6 +468,8 @@ class SolverElastic std::shared_ptr force; // Optional modules. std::shared_ptr heat_transfer; + std::shared_ptr contact; + contact_model_type contact_model; // Output files. std::string output_file; @@ -437,41 +485,40 @@ class SolverElastic }; template + class ForceModelType, class ContactModelType = NoContact> class SolverFracture - : public SolverElastic + : public SolverElastic { public: - using base_type = - SolverElastic; + using base_type = SolverElastic; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; using particle_type = typename base_type::particle_type; using integrator_type = typename base_type::integrator_type; using comm_type = typename base_type::comm_type; - using force_model_type = ForceModel; + using force_model_type = ForceModelType; using force_type = typename base_type::force_type; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = typename base_type::input_type; - template + using contact_model_type = ContactModelType; + SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, PrenotchType prenotch ) + force_model_type force_model ) : base_type( _inputs, _particles, force_model ) { init_mu(); - - // Create prenotch. - prenotch.create( exec_space{}, mu, *particles, force->_neigh_list ); - _init_time += prenotch.time(); } SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model ) - : base_type( _inputs, _particles, force_model ) + force_model_type force_model, + contact_model_type contact_model ) + : base_type( _inputs, _particles, force_model, contact_model ) { init_mu(); } @@ -488,6 +535,14 @@ class SolverFracture _init_timer.stop(); } + template + void init_prenotch( Prenotch prenotch ) + { + // Create prenotch. + prenotch.create( exec_space{}, mu, *particles, force->_neigh_list ); + _init_time += prenotch.time(); + } + void init( const bool initial_output = true ) { // Compute initial internal forces and energy. @@ -498,6 +553,14 @@ class SolverFracture particles->output( 0, 0.0, output_reference ); } + template + void init( Prenotch prenotch, + const bool initial_output = true ) + { + init_prenotch( prenotch ); + init( initial_output ); + } + template void init( BoundaryType boundary_condition, const bool initial_output = true ) @@ -522,6 +585,14 @@ class SolverFracture particles->output( 0, 0.0, output_reference ); } + template + void init( BoundaryType boundary_condition, Prenotch prenotch, + const bool initial_output = true ) + { + init_prenotch( prenotch ); + init( boundary_condition, initial_output ); + } + template void run( BoundaryType boundary_condition ) { @@ -549,6 +620,9 @@ class SolverFracture // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeForce( *contact, *particles, neigh_iter_tag{}, false ); + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space{}, *particles, step * dt ); @@ -585,6 +659,9 @@ class SolverFracture // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeForce( *contact, *particles, neigh_iter_tag{}, false ); + // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); @@ -637,6 +714,7 @@ class SolverFracture protected: using base_type::comm; + using base_type::contact; using base_type::force; using base_type::inputs; using base_type::integrator; @@ -651,37 +729,52 @@ class SolverFracture using base_type::print; }; +// =============================================================== + template + class ForceModelType> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model ) + ForceModelType model ) { return std::make_shared< - SolverElastic>( + SolverElastic>( inputs, particles, model ); } template + class ForceModelType, class ContactModelType> +auto createSolverElastic( InputsType inputs, + std::shared_ptr particles, + ForceModelType model, ContactModelType contact_model ) +{ + return std::make_shared>( + inputs, particles, model, contact_model ); +} + +template auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model ) + ForceModelType model ) { return std::make_shared< - SolverFracture>( + SolverFracture>( inputs, particles, model ); } template + class ForceModelType, class ContactModelType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, PrenotchType prenotch ) + ForceModelType model, + ContactModelType contact_model ) { return std::make_shared< - SolverFracture>( - inputs, particles, model, prenotch ); + SolverFracture>( inputs, particles, model, + contact_model ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 58894e93..4ca10d35 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -14,7 +14,7 @@ namespace CabanaPD { -// Mechanics types. +// Fracture tags. struct Elastic { }; @@ -22,7 +22,17 @@ struct Fracture { }; -// Thermal types. +// Contact and DEM (contact without PD) tags. + +struct NoContact +{ +}; +template +struct is_contact : public std::false_type +{ +}; + +// Thermal tags. struct TemperatureIndependent { using base_type = TemperatureIndependent; @@ -58,7 +68,7 @@ struct is_heat_transfer : public std::true_type { }; -// Model types. +// Force model tags. struct PMB { }; diff --git a/src/force/CabanaPD_ContactModels.hpp b/src/force/CabanaPD_ContactModels.hpp new file mode 100644 index 00000000..b588e52d --- /dev/null +++ b/src/force/CabanaPD_ContactModels.hpp @@ -0,0 +1,79 @@ +/**************************************************************************** + * Copyright (c) 2022 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef CONTACTMODELS_H +#define CONTACTMODELS_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Contact model +******************************************************************************/ +struct ContactModel +{ + double delta; + double Rc; + + ContactModel(){}; + // PD horizon + // Contact radius + ContactModel( const double _delta, const double _Rc ) + : delta( _delta ) + , Rc( _Rc ){}; +}; + +/* Normal repulsion */ + +struct NormalRepulsionModel : public ContactModel +{ + // FIXME: This is for use as the primary force model. + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; + + using ContactModel::delta; + using ContactModel::Rc; + + double c; + double K; + + NormalRepulsionModel(){}; + NormalRepulsionModel( const double delta, const double Rc, const double _K ) + : ContactModel( delta, Rc ) + , K( _K ) + { + set_param( delta, Rc, K ); + } + + void set_param( const double _delta, const double _Rc, const double _K ) + { + delta = _delta; + Rc = _Rc; + K = _K; + // This could inherit from PMB (same c) + c = 18.0 * K / ( pi * delta * delta * delta * delta ); + } +}; + +template <> +struct is_contact : public std::true_type +{ +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp new file mode 100644 index 00000000..f3b389ef --- /dev/null +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -0,0 +1,123 @@ +/**************************************************************************** + * Copyright (c) 2022 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef CONTACT_H +#define CONTACT_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Normal repulsion forces +******************************************************************************/ +template +class Force + : public Force +{ + public: + using base_type = Force; + using neighbor_list_type = typename base_type::neighbor_list_type; + + template + Force( const bool half_neigh, const ParticleType& particles, + const NormalRepulsionModel model ) + : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(), + particles.n_local, particles.ghost_mesh_lo, + particles.ghost_mesh_hi ) + , _model( model ) + { + for ( int d = 0; d < particles.dim; d++ ) + { + mesh_min[d] = particles.ghost_mesh_lo[d]; + mesh_max[d] = particles.ghost_mesh_hi[d]; + } + } + + template + void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, + const ParticleType& particles, const int n_local, + ParallelType& neigh_op_tag ) + { + auto delta = _model.delta; + auto Rc = _model.Rc; + auto c = _model.c; + const auto vol = particles.sliceVolume(); + const auto y = particles.sliceCurrentPosition(); + + _neigh_timer.start(); + _neigh_list.build( y, 0, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_timer.stop(); + + auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fcx_i = 0.0; + double fcy_i = 0.0; + double fcz_i = 0.0; + + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Contact "stretch" + const double sc = ( r - Rc ) / delta; + + // Normal repulsion uses a 15 factor compared to the PMB force + const double coeff = 15 * c * sc * vol( j ); + fcx_i = coeff * rx / r; + fcy_i = coeff * ry / r; + fcz_i = coeff * rz / r; + + fc( i, 0 ) += fcx_i; + fc( i, 1 ) += fcy_i; + fc( i, 2 ) += fcz_i; + }; + + _timer.start(); + + // FIXME: using default space for now. + using exec_space = typename MemorySpace::execution_space; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::Contact::compute_full" ); + + _timer.stop(); + } + + // FIXME: implement energy + template + double computeEnergyFull( WType&, const PosType&, const PosType&, + ParticleType&, const int, ParallelType& ) + { + return 0.0; + } + + protected: + NormalRepulsionModel _model; + using base_type::_half_neigh; + using base_type::_neigh_list; + using base_type::_timer; + Timer _neigh_timer; + + double mesh_max[3]; + double mesh_min[3]; +}; + +} // namespace CabanaPD + +#endif