ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3725 from akohlmey/collected-small-changes

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2023-04-10 13:53:16 -04:00
committed by GitHub
parent 96f3c88f0f 415be03f6c
commit ff5b4e971d
21 changed files with 1091 additions and 892 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/LAMMPSUtils.cmake
View File

@ -152,7 +152,7 @@ endfunction(FetchPotentials)
# set CMAKE_LINUX_DISTRO and CMAKE_DISTRO_VERSION on Linux # set CMAKE_LINUX_DISTRO and CMAKE_DISTRO_VERSION on Linux
if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (EXISTS /etc/os-release)) if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (EXISTS /etc/os-release))
file(STRINGS /etc/os-release distro REGEX "^NAME=") file(STRINGS /etc/os-release distro REGEX "^NAME=")
string(REGEX REPLACE "NAME=\"?([^\"]*)\"?" "\\1" distro "${distro}") string(REGEX REPLACE "NAME=\"?([^ ]+).*\"?" "\\1" distro "${distro}")
file(STRINGS /etc/os-release disversion REGEX "^VERSION_ID=") file(STRINGS /etc/os-release disversion REGEX "^VERSION_ID=")
string(REGEX REPLACE "VERSION_ID=\"?([^\"]*)\"?" "\\1" disversion "${disversion}") string(REGEX REPLACE "VERSION_ID=\"?([^\"]*)\"?" "\\1" disversion "${disversion}")
set(CMAKE_LINUX_DISTRO ${distro}) set(CMAKE_LINUX_DISTRO ${distro})

235
cmake/Modules/Packages/PLUMED.cmake
View File

@ -1,106 +1,169 @@
set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library") # Plumed2 support for PLUMED package
set(PLUMED_MODE_VALUES static shared runtime)
set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES})
validate_option(PLUMED_MODE PLUMED_MODE_VALUES)
string(TOUPPER ${PLUMED_MODE} PLUMED_MODE)
set(PLUMED_LINK_LIBS) if(BUILD_MPI)
if(PLUMED_MODE STREQUAL "STATIC") set(PLUMED_CONFIG_MPI "--enable-mpi")
find_package(LAPACK REQUIRED) set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER})
find_package(BLAS REQUIRED) set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER})
find_package(GSL REQUIRED) set(PLUMED_CONFIG_CPP "-I ${MPI_CXX_INCLUDE_PATH}")
list(APPEND PLUMED_LINK_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} GSL::gsl) set(PLUMED_CONFIG_LIB "${MPI_CXX_LIBRARIES}")
find_package(ZLIB QUIET) set(PLUMED_CONFIG_DEP "mpi4win_build")
if(ZLIB_FOUND) else()
list(APPEND PLUMED_LINK_LIBS ZLIB::ZLIB) set(PLUMED_CONFIG_MPI "--disable-mpi")
endif() set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER})
find_package(FFTW3 QUIET) set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER})
if(FFTW3_FOUND) set(PLUMED_CONFIG_CPP "")
list(APPEND PLUMED_LINK_LIBS FFTW3::FFTW3) set(PLUMED_CONFIG_LIB "")
endif() set(PLUMED_CONFIG_DEP "")
endif()
if(BUILD_OMP)
set(PLUMED_CONFIG_OMP "--enable-openmp")
else()
set(PLUMED_CONFIG_OMP "--disable-openmp")
endif() endif()
find_package(PkgConfig QUIET) set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz"
set(DOWNLOAD_PLUMED_DEFAULT ON) CACHE STRING "URL for PLUMED tarball")
if(PKG_CONFIG_FOUND) set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball")
pkg_check_modules(PLUMED QUIET plumed)
if(PLUMED_FOUND)
set(DOWNLOAD_PLUMED_DEFAULT OFF)
endif()
endif()
option(DOWNLOAD_PLUMED "Download Plumed package instead of using an already installed one" ${DOWNLOAD_PLUMED_DEFAULT}) mark_as_advanced(PLUMED_URL)
if(DOWNLOAD_PLUMED) mark_as_advanced(PLUMED_MD5)
if(BUILD_MPI) GetFallbackURL(PLUMED_URL PLUMED_FALLBACK)
set(PLUMED_CONFIG_MPI "--enable-mpi")
set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER}) if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING))
set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER}) if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64")
set(CROSS_CONFIGURE mingw64-configure)
elseif(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86")
set(CROSS_CONFIGURE mingw32-configure)
else() else()
set(PLUMED_CONFIG_MPI "--disable-mpi") message(FATAL_ERROR "Unsupported target system: ${CMAKE_SYSTEM_NAME}/${CMAKE_SYSTEM_PROCESSOR}")
set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER})
set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER})
endif() endif()
if(BUILD_OMP) message(STATUS "Downloading and cross-compiling Plumed2 for ${CMAKE_SYSTEM_NAME}/${CMAKE_SYSTEM_PROCESSOR} with ${CROSS_CONFIGURE}")
set(PLUMED_CONFIG_OMP "--enable-openmp")
else()
set(PLUMED_CONFIG_OMP "--disable-openmp")
endif()
message(STATUS "PLUMED download requested - we will build our own")
if(PLUMED_MODE STREQUAL "STATIC")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}")
elseif(PLUMED_MODE STREQUAL "SHARED")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX};<INSTALL_DIR>/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
elseif(PLUMED_MODE STREQUAL "RUNTIME")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_PREFIX}")
endif()
set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball")
set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball")
mark_as_advanced(PLUMED_URL)
mark_as_advanced(PLUMED_MD5)
GetFallbackURL(PLUMED_URL PLUMED_FALLBACK)
include(ExternalProject) include(ExternalProject)
ExternalProject_Add(plumed_build ExternalProject_Add(plumed_build
URL ${PLUMED_URL} ${PLUMED_FALLBACK} URL ${PLUMED_URL} ${PLUMED_FALLBACK}
URL_MD5 ${PLUMED_MD5} URL_MD5 ${PLUMED_MD5}
BUILD_IN_SOURCE 1 BUILD_IN_SOURCE 1
CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR> CONFIGURE_COMMAND ${CROSS_CONFIGURE} --disable-shared --disable-bsymbolic
--disable-python --enable-cxx=11
--enable-modules=-adjmat:+crystallization:-dimred:+drr:+eds:-fisst:+funnel:+logmfd:+manyrestraints:+maze:+opes:+multicolvar:-pamm:-piv:+s2cm:-sasa:-ves
${PLUMED_CONFIG_OMP}
${PLUMED_CONFIG_MPI}
CXX=${PLUMED_CONFIG_CXX}
CC=${PLUMED_CONFIG_CC}
CPPFLAGS=${PLUMED_CONFIG_CPP}
LIBS=${PLUMED_CONFIG_LIB}
INSTALL_COMMAND ""
BUILD_BYPRODUCTS "<SOURCE_DIR>/src/lib/install/libplumed.a" "<SOURCE_DIR>/src/lib/install/plumed.exe"
DEPENDS "${PLUMED_MPI_CONFIG_DEP}"
)
ExternalProject_Get_Property(plumed_build SOURCE_DIR)
set(PLUMED_BUILD_DIR ${SOURCE_DIR})
set(PLUMED_INSTALL_DIR ${PLUMED_BUILD_DIR}/src/lib/install)
file(MAKE_DIRECTORY ${PLUMED_BUILD_DIR}/src/include)
add_library(LAMMPS::PLUMED UNKNOWN IMPORTED)
add_dependencies(LAMMPS::PLUMED plumed_build)
set_target_properties(LAMMPS::PLUMED PROPERTIES
IMPORTED_LOCATION "${PLUMED_INSTALL_DIR}/libplumed.a"
INTERFACE_LINK_LIBRARIES "-Wl,--image-base -Wl,0x10000000 -lfftw3 -lz -fstack-protector -lssp -fopenmp"
INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_BUILD_DIR}/src/include")
add_custom_target(plumed_copy ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches
COMMAND ${CMAKE_COMMAND} -E copy ${PLUMED_INSTALL_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed.exe
COMMAND ${CMAKE_COMMAND} -E copy_directory ${PLUMED_BUILD_DIR}/patches ${CMAKE_BINARY_DIR}/patches
BYPRODUCTS ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/patches
DEPENDS plumed_build
COMMENT "Copying Plumed files"
)
else()
set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library")
set(PLUMED_MODE_VALUES static shared runtime)
set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES})
validate_option(PLUMED_MODE PLUMED_MODE_VALUES)
string(TOUPPER ${PLUMED_MODE} PLUMED_MODE)
set(PLUMED_LINK_LIBS)
if(PLUMED_MODE STREQUAL "STATIC")
find_package(LAPACK REQUIRED)
find_package(BLAS REQUIRED)
find_package(GSL REQUIRED)
list(APPEND PLUMED_LINK_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} GSL::gsl)
find_package(ZLIB QUIET)
if(ZLIB_FOUND)
list(APPEND PLUMED_LINK_LIBS ZLIB::ZLIB)
endif()
find_package(FFTW3 QUIET)
if(FFTW3_FOUND)
list(APPEND PLUMED_LINK_LIBS FFTW3::FFTW3)
endif()
endif()
find_package(PkgConfig QUIET)
set(DOWNLOAD_PLUMED_DEFAULT ON)
if(PKG_CONFIG_FOUND)
pkg_check_modules(PLUMED QUIET plumed)
if(PLUMED_FOUND)
set(DOWNLOAD_PLUMED_DEFAULT OFF)
endif()
endif()
option(DOWNLOAD_PLUMED "Download Plumed package instead of using an already installed one" ${DOWNLOAD_PLUMED_DEFAULT})
if(DOWNLOAD_PLUMED)
message(STATUS "PLUMED download requested - we will build our own")
if(PLUMED_MODE STREQUAL "STATIC")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}")
elseif(PLUMED_MODE STREQUAL "SHARED")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX};<INSTALL_DIR>/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
elseif(PLUMED_MODE STREQUAL "RUNTIME")
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_PREFIX}")
endif()
include(ExternalProject)
ExternalProject_Add(plumed_build
URL ${PLUMED_URL} ${PLUMED_FALLBACK}
URL_MD5 ${PLUMED_MD5}
BUILD_IN_SOURCE 1
CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR>
${CONFIGURE_REQUEST_PIC} ${CONFIGURE_REQUEST_PIC}
--enable-modules=all --enable-modules=all
--enable-cxx=11
--disable-python
${PLUMED_CONFIG_MPI} ${PLUMED_CONFIG_MPI}
${PLUMED_CONFIG_OMP} ${PLUMED_CONFIG_OMP}
CXX=${PLUMED_CONFIG_CXX} CXX=${PLUMED_CONFIG_CXX}
CC=${PLUMED_CONFIG_CC} CC=${PLUMED_CONFIG_CC}
BUILD_BYPRODUCTS ${PLUMED_BUILD_BYPRODUCTS} BUILD_BYPRODUCTS ${PLUMED_BUILD_BYPRODUCTS}
) )
ExternalProject_get_property(plumed_build INSTALL_DIR) ExternalProject_get_property(plumed_build INSTALL_DIR)
add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) add_library(LAMMPS::PLUMED UNKNOWN IMPORTED)
add_dependencies(LAMMPS::PLUMED plumed_build) add_dependencies(LAMMPS::PLUMED plumed_build)
if(PLUMED_MODE STREQUAL "STATIC") if(PLUMED_MODE STREQUAL "STATIC")
set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${PLUMED_LINK_LIBS};${CMAKE_DL_LIBS}") set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${PLUMED_LINK_LIBS};${CMAKE_DL_LIBS}")
elseif(PLUMED_MODE STREQUAL "SHARED") elseif(PLUMED_MODE STREQUAL "SHARED")
set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX};${CMAKE_DL_LIBS}") set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX};${CMAKE_DL_LIBS}")
elseif(PLUMED_MODE STREQUAL "RUNTIME") elseif(PLUMED_MODE STREQUAL "RUNTIME")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${CMAKE_DL_LIBS}") set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${CMAKE_DL_LIBS}")
endif()
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include)
file(MAKE_DIRECTORY ${INSTALL_DIR}/include)
else()
find_package(PkgConfig REQUIRED)
pkg_check_modules(PLUMED REQUIRED plumed)
add_library(LAMMPS::PLUMED INTERFACE IMPORTED)
if(PLUMED_MODE STREQUAL "STATIC")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static)
elseif(PLUMED_MODE STREQUAL "SHARED")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared)
elseif(PLUMED_MODE STREQUAL "RUNTIME")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime)
endif()
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}")
endif() endif()
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include)
file(MAKE_DIRECTORY ${INSTALL_DIR}/include)
else()
find_package(PkgConfig REQUIRED)
pkg_check_modules(PLUMED REQUIRED plumed)
add_library(LAMMPS::PLUMED INTERFACE IMPORTED)
if(PLUMED_MODE STREQUAL "STATIC")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static)
elseif(PLUMED_MODE STREQUAL "SHARED")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared)
elseif(PLUMED_MODE STREQUAL "RUNTIME")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime)
endif()
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}")
endif() endif()
target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) target_link_libraries(lammps PRIVATE LAMMPS::PLUMED)

18
cmake/Modules/Testing.cmake
View File

@ -18,29 +18,33 @@ if(ENABLE_TESTING)
# we need to build and link a LOT of tester executables, so it is worth checking if # we need to build and link a LOT of tester executables, so it is worth checking if
# a faster linker is available. requires GNU or Clang compiler, newer CMake. # a faster linker is available. requires GNU or Clang compiler, newer CMake.
# also only verified with Fedora Linux > 30 and Ubuntu <= 18.04 (Ubuntu 20.04 fails) # also only verified with Fedora Linux > 30 and Ubuntu 18.04 or 22.04+(Ubuntu 20.04 fails)
if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (CMAKE_VERSION VERSION_GREATER_EQUAL 3.13) if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (CMAKE_VERSION VERSION_GREATER_EQUAL 3.13)
AND ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") AND ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") OR (CMAKE_CXX_COMPILER_ID STREQUAL "Clang")))
OR (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))) if(((CMAKE_LINUX_DISTRO STREQUAL "Ubuntu") AND
if(((CMAKE_LINUX_DISTRO STREQUAL "Ubuntu") AND (CMAKE_DISTRO_VERSION VERSION_LESS_EQUAL 18.04)) ((CMAKE_DISTRO_VERSION VERSION_LESS_EQUAL 18.04) OR (CMAKE_DISTRO_VERSION VERSION_GREATER_EQUAL 22.04)))
OR ((CMAKE_LINUX_DISTRO STREQUAL "Fedora") AND (CMAKE_DISTRO_VERSION VERSION_GREATER 30))) OR ((CMAKE_LINUX_DISTRO STREQUAL "Fedora") AND (CMAKE_DISTRO_VERSION VERSION_GREATER 30)))
include(CheckCXXCompilerFlag) include(CheckCXXCompilerFlag)
set(CMAKE_CUSTOM_LINKER_DEFAULT default) set(CMAKE_CUSTOM_LINKER_DEFAULT default)
check_cxx_compiler_flag(-fuse-ld=mold HAVE_MOLD_LINKER_FLAG)
check_cxx_compiler_flag(-fuse-ld=lld HAVE_LLD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=lld HAVE_LLD_LINKER_FLAG)
check_cxx_compiler_flag(-fuse-ld=gold HAVE_GOLD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=gold HAVE_GOLD_LINKER_FLAG)
check_cxx_compiler_flag(-fuse-ld=bfd HAVE_BFD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=bfd HAVE_BFD_LINKER_FLAG)
find_program(HAVE_MOLD_LINKER_BIN ld.mold)
find_program(HAVE_LLD_LINKER_BIN lld ld.lld) find_program(HAVE_LLD_LINKER_BIN lld ld.lld)
find_program(HAVE_GOLD_LINKER_BIN ld.gold) find_program(HAVE_GOLD_LINKER_BIN ld.gold)
find_program(HAVE_BFD_LINKER_BIN ld.bfd) find_program(HAVE_BFD_LINKER_BIN ld.bfd)
if(HAVE_LLD_LINKER_FLAG AND HAVE_LLD_LINKER_BIN) if(HAVE_MOLD_LINKER_FLAG AND HAVE_MOLD_LINKER_BIN)
set(CMAKE_CUSTOM_LINKER_DEFAULT mold)
elseif(HAVE_LLD_LINKER_FLAG AND HAVE_LLD_LINKER_BIN)
set(CMAKE_CUSTOM_LINKER_DEFAULT lld) set(CMAKE_CUSTOM_LINKER_DEFAULT lld)
elseif(HAVE_GOLD_LINKER_FLAG AND HAVE_GOLD_LINKER_BIN) elseif(HAVE_GOLD_LINKER_FLAG AND HAVE_GOLD_LINKER_BIN)
set(CMAKE_CUSTOM_LINKER_DEFAULT gold) set(CMAKE_CUSTOM_LINKER_DEFAULT gold)
elseif(HAVE_BFD_LINKER_FLAG AND HAVE_BFD_LINKER_BIN) elseif(HAVE_BFD_LINKER_FLAG AND HAVE_BFD_LINKER_BIN)
set(CMAKE_CUSTOM_LINKER_DEFAULT bfd) set(CMAKE_CUSTOM_LINKER_DEFAULT bfd)
endif() endif()
set(CMAKE_CUSTOM_LINKER_VALUES lld gold bfd default) set(CMAKE_CUSTOM_LINKER_VALUES mold lld gold bfd default)
set(CMAKE_CUSTOM_LINKER ${CMAKE_CUSTOM_LINKER_DEFAULT} CACHE STRING "Choose a custom linker for faster linking (lld, gold, bfd, default)") set(CMAKE_CUSTOM_LINKER ${CMAKE_CUSTOM_LINKER_DEFAULT} CACHE STRING "Choose a custom linker for faster linking (mold, lld, gold, bfd, default)")
validate_option(CMAKE_CUSTOM_LINKER CMAKE_CUSTOM_LINKER_VALUES) validate_option(CMAKE_CUSTOM_LINKER CMAKE_CUSTOM_LINKER_VALUES)
mark_as_advanced(CMAKE_CUSTOM_LINKER) mark_as_advanced(CMAKE_CUSTOM_LINKER)
if(NOT "${CMAKE_CUSTOM_LINKER}" STREQUAL "default") if(NOT "${CMAKE_CUSTOM_LINKER}" STREQUAL "default")

11
doc/src/Build_extras.rst
View File

@ -683,20 +683,11 @@ This list was last updated for version 3.7.1 of the Kokkos library.
-D Kokkos_ARCH_GPUARCH=yes # GPUARCH = GPU from list above -D Kokkos_ARCH_GPUARCH=yes # GPUARCH = GPU from list above
-D Kokkos_ENABLE_CUDA=yes -D Kokkos_ENABLE_CUDA=yes
-D Kokkos_ENABLE_OPENMP=yes -D Kokkos_ENABLE_OPENMP=yes
-D CMAKE_CXX_COMPILER=wrapper # wrapper = full path to Cuda nvcc wrapper
This will also enable executing FFTs on the GPU, either via the This will also enable executing FFTs on the GPU, either via the
internal KISSFFT library, or - by preference - with the cuFFT internal KISSFFT library, or - by preference - with the cuFFT
library bundled with the CUDA toolkit, depending on whether CMake library bundled with the CUDA toolkit, depending on whether CMake
can identify its location. The *wrapper* value for can identify its location.
``CMAKE_CXX_COMPILER`` variable is the path to the CUDA nvcc
compiler wrapper provided in the Kokkos library:
``lib/kokkos/bin/nvcc_wrapper``\ . The setting should include the
full path name to the wrapper, e.g.
.. code-block:: bash
-D CMAKE_CXX_COMPILER=${HOME}/lammps/lib/kokkos/bin/nvcc_wrapper
For AMD or NVIDIA GPUs using HIP, set these variables: For AMD or NVIDIA GPUs using HIP, set these variables:

6
doc/src/fix_wall_gran.rst
View File

@ -71,9 +71,9 @@ Examples
fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0
fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL
fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0
fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox
fix 4 all wall/gran/region granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone fix 4 all wall/gran granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone
fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0 fix 5 all wall/gran granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0
fix 6 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts fix 6 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts
Description Description

2
doc/src/pair_colloid.rst
View File

@ -93,7 +93,7 @@ with :math:`A_{ss}` set appropriately, which results from letting both
particle sizes go to zero. particle sizes go to zero.
When used in combination with :doc:`pair_style yukawa/colloid When used in combination with :doc:`pair_style yukawa/colloid
<pair_colloid>`, the two terms become the so-called DLVO potential, <pair_yukawa_colloid>`, the two terms become the so-called DLVO potential,
which combines electrostatic repulsion and van der Waals attraction. which combines electrostatic repulsion and van der Waals attraction.
The following coefficients must be defined for each pair of atoms The following coefficients must be defined for each pair of atoms

6
examples/mliap/mliap_numpy_Ta06A.py
View File

@ -85,8 +85,7 @@ class LinearModel():
def __call__(self,elems,bispectrum,beta,energy): def __call__(self,elems,bispectrum,beta,energy):
energy[:] = bispectrum @ self.weights + self.bias energy[:] = bispectrum @ self.weights + self.bias
beta[:] = self.weights beta[:] = self.weights
mymodel = LinearModel("Ta06A.mliap.model") mymodel = LinearModel("Ta06A.mliap.model")
import lammps import lammps
@ -98,4 +97,5 @@ lammps.mliap.activate_mliappy(lmp)
lmp.commands_string(before_loading) lmp.commands_string(before_loading)
lammps.mliap.load_model(mymodel) lammps.mliap.load_model(mymodel)
lmp.commands_string(after_loading) lmp.commands_string(after_loading)
lmp.close()
lmp.finalize()

4
examples/mliap/mliap_pytorch_Ta06A.py
View File

@ -103,6 +103,8 @@ model = torch.load(torch_model)
# Connect the PyTorch model to the mliap pair style. # Connect the PyTorch model to the mliap pair style.
lammps.mliap.load_model(model) lammps.mliap.load_model(model)
# run the simulation with the mliap pair style # run the simulation with the mliap pair style
lmp.commands_string(after_loading) lmp.commands_string(after_loading)
lmp.close()
lmp.finalize()

4
examples/mliap/mliap_pytorch_Ta06A_kokkos.py
View File

@ -103,6 +103,8 @@ model = torch.load(torch_model)
# Connect the PyTorch model to the mliap pair style. # Connect the PyTorch model to the mliap pair style.
lammps.mliap.load_model_kokkos(model) lammps.mliap.load_model_kokkos(model)
# run the simulation with the mliap pair style # run the simulation with the mliap pair style
lmp.commands_string(after_loading) lmp.commands_string(after_loading)
lmp.close()
lmp.finalize()

32
examples/mliap/mliap_unified_lj_Ar.py
View File

@ -6,30 +6,30 @@
before_loading =\ before_loading =\
"""# 3d Lennard-Jones melt """# 3d Lennard-Jones melt
units lj units lj
atom_style atomic atom_style atomic
lattice fcc 0.8442 lattice fcc 0.8442
region box block 0 10 0 10 0 10 region box block 0 10 0 10 0 10
create_box 1 box create_box 1 box
create_atoms 1 box create_atoms 1 box
mass 1 1.0 mass 1 1.0
velocity all create 3.0 87287 loop geom velocity all create 3.0 87287 loop geom
""" """
after_loading =\ after_loading =\
""" """
pair_style mliap unified EXISTS pair_style mliap unified EXISTS
pair_coeff * * Ar pair_coeff * * Ar
neighbor 0.3 bin neighbor 0.3 bin
neigh_modify every 20 delay 0 check no neigh_modify every 20 delay 0 check no
fix 1 all nve fix 1 all nve
thermo 50 thermo 50
run 250 run 250
""" """
import lammps import lammps
@ -63,3 +63,5 @@ lammps.mliap.load_unified(unified)
# Run the simulation with the mliap unified pair style # Run the simulation with the mliap unified pair style
# Use pre-loaded model by specifying model filename as "EXISTS" # Use pre-loaded model by specifying model filename as "EXISTS"
lmp.commands_string(after_loading) lmp.commands_string(after_loading)
lmp.close()
lmp.finalize()

2
src/EXTRA-PAIR/pair_born_gauss.cpp
View File

@ -38,6 +38,8 @@ PairBornGauss::PairBornGauss(LAMMPS *lmp) :
single_enable = 1; single_enable = 1;
respa_enable = 0; respa_enable = 0;
writedata = 1; writedata = 1;
cut_global = 0.0;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

3
src/MANYBODY/pair_bop.cpp
View File

@ -1873,8 +1873,7 @@ void PairBOP::read_table(char *filename)
reader = new PotentialFileReader(lmp, filename, "BOP"); reader = new PotentialFileReader(lmp, filename, "BOP");
bop_types = reader->next_int(); bop_types = reader->next_int();
if (bop_types <= 0) if (bop_types <= 0)
error->one(FLERR,fmt::format("BOP potential file with {} " error->one(FLERR,fmt::format("BOP potential file with {} elements",bop_types));
"elements",bop_types));
bop_elements = new char*[bop_types]; bop_elements = new char*[bop_types];
bop_masses = new double[bop_types]; bop_masses = new double[bop_types];

7
src/MC/fix_charge_regulation.cpp
View File

@ -160,7 +160,12 @@ FixChargeRegulation::~FixChargeRegulation() {
if (exclusion_group_bit && group) { if (exclusion_group_bit && group) {
auto group_id = std::string("FixChargeRegulation:gcmc_exclusion_group:") + id; auto group_id = std::string("FixChargeRegulation:gcmc_exclusion_group:") + id;
group->assign(group_id + " delete"); try {
group->assign(group_id + " delete");
} catch (std::exception &e) {
if (comm->me == 0)
fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what());
}
} }
if (group) { if (group) {

14
src/MC/fix_gcmc.cpp
View File

@ -427,12 +427,22 @@ FixGCMC::~FixGCMC()
if (exclusion_group_bit && group) { if (exclusion_group_bit && group) {
auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id; auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id;
group->assign(group_id + " delete"); try {
group->assign(group_id + " delete");
} catch (std::exception &e) {
if (comm->me == 0)
fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what());
}
} }
if (molecule_group_bit && group) { if (molecule_group_bit && group) {
auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id; auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id;
group->assign(group_id + " delete"); try {
group->assign(group_id + " delete");
} catch (std::exception &e) {
if (comm->me == 0)
fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what());
}
} }
if (full_flag && group && neighbor) { if (full_flag && group && neighbor) {

14
src/MC/fix_widom.cpp
View File

@ -255,12 +255,22 @@ FixWidom::~FixWidom()
if (exclusion_group_bit && group) { if (exclusion_group_bit && group) {
auto group_id = std::string("FixWidom:widom_exclusion_group:") + id; auto group_id = std::string("FixWidom:widom_exclusion_group:") + id;
group->assign(group_id + " delete"); try {
group->assign(group_id + " delete");
} catch (std::exception &e) {
if (comm->me == 0)
fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what());
}
} }
if (molecule_group_bit && group) { if (molecule_group_bit && group) {
auto group_id = std::string("FixWidom:rotation_gas_atoms:") + id; auto group_id = std::string("FixWidom:rotation_gas_atoms:") + id;
group->assign(group_id + " delete"); try {
group->assign(group_id + " delete");
} catch (std::exception &e) {
if (comm->me == 0)
fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what());
}
} }
if (full_flag && group && neighbor) { if (full_flag && group && neighbor) {

40
src/RIGID/fix_shake.cpp
View File

@ -31,15 +31,15 @@
#include "respa.h" #include "respa.h"
#include "update.h" #include "update.h"
#include <cmath>
#include <cctype> #include <cctype>
#include <cmath>
#include <cstring> #include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
#define RVOUS 1 // 0 for irregular, 1 for all2all #define RVOUS 1 // 0 for irregular, 1 for all2all
static constexpr double BIG = 1.0e20; static constexpr double BIG = 1.0e20;
static constexpr double MASSDELTA = 0.1; static constexpr double MASSDELTA = 0.1;
@ -47,15 +47,15 @@ static constexpr double MASSDELTA = 0.1;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), bond_flag(nullptr), angle_flag(nullptr), Fix(lmp, narg, arg), bond_flag(nullptr), angle_flag(nullptr), type_flag(nullptr),
type_flag(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), loop_respa(nullptr),
loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr),
vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), shake_atom(nullptr),
shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), closest_list(nullptr),
b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr), b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr),
b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr), a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr), closest_list(nullptr) a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
{ {
energy_global_flag = energy_peratom_flag = 1; energy_global_flag = energy_peratom_flag = 1;
virial_global_flag = virial_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1;
@ -76,7 +76,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
molecular = atom->molecular; molecular = atom->molecular;
if (molecular == Atom::ATOMIC) if (molecular == Atom::ATOMIC)
error->all(FLERR,"Cannot use fix {} with non-molecular system", style); error->all(FLERR, "Cannot use fix {} with non-molecular system", style);
// perform initial allocation of atom-based arrays // perform initial allocation of atom-based arrays
// register with Atom class // register with Atom class
@ -97,13 +97,13 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 3; comm_forward = 3;
// parse SHAKE args // parse SHAKE args
auto mystyle = fmt::format("fix {}",style); auto mystyle = fmt::format("fix {}", style);
if (narg < 8) utils::missing_cmd_args(FLERR,mystyle, error); if (narg < 8) utils::missing_cmd_args(FLERR, mystyle, error);
tolerance = utils::numeric(FLERR,arg[3],false,lmp); tolerance = utils::numeric(FLERR, arg[3], false, lmp);
max_iter = utils::inumeric(FLERR,arg[4],false,lmp); max_iter = utils::inumeric(FLERR, arg[4], false, lmp);
output_every = utils::inumeric(FLERR,arg[5],false,lmp); output_every = utils::inumeric(FLERR, arg[5], false, lmp);
// parse SHAKE args for bond and angle types // parse SHAKE args for bond and angle types
// will be used by find_clusters // will be used by find_clusters
@ -111,11 +111,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
// store args for "m" in list of length nmass for looping over // store args for "m" in list of length nmass for looping over
// for "m" verify that atom masses have been set // for "m" verify that atom masses have been set
bond_flag = new int[atom->nbondtypes+1]; bond_flag = new int[atom->nbondtypes + 1];
for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0; for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0;
angle_flag = new int[atom->nangletypes+1]; angle_flag = new int[atom->nangletypes + 1];
for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0; for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0;
type_flag = new int[atom->ntypes+1]; type_flag = new int[atom->ntypes + 1];
for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0; for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0;
mass_list = new double[atom->ntypes]; mass_list = new double[atom->ntypes];
nmass = 0; nmass = 0;

362
src/region_block.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -24,100 +23,127 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE}; enum { CONSTANT, VARIABLE };
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), zlostr(nullptr), zhistr(nullptr) Region(lmp, narg, arg), xlostr(nullptr), ylostr(nullptr), zlostr(nullptr), xhistr(nullptr),
yhistr(nullptr), zhistr(nullptr)
{ {
options(narg-8,&arg[8]); xlovar = xhivar = ylovar = yhivar = zlovar = zhivar = -1;
options(narg - 8, &arg[8]);
xlostyle = CONSTANT; xlostyle = CONSTANT;
if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[2],"INF") == 0) xlo = -BIG; if (strcmp(arg[2], "INF") == 0)
else if (domain->triclinic == 0) xlo = domain->boxlo[0]; xlo = -BIG;
else xlo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[2],"^v_")) { xlo = domain->boxlo[0];
xlostr = utils::strdup(arg[2]+2); else
xlo = 0.0; xlo = domain->boxlo_bound[0];
xlostyle = VARIABLE; } else if (utils::strmatch(arg[2], "^v_")) {
varshape = 1; xlostr = utils::strdup(arg[2] + 2);
} else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); xlo = 0.0;
xlostyle = VARIABLE;
varshape = 1;
} else
xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp);
xhistyle = CONSTANT; xhistyle = CONSTANT;
if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[3],"INF") == 0) xhi = BIG; if (strcmp(arg[3], "INF") == 0)
else if (domain->triclinic == 0) xhi = domain->boxhi[0]; xhi = BIG;
else xhi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[3],"^v_")) { xhi = domain->boxhi[0];
xhistr = utils::strdup(arg[3]+2); else
xhi = 0.0; xhi = domain->boxhi_bound[0];
xhistyle = VARIABLE; } else if (utils::strmatch(arg[3], "^v_")) {
varshape = 1; xhistr = utils::strdup(arg[3] + 2);
} else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); xhi = 0.0;
xhistyle = VARIABLE;
varshape = 1;
} else
xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp);
ylostyle = CONSTANT; ylostyle = CONSTANT;
if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[4],"INF") == 0) ylo = -BIG; if (strcmp(arg[4], "INF") == 0)
else if (domain->triclinic == 0) ylo = domain->boxlo[1]; ylo = -BIG;
else ylo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[4],"^v_")) { ylo = domain->boxlo[1];
ylostr = utils::strdup(arg[4]+2); else
ylo = 0.0; ylo = domain->boxlo_bound[1];
ylostyle = VARIABLE; } else if (utils::strmatch(arg[4], "^v_")) {
varshape = 1; ylostr = utils::strdup(arg[4] + 2);
} else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); ylo = 0.0;
ylostyle = VARIABLE;
varshape = 1;
} else
ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp);
yhistyle = CONSTANT; yhistyle = CONSTANT;
if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[5],"INF") == 0) yhi = BIG; if (strcmp(arg[5], "INF") == 0)
else if (domain->triclinic == 0) yhi = domain->boxhi[1]; yhi = BIG;
else yhi = domain->boxhi_bound[1]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[5],"^v_")) { yhi = domain->boxhi[1];
yhistr = utils::strdup(arg[5]+2); else
yhi = 0.0; yhi = domain->boxhi_bound[1];
yhistyle = VARIABLE; } else if (utils::strmatch(arg[5], "^v_")) {
varshape = 1; yhistr = utils::strdup(arg[5] + 2);
} else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); yhi = 0.0;
yhistyle = VARIABLE;
varshape = 1;
} else
yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp);
zlostyle = CONSTANT; zlostyle = CONSTANT;
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[6],"INF") == 0) zlo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) zlo = domain->boxlo[2]; zlo = -BIG;
else zlo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[6],"^v_")) { zlo = domain->boxlo[2];
zlostr = utils::strdup(arg[6]+2); else
zlo = 0.0; zlo = domain->boxlo_bound[2];
zlostyle = VARIABLE; } else if (utils::strmatch(arg[6], "^v_")) {
varshape = 1; zlostr = utils::strdup(arg[6] + 2);
} else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); zlo = 0.0;
zlostyle = VARIABLE;
varshape = 1;
} else
zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
zhistyle = CONSTANT; zhistyle = CONSTANT;
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[7],"INF") == 0) zhi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) zhi = domain->boxhi[2]; zhi = BIG;
else zhi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[7],"^v_")) { zhi = domain->boxhi[2];
zhistr = utils::strdup(arg[7]+2); else
zhi = 0.0; zhi = domain->boxhi_bound[2];
zhistyle = VARIABLE; } else if (utils::strmatch(arg[7], "^v_")) {
varshape = 1; zhistr = utils::strdup(arg[7] + 2);
} else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); zhi = 0.0;
zhistyle = VARIABLE;
varshape = 1;
} else
zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
if (varshape) { if (varshape) {
variable_check(); variable_check();
@ -126,9 +152,9 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
// error check // error check
if (xlo > xhi) error->all(FLERR,"Illegal region block xlo: {} >= xhi: {}", xlo, xhi); if (xlo > xhi) error->all(FLERR, "Illegal region block xlo: {} >= xhi: {}", xlo, xhi);
if (ylo > yhi) error->all(FLERR,"Illegal region block ylo: {} >= yhi: {}", ylo, yhi); if (ylo > yhi) error->all(FLERR, "Illegal region block ylo: {} >= yhi: {}", ylo, yhi);
if (zlo > zhi) error->all(FLERR,"Illegal region block zlo: {} >= zhi: {}", zlo, zhi); if (zlo > zhi) error->all(FLERR, "Illegal region block zlo: {} >= zhi: {}", zlo, zhi);
// extent of block // extent of block
@ -140,15 +166,18 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
extent_yhi = yhi; extent_yhi = yhi;
extent_zlo = zlo; extent_zlo = zlo;
extent_zhi = zhi; extent_zhi = zhi;
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to all 6 planes // particle could be close to all 6 planes
// particle can only touch 3 planes // particle can only touch 3 planes
cmax = 6; cmax = 6;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 3; if (interior)
else tmax = 1; tmax = 3;
else
tmax = 1;
// open face data structs // open face data structs
@ -235,13 +264,13 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
RegBlock::~RegBlock() RegBlock::~RegBlock()
{ {
if (copymode) return; if (copymode) return;
delete [] xlostr; delete[] xlostr;
delete [] xhistr; delete[] xhistr;
delete [] ylostr; delete[] ylostr;
delete [] yhistr; delete[] yhistr;
delete [] zlostr; delete[] zlostr;
delete [] zhistr; delete[] zhistr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -259,8 +288,7 @@ void RegBlock::init()
int RegBlock::inside(double x, double y, double z) int RegBlock::inside(double x, double y, double z)
{ {
if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) return 1;
return 1;
return 0; return 0;
} }
@ -277,8 +305,7 @@ int RegBlock::surface_interior(double *x, double cutoff)
// x is exterior to block // x is exterior to block
if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0;
x[2] < zlo || x[2] > zhi) return 0;
// x is interior to block or on its surface // x is interior to block or on its surface
@ -352,17 +379,16 @@ int RegBlock::surface_interior(double *x, double cutoff)
int RegBlock::surface_exterior(double *x, double cutoff) int RegBlock::surface_exterior(double *x, double cutoff)
{ {
double xp,yp,zp; double xp, yp, zp;
double xc,yc,zc,dist,mindist; double xc, yc, zc, dist, mindist;
// x is far enough from block that there is no contact // x is far enough from block that there is no contact
// x is interior to block // x is interior to block
if (x[0] <= xlo-cutoff || x[0] >= xhi+cutoff || if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff ||
x[1] <= ylo-cutoff || x[1] >= yhi+cutoff || x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff)
x[2] <= zlo-cutoff || x[2] >= zhi+cutoff) return 0; return 0;
if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0;
x[2] > zlo && x[2] < zhi) return 0;
// x is exterior to block or on its surface // x is exterior to block or on its surface
// xp,yp,zp = point on surface of block that x is closest to // xp,yp,zp = point on surface of block that x is closest to
@ -370,20 +396,29 @@ int RegBlock::surface_exterior(double *x, double cutoff)
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
if (!openflag) { if (!openflag) {
if (x[0] < xlo) xp = xlo; if (x[0] < xlo)
else if (x[0] > xhi) xp = xhi; xp = xlo;
else xp = x[0]; else if (x[0] > xhi)
if (x[1] < ylo) yp = ylo; xp = xhi;
else if (x[1] > yhi) yp = yhi; else
else yp = x[1]; xp = x[0];
if (x[2] < zlo) zp = zlo; if (x[1] < ylo)
else if (x[2] > zhi) zp = zhi; yp = ylo;
else zp = x[2]; else if (x[1] > yhi)
yp = yhi;
else
yp = x[1];
if (x[2] < zlo)
zp = zlo;
else if (x[2] > zhi)
zp = zhi;
else
zp = x[2];
} else { } else {
mindist = BIG; mindist = BIG;
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
if (open_faces[i]) continue; if (open_faces[i]) continue;
dist = find_closest_point(i,x,xc,yc,zc); dist = find_closest_point(i, x, xc, yc, zc);
if (dist < mindist) { if (dist < mindist) {
xp = xc; xp = xc;
yp = yc; yp = yc;
@ -393,7 +428,7 @@ int RegBlock::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
return 0; return 0;
@ -403,28 +438,17 @@ int RegBlock::surface_exterior(double *x, double cutoff)
change region shape via variable evaluation change region shape via variable evaluation
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegBlock::shape_update() // addition void RegBlock::shape_update() // addition
{ {
if (xlostyle == VARIABLE) if (xlostyle == VARIABLE) xlo = xscale * input->variable->compute_equal(xlovar);
xlo = xscale * input->variable->compute_equal(xlovar); if (xhistyle == VARIABLE) xhi = xscale * input->variable->compute_equal(xhivar);
if (ylostyle == VARIABLE) ylo = yscale * input->variable->compute_equal(ylovar);
if (xhistyle == VARIABLE) if (yhistyle == VARIABLE) yhi = yscale * input->variable->compute_equal(yhivar);
xhi = xscale * input->variable->compute_equal(xhivar); if (zlostyle == VARIABLE) zlo = zscale * input->variable->compute_equal(zlovar);
if (zhistyle == VARIABLE) zhi = zscale * input->variable->compute_equal(zhivar);
if (ylostyle == VARIABLE)
ylo = yscale * input->variable->compute_equal(ylovar);
if (yhistyle == VARIABLE)
yhi = yscale * input->variable->compute_equal(yhivar);
if (zlostyle == VARIABLE)
zlo = zscale * input->variable->compute_equal(zlovar);
if (zhistyle == VARIABLE)
zhi = zscale * input->variable->compute_equal(zhivar);
if (xlo > xhi || ylo > yhi || zlo > zhi) if (xlo > xhi || ylo > yhi || zlo > zhi)
error->one(FLERR,"Variable evaluation in region gave bad value"); error->one(FLERR, "Variable evaluation in region gave bad value");
// face[0] // face[0]
@ -489,54 +513,45 @@ void RegBlock::shape_update() // addition
error check on existence of variable error check on existence of variable
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegBlock::variable_check() // addition void RegBlock::variable_check() // addition
{ {
if (xlostyle == VARIABLE) { if (xlostyle == VARIABLE) {
xlovar = input->variable->find(xlostr); xlovar = input->variable->find(xlostr);
if (xlovar < 0) if (xlovar < 0) error->all(FLERR, "Variable {} for region block does not exist", xlostr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xlovar)) if (!input->variable->equalstyle(xlovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", xlostr);
} }
if (xhistyle == VARIABLE) { if (xhistyle == VARIABLE) {
xhivar = input->variable->find(xhistr); xhivar = input->variable->find(xhistr);
if (xhivar < 0) if (xhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", xhistr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xhivar)) if (!input->variable->equalstyle(xhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", xhistr);
} }
if (ylostyle == VARIABLE) { if (ylostyle == VARIABLE) {
ylovar = input->variable->find(ylostr); ylovar = input->variable->find(ylostr);
if (ylovar < 0) if (ylovar < 0) error->all(FLERR, "Variable {} for region block does not exist", ylostr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(ylovar)) if (!input->variable->equalstyle(ylovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", ylostr);
} }
if (yhistyle == VARIABLE) { if (yhistyle == VARIABLE) {
yhivar = input->variable->find(yhistr); yhivar = input->variable->find(yhistr);
if (yhivar < 0) if (yhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", yhistr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(yhivar)) if (!input->variable->equalstyle(yhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", yhistr);
} }
if (zlostyle == VARIABLE) { if (zlostyle == VARIABLE) {
zlovar = input->variable->find(zlostr); zlovar = input->variable->find(zlostr);
if (zlovar < 0) if (zlovar < 0) error->all(FLERR, "Variable {} for region block does not exist", zlostr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zlovar)) if (!input->variable->equalstyle(zlovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", zlostr);
} }
if (zhistyle == VARIABLE) { if (zhistyle == VARIABLE) {
zhivar = input->variable->find(zhistr); zhivar = input->variable->find(zhistr);
if (zhivar < 0) if (zhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", zhistr);
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zhivar)) if (!input->variable->equalstyle(zhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable {} for region block is invalid style", zhistr);
} }
} }
@ -545,36 +560,35 @@ void RegBlock::variable_check() // addition
store closest point in xc,yc,zc store closest point in xc,yc,zc
--------------------------------------------------------------------------*/ --------------------------------------------------------------------------*/
double RegBlock::find_closest_point(int i, double *x, double RegBlock::find_closest_point(int i, double *x, double &xc, double &yc, double &zc)
double &xc, double &yc, double &zc)
{ {
double dot,d2,d2min; double dot, d2, d2min;
double xr[3],xproj[3],p[3]; double xr[3], xproj[3], p[3];
xr[0] = x[0] - corners[i][0][0]; xr[0] = x[0] - corners[i][0][0];
xr[1] = x[1] - corners[i][0][1]; xr[1] = x[1] - corners[i][0][1];
xr[2] = x[2] - corners[i][0][2]; xr[2] = x[2] - corners[i][0][2];
dot = face[i][0]*xr[0] + face[i][1]*xr[1] + face[i][2]*xr[2]; dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2];
xproj[0] = xr[0] - dot*face[i][0]; xproj[0] = xr[0] - dot * face[i][0];
xproj[1] = xr[1] - dot*face[i][1]; xproj[1] = xr[1] - dot * face[i][1];
xproj[2] = xr[2] - dot*face[i][2]; xproj[2] = xr[2] - dot * face[i][2];
d2min = BIG; d2min = BIG;
// check if point projects inside of face // check if point projects inside of face
if (inside_face(xproj, i)) { if (inside_face(xproj, i)) {
d2 = d2min = dot*dot; d2 = d2min = dot * dot;
xc = xproj[0] + corners[i][0][0]; xc = xproj[0] + corners[i][0][0];
yc = xproj[1] + corners[i][0][1]; yc = xproj[1] + corners[i][0][1];
zc = xproj[2] + corners[i][0][2]; zc = xproj[2] + corners[i][0][2];
// check each edge // check each edge
} else { } else {
point_on_line_segment(corners[i][0],corners[i][1],x,p); point_on_line_segment(corners[i][0], corners[i][1], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -582,9 +596,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][1],corners[i][2],x,p); point_on_line_segment(corners[i][1], corners[i][2], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -592,9 +606,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][2],corners[i][3],x,p); point_on_line_segment(corners[i][2], corners[i][3], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -602,9 +616,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][3],corners[i][0],x,p); point_on_line_segment(corners[i][3], corners[i][0], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -623,14 +637,12 @@ double RegBlock::find_closest_point(int i, double *x,
int RegBlock::inside_face(double *xproj, int iface) int RegBlock::inside_face(double *xproj, int iface)
{ {
if (iface < 2) { if (iface < 2) {
if (xproj[1] > 0 && (xproj[1] < yhi-ylo) && if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1;
xproj[2] > 0 && (xproj[2] < zhi-zlo)) return 1;
} else if (iface < 4) { } else if (iface < 4) {
if (xproj[0] > 0 && (xproj[0] < (xhi-xlo)) && if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo)))
xproj[2] > 0 && (xproj[2] < (zhi-zlo))) return 1; return 1;
} else { } else {
if (xproj[0] > 0 && xproj[0] < (xhi-xlo) && if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1;
xproj[1] > 0 && xproj[1] < (yhi-ylo)) return 1;
} }
return 0; return 0;

335
src/region_cone.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -26,102 +25,120 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0)
Region(lmp, narg, arg), lo(0.0), hi(0.0)
{ {
options(narg-9,&arg[9]); options(narg - 9, &arg[9]);
// check open face settings // check open face settings
if (openflag) if (openflag)
for (int i=3; i<6; i++) for (int i = 3; i < 6; i++)
if (open_faces[i]) error->all(FLERR,"Illegal region cone open face: {}", i+1); if (open_faces[i]) error->all(FLERR, "Illegal region cone open face: {}", i + 1);
if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0)
error->all(FLERR,"Illegal region cone axis: {}", arg[2]); error->all(FLERR, "Illegal region cone axis: {}", arg[2]);
axis = arg[2][0]; axis = arg[2][0];
if (axis == 'x') { if (axis == 'x') {
c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = yscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = yscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp);
} else if (axis == 'y') { } else if (axis == 'y') {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
} else if (axis == 'z') { } else if (axis == 'z') {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
} }
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[0]; lo = -BIG;
else lo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
lo = domain->boxlo[0];
else
lo = domain->boxlo_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[1]; lo = -BIG;
else lo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
lo = domain->boxlo[1];
else
lo = domain->boxlo_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[2]; lo = -BIG;
else lo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
lo = domain->boxlo[2];
else
lo = domain->boxlo_bound[2];
} }
} else { } else {
if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} }
if (strcmp(arg[8],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[8], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[0]; hi = BIG;
else hi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
hi = domain->boxhi[0];
else
hi = domain->boxhi_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0) hi = BIG;
if (domain->triclinic == 0) hi = domain->boxhi[1]; if (domain->triclinic == 0)
else hi = domain->boxhi_bound[1]; hi = domain->boxhi[1];
else
hi = domain->boxhi_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[2]; hi = BIG;
else hi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
hi = domain->boxhi[2];
else
hi = domain->boxhi_bound[2];
} }
} else { } else {
if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[8], false, lmp);
} }
// error check // error check
if (radiuslo < 0.0) error->all(FLERR,"Illegal radius in region cone command"); if (radiuslo < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiushi < 0.0) error->all(FLERR,"Illegal radius in region cone command"); if (radiushi < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiuslo == 0.0 && radiushi == 0.0) if (radiuslo == 0.0 && radiushi == 0.0)
error->all(FLERR,"Illegal radius in region cone command"); error->all(FLERR, "Illegal radius in region cone command");
if (hi == lo) error->all(FLERR,"Illegal cone length in region cone command"); if (hi == lo) error->all(FLERR, "Illegal cone length in region cone command");
// extent of cone // extent of cone
if (radiuslo > radiushi) maxradius = radiuslo; if (radiuslo > radiushi)
else maxradius = radiushi; maxradius = radiuslo;
else
maxradius = radiushi;
if (interior) { if (interior) {
bboxflag = 1; bboxflag = 1;
@ -149,22 +166,25 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
extent_zlo = lo; extent_zlo = lo;
extent_zhi = hi; extent_zhi = hi;
} }
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to cone surface and 2 ends // particle could be close to cone surface and 2 ends
// particle can only touch surface and 1 end // particle can only touch surface and 1 end
cmax = 3; cmax = 3;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 2; if (interior)
else tmax = 1; tmax = 2;
else
tmax = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCone::~RegCone() RegCone::~RegCone()
{ {
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -174,30 +194,36 @@ RegCone::~RegCone()
int RegCone::inside(double x, double y, double z) int RegCone::inside(double x, double y, double z)
{ {
double del1,del2,dist; double del1, del2, dist;
double currentradius; double currentradius;
if (axis == 'x') { if (axis == 'x') {
del1 = y - c1; del1 = y - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && x >= lo && x <= hi) return 1; if (dist <= currentradius && x >= lo && x <= hi)
else return 0; return 1;
else
return 0;
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x - c1; del1 = x - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (y-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (y - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && y >= lo && y <= hi) return 1; if (dist <= currentradius && y >= lo && y <= hi)
else return 0; return 1;
else
return 0;
} else if (axis == 'z') { } else if (axis == 'z') {
del1 = x - c1; del1 = x - c1;
del2 = y - c2; del2 = y - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (z-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (z - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && z >= lo && z <= hi) return 1; if (dist <= currentradius && z >= lo && z <= hi)
else return 0; return 1;
else
return 0;
} }
return 0; return 0;
@ -213,16 +239,16 @@ int RegCone::inside(double x, double y, double z)
int RegCone::surface_interior(double *x, double cutoff) int RegCone::surface_interior(double *x, double cutoff)
{ {
double del1,del2,r,currentradius,delx,dely,delz,dist,delta; double del1, del2, r, currentradius, delx, dely, delz, dist, delta;
double surflo[3],surfhi[3],xs[3]; double surflo[3], surfhi[3], xs[3];
int n = 0; int n = 0;
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo);
// x is exterior to cone // x is exterior to cone
@ -234,23 +260,22 @@ int RegCone::surface_interior(double *x, double cutoff)
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = lo; surflo[0] = lo;
surflo[1] = c1 + del1*radiuslo/r; surflo[1] = c1 + del1 * radiuslo / r;
surflo[2] = c2 + del2*radiuslo/r; surflo[2] = c2 + del2 * radiuslo / r;
surfhi[0] = hi; surfhi[0] = hi;
surfhi[1] = c1 + del1*radiushi/r; surfhi[1] = c1 + del1 * radiushi / r;
surfhi[2] = c2 + del2*radiushi/r; surfhi[2] = c2 + del2 * radiushi / r;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].radius = -2.0*(radiuslo + (xs[0]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[0] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
contact[n].iwall = 2; contact[n].iwall = 2;
n++; n++;
} }
@ -278,9 +303,8 @@ int RegCone::surface_interior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[1]-lo)* currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo);
(radiushi-radiuslo)/(hi-lo);
// y is exterior to cone // y is exterior to cone
@ -291,25 +315,24 @@ int RegCone::surface_interior(double *x, double cutoff)
// surfhi = pt on outer circle of top end plane, same dir as y vs axis // surfhi = pt on outer circle of top end plane, same dir as y vs axis
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r; surflo[0] = c1 + del1 * radiuslo / r;
surflo[1] = lo; surflo[1] = lo;
surflo[2] = c2 + del2*radiuslo/r; surflo[2] = c2 + del2 * radiuslo / r;
surfhi[0] = c1 + del1*radiushi/r; surfhi[0] = c1 + del1 * radiushi / r;
surfhi[1] = hi; surfhi[1] = hi;
surfhi[2] = c2 + del2*radiushi/r; surfhi[2] = c2 + del2 * radiushi / r;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[1]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[1] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
n++; n++;
} }
} }
@ -336,8 +359,8 @@ int RegCone::surface_interior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo);
// z is exterior to cone // z is exterior to cone
@ -348,31 +371,30 @@ int RegCone::surface_interior(double *x, double cutoff)
// surfhi = pt on outer circle of top end plane, same dir as z vs axis // surfhi = pt on outer circle of top end plane, same dir as z vs axis
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r; surflo[0] = c1 + del1 * radiuslo / r;
surflo[1] = c2 + del2*radiuslo/r; surflo[1] = c2 + del2 * radiuslo / r;
surflo[2] = lo; surflo[2] = lo;
surfhi[0] = c1 + del1*radiushi/r; surfhi[0] = c1 + del1 * radiushi / r;
surfhi[1] = c2 + del2*radiushi/r; surfhi[1] = c2 + del2 * radiushi / r;
surfhi[2] = hi; surfhi[2] = hi;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[2]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[2] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
n++; n++;
} }
} }
delta = x[2] - lo; delta = x[2] - lo;
if (delta < cutoff && !open_faces[0]) { if (delta < cutoff && !open_faces[0]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delz = delta; contact[n].delz = delta;
contact[n].delx = contact[n].dely = 0.0; contact[n].delx = contact[n].dely = 0.0;
@ -381,7 +403,7 @@ int RegCone::surface_interior(double *x, double cutoff)
n++; n++;
} }
delta = hi - x[2]; delta = hi - x[2];
if (delta < cutoff && !open_faces[1]) { if (delta < cutoff && !open_faces[1]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delz = -delta; contact[n].delz = -delta;
contact[n].delx = contact[n].dely = 0.0; contact[n].delx = contact[n].dely = 0.0;
@ -402,14 +424,14 @@ int RegCone::surface_interior(double *x, double cutoff)
int RegCone::surface_exterior(double *x, double cutoff) int RegCone::surface_exterior(double *x, double cutoff)
{ {
double del1,del2,r,currentradius,distsq,distsqprev,crad; double del1, del2, r, currentradius, distsq, distsqprev, crad;
double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3]; double corner1[3], corner2[3], corner3[3], corner4[3], xp[3], nearest[3];
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -418,8 +440,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// x is far enough from cone that there is no contact // x is far enough from cone that there is no contact
// x is interior to cone // x is interior to cone
if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) if (r >= maxradius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0;
return 0;
if (r < currentradius && x[0] > lo && x[0] < hi) return 0; if (r < currentradius && x[0] > lo && x[0] < hi) return 0;
// x is exterior to cone or on its surface // x is exterior to cone or on its surface
@ -430,11 +451,11 @@ int RegCone::surface_exterior(double *x, double cutoff)
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = lo; corner1[0] = lo;
corner1[1] = c1 + del1*radiuslo/r; corner1[1] = c1 + del1 * radiuslo / r;
corner1[2] = c2 + del2*radiuslo/r; corner1[2] = c2 + del2 * radiuslo / r;
corner2[0] = hi; corner2[0] = hi;
corner2[1] = c1 + del1*radiushi/r; corner2[1] = c1 + del1 * radiushi / r;
corner2[2] = c2 + del2*radiushi/r; corner2[2] = c2 + del2 * radiushi / r;
corner3[0] = lo; corner3[0] = lo;
corner3[1] = c1; corner3[1] = c1;
corner3[2] = c2; corner3[2] = c2;
@ -445,27 +466,27 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[0]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[0] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0.0; if (distsq < distsqprev) crad = 0.0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0.0; if (distsq < distsqprev) crad = 0.0;
} }
if (distsq == BIG) return 0; if (distsq == BIG) return 0;
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -474,8 +495,8 @@ int RegCone::surface_exterior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -484,8 +505,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// y is far enough from cone that there is no contact // y is far enough from cone that there is no contact
// y is interior to cone // y is interior to cone
if (r >= maxradius+cutoff || if (r >= maxradius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0;
x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
if (r < currentradius && x[1] > lo && x[1] < hi) return 0; if (r < currentradius && x[1] > lo && x[1] < hi) return 0;
// y is exterior to cone or on its surface // y is exterior to cone or on its surface
@ -495,12 +515,12 @@ int RegCone::surface_exterior(double *x, double cutoff)
// could be edge of cone // could be edge of cone
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = c1 + del1*radiuslo/r; corner1[0] = c1 + del1 * radiuslo / r;
corner1[1] = lo; corner1[1] = lo;
corner1[2] = c2 + del2*radiuslo/r; corner1[2] = c2 + del2 * radiuslo / r;
corner2[0] = c1 + del1*radiushi/r; corner2[0] = c1 + del1 * radiushi / r;
corner2[1] = hi; corner2[1] = hi;
corner2[2] = c2 + del2*radiushi/r; corner2[2] = c2 + del2 * radiushi / r;
corner3[0] = c1; corner3[0] = c1;
corner3[1] = lo; corner3[1] = lo;
corner3[2] = c2; corner3[2] = c2;
@ -511,26 +531,26 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[1]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[1] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -539,8 +559,8 @@ int RegCone::surface_exterior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -549,8 +569,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// z is far enough from cone that there is no contact // z is far enough from cone that there is no contact
// z is interior to cone // z is interior to cone
if (r >= maxradius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) if (r >= maxradius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0;
return 0;
if (r < currentradius && x[2] > lo && x[2] < hi) return 0; if (r < currentradius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cone or on its surface // z is exterior to cone or on its surface
@ -560,11 +579,11 @@ int RegCone::surface_exterior(double *x, double cutoff)
// could be edge of cone // could be edge of cone
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = c1 + del1*radiuslo/r; corner1[0] = c1 + del1 * radiuslo / r;
corner1[1] = c2 + del2*radiuslo/r; corner1[1] = c2 + del2 * radiuslo / r;
corner1[2] = lo; corner1[2] = lo;
corner2[0] = c1 + del1*radiushi/r; corner2[0] = c1 + del1 * radiushi / r;
corner2[1] = c2 + del2*radiushi/r; corner2[1] = c2 + del2 * radiushi / r;
corner2[2] = hi; corner2[2] = hi;
corner3[0] = c1; corner3[0] = c1;
corner3[1] = c2; corner3[1] = c2;
@ -576,26 +595,26 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[2]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[2] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -610,7 +629,7 @@ double RegCone::closest(double *x, double *near, double *nearest, double dsq)
double delx = x[0] - near[0]; double delx = x[0] - near[0];
double dely = x[1] - near[1]; double dely = x[1] - near[1];
double delz = x[2] - near[2]; double delz = x[2] - near[2];
double rsq = delx*delx + dely*dely + delz*delz; double rsq = delx * delx + dely * dely + delz * delz;
if (rsq >= dsq) return dsq; if (rsq >= dsq) return dsq;
nearest[0] = near[0]; nearest[0] = near[0];

428
src/region_cylinder.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -25,94 +24,97 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
enum{CONSTANT,VARIABLE};
enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr) Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr)
{ {
options(narg-8,&arg[8]); options(narg - 8, &arg[8]);
// check open face settings // check open face settings
if (openflag) if (openflag)
for (int i=3; i<6; i++) for (int i = 3; i < 6; i++)
if (open_faces[i]) error->all(FLERR,"Illegal region cylinder open face: {}", i+1); if (open_faces[i]) error->all(FLERR, "Illegal region cylinder open face: {}", i + 1);
if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0)
error->all(FLERR,"Illegal region cylinder axis: {}", arg[2]); error->all(FLERR, "Illegal region cylinder axis: {}", arg[2]);
axis = arg[2][0]; axis = arg[2][0];
if (axis == 'x') { if (axis == 'x') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} else if (axis == 'y') { } else if (axis == 'y') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} else if (axis == 'z') { } else if (axis == 'z') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} }
if (utils::strmatch(arg[5],"^v_")) { if (utils::strmatch(arg[5], "^v_")) {
rstr = utils::strdup(arg[5]+2); rstr = utils::strdup(arg[5] + 2);
radius = 0.0; radius = 0.0;
rstyle = VARIABLE; rstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
radius = utils::numeric(FLERR,arg[5],false,lmp); radius = utils::numeric(FLERR, arg[5], false, lmp);
if (axis == 'x') radius *= yscale; if (axis == 'x')
else radius *= xscale; radius *= yscale;
else
radius *= xscale;
rstyle = CONSTANT; rstyle = CONSTANT;
} }
@ -121,57 +123,75 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
RegCylinder::shape_update(); RegCylinder::shape_update();
} }
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[0]; lo = -BIG;
else lo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
lo = domain->boxlo[0];
else
lo = domain->boxlo_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[1]; lo = -BIG;
else lo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
lo = domain->boxlo[1];
else
lo = domain->boxlo_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[2]; lo = -BIG;
else lo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
lo = domain->boxlo[2];
else
lo = domain->boxlo_bound[2];
} }
} else { } else {
if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[6], false, lmp);
if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[6], false, lmp);
if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
} }
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[0]; hi = BIG;
else hi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
hi = domain->boxhi[0];
else
hi = domain->boxhi_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[1]; hi = BIG;
else hi = domain->boxhi_bound[1]; else if (domain->triclinic == 0)
hi = domain->boxhi[1];
else
hi = domain->boxhi_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[2]; hi = BIG;
else hi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
hi = domain->boxhi[2];
else
hi = domain->boxhi_bound[2];
} }
} else { } else {
if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} }
// error check // error check
if (radius <= 0.0) error->all(FLERR,"Illegal radius {} in region cylinder command", radius); if (radius <= 0.0) error->all(FLERR, "Illegal radius {} in region cylinder command", radius);
// extent of cylinder // extent of cylinder
// for variable radius, uses initial radius // for variable radius, uses initial radius
@ -202,25 +222,28 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
extent_zlo = lo; extent_zlo = lo;
extent_zhi = hi; extent_zhi = hi;
} }
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to cylinder surface and 2 ends // particle could be close to cylinder surface and 2 ends
// particle can only touch surface and 1 end // particle can only touch surface and 1 end
cmax = 3; cmax = 3;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 2; if (interior)
else tmax = 1; tmax = 2;
else
tmax = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCylinder::~RegCylinder() RegCylinder::~RegCylinder()
{ {
delete [] c1str; delete[] c1str;
delete [] c2str; delete[] c2str;
delete [] rstr; delete[] rstr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -238,27 +261,33 @@ void RegCylinder::init()
int RegCylinder::inside(double x, double y, double z) int RegCylinder::inside(double x, double y, double z)
{ {
double del1,del2,dist; double del1, del2, dist;
int inside; int inside;
if (axis == 'x') { if (axis == 'x') {
del1 = y - c1; del1 = y - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && x >= lo && x <= hi) inside = 1; if (dist <= radius && x >= lo && x <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x - c1; del1 = x - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && y >= lo && y <= hi) inside = 1; if (dist <= radius && y >= lo && y <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} else { } else {
del1 = x - c1; del1 = x - c1;
del2 = y - c2; del2 = y - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && z >= lo && z <= hi) inside = 1; if (dist <= radius && z >= lo && z <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} }
return inside; return inside;
@ -274,14 +303,14 @@ int RegCylinder::inside(double x, double y, double z)
int RegCylinder::surface_interior(double *x, double cutoff) int RegCylinder::surface_interior(double *x, double cutoff)
{ {
double del1,del2,r,delta; double del1, del2, r, delta;
int n = 0; int n = 0;
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// x is exterior to cylinder // x is exterior to cylinder
@ -293,9 +322,9 @@ int RegCylinder::surface_interior(double *x, double cutoff)
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = 0.0; contact[n].delx = 0.0;
contact[n].dely = del1*(1.0-radius/r); contact[n].dely = del1 * (1.0 - radius / r);
contact[n].delz = del2*(1.0-radius/r); contact[n].delz = del2 * (1.0 - radius / r);
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -324,7 +353,7 @@ int RegCylinder::surface_interior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// y is exterior to cylinder // y is exterior to cylinder
@ -335,10 +364,10 @@ int RegCylinder::surface_interior(double *x, double cutoff)
delta = radius - r; delta = radius - r;
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = del1*(1.0-radius/r); contact[n].delx = del1 * (1.0 - radius / r);
contact[n].dely = 0.0; contact[n].dely = 0.0;
contact[n].delz = del2*(1.0-radius/r); contact[n].delz = del2 * (1.0 - radius / r);
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -367,7 +396,7 @@ int RegCylinder::surface_interior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// z is exterior to cylinder // z is exterior to cylinder
@ -378,10 +407,10 @@ int RegCylinder::surface_interior(double *x, double cutoff)
delta = radius - r; delta = radius - r;
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = del1*(1.0-radius/r); contact[n].delx = del1 * (1.0 - radius / r);
contact[n].dely = del2*(1.0-radius/r); contact[n].dely = del2 * (1.0 - radius / r);
contact[n].delz = 0.0; contact[n].delz = 0.0;
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -419,12 +448,12 @@ int RegCylinder::surface_interior(double *x, double cutoff)
int RegCylinder::surface_exterior(double *x, double cutoff) int RegCylinder::surface_exterior(double *x, double cutoff)
{ {
double del1,del2,r; double del1, del2, r;
double xp,yp,zp; double xp, yp, zp;
double dx, dr, dr2, d2, d2prev; double dx, dr, dr2, d2, d2prev;
// radius of curvature for granular // radius of curvature for granular
// 0 for flat surfaces (infinite case), 2*radius for curved portion // 0 for flat surfaces (infinite case), 2*radius for curved portion
double crad = 0.0; double crad = 0.0;
int varflag = 0; int varflag = 0;
@ -432,12 +461,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// x is far enough from cylinder that there is no contact // x is far enough from cylinder that there is no contact
// x is interior to cylinder // x is interior to cylinder
if (r >= radius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0;
if (r < radius && x[0] > lo && x[0] < hi) return 0; if (r < radius && x[0] > lo && x[0] < hi) return 0;
// x is exterior to cylinder or on its surface // x is exterior to cylinder or on its surface
@ -448,48 +477,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
yp = c1 + del1*radius/r; yp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
yp = x[1]; yp = x[1];
zp = x[2]; zp = x[2];
} }
if (x[0] < lo) xp = lo; if (x[0] < lo)
else if (x[0] > hi) xp = hi; xp = lo;
else xp = x[0]; else if (x[0] > hi)
xp = hi;
else
xp = x[0];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
yp = c1 + del1*radius/r; yp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
if (x[0] < lo) { if (x[0] < lo) {
dx = lo-x[0]; dx = lo - x[0];
xp = lo; xp = lo;
} } else if (x[0] > hi) {
else if (x[0] > hi) { dx = x[0] - hi;
dx = x[0]-hi;
xp = hi; xp = hi;
} } else {
else {
dx = 0; dx = 0;
xp = x[0]; xp = x[0];
} }
d2 = d2prev = dr2 + dx*dx; d2 = d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[0]; dx = lo - x[0];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
xp = lo; xp = lo;
if (r < radius) { if (r < radius) {
@ -504,8 +536,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[0]; dx = hi - x[0];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
xp = hi; xp = hi;
if (r < radius) { if (r < radius) {
@ -516,7 +550,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -526,12 +560,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// y is far enough from cylinder that there is no contact // y is far enough from cylinder that there is no contact
// y is interior to cylinder // y is interior to cylinder
if (r >= radius+cutoff || x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0;
if (r < radius && x[1] > lo && x[1] < hi) return 0; if (r < radius && x[1] > lo && x[1] < hi) return 0;
// y is exterior to cylinder or on its surface // y is exterior to cylinder or on its surface
@ -542,48 +576,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
xp = x[0]; xp = x[0];
zp = x[2]; zp = x[2];
} }
if (x[1] < lo) yp = lo; if (x[1] < lo)
else if (x[1] > hi) yp = hi; yp = lo;
else yp = x[1]; else if (x[1] > hi)
yp = hi;
else
yp = x[1];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
if (x[1] < lo) { if (x[1] < lo) {
dx = lo-x[1]; dx = lo - x[1];
yp = lo; yp = lo;
} } else if (x[1] > hi) {
else if (x[1] > hi) { dx = x[1] - hi;
dx = x[1]-hi;
yp = hi; yp = hi;
} } else {
else {
dx = 0; dx = 0;
yp = x[1]; yp = x[1];
} }
d2 = d2prev = dr2 + dx*dx; d2 = d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[1]; dx = lo - x[1];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
yp = lo; yp = lo;
if (r < radius) { if (r < radius) {
@ -598,8 +635,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[1]; dx = hi - x[1];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
yp = hi; yp = hi;
if (r < radius) { if (r < radius) {
@ -610,7 +649,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -620,12 +659,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// z is far enough from cylinder that there is no contact // z is far enough from cylinder that there is no contact
// z is interior to cylinder // z is interior to cylinder
if (r >= radius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0;
if (r < radius && x[2] > lo && x[2] < hi) return 0; if (r < radius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cylinder or on its surface // z is exterior to cylinder or on its surface
@ -636,46 +675,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
yp = c2 + del2*radius/r; yp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
xp = x[0]; xp = x[0];
yp = x[1]; yp = x[1];
} }
if (x[2] < lo) zp = lo; if (x[2] < lo)
else if (x[2] > hi) zp = hi; zp = lo;
else zp = x[2]; else if (x[2] > hi)
zp = hi;
else
zp = x[2];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
yp = c2 + del2*radius/r; yp = c2 + del2 * radius / r;
if (x[2] < lo) { if (x[2] < lo) {
dx = lo-x[2]; dx = lo - x[2];
zp = lo; zp = lo;
} else if (x[2] > hi) { } else if (x[2] > hi) {
dx = x[2]-hi; dx = x[2] - hi;
zp = hi; zp = hi;
} else { } else {
dx = 0; dx = 0;
zp = x[2]; zp = x[2];
} }
d2prev = dr2 + dx*dx; d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[2]; dx = lo - x[2];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
zp = lo; zp = lo;
if (r < radius) { if (r < radius) {
@ -690,8 +734,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[2]; dx = hi - x[2];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
zp = hi; zp = hi;
if (r < radius) { if (r < radius) {
@ -702,7 +748,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -721,22 +767,21 @@ void RegCylinder::shape_update()
if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var); if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var);
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
radius = input->variable->compute_equal(rvar); radius = input->variable->compute_equal(rvar);
if (radius < 0.0) if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
error->one(FLERR,"Variable evaluation in region gave bad value");
} }
if (axis == 'x') { if (axis == 'x') {
if (c1style == VARIABLE) c1 *= yscale; if (c1style == VARIABLE) c1 *= yscale;
if (c2style == VARIABLE) c2 *= zscale; if (c2style == VARIABLE) c2 *= zscale;
if (rstyle == VARIABLE) radius *= yscale; if (rstyle == VARIABLE) radius *= yscale;
} else if (axis == 'y') { } else if (axis == 'y') {
if (c1style == VARIABLE) c1 *= xscale; if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= zscale; if (c2style == VARIABLE) c2 *= zscale;
if (rstyle == VARIABLE) radius *= xscale; if (rstyle == VARIABLE) radius *= xscale;
} else { // axis == 'z' } else { // axis == 'z'
if (c1style == VARIABLE) c1 *= xscale; if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= yscale; if (c2style == VARIABLE) c2 *= yscale;
if (rstyle == VARIABLE) radius *= xscale; if (rstyle == VARIABLE) radius *= xscale;
} }
} }
@ -748,30 +793,26 @@ void RegCylinder::variable_check()
{ {
if (c1style == VARIABLE) { if (c1style == VARIABLE) {
c1var = input->variable->find(c1str); c1var = input->variable->find(c1str);
if (c1var < 0) if (c1var < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", c1str);
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(c1var)) if (!input->variable->equalstyle(c1var))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable {} for region cylinder is invalid style", c1str);
} }
if (c2style == VARIABLE) { if (c2style == VARIABLE) {
c2var = input->variable->find(c2str); c2var = input->variable->find(c2str);
if (c2var < 0) if (c2var < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", c2str);
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(c2var)) if (!input->variable->equalstyle(c2var))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable {} for region cylinder is invalid style", c2str);
} }
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
rvar = input->variable->find(rstr); rvar = input->variable->find(rstr);
if (rvar < 0) if (rvar < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", rstr);
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(rvar)) if (!input->variable->equalstyle(rvar))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable {} for region cylinder is invalid style", rstr);
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Set values needed to calculate velocity due to shape changes. Set values needed to calculate velocity due to shape changes.
These values do not depend on the contact, so this function is These values do not depend on the contact, so this function is
@ -795,35 +836,34 @@ void RegCylinder::set_velocity_shape()
xcenter[2] = 0; xcenter[2] = 0;
} }
forward_transform(xcenter[0], xcenter[1], xcenter[2]); forward_transform(xcenter[0], xcenter[1], xcenter[2]);
if (update->ntimestep > 0) rprev = prev[4]; if (update->ntimestep > 0)
else rprev = radius; rprev = prev[4];
else
rprev = radius;
prev[4] = radius; prev[4] = radius;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add velocity due to shape change to wall velocity add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegCylinder::velocity_contact_shape(double *vwall, double *xc) void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
{ {
double delx, dely, delz; // Displacement of contact point in x,y,z double delx, dely, delz; // Displacement of contact point in x,y,z
if (axis == 'x') { if (axis == 'x') {
delx = 0; delx = 0;
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
} else if (axis == 'y') { } else if (axis == 'y') {
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = 0; dely = 0;
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
} else { } else {
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = 0; delz = 0;
} }
vwall[0] += delx/update->dt; vwall[0] += delx / update->dt;
vwall[1] += dely/update->dt; vwall[1] += dely / update->dt;
vwall[2] += delz/update->dt; vwall[2] += delz / update->dt;
} }

339
src/region_prism.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -26,105 +25,126 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
{ {
options(narg-11,&arg[11]); options(narg - 11, &arg[11]);
if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[2],"INF") == 0) xlo = -BIG; if (strcmp(arg[2], "INF") == 0)
else xlo = domain->boxlo[0]; xlo = -BIG;
} else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); else
xlo = domain->boxlo[0];
} else
xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp);
if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[3],"INF") == 0) xhi = BIG; if (strcmp(arg[3], "INF") == 0)
else xhi = domain->boxhi[0]; xhi = BIG;
} else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); else
xhi = domain->boxhi[0];
} else
xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp);
if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[4],"INF") == 0) ylo = -BIG; if (strcmp(arg[4], "INF") == 0)
else ylo = domain->boxlo[1]; ylo = -BIG;
} else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); else
ylo = domain->boxlo[1];
} else
ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp);
if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[5],"INF") == 0) yhi = BIG; if (strcmp(arg[5], "INF") == 0)
else yhi = domain->boxhi[1]; yhi = BIG;
} else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); else
yhi = domain->boxhi[1];
} else
yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp);
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[6],"INF") == 0) zlo = -BIG; if (strcmp(arg[6], "INF") == 0)
else zlo = domain->boxlo[2]; zlo = -BIG;
} else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); else
zlo = domain->boxlo[2];
} else
zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[7],"INF") == 0) zhi = BIG; if (strcmp(arg[7], "INF") == 0)
else zhi = domain->boxhi[2]; zhi = BIG;
} else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); else
zhi = domain->boxhi[2];
} else
zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
xy = xscale*utils::numeric(FLERR,arg[8],false,lmp); xy = xscale * utils::numeric(FLERR, arg[8], false, lmp);
xz = xscale*utils::numeric(FLERR,arg[9],false,lmp); xz = xscale * utils::numeric(FLERR, arg[9], false, lmp);
yz = yscale*utils::numeric(FLERR,arg[10],false,lmp); yz = yscale * utils::numeric(FLERR, arg[10], false, lmp);
// error check // error check
// prism cannot be 0 thickness in any dim, else inverse blows up // prism cannot be 0 thickness in any dim, else inverse blows up
// non-zero tilt values cannot be used if either dim is INF on both ends // non-zero tilt values cannot be used if either dim is INF on both ends
if (xlo >= xhi) error->all(FLERR,"Illegal region prism xlo: {} >= xhi: {}", xlo, xhi); if (xlo >= xhi) error->all(FLERR, "Illegal region prism xlo: {} >= xhi: {}", xlo, xhi);
if (ylo >= yhi) error->all(FLERR,"Illegal region prism ylo: {} >= yhi: {}", ylo, yhi); if (ylo >= yhi) error->all(FLERR, "Illegal region prism ylo: {} >= yhi: {}", ylo, yhi);
if (zlo >= zhi) error->all(FLERR,"Illegal region prism zlo: {} >= zhi: {}", zlo ,zhi); if (zlo >= zhi) error->all(FLERR, "Illegal region prism zlo: {} >= zhi: {}", zlo, zhi);
if (xy != 0.0 && xlo == -BIG && xhi == BIG) if (xy != 0.0 && xlo == -BIG && xhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite x size"); error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite x size");
if (xy != 0.0 && ylo == -BIG && yhi == BIG) if (xy != 0.0 && ylo == -BIG && yhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite y size"); error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite y size");
if (xz != 0.0 && xlo == -BIG && xhi == BIG) if (xz != 0.0 && xlo == -BIG && xhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite x size"); error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite x size");
if (xz != 0.0 && zlo == -BIG && zhi == BIG) if (xz != 0.0 && zlo == -BIG && zhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite z size"); error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite z size");
if (yz != 0.0 && ylo == -BIG && yhi == BIG) if (yz != 0.0 && ylo == -BIG && yhi == BIG)
error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite y size"); error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite y size");
if (yz != 0.0 && zlo == -BIG && zhi == BIG) if (yz != 0.0 && zlo == -BIG && zhi == BIG)
error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite z size"); error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite z size");
// extent of prism // extent of prism
if (interior) { if (interior) {
bboxflag = 1; bboxflag = 1;
extent_xlo = MIN(xlo,xlo+xy); extent_xlo = MIN(xlo, xlo + xy);
extent_xlo = MIN(extent_xlo,extent_xlo+xz); extent_xlo = MIN(extent_xlo, extent_xlo + xz);
extent_ylo = MIN(ylo,ylo+yz); extent_ylo = MIN(ylo, ylo + yz);
extent_zlo = zlo; extent_zlo = zlo;
extent_xhi = MAX(xhi,xhi+xy); extent_xhi = MAX(xhi, xhi + xy);
extent_xhi = MAX(extent_xhi,extent_xhi+xz); extent_xhi = MAX(extent_xhi, extent_xhi + xz);
extent_yhi = MAX(yhi,yhi+yz); extent_yhi = MAX(yhi, yhi + yz);
extent_zhi = zhi; extent_zhi = zhi;
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to all 6 planes // particle could be close to all 6 planes
// particle can only touch 3 planes // particle can only touch 3 planes
cmax = 6; cmax = 6;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 3; if (interior)
else tmax = 1; tmax = 3;
else
tmax = 1;
// h = transformation matrix from tilt coords (0-1) to box coords (xyz) // h = transformation matrix from tilt coords (0-1) to box coords (xyz)
// columns of h are edge vectors of tilted box // columns of h are edge vectors of tilted box
@ -140,26 +160,26 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
h[1][2] = yz; h[1][2] = yz;
h[2][2] = zhi - zlo; h[2][2] = zhi - zlo;
hinv[0][0] = 1.0/h[0][0]; hinv[0][0] = 1.0 / h[0][0];
hinv[0][1] = -h[0][1] / (h[0][0]*h[1][1]); hinv[0][1] = -h[0][1] / (h[0][0] * h[1][1]);
hinv[0][2] = (h[0][1]*h[1][2] - h[0][2]*h[1][1]) / (h[0][0]*h[1][1]*h[2][2]); hinv[0][2] = (h[0][1] * h[1][2] - h[0][2] * h[1][1]) / (h[0][0] * h[1][1] * h[2][2]);
hinv[1][1] = 1.0/h[1][1]; hinv[1][1] = 1.0 / h[1][1];
hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]); hinv[1][2] = -h[1][2] / (h[1][1] * h[2][2]);
hinv[2][2] = 1.0/h[2][2]; hinv[2][2] = 1.0 / h[2][2];
// corners = 8 corner points of prism // corners = 8 corner points of prism
// order = x varies fastest, then y, finally z // order = x varies fastest, then y, finally z
// clo/chi = lo and hi corner pts of prism // clo/chi = lo and hi corner pts of prism
a[0] = xhi-xlo; a[0] = xhi - xlo;
a[1] = 0.0; a[1] = 0.0;
a[2] = 0.0; a[2] = 0.0;
b[0] = xy; b[0] = xy;
b[1] = yhi-ylo; b[1] = yhi - ylo;
b[2] = 0.0; b[2] = 0.0;
c[0] = xz; c[0] = xz;
c[1] = yz; c[1] = yz;
c[2] = zhi-zlo; c[2] = zhi - zlo;
clo[0] = corners[0][0] = xlo; clo[0] = corners[0][0] = xlo;
clo[1] = corners[0][1] = ylo; clo[1] = corners[0][1] = ylo;
@ -196,19 +216,18 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// face = 6 inward-facing unit normals to prism faces // face = 6 inward-facing unit normals to prism faces
// order = xy plane, xz plane, yz plane // order = xy plane, xz plane, yz plane
MathExtra::cross3(a,b,face[0]); MathExtra::cross3(a, b, face[0]);
MathExtra::cross3(b,a,face[1]); MathExtra::cross3(b, a, face[1]);
MathExtra::cross3(c,a,face[2]); MathExtra::cross3(c, a, face[2]);
MathExtra::cross3(a,c,face[3]); MathExtra::cross3(a, c, face[3]);
MathExtra::cross3(b,c,face[4]); MathExtra::cross3(b, c, face[4]);
MathExtra::cross3(c,b,face[5]); MathExtra::cross3(c, b, face[5]);
// remap open face indices to be consistent // remap open face indices to be consistent
if (openflag) { if (openflag) {
int temp[6]; int temp[6];
for (int i = 0; i < 6; i++) for (int i = 0; i < 6; i++) temp[i] = open_faces[i];
temp[i] = open_faces[i];
open_faces[0] = temp[4]; open_faces[0] = temp[4];
open_faces[1] = temp[5]; open_faces[1] = temp[5];
open_faces[2] = temp[2]; open_faces[2] = temp[2];
@ -223,27 +242,51 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// verts in each tri are ordered so that right-hand rule gives inward norm // verts in each tri are ordered so that right-hand rule gives inward norm
// order = xy plane, xz plane, yz plane // order = xy plane, xz plane, yz plane
tri[0][0] = 0; tri[0][1] = 1; tri[0][2] = 3; tri[0][0] = 0;
tri[1][0] = 0; tri[1][1] = 3; tri[1][2] = 2; tri[0][1] = 1;
tri[2][0] = 4; tri[2][1] = 7; tri[2][2] = 5; tri[0][2] = 3;
tri[3][0] = 4; tri[3][1] = 6; tri[3][2] = 7; tri[1][0] = 0;
tri[1][1] = 3;
tri[1][2] = 2;
tri[2][0] = 4;
tri[2][1] = 7;
tri[2][2] = 5;
tri[3][0] = 4;
tri[3][1] = 6;
tri[3][2] = 7;
tri[4][0] = 0; tri[4][1] = 4; tri[4][2] = 5; tri[4][0] = 0;
tri[5][0] = 0; tri[5][1] = 5; tri[5][2] = 1; tri[4][1] = 4;
tri[6][0] = 2; tri[6][1] = 7; tri[6][2] = 6; tri[4][2] = 5;
tri[7][0] = 2; tri[7][1] = 3; tri[7][2] = 7; tri[5][0] = 0;
tri[5][1] = 5;
tri[5][2] = 1;
tri[6][0] = 2;
tri[6][1] = 7;
tri[6][2] = 6;
tri[7][0] = 2;
tri[7][1] = 3;
tri[7][2] = 7;
tri[8][0] = 2; tri[8][1] = 6; tri[8][2] = 4; tri[8][0] = 2;
tri[9][0] = 2; tri[9][1] = 4; tri[9][2] = 0; tri[8][1] = 6;
tri[10][0] = 1; tri[10][1] = 5; tri[10][2] = 7; tri[8][2] = 4;
tri[11][0] = 1; tri[11][1] = 7; tri[11][2] = 3; tri[9][0] = 2;
tri[9][1] = 4;
tri[9][2] = 0;
tri[10][0] = 1;
tri[10][1] = 5;
tri[10][2] = 7;
tri[11][0] = 1;
tri[11][1] = 7;
tri[11][2] = 3;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegPrism::~RegPrism() RegPrism::~RegPrism()
{ {
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -258,12 +301,11 @@ RegPrism::~RegPrism()
int RegPrism::inside(double x, double y, double z) int RegPrism::inside(double x, double y, double z)
{ {
double a = hinv[0][0]*(x-xlo) + hinv[0][1]*(y-ylo) + hinv[0][2]*(z-zlo); double a = hinv[0][0] * (x - xlo) + hinv[0][1] * (y - ylo) + hinv[0][2] * (z - zlo);
double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo); double b = hinv[1][1] * (y - ylo) + hinv[1][2] * (z - zlo);
double c = hinv[2][2]*(z-zlo); double c = hinv[2][2] * (z - zlo);
if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) return 1;
return 1;
return 0; return 0;
} }
@ -283,10 +325,12 @@ int RegPrism::surface_interior(double *x, double cutoff)
// x is exterior to prism // x is exterior to prism
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot < 0.0) return 0; if (dot < 0.0) return 0;
} }
@ -296,15 +340,17 @@ int RegPrism::surface_interior(double *x, double cutoff)
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (open_faces[i]) continue; if (open_faces[i]) continue;
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot < cutoff) { if (dot < cutoff) {
contact[n].r = dot; contact[n].r = dot;
contact[n].delx = dot*face[i][0]; contact[n].delx = dot * face[i][0];
contact[n].dely = dot*face[i][1]; contact[n].dely = dot * face[i][1];
contact[n].delz = dot*face[i][2]; contact[n].delz = dot * face[i][2];
contact[n].radius = 0; contact[n].radius = 0;
contact[n].iwall = i; contact[n].iwall = i;
n++; n++;
@ -325,25 +371,29 @@ int RegPrism::surface_exterior(double *x, double cutoff)
int i; int i;
double dot; double dot;
double *corner; double *corner;
double xp,yp,zp; double xp, yp, zp;
// x is far enough from prism that there is no contact // x is far enough from prism that there is no contact
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot <= -cutoff) return 0; if (dot <= -cutoff) return 0;
} }
// x is interior to prism // x is interior to prism
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot <= 0.0) break; if (dot <= 0.0) break;
} }
@ -354,8 +404,8 @@ int RegPrism::surface_exterior(double *x, double cutoff)
// could be edge or corner pt of prism // could be edge or corner pt of prism
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
find_nearest(x,xp,yp,zp); find_nearest(x, xp, yp, zp);
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = 0; contact[0].radius = 0;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -369,8 +419,8 @@ int RegPrism::surface_exterior(double *x, double cutoff)
void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp) void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
{ {
int i,j,k,iface; int i, j, k, iface;
double xproj[3],xline[3],nearest[3]; double xproj[3], xline[3], nearest[3];
double dot; double dot;
// generate successive xnear points, one nearest to x is (xp,yp,zp) // generate successive xnear points, one nearest to x is (xp,yp,zp)
@ -384,27 +434,25 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
double distsq = BIG; double distsq = BIG;
for (int itri = 0; itri < 12; itri++) { for (int itri = 0; itri < 12; itri++) {
iface = itri/2; iface = itri / 2;
if (open_faces[iface]) continue; if (open_faces[iface]) continue;
i = tri[itri][0]; i = tri[itri][0];
j = tri[itri][1]; j = tri[itri][1];
k = tri[itri][2]; k = tri[itri][2];
dot = (x[0]-corners[i][0])*face[iface][0] + dot = (x[0] - corners[i][0]) * face[iface][0] + (x[1] - corners[i][1]) * face[iface][1] +
(x[1]-corners[i][1])*face[iface][1] + (x[2] - corners[i][2]) * face[iface][2];
(x[2]-corners[i][2])*face[iface][2]; xproj[0] = x[0] - dot * face[iface][0];
xproj[0] = x[0] - dot*face[iface][0]; xproj[1] = x[1] - dot * face[iface][1];
xproj[1] = x[1] - dot*face[iface][1]; xproj[2] = x[2] - dot * face[iface][2];
xproj[2] = x[2] - dot*face[iface][2]; if (inside_tri(xproj, corners[i], corners[j], corners[k], face[iface])) {
if (inside_tri(xproj,corners[i],corners[j],corners[k],face[iface])) { distsq = closest(x, xproj, nearest, distsq);
distsq = closest(x,xproj,nearest,distsq); } else {
} point_on_line_segment(corners[i], corners[j], xproj, xline);
else { distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[i],corners[j],xproj,xline); point_on_line_segment(corners[j], corners[k], xproj, xline);
distsq = closest(x,xline,nearest,distsq); distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[j],corners[k],xproj,xline); point_on_line_segment(corners[i], corners[k], xproj, xline);
distsq = closest(x,xline,nearest,distsq); distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[i],corners[k],xproj,xline);
distsq = closest(x,xline,nearest,distsq);
} }
} }
@ -422,25 +470,24 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, double *norm)
double *norm)
{ {
double edge[3],pvec[3],xproduct[3]; double edge[3], pvec[3], xproduct[3];
MathExtra::sub3(v2,v1,edge); MathExtra::sub3(v2, v1, edge);
MathExtra::sub3(x,v1,pvec); MathExtra::sub3(x, v1, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
MathExtra::sub3(v3,v2,edge); MathExtra::sub3(v3, v2, edge);
MathExtra::sub3(x,v2,pvec); MathExtra::sub3(x, v2, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
MathExtra::sub3(v1,v3,edge); MathExtra::sub3(v1, v3, edge);
MathExtra::sub3(x,v3,pvec); MathExtra::sub3(x, v3, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
return 1; return 1;
} }
@ -452,7 +499,7 @@ double RegPrism::closest(double *x, double *near, double *nearest, double dsq)
double delx = x[0] - near[0]; double delx = x[0] - near[0];
double dely = x[1] - near[1]; double dely = x[1] - near[1];
double delz = x[2] - near[2]; double delz = x[2] - near[2];
double rsq = delx*delx + dely*dely + delz*delz; double rsq = delx * delx + dely * dely + delz * delz;
if (rsq >= dsq) return dsq; if (rsq >= dsq) return dsq;
nearest[0] = near[0]; nearest[0] = near[0];

119
src/region_sphere.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -23,52 +22,52 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE}; enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr) Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr)
{ {
options(narg-6,&arg[6]); options(narg - 6, &arg[6]);
if (utils::strmatch(arg[2],"^v_")) { if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2]+2); xstr = utils::strdup(arg[2] + 2);
xc = 0.0; xc = 0.0;
xstyle = VARIABLE; xstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
xc = xscale*utils::numeric(FLERR,arg[2],false,lmp); xc = xscale * utils::numeric(FLERR, arg[2], false, lmp);
xstyle = CONSTANT; xstyle = CONSTANT;
} }
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
ystr = utils::strdup(arg[3]+2); ystr = utils::strdup(arg[3] + 2);
yc = 0.0; yc = 0.0;
ystyle = VARIABLE; ystyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
yc = yscale*utils::numeric(FLERR,arg[3],false,lmp); yc = yscale * utils::numeric(FLERR, arg[3], false, lmp);
ystyle = CONSTANT; ystyle = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
zstr = utils::strdup(arg[4]+2); zstr = utils::strdup(arg[4] + 2);
zc = 0.0; zc = 0.0;
zstyle = VARIABLE; zstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
zc = zscale*utils::numeric(FLERR,arg[4],false,lmp); zc = zscale * utils::numeric(FLERR, arg[4], false, lmp);
zstyle = CONSTANT; zstyle = CONSTANT;
} }
if (utils::strmatch(arg[5],"^v_")) { if (utils::strmatch(arg[5], "^v_")) {
rstr = utils::strdup(arg[5]+2); rstr = utils::strdup(arg[5] + 2);
radius = 0.0; radius = 0.0;
rstyle = VARIABLE; rstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
radius = xscale*utils::numeric(FLERR,arg[5],false,lmp); radius = xscale * utils::numeric(FLERR, arg[5], false, lmp);
rstyle = CONSTANT; rstyle = CONSTANT;
} }
@ -79,7 +78,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
// error check // error check
if (radius < 0.0) error->all(FLERR,"Illegal region sphere radius: {}", radius); if (radius < 0.0) error->all(FLERR, "Illegal region sphere radius: {}", radius);
// extent of sphere // extent of sphere
// for variable radius, uses initial radius and origin for variable center // for variable radius, uses initial radius and origin for variable center
@ -92,7 +91,8 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
extent_yhi = yc + radius; extent_yhi = yc + radius;
extent_zlo = zc - radius; extent_zlo = zc - radius;
extent_zhi = zc + radius; extent_zhi = zc + radius;
} else bboxflag = 0; } else
bboxflag = 0;
cmax = 1; cmax = 1;
contact = new Contact[cmax]; contact = new Contact[cmax];
@ -103,11 +103,11 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
RegSphere::~RegSphere() RegSphere::~RegSphere()
{ {
delete [] xstr; delete[] xstr;
delete [] ystr; delete[] ystr;
delete [] zstr; delete[] zstr;
delete [] rstr; delete[] rstr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -128,7 +128,7 @@ int RegSphere::inside(double x, double y, double z)
double delx = x - xc; double delx = x - xc;
double dely = y - yc; double dely = y - yc;
double delz = z - zc; double delz = z - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r <= radius) return 1; if (r <= radius) return 1;
return 0; return 0;
@ -146,15 +146,15 @@ int RegSphere::surface_interior(double *x, double cutoff)
double delx = x[0] - xc; double delx = x[0] - xc;
double dely = x[1] - yc; double dely = x[1] - yc;
double delz = x[2] - zc; double delz = x[2] - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r > radius || r == 0.0) return 0; if (r > radius || r == 0.0) return 0;
double delta = radius - r; double delta = radius - r;
if (delta < cutoff) { if (delta < cutoff) {
contact[0].r = delta; contact[0].r = delta;
contact[0].delx = delx*(1.0-radius/r); contact[0].delx = delx * (1.0 - radius / r);
contact[0].dely = dely*(1.0-radius/r); contact[0].dely = dely * (1.0 - radius / r);
contact[0].delz = delz*(1.0-radius/r); contact[0].delz = delz * (1.0 - radius / r);
contact[0].radius = -radius; contact[0].radius = -radius;
contact[0].iwall = 0; contact[0].iwall = 0;
contact[0].varflag = 1; contact[0].varflag = 1;
@ -174,15 +174,15 @@ int RegSphere::surface_exterior(double *x, double cutoff)
double delx = x[0] - xc; double delx = x[0] - xc;
double dely = x[1] - yc; double dely = x[1] - yc;
double delz = x[2] - zc; double delz = x[2] - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r < radius) return 0; if (r < radius) return 0;
double delta = r - radius; double delta = r - radius;
if (delta < cutoff) { if (delta < cutoff) {
contact[0].r = delta; contact[0].r = delta;
contact[0].delx = delx*(1.0-radius/r); contact[0].delx = delx * (1.0 - radius / r);
contact[0].dely = dely*(1.0-radius/r); contact[0].dely = dely * (1.0 - radius / r);
contact[0].delz = delz*(1.0-radius/r); contact[0].delz = delz * (1.0 - radius / r);
contact[0].radius = radius; contact[0].radius = radius;
contact[0].iwall = 0; contact[0].iwall = 0;
contact[0].varflag = 1; contact[0].varflag = 1;
@ -197,19 +197,15 @@ int RegSphere::surface_exterior(double *x, double cutoff)
void RegSphere::shape_update() void RegSphere::shape_update()
{ {
if (xstyle == VARIABLE) if (xstyle == VARIABLE) xc = xscale * input->variable->compute_equal(xvar);
xc = xscale * input->variable->compute_equal(xvar);
if (ystyle == VARIABLE) if (ystyle == VARIABLE) yc = yscale * input->variable->compute_equal(yvar);
yc = yscale * input->variable->compute_equal(yvar);
if (zstyle == VARIABLE) if (zstyle == VARIABLE) zc = zscale * input->variable->compute_equal(zvar);
zc = zscale * input->variable->compute_equal(zvar);
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
radius = xscale * input->variable->compute_equal(rvar); radius = xscale * input->variable->compute_equal(rvar);
if (radius < 0.0) if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
error->one(FLERR,"Variable evaluation in region gave bad value");
} }
} }
@ -221,34 +217,30 @@ void RegSphere::variable_check()
{ {
if (xstyle == VARIABLE) { if (xstyle == VARIABLE) {
xvar = input->variable->find(xstr); xvar = input->variable->find(xstr);
if (xvar < 0) if (xvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", xstr);
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(xvar)) if (!input->variable->equalstyle(xvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable {} for region sphere is invalid style", xstr);
} }
if (ystyle == VARIABLE) { if (ystyle == VARIABLE) {
yvar = input->variable->find(ystr); yvar = input->variable->find(ystr);
if (yvar < 0) if (yvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", ystr);
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(yvar)) if (!input->variable->equalstyle(yvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable {} for region sphere is invalid style", ystr);
} }
if (zstyle == VARIABLE) { if (zstyle == VARIABLE) {
zvar = input->variable->find(zstr); zvar = input->variable->find(zstr);
if (zvar < 0) if (zvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", zstr);
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(zvar)) if (!input->variable->equalstyle(zvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable {} for region sphere is invalid style", zstr);
} }
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
rvar = input->variable->find(rstr); rvar = input->variable->find(rstr);
if (rvar < 0) if (rvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", rstr);
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(rvar)) if (!input->variable->equalstyle(rvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable {} for region sphere is invalid style", rstr);
} }
} }
@ -265,27 +257,26 @@ void RegSphere::set_velocity_shape()
xcenter[1] = yc; xcenter[1] = yc;
xcenter[2] = zc; xcenter[2] = zc;
forward_transform(xcenter[0], xcenter[1], xcenter[2]); forward_transform(xcenter[0], xcenter[1], xcenter[2]);
if (update->ntimestep > 0) rprev = prev[4]; if (update->ntimestep > 0)
else rprev = radius; rprev = prev[4];
else
rprev = radius;
prev[4] = radius; prev[4] = radius;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add velocity due to shape change to wall velocity add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegSphere::velocity_contact_shape(double *vwall, double *xc) void RegSphere::velocity_contact_shape(double *vwall, double *xc)
{ {
double delx, dely, delz; // Displacement of contact point in x,y,z double delx, dely, delz; // Displacement of contact point in x,y,z
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
vwall[0] += delx/update->dt; vwall[0] += delx / update->dt;
vwall[1] += dely/update->dt; vwall[1] += dely / update->dt;
vwall[2] += delz/update->dt; vwall[2] += delz / update->dt;
} }