diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b21b7f68..a7ca19370 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -185,6 +185,7 @@ target_link_libraries(mpm git) if(MPM_BUILD_TESTING) SET(test_src ${mpm_SOURCE_DIR}/tests/test_main.cc + ${mpm_SOURCE_DIR}/tests/acceleration_constraint_test.cc ${mpm_SOURCE_DIR}/tests/cell_test.cc ${mpm_SOURCE_DIR}/tests/cell_vector_test.cc ${mpm_SOURCE_DIR}/tests/contact_test.cc diff --git a/external/csv/csv.h b/external/csv/csv.h new file mode 100644 index 000000000..1f697e17b --- /dev/null +++ b/external/csv/csv.h @@ -0,0 +1,1272 @@ +// Copyright: (2012-2015) Ben Strasser +// License: BSD-3 +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#ifndef CSV_H +#define CSV_H + +#include +#include +#include +#include +#include +#include +#include +#ifndef CSV_IO_NO_THREAD +#include +#include +#include +#endif +#include +#include +#include +#include + +namespace io{ + //////////////////////////////////////////////////////////////////////////// + // LineReader // + //////////////////////////////////////////////////////////////////////////// + + namespace error{ + struct base : std::exception{ + virtual void format_error_message()const = 0; + + const char*what()const noexcept override{ + format_error_message(); + return error_message_buffer; + } + + mutable char error_message_buffer[512]; + }; + + const int max_file_name_length = 255; + + struct with_file_name{ + with_file_name(){ + std::memset(file_name, 0, sizeof(file_name)); + } + + void set_file_name(const char*file_name){ + if(file_name != nullptr){ + // This call to strncpy has parenthesis around it + // to silence the GCC -Wstringop-truncation warning + (strncpy(this->file_name, file_name, sizeof(this->file_name))); + this->file_name[sizeof(this->file_name)-1] = '\0'; + }else{ + this->file_name[0] = '\0'; + } + } + + char file_name[max_file_name_length+1]; + }; + + struct with_file_line{ + with_file_line(){ + file_line = -1; + } + + void set_file_line(int file_line){ + this->file_line = file_line; + } + + int file_line; + }; + + struct with_errno{ + with_errno(){ + errno_value = 0; + } + + void set_errno(int errno_value){ + this->errno_value = errno_value; + } + + int errno_value; + }; + + struct can_not_open_file : + base, + with_file_name, + with_errno{ + void format_error_message()const override{ + if(errno_value != 0) + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Can not open file \"%s\" because \"%s\"." + , file_name, std::strerror(errno_value)); + else + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Can not open file \"%s\"." + , file_name); + } + }; + + struct line_length_limit_exceeded : + base, + with_file_name, + with_file_line{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Line number %d in file \"%s\" exceeds the maximum length of 2^24-1." + , file_line, file_name); + } + }; + } + + class ByteSourceBase{ + public: + virtual int read(char*buffer, int size)=0; + virtual ~ByteSourceBase(){} + }; + + namespace detail{ + + class OwningStdIOByteSourceBase : public ByteSourceBase{ + public: + explicit OwningStdIOByteSourceBase(FILE*file):file(file){ + // Tell the std library that we want to do the buffering ourself. + std::setvbuf(file, 0, _IONBF, 0); + } + + int read(char*buffer, int size){ + return std::fread(buffer, 1, size, file); + } + + ~OwningStdIOByteSourceBase(){ + std::fclose(file); + } + + private: + FILE*file; + }; + + class NonOwningIStreamByteSource : public ByteSourceBase{ + public: + explicit NonOwningIStreamByteSource(std::istream&in):in(in){} + + int read(char*buffer, int size){ + in.read(buffer, size); + return in.gcount(); + } + + ~NonOwningIStreamByteSource(){} + + private: + std::istream∈ + }; + + class NonOwningStringByteSource : public ByteSourceBase{ + public: + NonOwningStringByteSource(const char*str, long long size):str(str), remaining_byte_count(size){} + + int read(char*buffer, int desired_byte_count){ + int to_copy_byte_count = desired_byte_count; + if(remaining_byte_count < to_copy_byte_count) + to_copy_byte_count = remaining_byte_count; + std::memcpy(buffer, str, to_copy_byte_count); + remaining_byte_count -= to_copy_byte_count; + str += to_copy_byte_count; + return to_copy_byte_count; + } + + ~NonOwningStringByteSource(){} + + private: + const char*str; + long long remaining_byte_count; + }; + + #ifndef CSV_IO_NO_THREAD + class AsynchronousReader{ + public: + void init(std::unique_ptrarg_byte_source){ + std::unique_lockguard(lock); + byte_source = std::move(arg_byte_source); + desired_byte_count = -1; + termination_requested = false; + worker = std::thread( + [&]{ + std::unique_lockguard(lock); + try{ + for(;;){ + read_requested_condition.wait( + guard, + [&]{ + return desired_byte_count != -1 || termination_requested; + } + ); + if(termination_requested) + return; + + read_byte_count = byte_source->read(buffer, desired_byte_count); + desired_byte_count = -1; + if(read_byte_count == 0) + break; + read_finished_condition.notify_one(); + } + }catch(...){ + read_error = std::current_exception(); + } + read_finished_condition.notify_one(); + } + ); + } + + bool is_valid()const{ + return byte_source != nullptr; + } + + void start_read(char*arg_buffer, int arg_desired_byte_count){ + std::unique_lockguard(lock); + buffer = arg_buffer; + desired_byte_count = arg_desired_byte_count; + read_byte_count = -1; + read_requested_condition.notify_one(); + } + + int finish_read(){ + std::unique_lockguard(lock); + read_finished_condition.wait( + guard, + [&]{ + return read_byte_count != -1 || read_error; + } + ); + if(read_error) + std::rethrow_exception(read_error); + else + return read_byte_count; + } + + ~AsynchronousReader(){ + if(byte_source != nullptr){ + { + std::unique_lockguard(lock); + termination_requested = true; + } + read_requested_condition.notify_one(); + worker.join(); + } + } + + private: + std::unique_ptrbyte_source; + + std::thread worker; + + bool termination_requested; + std::exception_ptr read_error; + char*buffer; + int desired_byte_count; + int read_byte_count; + + std::mutex lock; + std::condition_variable read_finished_condition; + std::condition_variable read_requested_condition; + }; + #endif + + class SynchronousReader{ + public: + void init(std::unique_ptrarg_byte_source){ + byte_source = std::move(arg_byte_source); + } + + bool is_valid()const{ + return byte_source != nullptr; + } + + void start_read(char*arg_buffer, int arg_desired_byte_count){ + buffer = arg_buffer; + desired_byte_count = arg_desired_byte_count; + } + + int finish_read(){ + return byte_source->read(buffer, desired_byte_count); + } + private: + std::unique_ptrbyte_source; + char*buffer; + int desired_byte_count; + }; + } + + class LineReader{ + private: + static const int block_len = 1<<20; + std::unique_ptrbuffer; // must be constructed before (and thus destructed after) the reader! + #ifdef CSV_IO_NO_THREAD + detail::SynchronousReader reader; + #else + detail::AsynchronousReader reader; + #endif + int data_begin; + int data_end; + + char file_name[error::max_file_name_length+1]; + unsigned file_line; + + static std::unique_ptr open_file(const char*file_name){ + // We open the file in binary mode as it makes no difference under *nix + // and under Windows we handle \r\n newlines ourself. + FILE*file = std::fopen(file_name, "rb"); + if(file == 0){ + int x = errno; // store errno as soon as possible, doing it after constructor call can fail. + error::can_not_open_file err; + err.set_errno(x); + err.set_file_name(file_name); + throw err; + } + return std::unique_ptr(new detail::OwningStdIOByteSourceBase(file)); + } + + void init(std::unique_ptrbyte_source){ + file_line = 0; + + buffer = std::unique_ptr(new char[3*block_len]); + data_begin = 0; + data_end = byte_source->read(buffer.get(), 2*block_len); + + // Ignore UTF-8 BOM + if(data_end >= 3 && buffer[0] == '\xEF' && buffer[1] == '\xBB' && buffer[2] == '\xBF') + data_begin = 3; + + if(data_end == 2*block_len){ + reader.init(std::move(byte_source)); + reader.start_read(buffer.get() + 2*block_len, block_len); + } + } + + public: + LineReader() = delete; + LineReader(const LineReader&) = delete; + LineReader&operator=(const LineReader&) = delete; + + explicit LineReader(const char*file_name){ + set_file_name(file_name); + init(open_file(file_name)); + } + + explicit LineReader(const std::string&file_name){ + set_file_name(file_name.c_str()); + init(open_file(file_name.c_str())); + } + + LineReader(const char*file_name, std::unique_ptrbyte_source){ + set_file_name(file_name); + init(std::move(byte_source)); + } + + LineReader(const std::string&file_name, std::unique_ptrbyte_source){ + set_file_name(file_name.c_str()); + init(std::move(byte_source)); + } + + LineReader(const char*file_name, const char*data_begin, const char*data_end){ + set_file_name(file_name); + init(std::unique_ptr(new detail::NonOwningStringByteSource(data_begin, data_end-data_begin))); + } + + LineReader(const std::string&file_name, const char*data_begin, const char*data_end){ + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::NonOwningStringByteSource(data_begin, data_end-data_begin))); + } + + LineReader(const char*file_name, FILE*file){ + set_file_name(file_name); + init(std::unique_ptr(new detail::OwningStdIOByteSourceBase(file))); + } + + LineReader(const std::string&file_name, FILE*file){ + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::OwningStdIOByteSourceBase(file))); + } + + LineReader(const char*file_name, std::istream&in){ + set_file_name(file_name); + init(std::unique_ptr(new detail::NonOwningIStreamByteSource(in))); + } + + LineReader(const std::string&file_name, std::istream&in){ + set_file_name(file_name.c_str()); + init(std::unique_ptr(new detail::NonOwningIStreamByteSource(in))); + } + + void set_file_name(const std::string&file_name){ + set_file_name(file_name.c_str()); + } + + void set_file_name(const char*file_name){ + if(file_name != nullptr){ + strncpy(this->file_name, file_name, sizeof(this->file_name)); + this->file_name[sizeof(this->file_name)-1] = '\0'; + }else{ + this->file_name[0] = '\0'; + } + } + + const char*get_truncated_file_name()const{ + return file_name; + } + + void set_file_line(unsigned file_line){ + this->file_line = file_line; + } + + unsigned get_file_line()const{ + return file_line; + } + + char*next_line(){ + if(data_begin == data_end) + return nullptr; + + ++file_line; + + assert(data_begin < data_end); + assert(data_end <= block_len*2); + + if(data_begin >= block_len){ + std::memcpy(buffer.get(), buffer.get()+block_len, block_len); + data_begin -= block_len; + data_end -= block_len; + if(reader.is_valid()) + { + data_end += reader.finish_read(); + std::memcpy(buffer.get()+block_len, buffer.get()+2*block_len, block_len); + reader.start_read(buffer.get() + 2*block_len, block_len); + } + } + + int line_end = data_begin; + while(buffer[line_end] != '\n' && line_end != data_end){ + ++line_end; + } + + if(line_end - data_begin + 1 > block_len){ + error::line_length_limit_exceeded err; + err.set_file_name(file_name); + err.set_file_line(file_line); + throw err; + } + + if(buffer[line_end] == '\n' && line_end != data_end){ + buffer[line_end] = '\0'; + }else{ + // some files are missing the newline at the end of the + // last line + ++data_end; + buffer[line_end] = '\0'; + } + + // handle windows \r\n-line breaks + if(line_end != data_begin && buffer[line_end-1] == '\r') + buffer[line_end-1] = '\0'; + + char*ret = buffer.get() + data_begin; + data_begin = line_end+1; + return ret; + } + }; + + + //////////////////////////////////////////////////////////////////////////// + // CSV // + //////////////////////////////////////////////////////////////////////////// + + namespace error{ + const int max_column_name_length = 63; + struct with_column_name{ + with_column_name(){ + std::memset(column_name, 0, max_column_name_length+1); + } + + void set_column_name(const char*column_name){ + if(column_name != nullptr){ + std::strncpy(this->column_name, column_name, max_column_name_length); + this->column_name[max_column_name_length] = '\0'; + }else{ + this->column_name[0] = '\0'; + } + } + + char column_name[max_column_name_length+1]; + }; + + + const int max_column_content_length = 63; + + struct with_column_content{ + with_column_content(){ + std::memset(column_content, 0, max_column_content_length+1); + } + + void set_column_content(const char*column_content){ + if(column_content != nullptr){ + std::strncpy(this->column_content, column_content, max_column_content_length); + this->column_content[max_column_content_length] = '\0'; + }else{ + this->column_content[0] = '\0'; + } + } + + char column_content[max_column_content_length+1]; + }; + + + struct extra_column_in_header : + base, + with_file_name, + with_column_name{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Extra column "%s" in header of file "%s".)" + , column_name, file_name); + } + }; + + struct missing_column_in_header : + base, + with_file_name, + with_column_name{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Missing column "%s" in header of file "%s".)" + , column_name, file_name); + } + }; + + struct duplicated_column_in_header : + base, + with_file_name, + with_column_name{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(Duplicated column "%s" in header of file "%s".)" + , column_name, file_name); + } + }; + + struct header_missing : + base, + with_file_name{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Header missing in file \"%s\"." + , file_name); + } + }; + + struct too_few_columns : + base, + with_file_name, + with_file_line{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Too few columns in line %d in file \"%s\"." + , file_line, file_name); + } + }; + + struct too_many_columns : + base, + with_file_name, + with_file_line{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Too many columns in line %d in file \"%s\"." + , file_line, file_name); + } + }; + + struct escaped_string_not_closed : + base, + with_file_name, + with_file_line{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + "Escaped string was not closed in line %d in file \"%s\"." + , file_line, file_name); + } + }; + + struct integer_must_be_positive : + base, + with_file_name, + with_file_line, + with_column_name, + with_column_content{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" must be positive or 0 in column "%s" in file "%s" in line "%d".)" + , column_content, column_name, file_name, file_line); + } + }; + + struct no_digit : + base, + with_file_name, + with_file_line, + with_column_name, + with_column_content{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" contains an invalid digit in column "%s" in file "%s" in line "%d".)" + , column_content, column_name, file_name, file_line); + } + }; + + struct integer_overflow : + base, + with_file_name, + with_file_line, + with_column_name, + with_column_content{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" overflows in column "%s" in file "%s" in line "%d".)" + , column_content, column_name, file_name, file_line); + } + }; + + struct integer_underflow : + base, + with_file_name, + with_file_line, + with_column_name, + with_column_content{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The integer "%s" underflows in column "%s" in file "%s" in line "%d".)" + , column_content, column_name, file_name, file_line); + } + }; + + struct invalid_single_character : + base, + with_file_name, + with_file_line, + with_column_name, + with_column_content{ + void format_error_message()const override{ + std::snprintf(error_message_buffer, sizeof(error_message_buffer), + R"(The content "%s" of column "%s" in file "%s" in line "%d" is not a single character.)" + , column_content, column_name, file_name, file_line); + } + }; + } + + using ignore_column = unsigned int; + static const ignore_column ignore_no_column = 0; + static const ignore_column ignore_extra_column = 1; + static const ignore_column ignore_missing_column = 2; + + template + struct trim_chars{ + private: + constexpr static bool is_trim_char(char){ + return false; + } + + template + constexpr static bool is_trim_char(char c, char trim_char, OtherTrimChars...other_trim_chars){ + return c == trim_char || is_trim_char(c, other_trim_chars...); + } + + public: + static void trim(char*&str_begin, char*&str_end){ + while(str_begin != str_end && is_trim_char(*str_begin, trim_char_list...)) + ++str_begin; + while(str_begin != str_end && is_trim_char(*(str_end-1), trim_char_list...)) + --str_end; + *str_end = '\0'; + } + }; + + + struct no_comment{ + static bool is_comment(const char*){ + return false; + } + }; + + template + struct single_line_comment{ + private: + constexpr static bool is_comment_start_char(char){ + return false; + } + + template + constexpr static bool is_comment_start_char(char c, char comment_start_char, OtherCommentStartChars...other_comment_start_chars){ + return c == comment_start_char || is_comment_start_char(c, other_comment_start_chars...); + } + + public: + + static bool is_comment(const char*line){ + return is_comment_start_char(*line, comment_start_char_list...); + } + }; + + struct empty_line_comment{ + static bool is_comment(const char*line){ + if(*line == '\0') + return true; + while(*line == ' ' || *line == '\t'){ + ++line; + if(*line == 0) + return true; + } + return false; + } + }; + + template + struct single_and_empty_line_comment{ + static bool is_comment(const char*line){ + return single_line_comment::is_comment(line) || empty_line_comment::is_comment(line); + } + }; + + template + struct no_quote_escape{ + static const char*find_next_column_end(const char*col_begin){ + while(*col_begin != sep && *col_begin != '\0') + ++col_begin; + return col_begin; + } + + static void unescape(char*&, char*&){ + + } + }; + + template + struct double_quote_escape{ + static const char*find_next_column_end(const char*col_begin){ + while(*col_begin != sep && *col_begin != '\0') + if(*col_begin != quote) + ++col_begin; + else{ + do{ + ++col_begin; + while(*col_begin != quote){ + if(*col_begin == '\0') + throw error::escaped_string_not_closed(); + ++col_begin; + } + ++col_begin; + }while(*col_begin == quote); + } + return col_begin; + } + + static void unescape(char*&col_begin, char*&col_end){ + if(col_end - col_begin >= 2){ + if(*col_begin == quote && *(col_end-1) == quote){ + ++col_begin; + --col_end; + char*out = col_begin; + for(char*in = col_begin; in!=col_end; ++in){ + if(*in == quote && (in+1) != col_end && *(in+1) == quote){ + ++in; + } + *out = *in; + ++out; + } + col_end = out; + *col_end = '\0'; + } + } + + } + }; + + struct throw_on_overflow{ + template + static void on_overflow(T&){ + throw error::integer_overflow(); + } + + template + static void on_underflow(T&){ + throw error::integer_underflow(); + } + }; + + struct ignore_overflow{ + template + static void on_overflow(T&){} + + template + static void on_underflow(T&){} + }; + + struct set_to_max_on_overflow{ + template + static void on_overflow(T&x){ + // using (std::numeric_limits::max) instead of std::numeric_limits::max + // to make code including windows.h with its max macro happy + x = (std::numeric_limits::max)(); + } + + template + static void on_underflow(T&x){ + x = (std::numeric_limits::min)(); + } + }; + + + namespace detail{ + template + void chop_next_column( + char*&line, char*&col_begin, char*&col_end + ){ + assert(line != nullptr); + + col_begin = line; + // the col_begin + (... - col_begin) removes the constness + col_end = col_begin + (quote_policy::find_next_column_end(col_begin) - col_begin); + + if(*col_end == '\0'){ + line = nullptr; + }else{ + *col_end = '\0'; + line = col_end + 1; + } + } + + template + void parse_line( + char*line, + char**sorted_col, + const std::vector&col_order + ){ + for (int i : col_order) { + if(line == nullptr) + throw ::io::error::too_few_columns(); + char*col_begin, *col_end; + chop_next_column(line, col_begin, col_end); + + if (i != -1) { + trim_policy::trim(col_begin, col_end); + quote_policy::unescape(col_begin, col_end); + + sorted_col[i] = col_begin; + } + } + if(line != nullptr) + throw ::io::error::too_many_columns(); + } + + template + void parse_header_line( + char*line, + std::vector&col_order, + const std::string*col_name, + ignore_column ignore_policy + ){ + col_order.clear(); + + bool found[column_count]; + std::fill(found, found + column_count, false); + while(line){ + char*col_begin,*col_end; + chop_next_column(line, col_begin, col_end); + + trim_policy::trim(col_begin, col_end); + quote_policy::unescape(col_begin, col_end); + + for(unsigned i=0; i + void parse(char*col, char &x){ + if(!*col) + throw error::invalid_single_character(); + x = *col; + ++col; + if(*col) + throw error::invalid_single_character(); + } + + template + void parse(char*col, std::string&x){ + x = col; + } + + template + void parse(char*col, const char*&x){ + x = col; + } + + template + void parse(char*col, char*&x){ + x = col; + } + + template + void parse_unsigned_integer(const char*col, T&x){ + x = 0; + while(*col != '\0'){ + if('0' <= *col && *col <= '9'){ + T y = *col - '0'; + if(x > ((std::numeric_limits::max)()-y)/10){ + overflow_policy::on_overflow(x); + return; + } + x = 10*x+y; + }else + throw error::no_digit(); + ++col; + } + } + + templatevoid parse(char*col, unsigned char &x) + {parse_unsigned_integer(col, x);} + templatevoid parse(char*col, unsigned short &x) + {parse_unsigned_integer(col, x);} + templatevoid parse(char*col, unsigned int &x) + {parse_unsigned_integer(col, x);} + templatevoid parse(char*col, unsigned long &x) + {parse_unsigned_integer(col, x);} + templatevoid parse(char*col, unsigned long long &x) + {parse_unsigned_integer(col, x);} + + template + void parse_signed_integer(const char*col, T&x){ + if(*col == '-'){ + ++col; + + x = 0; + while(*col != '\0'){ + if('0' <= *col && *col <= '9'){ + T y = *col - '0'; + if(x < ((std::numeric_limits::min)()+y)/10){ + overflow_policy::on_underflow(x); + return; + } + x = 10*x-y; + }else + throw error::no_digit(); + ++col; + } + return; + }else if(*col == '+') + ++col; + parse_unsigned_integer(col, x); + } + + templatevoid parse(char*col, signed char &x) + {parse_signed_integer(col, x);} + templatevoid parse(char*col, signed short &x) + {parse_signed_integer(col, x);} + templatevoid parse(char*col, signed int &x) + {parse_signed_integer(col, x);} + templatevoid parse(char*col, signed long &x) + {parse_signed_integer(col, x);} + templatevoid parse(char*col, signed long long &x) + {parse_signed_integer(col, x);} + + template + void parse_float(const char*col, T&x){ + bool is_neg = false; + if(*col == '-'){ + is_neg = true; + ++col; + }else if(*col == '+') + ++col; + + x = 0; + while('0' <= *col && *col <= '9'){ + int y = *col - '0'; + x *= 10; + x += y; + ++col; + } + + if(*col == '.'|| *col == ','){ + ++col; + T pos = 1; + while('0' <= *col && *col <= '9'){ + pos /= 10; + int y = *col - '0'; + ++col; + x += y*pos; + } + } + + if(*col == 'e' || *col == 'E'){ + ++col; + int e; + + parse_signed_integer(col, e); + + if(e != 0){ + T base; + if(e < 0){ + base = T(0.1); + e = -e; + }else{ + base = T(10); + } + + while(e != 1){ + if((e & 1) == 0){ + base = base*base; + e >>= 1; + }else{ + x *= base; + --e; + } + } + x *= base; + } + }else{ + if(*col != '\0') + throw error::no_digit(); + } + + if(is_neg) + x = -x; + } + + template void parse(char*col, float&x) { parse_float(col, x); } + template void parse(char*col, double&x) { parse_float(col, x); } + template void parse(char*col, long double&x) { parse_float(col, x); } + + template + void parse(char*col, T&x){ + // Mute unused variable compiler warning + (void)col; + (void)x; + // GCC evalutes "false" when reading the template and + // "sizeof(T)!=sizeof(T)" only when instantiating it. This is why + // this strange construct is used. + static_assert(sizeof(T)!=sizeof(T), + "Can not parse this type. Only buildin integrals, floats, char, char*, const char* and std::string are supported"); + } + + } + + template, + class quote_policy = no_quote_escape<','>, + class overflow_policy = throw_on_overflow, + class comment_policy = no_comment + > + class CSVReader{ + private: + LineReader in; + + char*row[column_count]; + std::string column_names[column_count]; + + std::vectorcol_order; + + template + void set_column_names(std::string s, ColNames...cols){ + column_names[column_count-sizeof...(ColNames)-1] = std::move(s); + set_column_names(std::forward(cols)...); + } + + void set_column_names(){} + + + public: + CSVReader() = delete; + CSVReader(const CSVReader&) = delete; + CSVReader&operator=(const CSVReader&); + + template + explicit CSVReader(Args&&...args):in(std::forward(args)...){ + std::fill(row, row+column_count, nullptr); + col_order.resize(column_count); + for(unsigned i=0; i + void read_header(ignore_column ignore_policy, ColNames...cols){ + static_assert(sizeof...(ColNames)>=column_count, "not enough column names specified"); + static_assert(sizeof...(ColNames)<=column_count, "too many column names specified"); + try{ + set_column_names(std::forward(cols)...); + + char*line; + do{ + line = in.next_line(); + if(!line) + throw error::header_missing(); + }while(comment_policy::is_comment(line)); + + detail::parse_header_line + + (line, col_order, column_names, ignore_policy); + }catch(error::with_file_name&err){ + err.set_file_name(in.get_truncated_file_name()); + throw; + } + } + + template + void set_header(ColNames...cols){ + static_assert(sizeof...(ColNames)>=column_count, + "not enough column names specified"); + static_assert(sizeof...(ColNames)<=column_count, + "too many column names specified"); + set_column_names(std::forward(cols)...); + std::fill(row, row+column_count, nullptr); + col_order.resize(column_count); + for(unsigned i=0; i + void parse_helper(std::size_t r, T&t, ColType&...cols){ + if(row[r]){ + try{ + try{ + ::io::detail::parse(row[r], t); + }catch(error::with_column_content&err){ + err.set_column_content(row[r]); + throw; + } + }catch(error::with_column_name&err){ + err.set_column_name(column_names[r].c_str()); + throw; + } + } + parse_helper(r+1, cols...); + } + + + public: + template + bool read_row(ColType& ...cols){ + static_assert(sizeof...(ColType)>=column_count, + "not enough columns specified"); + static_assert(sizeof...(ColType)<=column_count, + "too many columns specified"); + try{ + try{ + + char*line; + do{ + line = in.next_line(); + if(!line) + return false; + }while(comment_policy::is_comment(line)); + + detail::parse_line + (line, row, col_order); + + parse_helper(0, cols...); + }catch(error::with_file_name&err){ + err.set_file_name(in.get_truncated_file_name()); + throw; + } + }catch(error::with_file_line&err){ + err.set_file_line(in.get_file_line()); + throw; + } + + return true; + } + }; +} +#endif + diff --git a/external/csv_license.md b/external/csv_license.md new file mode 100644 index 000000000..82d53e1ec --- /dev/null +++ b/external/csv_license.md @@ -0,0 +1,27 @@ +Copyright (c) 2015, ben-strasser +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +* Neither the name of fast-cpp-csv-parser nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h index 4150af5b0..183d51c91 100644 --- a/include/io/io_mesh.h +++ b/include/io/io_mesh.h @@ -12,6 +12,7 @@ #include "Eigen/Dense" #include +#include "csv/csv.h" #include "spdlog/spdlog.h" #include "data_types.h" @@ -87,22 +88,33 @@ class IOMesh { const std::vector>& particles_cells) = 0; //! Read velocity constraints file - //! \param[in] velocity_constraints_files file name with constraints + //! \param[in] velocity_constraints_file file name with constraints virtual std::vector> read_velocity_constraints( const std::string& velocity_constraints_file) = 0; + //! Read acceleration constraints file + //! \param[in] acceleration_constraints_file file name with constraints + virtual std::vector> + read_acceleration_constraints( + const std::string& acceleration_constraints_file) = 0; + //! Read friction constraints file - //! \param[in] friction_constraints_files file name with frictions + //! \param[in] friction_constraints_file file name with frictions virtual std::vector> read_friction_constraints( const std::string& friction_constraints_file) = 0; //! Read forces file - //! \param[in] forces_files file name with nodal concentrated force + //! \param[in] forces_file file name with nodal concentrated force virtual std::vector> read_forces( const std::string& forces_file) = 0; + //! Read math function file + //! \param[in] function_file file name with linear math function entries + virtual std::array, 2> read_math_functions( + const std::string& math_file) = 0; + }; // IOMesh class } // namespace mpm diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h index 433661546..ce82674f7 100644 --- a/include/io/io_mesh_ascii.h +++ b/include/io/io_mesh_ascii.h @@ -74,14 +74,20 @@ class IOMeshAscii : public IOMesh { const std::string& particles_cells_file, const std::vector>& particles_cells) override; - //! Read constraints file - //! \param[in] velocity_constraints_files file name with constraints + //! Read velocity constraints file + //! \param[in] velocity_constraints_file file name with constraints std::vector> read_velocity_constraints( const std::string& velocity_constraints_file) override; + //! Read acceleration constraints file + //! \param[in] acceleration_constraints_file file name with constraints + std::vector> + read_acceleration_constraints( + const std::string& acceleration_constraints_file) override; + //! Read friction constraints file - //! \param[in] friction_constraints_files file name with frictions + //! \param[in] friction_constraints_file file name with frictions std::vector> read_friction_constraints( const std::string& friction_constraints_file) override; @@ -91,6 +97,11 @@ class IOMeshAscii : public IOMesh { std::vector> read_forces( const std::string& forces_file) override; + //! Read math function file + //! \param[in] function_file file name with linear math function entries + std::array, 2> read_math_functions( + const std::string& math_file) override; + private: //! Logger std::shared_ptr console_; diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc index 0f09fdc1a..485868896 100644 --- a/include/io/io_mesh_ascii.tcc +++ b/include/io/io_mesh_ascii.tcc @@ -270,12 +270,12 @@ std::map> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID and read stream + mpm::Index id; + // Angles and ream stream + Eigen::Matrix angles; while (istream.good()) { - // ID and read stream - mpm::Index id; istream >> id; - // Angles and ream stream - Eigen::Matrix angles; for (unsigned i = 0; i < Tdim; ++i) istream >> angles[i]; euler_angles.emplace(std::make_pair(id, angles)); } @@ -314,11 +314,11 @@ std::vector> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Volume + double volume; while (istream.good()) { - // ID - mpm::Index id; - // Volume - double volume; // Read stream istream >> id >> volume; volumes.emplace_back(std::make_tuple(id, volume)); @@ -358,9 +358,9 @@ std::vector> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index pid, cid; while (istream.good()) { - // ID - mpm::Index pid, cid; // Read stream istream >> pid >> cid; particles_cells.emplace_back(std::array({pid, cid})); @@ -416,13 +416,13 @@ std::vector> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Direction + unsigned dir; + // Velocity + double velocity; while (istream.good()) { - // ID - mpm::Index id; - // Direction - unsigned dir; - // Velocity - double velocity; // Read stream istream >> id >> dir >> velocity; constraints.emplace_back(std::make_tuple(id, dir, velocity)); @@ -438,6 +438,52 @@ std::vector> return constraints; } +//! Return acceleration constraints of nodes or particles +template +std::vector> + mpm::IOMeshAscii::read_acceleration_constraints( + const std::string& acceleration_constraints_file) { + + // Nodal or particle acceleration constraints + std::vector> constraints; + constraints.clear(); + + // input file stream + std::fstream file; + file.open(acceleration_constraints_file.c_str(), std::ios::in); + + try { + if (file.is_open() && file.good()) { + // Line + std::string line; + while (std::getline(file, line)) { + boost::algorithm::trim(line); + std::istringstream istream(line); + // ignore comment lines (# or !) or blank lines + if ((line.find('#') == std::string::npos) && + (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Direction + unsigned dir; + // Velocity + double acceleration; + while (istream.good()) { + // Read stream + istream >> id >> dir >> acceleration; + constraints.emplace_back(std::make_tuple(id, dir, acceleration)); + } + } + } + } + file.close(); + } catch (std::exception& exception) { + console_->error("Read acceleration constraints: {}", exception.what()); + file.close(); + } + return constraints; +} + //! Return friction constraints of particles template std::vector> @@ -462,15 +508,15 @@ std::vector> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Direction + unsigned dir; + // Sign + int sign; + // Friction + double friction; while (istream.good()) { - // ID - mpm::Index id; - // Direction - unsigned dir; - // Sign - int sign; - // Friction - double friction; // Read stream istream >> id >> dir >> sign >> friction; constraints.emplace_back(std::make_tuple(id, dir, sign, friction)); @@ -509,13 +555,13 @@ std::vector> // ignore comment lines (# or !) or blank lines if ((line.find('#') == std::string::npos) && (line.find('!') == std::string::npos) && (line != "")) { + // ID + mpm::Index id; + // Direction + unsigned dir; + // Force + double force; while (istream.good()) { - // ID - mpm::Index id; - // Direction - unsigned dir; - // Force - double force; // Read stream istream >> id >> dir >> force; forces.emplace_back(std::make_tuple(id, dir, force)); @@ -530,3 +576,25 @@ std::vector> } return forces; } + +// Return array with math function entries +template +std::array, 2> mpm::IOMeshAscii::read_math_functions( + const std::string& math_file) { + // Initialise vector with 2 empty vectors + std::array, 2> xfx_values; + + // Read from csv file + try { + io::CSVReader<2> in(math_file.c_str()); + double x_value, fx_value; + while (in.read_row(x_value, fx_value)) { + xfx_values[0].push_back(x_value); + xfx_values[1].push_back(fx_value); + } + } catch (std::exception& exception) { + console_->error("Read math functions: {}", exception.what()); + } + + return xfx_values; +} diff --git a/include/loads_bcs/acceleration_constraint.h b/include/loads_bcs/acceleration_constraint.h new file mode 100644 index 000000000..a1c4fbf8c --- /dev/null +++ b/include/loads_bcs/acceleration_constraint.h @@ -0,0 +1,56 @@ +#ifndef MPM_ACCELERATION_CONSTRAINT_H_ +#define MPM_ACCELERATION_CONSTRAINT_H_ + +//! Alias for JSON +#include "json.hpp" +using Json = nlohmann::json; + +#include "function_base.h" + +namespace mpm { + +//! AccelerationConstraint class to store acceleration constraint on a set +//! \brief AccelerationConstraint class to store a constraint on a set +//! \details AccelerationConstraint stores the constraint as a static value +class AccelerationConstraint { + public: + // Constructor + //! \param[in] setid set id + //! \param[in] acceleration_fn Math function if defined + //! \param[in] dir Direction of constraint load + //! \param[in] acceleration Constraint acceleration + AccelerationConstraint( + int setid, const std::shared_ptr& acceleration_fn, + unsigned dir, double acceleration) + : setid_{setid}, + acceleration_fn_{acceleration_fn}, + dir_{dir}, + acceleration_{acceleration} {}; + + // Set id + int setid() const { return setid_; } + + // Direction + unsigned dir() const { return dir_; } + + // Return acceleration + double acceleration(double current_time) const { + // Static load when no math function is defined + double scalar = (this->acceleration_fn_ != nullptr) + ? (this->acceleration_fn_)->value(current_time) + : 1.0; + return acceleration_ * scalar; + } + + private: + // ID + int setid_; + // Math function + std::shared_ptr acceleration_fn_; + // Direction + unsigned dir_; + // Acceleration + double acceleration_; +}; +} // namespace mpm +#endif // MPM_ACCELERATION_CONSTRAINT_H_ diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 0efdef6f6..5354b976e 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -3,6 +3,7 @@ #include +#include "acceleration_constraint.h" #include "friction_constraint.h" #include "logger.h" #include "mesh.h" @@ -34,6 +35,21 @@ class Constraints { const std::vector>& velocity_constraints); + //! Assign nodal acceleration constraints + //! \param[in] setid Node set id + //! \param[in] acceleration_constraints Accelerartion constraint at node, dir, + //! acceleration + bool assign_nodal_acceleration_constraint( + int set_id, + const std::shared_ptr& constraint); + + //! Assign acceleartion constraints to nodes + //! \param[in] acceleration_constraints Constraint at node, dir, and + //! acceleration + bool assign_nodal_acceleration_constraints( + const std::vector>& + acceleration_constraints); + //! Assign nodal frictional constraints //! \param[in] setid Node set id //! \param[in] friction_constraints Constraint at node, dir, sign, friction diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 5813e609f..c12d554d2 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -51,6 +51,60 @@ bool mpm::Constraints::assign_nodal_velocity_constraints( return status; } +//! Assign nodal acceleration constraints +template +bool mpm::Constraints::assign_nodal_acceleration_constraint( + int set_id, + const std::shared_ptr& constraint) { + bool status = true; + try { + int set_id = constraint->setid(); + auto nset = mesh_->nodes(set_id); + if (nset.size() == 0) + throw std::runtime_error( + "Node set is empty for assignment of acceleration constraints"); + + unsigned dir = constraint->dir(); + double acceleration = constraint->acceleration(0); + for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { + if (!(*nitr)->assign_acceleration_constraint(dir, acceleration)) + throw std::runtime_error( + "Failed to initialise acceleration constraint at node"); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assign acceleration constraints to nodes +template +bool mpm::Constraints::assign_nodal_acceleration_constraints( + const std::vector>& + acceleration_constraints) { + bool status = true; + try { + for (const auto& acceleration_constraint : acceleration_constraints) { + // Node id + mpm::Index nid = std::get<0>(acceleration_constraint); + // Direction + unsigned dir = std::get<1>(acceleration_constraint); + // Acceleration + double acceleration = std::get<2>(acceleration_constraint); + + // Apply constraint + if (!mesh_->node(nid)->assign_acceleration_constraint(dir, acceleration)) + throw std::runtime_error( + "Nodal acceleration constraints assignment failed"); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + //! Assign friction constraints to nodes template bool mpm::Constraints::assign_nodal_frictional_constraint( @@ -61,14 +115,14 @@ bool mpm::Constraints::assign_nodal_frictional_constraint( auto nset = mesh_->nodes(set_id); if (nset.size() == 0) throw std::runtime_error( - "Node set is empty for assignment of velocity constraints"); + "Node set is empty for assignment of friction constraints"); unsigned dir = fconstraint->dir(); int nsign_n = fconstraint->sign_n(); double friction = fconstraint->friction(); for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { if (!(*nitr)->assign_friction_constraint(dir, nsign_n, friction)) throw std::runtime_error( - "Failed to initialise velocity constraint at node"); + "Failed to initialise friction constraint at node"); } } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); diff --git a/include/mesh.h b/include/mesh.h index bb51f161a..40735b0ab 100644 --- a/include/mesh.h +++ b/include/mesh.h @@ -24,6 +24,7 @@ #include "json.hpp" using Json = nlohmann::json; +#include "acceleration_constraint.h" #include "cell.h" #include "factory.h" #include "friction_constraint.h" @@ -114,6 +115,12 @@ class Mesh { template void iterate_over_nodes_predicate(Toper oper, Tpred pred); + //! Iterate over node set + //! \tparam Toper Callable object typically a baseclass functor + //! \param[in] set_id particle set id + template + void iterate_over_node_set(int set_id, Toper oper); + //! Return a vector of nodes //! \param[in] set_id Set of id of nodes (-1 for all nodes) Vector> nodes(unsigned set_id) const { @@ -308,14 +315,26 @@ class Mesh { //! \param[in] current_time Current time void apply_traction_on_particles(double current_time); - //! Create particle velocity constraints tractions + //! Create nodal acceleration constraints //! \param[in] setid Node set id + //! \param[in] constraint Acceleration constraint + bool create_nodal_acceleration_constraint( + int set_id, + const std::shared_ptr& constraint); + + //! Create particle velocity constraints + //! \param[in] setid Node set id + //! \param[in] constraint Velocity constraint bool create_particle_velocity_constraint( int set_id, const std::shared_ptr& constraint); //! Apply particles velocity constraints void apply_particle_velocity_constraints(); + //! Update nodal acceleration constraints + //! \param[in] current_time Current time + void update_nodal_acceleration_constraints(double current_time); + //! Assign nodal concentrated force //! \param[in] nodal_forces Force at dir on nodes bool assign_nodal_concentrated_forces( @@ -519,6 +538,9 @@ class Mesh { std::map>> materials_; //! Loading (Particle tractions) std::vector> particle_tractions_; + //! Nodal acceleration constraints + std::vector> + nodal_acceleration_constraints_; //! Particle velocity constraints std::vector> particle_velocity_constraints_; diff --git a/include/mesh.tcc b/include/mesh.tcc index 76e8c161e..30536c70e 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -83,6 +83,23 @@ void mpm::Mesh::iterate_over_nodes_predicate(Toper oper, Tpred pred) { } } +//! Iterate over particle set +template +template +void mpm::Mesh::iterate_over_node_set(int set_id, Toper oper) { + // If set id is -1, use all nodes + if (set_id == -1) { + this->iterate_over_nodes(oper); + } else { + // Iterate over the node set + auto nodes = node_sets_.at(set_id); +#pragma omp parallel for schedule(runtime) + for (auto nitr = nodes.cbegin(); nitr != nodes.cend(); ++nitr) { + oper(*nitr); + } + } +} + //! Create a list of active nodes in mesh template void mpm::Mesh::find_active_nodes() { @@ -1326,6 +1343,31 @@ void mpm::Mesh::apply_traction_on_particles(double current_time) { } } +//! Create nodal acceleration constraints +template +bool mpm::Mesh::create_nodal_acceleration_constraint( + int set_id, + const std::shared_ptr& constraint) { + bool status = true; + try { + if (set_id == -1 || node_sets_.find(set_id) != node_sets_.end()) { + // Create a nodal acceleration constraint + if (constraint->dir() < Tdim) + nodal_acceleration_constraints_.emplace_back(constraint); + else + throw std::runtime_error( + "Invalid direction of nodal acceleration constraint"); + } else + throw std::runtime_error( + "No node set found to assign nodal acceleration constraint"); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + //! Create particle velocity constraints template bool mpm::Mesh::create_particle_velocity_constraint( @@ -1349,7 +1391,7 @@ bool mpm::Mesh::create_particle_velocity_constraint( return status; } -//! Apply particle tractions +//! Apply particle velocity constraints template void mpm::Mesh::apply_particle_velocity_constraints() { // Iterate over all particle velocity constraints @@ -1366,6 +1408,23 @@ void mpm::Mesh::apply_particle_velocity_constraints() { } } +//! Update nodal acceleration constraints +template +void mpm::Mesh::update_nodal_acceleration_constraints( + double current_time) { + // Iterate over all nodal acceleration constraints + for (const auto& nacceleration : nodal_acceleration_constraints_) { + // If set id is -1, use all particles + int set_id = nacceleration->setid(); + unsigned dir = nacceleration->dir(); + double acceleration = nacceleration->acceleration(current_time); + + this->iterate_over_node_set( + set_id, std::bind(&mpm::NodeBase::update_acceleration_constraint, + std::placeholders::_1, dir, acceleration)); + } +} + //! Assign node tractions template bool mpm::Mesh::assign_nodal_concentrated_forces( diff --git a/include/node.h b/include/node.h index 31a6b9a03..bd61ecf62 100644 --- a/include/node.h +++ b/include/node.h @@ -198,9 +198,26 @@ class Node : public NodeBase { //! \param[in] velocity Applied velocity constraint bool assign_velocity_constraint(unsigned dir, double velocity) override; + //! Assign acceleration constraint + //! Directions can take values between 0 and Dim * Nphases + //! \param[in] dir Direction of acceleration constraint + //! \param[in] acceleration Applied acceleration constraint + bool assign_acceleration_constraint(unsigned dir, + double acceleration) override; + + //! Update acceleration constraint + //! Directions can take values between 0 and Dim * Nphases + //! \param[in] dir Direction of acceleration constraint + //! \param[in] acceleration Applied acceleration constraint + bool update_acceleration_constraint(unsigned dir, + double acceleration) override; + //! Apply velocity constraints void apply_velocity_constraints() override; + //! Apply acceleration constraints + void apply_acceleration_constraints() override; + //! Assign friction constraint //! Directions can take values between 0 and Dim * Nphases //! \param[in] dir Direction of friction constraint @@ -299,6 +316,8 @@ class Node : public NodeBase { Eigen::Matrix acceleration_; //! Velocity constraints std::map velocity_constraints_; + //! Acceleration constraints + std::map acceleration_constraints_; //! Rotation matrix for general velocity constraints Eigen::Matrix rotation_matrix_; //! Material ids whose information was passed to this node diff --git a/include/node.tcc b/include/node.tcc index 7d3a9cfcb..503c17416 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -14,8 +14,9 @@ mpm::Node::Node( "node" + std::to_string(Tdim) + "d::" + std::to_string(id); console_ = std::make_unique(logger, mpm::stdout_sink); - // Clear any velocity constraints + // Clear any velocity and acceleration constraints velocity_constraints_.clear(); + acceleration_constraints_.clear(); concentrated_force_.setZero(); this->initialise(); } @@ -234,6 +235,9 @@ bool mpm::Node::compute_acceleration_velocity( // Apply friction constraints this->apply_friction_constraints(dt); + // Apply acceleration constraints + this->apply_acceleration_constraints(); + // Velocity += acceleration * dt this->velocity_.col(phase) += this->acceleration_.col(phase) * dt; // Apply velocity constraints, which also sets acceleration to 0, @@ -270,6 +274,9 @@ bool mpm::Node::compute_acceleration_velocity_cundall( // Apply friction constraints this->apply_friction_constraints(dt); + // Apply acceleration constraints + this->apply_acceleration_constraints(); + // Velocity += acceleration * dt this->velocity_.col(phase) += this->acceleration_.col(phase) * dt; // Apply velocity constraints, which also sets acceleration to 0, @@ -309,6 +316,46 @@ bool mpm::Node::assign_velocity_constraint( return status; } +//! Assign acceleration constraint +//! Constrain directions can take values between 0 and Dim * Nphases +template +bool mpm::Node::assign_acceleration_constraint( + unsigned dir, double acceleration) { + bool status = true; + try { + //! Constrain directions can take values between 0 and Dim * Nphases + if (dir < (Tdim * Tnphases)) + this->acceleration_constraints_.insert(std::make_pair( + static_cast(dir), static_cast(acceleration))); + else + throw std::runtime_error("Constraint direction is out of bounds"); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Update acceleration constraint +template +bool mpm::Node::update_acceleration_constraint( + unsigned dir, double acceleration) { + bool status = true; + try { + // Check if an acceleration constraint was assigned to dir + if (acceleration_constraints_.find(dir) != acceleration_constraints_.end()) + acceleration_constraints_.find(dir)->second = acceleration; + else + throw std::runtime_error("Acceleration constraint direction is invalid"); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + //! Apply velocity constraints template void mpm::Node::apply_velocity_constraints() { @@ -346,6 +393,37 @@ void mpm::Node::apply_velocity_constraints() { } } +//! Apply acceleration constraints +template +void mpm::Node::apply_acceleration_constraints() { + // Set acceleration constraint + for (const auto& constraint : this->acceleration_constraints_) { + // Direction value in the constraint (0, Dim * Nphases) + const unsigned dir = constraint.first; + // Direction: dir % Tdim (modulus) + const auto direction = static_cast(dir % Tdim); + // Phase: Integer value of division (dir / Tdim) + const auto phase = static_cast(dir / Tdim); + + if (!generic_boundary_constraints_) { + // Acceleration constraints are applied on Cartesian boundaries + this->acceleration_(direction, phase) = constraint.second; + } else { + // Acceleration constraints on general boundaries + // Compute inverse rotation matrix + const Eigen::Matrix inverse_rotation_matrix = + rotation_matrix_.inverse(); + // Transform to local coordinate + Eigen::Matrix local_acceleration = + inverse_rotation_matrix * this->acceleration_; + // Apply boundary condition in local coordinate + local_acceleration(direction, phase) = constraint.second; + // Transform back to global coordinate + this->acceleration_ = rotation_matrix_ * local_acceleration; + } + } +} + //! Assign friction constraint //! Constrain directions can take values between 0 and Dim * Nphases template diff --git a/include/node_base.h b/include/node_base.h index 17eddfdca..010d98869 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -190,9 +190,26 @@ class NodeBase { //! \param[in] velocity Applied velocity constraint virtual bool assign_velocity_constraint(unsigned dir, double velocity) = 0; + //! Assign acceleration constraint + //! Directions can take values between 0 and Dim * Nphases + //! \param[in] dir Direction of acceleration constraint + //! \param[in] acceleration Applied acceleration constraint + virtual bool assign_acceleration_constraint(unsigned dir, + double acceleration) = 0; + + //! Update acceleration constraint + //! Directions can take values between 0 and Dim * Nphases + //! \param[in] dir Direction of acceleration constraint + //! \param[in] acceleration Applied acceleration constraint + virtual bool update_acceleration_constraint(unsigned dir, + double acceleration) = 0; + //! Apply velocity constraints virtual void apply_velocity_constraints() = 0; + //! Apply acceleration constraints + virtual void apply_acceleration_constraints() = 0; + //! Assign friction constraint //! Directions can take values between 0 and Dim * Nphases //! \param[in] dir Direction of friction constraint diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index c84eaf1f8..d05f627b1 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -27,6 +27,9 @@ #include "particle.h" #include "vector.h" +// CSV-parser +#include "csv/csv.h" + namespace mpm { //! Variable type @@ -121,6 +124,12 @@ class MPMBase : public MPM { void nodal_velocity_constraints( const Json& mesh_prop, const std::shared_ptr>& mesh_io); + //! Nodal acceleration constraints + //! \param[in] mesh_prop Mesh properties + //! \param[in] mesh_io Mesh IO handle + void nodal_acceleration_constraints( + const Json& mesh_prop, const std::shared_ptr>& mesh_io); + //! Nodal frictional constraints //! \param[in] mesh_prop Mesh properties //! \param[in] mesh_io Mesh IO handle diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index f829b5f95..9e135e9c3 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -233,6 +233,9 @@ void mpm::MPMBase::initialise_mesh() { // Read and assign velocity constraints this->nodal_velocity_constraints(mesh_props, mesh_io); + // Read and assign acceleration constraints + this->nodal_acceleration_constraints(mesh_props, mesh_io); + // Read and assign friction constraints this->nodal_frictional_constraints(mesh_props, mesh_io); @@ -760,16 +763,35 @@ bool mpm::MPMBase::initialise_math_functions(const Json& math_functions) { const std::string function_type = function_props["type"].template get(); + // Initiate another function_prop to be passed + auto function_props_update = function_props; + + // Create a file reader + const std::string io_type = + io_->json_object("mesh")["io_type"].template get(); + auto reader = Factory>::instance()->create(io_type); + + // Math function is specified in a file, replace function_props_update + if (function_props.find("file") != function_props.end()) { + // Read file and store in array of vectors + std::string math_file = io_->file_name( + function_props.at("file").template get()); + auto xfx_values = reader->read_math_functions(math_file); + + function_props_update["xvalues"] = xfx_values[0]; + function_props_update["fxvalues"] = xfx_values[1]; + } + // Create a new function from JSON object auto function = Factory::instance()->create( - function_type, std::move(function_id), function_props); + function_type, std::move(function_id), function_props_update); - // Add material to list + // Add math function to list auto insert_status = math_functions_.insert(std::make_pair(function->id(), function)); - // If insert material failed + // If insert math function failed if (!insert_status.second) { status = false; throw std::runtime_error( @@ -883,6 +905,65 @@ void mpm::MPMBase::nodal_velocity_constraints( } } +// Nodal acceleration constraints +template +void mpm::MPMBase::nodal_acceleration_constraints( + const Json& mesh_props, const std::shared_ptr>& mesh_io) { + try { + // Read and assign acceleration constraints + if (mesh_props.find("boundary_conditions") != mesh_props.end() && + mesh_props["boundary_conditions"].find("acceleration_constraints") != + mesh_props["boundary_conditions"].end()) { + // Iterate over acceleration constraints + for (const auto& constraints : + mesh_props["boundary_conditions"]["acceleration_constraints"]) { + // Acceleration constraints are specified in a file + if (constraints.find("file") != constraints.end()) { + std::string acceleration_constraints_file = + constraints.at("file").template get(); + bool acceleration_constraints = + constraints_->assign_nodal_acceleration_constraints( + mesh_io->read_acceleration_constraints( + io_->file_name(acceleration_constraints_file))); + if (!acceleration_constraints) + throw std::runtime_error( + "Acceleration constraints are not properly assigned"); + + } else { + // Set id + int nset_id = constraints.at("nset_id").template get(); + // Direction + unsigned dir = constraints.at("dir").template get(); + // Acceleration + double acceleration = + constraints.at("acceleration").template get(); + // Get the math function + std::shared_ptr afunction = nullptr; + if (constraints.find("math_function_id") != constraints.end()) + afunction = math_functions_.at( + constraints.at("math_function_id").template get()); + // Add acceleration constraint to mesh + auto acceleration_constraint = + std::make_shared(nset_id, afunction, + dir, acceleration); + mesh_->create_nodal_acceleration_constraint(nset_id, + acceleration_constraint); + bool acceleration_constraints = + constraints_->assign_nodal_acceleration_constraint( + nset_id, acceleration_constraint); + if (!acceleration_constraints) + throw std::runtime_error( + "Nodal acceleration constraint is not properly assigned"); + } + } + } else + throw std::runtime_error("Acceleration constraints JSON not found"); + } catch (std::exception& exception) { + console_->warn("#{}: Acceleration constraints are undefined {} ", __LINE__, + exception.what()); + } +} + // Nodal frictional constraints template void mpm::MPMBase::nodal_frictional_constraints( @@ -892,7 +973,7 @@ void mpm::MPMBase::nodal_frictional_constraints( if (mesh_props.find("boundary_conditions") != mesh_props.end() && mesh_props["boundary_conditions"].find("friction_constraints") != mesh_props["boundary_conditions"].end()) { - // Iterate over velocity constraints + // Iterate over friction constraints for (const auto& constraints : mesh_props["boundary_conditions"]["friction_constraints"]) { // Friction constraints are specified in a file diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc index 098a1488f..bbaff9f24 100644 --- a/include/solvers/mpm_explicit.tcc +++ b/include/solvers/mpm_explicit.tcc @@ -147,7 +147,7 @@ bool mpm::MPMExplicit::solve() { // Particle kinematics mpm_scheme_->compute_particle_kinematics(velocity_update_, phase, "Cundall", - damping_factor_); + damping_factor_, step_); // Update Stress Last mpm_scheme_->postcompute_stress_strain(phase, pressure_smoothing_); diff --git a/include/solvers/mpm_scheme/mpm_scheme.h b/include/solvers/mpm_scheme/mpm_scheme.h index 683d149ad..31739938e 100644 --- a/include/solvers/mpm_scheme/mpm_scheme.h +++ b/include/solvers/mpm_scheme/mpm_scheme.h @@ -62,9 +62,10 @@ class MPMScheme { //! \param[in] phase Phase of particle //! \param[in] damping_type Type of damping //! \param[in] damping_factor Value of critical damping + //! \param[in] step Number of step in solver virtual inline void compute_particle_kinematics( bool velocity_update, unsigned phase, const std::string& damping_type, - double damping_factor); + double damping_factor, unsigned step); //! Compute particle location //! \param[in] locate_particles Flag to enable locate particles, if set to diff --git a/include/solvers/mpm_scheme/mpm_scheme.tcc b/include/solvers/mpm_scheme/mpm_scheme.tcc index e95c61795..f6c5fc09b 100644 --- a/include/solvers/mpm_scheme/mpm_scheme.tcc +++ b/include/solvers/mpm_scheme/mpm_scheme.tcc @@ -172,7 +172,10 @@ inline void mpm::MPMScheme::compute_forces( template inline void mpm::MPMScheme::compute_particle_kinematics( bool velocity_update, unsigned phase, const std::string& damping_type, - double damping_factor) { + double damping_factor, unsigned step) { + + // Update nodal acceleration constraints + mesh_->update_nodal_acceleration_constraints(step * dt_); // Check if damping has been specified and accordingly Iterate over // active nodes to compute acceleratation and velocity diff --git a/tests/acceleration_constraint_test.cc b/tests/acceleration_constraint_test.cc new file mode 100644 index 000000000..12764a371 --- /dev/null +++ b/tests/acceleration_constraint_test.cc @@ -0,0 +1,56 @@ +#include + +#include "catch.hpp" + +#include "acceleration_constraint.h" +#include "function_base.h" +#include "linear_function.h" + +//! Check acceleration constraint +TEST_CASE("Acceleration constraint is checked", "[acceleration][constraint]") { + // Tolerance + const double Tolerance = 1.E-9; + + SECTION("Check acceleration constraint") { + + // Json property + Json jfunctionproperties; + jfunctionproperties["id"] = 0; + std::vector x_values{{0.0, 0.5, 1.0, 1.5}}; + std::vector fx_values{{0.0, 1.0, 1.0, 0.0}}; + jfunctionproperties["xvalues"] = x_values; + jfunctionproperties["fxvalues"] = fx_values; + + // math function + std::shared_ptr mfunction = + std::make_shared(0, jfunctionproperties); + + int set_id = -1; + unsigned dir = 1; + double acceleration = 10; + auto nacceleration = std::make_shared( + set_id, mfunction, dir, acceleration); + + // Check particle set id + REQUIRE(nacceleration->setid() == set_id); + // Check direction + REQUIRE(nacceleration->dir() == dir); + // Check acceleration at x = 0 + REQUIRE(nacceleration->acceleration(0.00) == Approx(0.).epsilon(Tolerance)); + // Check acceleration at x = 0.25 + REQUIRE(nacceleration->acceleration(0.25) == Approx(5.).epsilon(Tolerance)); + // Check acceleration at x = 0.5 + REQUIRE(nacceleration->acceleration(0.50) == + Approx(10.).epsilon(Tolerance)); + // Check acceleration at x = 0.75 + REQUIRE(nacceleration->acceleration(0.75) == + Approx(10.).epsilon(Tolerance)); + // Check acceleration at x = 1.0 + REQUIRE(nacceleration->acceleration(1.00) == + Approx(10.).epsilon(Tolerance)); + // Check acceleration at x = 1.25 + REQUIRE(nacceleration->acceleration(1.25) == Approx(5.).epsilon(Tolerance)); + // Check acceleration at x = 5 + REQUIRE(nacceleration->acceleration(5.00) == Approx(0.).epsilon(Tolerance)); + } +} diff --git a/tests/io/io_mesh_ascii_test.cc b/tests/io/io_mesh_ascii_test.cc index 13779101b..79c349ab6 100644 --- a/tests/io/io_mesh_ascii_test.cc +++ b/tests/io/io_mesh_ascii_test.cc @@ -247,6 +247,63 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") { } } + SECTION("Check acceleration constraints file") { + // Vector of acceleration constraints + std::vector> + acceleration_constraints; + + // Constraint + acceleration_constraints.emplace_back(std::make_tuple(0, 0, 10.5)); + acceleration_constraints.emplace_back(std::make_tuple(1, 1, -10.5)); + acceleration_constraints.emplace_back(std::make_tuple(2, 0, -12.5)); + acceleration_constraints.emplace_back(std::make_tuple(3, 1, 0.0)); + + // Dump constraints as an input file to be read + std::ofstream file; + file.open("acceleration-constraints-2d.txt"); + // Write particle coordinates + for (const auto& acceleration_constraint : acceleration_constraints) { + file << std::get<0>(acceleration_constraint) << "\t"; + file << std::get<1>(acceleration_constraint) << "\t"; + file << std::get<2>(acceleration_constraint) << "\t"; + + file << "\n"; + } + + file.close(); + + // Check read acceleration constraints + SECTION("Check acceleration constraints") { + // Create a read_mesh object + auto read_mesh = std::make_unique>(); + + // Try to read constraints from a non-existant file + auto constraints = read_mesh->read_acceleration_constraints( + "acceleration-constraints-missing.txt"); + // Check number of constraints + REQUIRE(constraints.size() == 0); + + // Check constraints + constraints = read_mesh->read_acceleration_constraints( + "acceleration-constraints-2d.txt"); + // Check number of particles + REQUIRE(constraints.size() == acceleration_constraints.size()); + + // Check coordinates of nodes + for (unsigned i = 0; i < acceleration_constraints.size(); ++i) { + REQUIRE(std::get<0>(constraints.at(i)) == + Approx(std::get<0>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + REQUIRE(std::get<1>(constraints.at(i)) == + Approx(std::get<1>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + REQUIRE(std::get<2>(constraints.at(i)) == + Approx(std::get<2>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + } + } + } + SECTION("Check friction constraints file") { // Vector of friction constraints std::vector> @@ -337,7 +394,7 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") { // Create a read_mesh object auto read_mesh = std::make_unique>(); - // Try to read constrtaints from a non-existant file + // Try to read forces from a non-existant file auto forces = read_mesh->read_forces("forces-missing.txt"); // Check number of forces REQUIRE(forces.size() == 0); @@ -482,6 +539,54 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") { } } } + + SECTION("Check math function file") { + // Vector of math function entries + std::array, 2> entries; + + // Populate the math function entries + for (int i = 0; i < 3; ++i) { + entries[0].push_back(0.7 * i); + entries[1].push_back(2.2 * i); + } + + // Dump the math entries as an input file to be read + std::ofstream file; + file.open("math-function-2d.csv"); + // Write math entries for x and fx + for (int i = 0; i < 3; ++i) + file << entries[0][i] << "," << entries[1][i] << "\n"; + + file.close(); + + // Check read math funciton file + SECTION("Check math function entries") { + // Create an io_mesh object + auto io_mesh = std::make_unique>(); + + // Try to read math funciton entries from a non-existant file + auto math_function_values = + io_mesh->read_math_functions("math-function-missing.csv"); + // Check number of math function entries + REQUIRE(math_function_values[0].size() == 0); + REQUIRE(math_function_values[1].size() == 0); + + // Check math function reading + math_function_values = + io_mesh->read_math_functions("math-function-2d.csv"); + // Check number of entries from math function file + REQUIRE(math_function_values[0].size() == entries[0].size()); + REQUIRE(math_function_values[1].size() == entries[1].size()); + REQUIRE(math_function_values[0].size() == math_function_values[1].size()); + + // Check entry values + for (unsigned i = 0; i < 2; ++i) { + for (unsigned j = 0; j < entries[0].size(); ++j) + REQUIRE(math_function_values[i][j] == + Approx(entries[i][j]).epsilon(Tolerance)); + } + } + } } // Check IOMeshAscii @@ -740,7 +845,7 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") { // Create a read_mesh object auto read_mesh = std::make_unique>(); - // Try to read constrtaints from a non-existant file + // Try to read constraints from a non-existant file auto constraints = read_mesh->read_velocity_constraints( "velocity-constraints-missing.txt"); // Check number of constraints @@ -767,6 +872,63 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") { } } + SECTION("Check acceleration constraints file") { + // Vector of particle coordinates + std::vector> + acceleration_constraints; + + // Constraint + acceleration_constraints.emplace_back(std::make_tuple(0, 0, 10.5)); + acceleration_constraints.emplace_back(std::make_tuple(1, 1, -10.5)); + acceleration_constraints.emplace_back(std::make_tuple(2, 2, -12.5)); + acceleration_constraints.emplace_back(std::make_tuple(3, 0, 0.0)); + + // Dump constratints as an input file to be read + std::ofstream file; + file.open("acceleration-constraints-3d.txt"); + // Write particle coordinates + for (const auto& acceleration_constraint : acceleration_constraints) { + file << std::get<0>(acceleration_constraint) << "\t"; + file << std::get<1>(acceleration_constraint) << "\t"; + file << std::get<2>(acceleration_constraint) << "\t"; + + file << "\n"; + } + + file.close(); + + // Check read acceleration constraints + SECTION("Check acceleration constraints") { + // Create a read_mesh object + auto read_mesh = std::make_unique>(); + + // Try to read constraints from a non-existant file + auto constraints = read_mesh->read_acceleration_constraints( + "acceleration-constraints-missing.txt"); + // Check number of constraints + REQUIRE(constraints.size() == 0); + + // Check constraints + constraints = read_mesh->read_acceleration_constraints( + "acceleration-constraints-3d.txt"); + // Check number of particles + REQUIRE(constraints.size() == acceleration_constraints.size()); + + // Check coordinates of nodes + for (unsigned i = 0; i < acceleration_constraints.size(); ++i) { + REQUIRE(std::get<0>(constraints.at(i)) == + Approx(std::get<0>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + REQUIRE(std::get<1>(constraints.at(i)) == + Approx(std::get<1>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + REQUIRE(std::get<2>(constraints.at(i)) == + Approx(std::get<2>(acceleration_constraints.at(i))) + .epsilon(Tolerance)); + } + } + } + SECTION("Check friction constraints file") { // Vector of friction constraints std::vector> @@ -1003,4 +1165,52 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") { } } } + + SECTION("Check math function file") { + // Vector of math function entries + std::array, 2> entries; + + // Populate the math function entries + for (int i = 0; i < 4; ++i) { + entries[0].push_back(0.3 * i); + entries[1].push_back(2.6 * i); + } + + // Dump the math entries as an input file to be read + std::ofstream file; + file.open("math-function-3d.csv"); + // Write math entries for x and fx + for (int i = 0; i < 4; ++i) + file << entries[0][i] << "," << entries[1][i] << "\n"; + + file.close(); + + // Check read math funciton file + SECTION("Check math function entries") { + // Create an io_mesh object + auto io_mesh = std::make_unique>(); + + // Try to read math funciton entries from a non-existant file + auto math_function_values = + io_mesh->read_math_functions("math-function-missing.csv"); + // Check number of math function entries + REQUIRE(math_function_values[0].size() == 0); + REQUIRE(math_function_values[1].size() == 0); + + // Check math function reading + math_function_values = + io_mesh->read_math_functions("math-function-3d.csv"); + // Check number of entries from math function file + REQUIRE(math_function_values[0].size() == entries[0].size()); + REQUIRE(math_function_values[1].size() == entries[1].size()); + REQUIRE(math_function_values[0].size() == math_function_values[1].size()); + + // Check entry values + for (unsigned i = 0; i < 2; ++i) { + for (unsigned j = 0; j < entries[0].size(); ++j) + REQUIRE(math_function_values[i][j] == + Approx(entries[i][j]).epsilon(Tolerance)); + } + } + } } diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index eedc61044..2fd40135c 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1104,6 +1104,140 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { REQUIRE(constraints->assign_nodal_velocity_constraint( set_id, velocity_constraint) == false); } + // Test assign acceleration constraints to nodes + SECTION("Check assign constant acceleration constraints to nodes") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + int set_id = 0; + int dir = 0; + double constraint = 10.5; + // Add acceleration constraint to mesh + auto acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + set_id = 1; + dir = 1; + constraint = -12.5; + // Add acceleration constraint to mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Add acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, nullptr, dir, + constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // When constraints fail + dir = 2; + // Add acceleration constraint to mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + dir = 1; + // Create acceleration constraint from mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Create acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, nullptr, dir, + constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Failed attempt to create acceleration constraint with + // non-existing set id + set_id = 2; + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + // Failed attempt to create acceleration constraint with + // out of bounds direction + set_id = 0; + dir = 3; + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + } + + SECTION("Check assign varying acceleration constraints to nodes") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + int set_id = 0; + int dir = 0; + double constraint = 10.5; + + // Json property + Json jfunctionproperties; + jfunctionproperties["id"] = 0; + std::vector x_values{{0.0, 0.5, 1.0}}; + std::vector fx_values{{0.0, 1.0, 1.0}}; + jfunctionproperties["xvalues"] = x_values; + jfunctionproperties["fxvalues"] = fx_values; + + // math function + std::shared_ptr mfunction = + std::make_shared(0, jfunctionproperties); + + // Add acceleration constraint to mesh + auto acceleration_constraint = + std::make_shared(set_id, mfunction, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Create acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, mfunction, dir, + constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Failed attempt to create acceleration constraint with + // non-existing set id + set_id = 2; + acceleration_constraint = + std::make_shared(set_id, mfunction, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + // Update acceleration constraints + REQUIRE_NOTHROW(mesh->update_nodal_acceleration_constraints(0.5)); + } SECTION("Check assign friction constraints to nodes") { tsl::robin_map> node_sets; @@ -1171,6 +1305,27 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { velocity_constraints) == false); } + // Test assign acceleration constraints to nodes + SECTION("Check assign acceleration constraints to nodes") { + // Vector of particle coordinates + std::vector> + acceleration_constraints; + //! Constraints object + auto constraints = std::make_shared>(mesh); + // Constraint + acceleration_constraints.emplace_back(std::make_tuple(0, 0, 10.5)); + acceleration_constraints.emplace_back(std::make_tuple(1, 1, -10.5)); + acceleration_constraints.emplace_back(std::make_tuple(2, 0, -12.5)); + acceleration_constraints.emplace_back(std::make_tuple(3, 1, 0.0)); + + REQUIRE(constraints->assign_nodal_acceleration_constraints( + acceleration_constraints) == true); + // When constraints fail + acceleration_constraints.emplace_back(std::make_tuple(3, 2, 0.0)); + REQUIRE(constraints->assign_nodal_acceleration_constraints( + acceleration_constraints) == false); + } + // Test assign friction constraints to nodes SECTION("Check assign friction constraints to nodes") { // Vector of particle coordinates diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc index 06f4e0cf2..1bbfe3421 100644 --- a/tests/mesh_test_3d.cc +++ b/tests/mesh_test_3d.cc @@ -1215,6 +1215,140 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") { REQUIRE(constraints->assign_nodal_velocity_constraint( set_id, velocity_constraint) == false); } + // Test assign acceleration constraints to nodes + SECTION("Check assign constant acceleration constraints to nodes") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + int set_id = 0; + int dir = 0; + double constraint = 10.5; + // Add acceleration constraint to mesh + auto acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + set_id = 1; + dir = 1; + constraint = -12.5; + // Add acceleration constraint to mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Add acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, nullptr, dir, + constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // When constraints fail + dir = 3; + // Add acceleration constraint to mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(constraints->assign_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + dir = 0; + // Create acceleration constraint from mesh + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Create acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, nullptr, dir, + constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Failed attempt to create acceleration constraint with + // non-existing set id + set_id = 2; + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + // Failed attempt to create acceleration constraint with + // out of bounds direction + set_id = 0; + dir = 3; + acceleration_constraint = + std::make_shared(set_id, nullptr, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + } + + SECTION("Check assign varying acceleration constraints to nodes") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + int set_id = 0; + int dir = 0; + double constraint = 10.5; + + // Json property + Json jfunctionproperties; + jfunctionproperties["id"] = 0; + std::vector x_values{{0.0, 0.5, 1.0}}; + std::vector fx_values{{0.0, 1.0, 1.0}}; + jfunctionproperties["xvalues"] = x_values; + jfunctionproperties["fxvalues"] = fx_values; + + // math function + std::shared_ptr mfunction = + std::make_shared(0, jfunctionproperties); + + // Add acceleration constraint to mesh + auto acceleration_constraint = + std::make_shared(set_id, mfunction, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Create acceleration constraint to all nodes in mesh + acceleration_constraint = + std::make_shared(-1, mfunction, dir, + constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == true); + + // Failed attempt to create acceleration constraint with + // non-existing set id + set_id = 2; + acceleration_constraint = + std::make_shared(set_id, mfunction, + dir, constraint); + REQUIRE(mesh->create_nodal_acceleration_constraint( + set_id, acceleration_constraint) == false); + + // Update acceleration constraints + REQUIRE_NOTHROW(mesh->update_nodal_acceleration_constraints(0.5)); + } SECTION("Check assign friction constraints to nodes") { tsl::robin_map> node_sets; @@ -1282,6 +1416,27 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") { velocity_constraints) == false); } + // Test assign acceleration constraints to nodes + SECTION("Check assign acceleration constraints to nodes") { + // Vector of particle coordinates + std::vector> + acceleration_constraints; + //! Constraints object + auto constraints = std::make_shared>(mesh); + // Constraint + acceleration_constraints.emplace_back(std::make_tuple(0, 0, 10.5)); + acceleration_constraints.emplace_back(std::make_tuple(1, 1, -10.5)); + acceleration_constraints.emplace_back(std::make_tuple(2, 0, -12.5)); + acceleration_constraints.emplace_back(std::make_tuple(3, 1, 0.0)); + + REQUIRE(constraints->assign_nodal_acceleration_constraints( + acceleration_constraints) == true); + // When constraints fail + acceleration_constraints.emplace_back(std::make_tuple(3, 3, 0.0)); + REQUIRE(constraints->assign_nodal_acceleration_constraints( + acceleration_constraints) == false); + } + // Test assign friction constraints to nodes SECTION("Check assign friction constraints to nodes") { // Vector of particle coordinates diff --git a/tests/node_test.cc b/tests/node_test.cc index 3cc92d2a5..c2230a589 100644 --- a/tests/node_test.cc +++ b/tests/node_test.cc @@ -307,51 +307,74 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(node->velocity(Nphase)(i) == Approx(velocity(i)).epsilon(Tolerance)); + SECTION("Check with friction constraint") { + // Apply friction constraints + REQUIRE(node->assign_friction_constraint(0, 1., 0.5) == true); + // Apply friction constraints + REQUIRE(node->assign_friction_constraint(-1, 1., 0.5) == false); + REQUIRE(node->assign_friction_constraint(3, 1., 0.5) == false); - // Apply friction constraints - REQUIRE(node->assign_friction_constraint(0, 1., 0.5) == true); - // Apply friction constraints - REQUIRE(node->assign_friction_constraint(-1, 1., 0.5) == false); - REQUIRE(node->assign_friction_constraint(3, 1., 0.5) == false); + // Test acceleration with constraints + acceleration[0] = 0.5 * acceleration[0]; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); - // Test acceleration with constraints - acceleration[0] = 0.5 * acceleration[0]; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + // Apply cundall damping when calculating acceleration + REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == + true); - // Apply cundall damping when calculating acceleration - REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == - true); + // Test acceleration with cundall damping + acceleration[0] = 0.; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } - // Test acceleration with cundall damping - acceleration[0] = 0.; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + SECTION("Check with velocity constraint") { + // Apply velocity constraints + REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); - // Apply velocity constraints - REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + // Test velocity with constraints + velocity[0] = 10.5; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); - // Test velocity with constraints - velocity[0] = 10.5; - for (unsigned i = 0; i < velocity.size(); ++i) - REQUIRE(node->velocity(Nphase)(i) == - Approx(velocity(i)).epsilon(Tolerance)); + // Test acceleration with constraints + acceleration[0] = 0.; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } - // Test acceleration with constraints - acceleration[0] = 0.; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + SECTION("Check with acceleration constraint") { + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, 5.0) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); - // Exception check when mass is zero - mass = 0.; - // Update mass to 0. - REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); - REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + // Test velocity with constraints + velocity[0] = 5.0 * dt; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); + + // Test acceleration with constraints + acceleration[0] = 5.0; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } + + SECTION("Check with no mass") { + // Exception check when mass is zero + mass = 0.; + // Update mass to 0. + REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); + REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + } } SECTION("Check momentum and velocity") { @@ -453,6 +476,36 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { for (unsigned i = 0; i < acceleration.size(); ++i) REQUIRE(node->acceleration(Nphase)(i) == Approx(acceleration(i)).epsilon(Tolerance)); + + // Update acceleration constraints before assignement + REQUIRE(node->update_acceleration_constraint(0, 1.5) == false); + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, 10.5) == true); + // Check out of bounds condition + REQUIRE(node->assign_acceleration_constraint(1, 12.5) == false); + + // Apply acceleration constraints + node->apply_acceleration_constraints(); + + // Check apply constraints + acceleration << 10.5; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Update acceleration constraints + REQUIRE(node->update_acceleration_constraint(0, 3.5) == true); + // Check out of bounds condition + REQUIRE(node->update_acceleration_constraint(1, 12.5) == false); + + // Apply acceleration constraints + node->apply_acceleration_constraints(); + + // Check updated constraints + acceleration << 3.5; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); } SECTION("Check node material ids") { @@ -803,38 +856,60 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE(node->velocity(Nphase)(i) == Approx(velocity(i)).epsilon(Tolerance)); - // Apply velocity constraints - REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + SECTION("Check with velocity constraint") { + // Apply velocity constraints + REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); - // Test velocity with constraints - velocity << 10.5, 0.03; - for (unsigned i = 0; i < velocity.size(); ++i) - REQUIRE(node->velocity(Nphase)(i) == - Approx(velocity(i)).epsilon(Tolerance)); + // Test velocity with constraints + velocity << 10.5, 0.03; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); - // Test acceleration with constraints - acceleration[0] = 0.; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + // Test acceleration with constraints + acceleration[0] = 0.; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); - // Apply cundall damping when calculating acceleration - REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == - true); + // Apply cundall damping when calculating acceleration + REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == + true); - // Test acceleration with cundall damping - acceleration << 0., 0.1425; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + // Test acceleration with cundall damping + acceleration << 0., 0.1425; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } - // Exception check when mass is zero - mass = 0.; - // Update mass to 0. - REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); - REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + SECTION("Check with acceleration constraint") { + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, 5.0) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + + // Test acceleration with constraints + acceleration << 5.0, 0.15; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Test velocity with constraints + velocity << 0.5, 0.03; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); + } + + SECTION("Check with no mass") { + // Exception check when mass is zero + mass = 0.; + // Update mass to 0. + REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); + REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + } } SECTION("Check momentum, velocity and acceleration") { @@ -1056,6 +1131,70 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(i) == Approx(acceleration(i)).epsilon(Tolerance)); } + + SECTION("Check Cartesian acceleration constraints") { + // Update acceleration constraints before any assignement + REQUIRE(node->update_acceleration_constraint(0, 4.5) == false); + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, -12.5) == true); + REQUIRE(node->assign_acceleration_constraint(1, 4.1) == true); + // Check out of bounds condition + REQUIRE(node->assign_acceleration_constraint(2, 0.) == false); + + // Apply constraints + node->apply_acceleration_constraints(); + + // Check apply constraints + acceleration << -12.5, 4.1; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Update acceleration constraints + REQUIRE(node->update_acceleration_constraint(0, 4.5) == true); + REQUIRE(node->update_acceleration_constraint(1, 1.9) == true); + // Check out of bounds condition + REQUIRE(node->update_acceleration_constraint(2, 0.) == false); + + // Apply constraints + node->apply_acceleration_constraints(); + + // Check updated constraints + acceleration << 4.5, 1.9; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } + + SECTION("Check general acceleration constraints in all directions") { + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, -12.5) == true); + REQUIRE(node->assign_acceleration_constraint(1, 7.5) == true); + + // Apply rotation matrix with Euler angles alpha = -10 deg, beta = 30 + // deg + Eigen::Matrix euler_angles; + euler_angles << -10. * M_PI / 180, 30. * M_PI / 180; + const auto rotation_matrix = + mpm::geometry::rotation_matrix(euler_angles); + node->assign_rotation_matrix(rotation_matrix); + const auto inverse_rotation_matrix = rotation_matrix.inverse(); + + // Apply inclined acceleration constraints + node->apply_acceleration_constraints(); + + // Check applied acceleration constraints in the global coordinates + acceleration << -14.311308834766370, 2.772442864323454; + for (unsigned i = 0; i < Dim; ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Check that the acceleration is as specified in local coordinate + REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(0) == + Approx(-12.5).epsilon(Tolerance)); + REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(1) == + Approx(7.5).epsilon(Tolerance)); + } } SECTION("Check node material ids") { @@ -1378,42 +1517,64 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { REQUIRE(node->velocity(Nphase)(i) == Approx(velocity(i)).epsilon(Tolerance)); - // Apply velocity constraints - REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + SECTION("Check with velocity constraint") { + // Apply velocity constraints + REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); - // Test velocity with constraints - velocity << 10.5, 0.03, 0.06; - for (unsigned i = 0; i < velocity.size(); ++i) - REQUIRE(node->velocity(Nphase)(i) == - Approx(velocity(i)).epsilon(Tolerance)); + // Test velocity with constraints + velocity << 10.5, 0.03, 0.06; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); - // Test acceleration with constraints - acceleration[0] = 0.; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + // Test acceleration with constraints + acceleration[0] = 0.; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); - // Apply cundall damping when calculating acceleration - REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == - true); + // Apply cundall damping when calculating acceleration + REQUIRE(node->compute_acceleration_velocity_cundall(Nphase, dt, 0.05) == + true); - // Test acceleration with cundall damping - acceleration << 0.0, 0.13322949016875, 0.28322949016875; - for (unsigned i = 0; i < acceleration.size(); ++i) - REQUIRE(node->acceleration(Nphase)(i) == - Approx(acceleration(i)).epsilon(Tolerance)); + // Test acceleration with cundall damping + acceleration << 0.0, 0.13322949016875, 0.28322949016875; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); - // Apply velocity constraints - REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + // Apply velocity constraints + REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + } - // Exception check when mass is zero - mass = 0.; - // Update mass to 0. - REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); - REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + SECTION("Check with acceleration constraint") { + // Apply velocity constraints + REQUIRE(node->assign_acceleration_constraint(0, 5.0) == true); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + + // Test velocity with constraints + velocity << 0.5, 0.03, 0.06; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); + + // Test acceleration with constraints + acceleration << 5.0, 0.15, 0.3; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } + + SECTION("Check with no mass") { + // Exception check when mass is zero + mass = 0.; + // Update mass to 0. + REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); + REQUIRE(node->mass(Nphase) == Approx(mass).epsilon(Tolerance)); + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == false); + } } SECTION("Check momentum, velocity and acceleration") { @@ -1645,6 +1806,76 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(i) == Approx(acceleration(i)).epsilon(Tolerance)); } + + SECTION("Check Cartesian acceleration constraints") { + // Update acceleration constraints before assignement + REQUIRE(node->update_acceleration_constraint(0, 2.1) == false); + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, 10.5) == true); + REQUIRE(node->assign_acceleration_constraint(1, -12.5) == true); + REQUIRE(node->assign_acceleration_constraint(2, -7.5) == true); + // Check out of bounds condition + REQUIRE(node->assign_acceleration_constraint(4, 0.) == false); + + // Apply constraints + node->apply_acceleration_constraints(); + + // Check apply constraints + acceleration << 10.5, -12.5, -7.5; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Update acceleration constraints + REQUIRE(node->update_acceleration_constraint(0, 2.1) == true); + REQUIRE(node->update_acceleration_constraint(1, -3.0) == true); + REQUIRE(node->update_acceleration_constraint(2, 4.6) == true); + // Check out of bounds condition + REQUIRE(node->update_acceleration_constraint(4, 0.) == false); + + // Apply constraints + node->apply_acceleration_constraints(); + + // Check updated constraints + acceleration << 2.1, -3.0, 4.6; + for (unsigned i = 0; i < acceleration.size(); ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + } + + SECTION("Check general acceleration constraints in all directions") { + // Apply acceleration constraints + REQUIRE(node->assign_acceleration_constraint(0, 10.5) == true); + REQUIRE(node->assign_acceleration_constraint(1, -12.5) == true); + REQUIRE(node->assign_acceleration_constraint(2, 7.5) == true); + + // Apply rotation matrix with Euler angles alpha = -10 deg, beta = 20 + // deg and gamma = -30 deg + Eigen::Matrix euler_angles; + euler_angles << -10. * M_PI / 180, 20. * M_PI / 180, -30. * M_PI / 180; + const auto rotation_matrix = + mpm::geometry::rotation_matrix(euler_angles); + node->assign_rotation_matrix(rotation_matrix); + const auto inverse_rotation_matrix = rotation_matrix.inverse(); + + // Apply constraints + node->apply_acceleration_constraints(); + + // Check apply constraints + acceleration << 13.351984588153375, -5.717804716697730, + 10.572663655835457; + for (unsigned i = 0; i < Dim; ++i) + REQUIRE(node->acceleration(Nphase)(i) == + Approx(acceleration(i)).epsilon(Tolerance)); + + // Check that the acceleration is as specified in local coordinate + REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(0) == + Approx(10.5).epsilon(Tolerance)); + REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(1) == + Approx(-12.5).epsilon(Tolerance)); + REQUIRE((inverse_rotation_matrix * node->acceleration(Nphase))(2) == + Approx(7.5).epsilon(Tolerance)); + } } SECTION("Check node material ids") { diff --git a/tests/solvers/mpm_scheme_test.cc b/tests/solvers/mpm_scheme_test.cc index f0f351aa1..dab90a997 100644 --- a/tests/solvers/mpm_scheme_test.cc +++ b/tests/solvers/mpm_scheme_test.cc @@ -194,14 +194,14 @@ TEST_CASE("Stress update is checked for USF and USL", REQUIRE_NOTHROW(mpm_scheme->compute_forces(gravity, phase, step, true)); // Particle kinematics - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(true, phase, "Cundall", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(false, phase, "Cundall", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(true, phase, "None", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(false, phase, "None", 0.02)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + true, phase, "Cundall", 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + false, phase, "Cundall", 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics(true, phase, "None", + 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + false, phase, "None", 0.02, step)); // Update Stress Last REQUIRE_NOTHROW(mpm_scheme->postcompute_stress_strain(phase, true)); @@ -235,14 +235,14 @@ TEST_CASE("Stress update is checked for USF and USL", REQUIRE_NOTHROW(mpm_scheme->compute_forces(gravity, phase, step, true)); // Particle kinematics - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(true, phase, "Cundall", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(false, phase, "Cundall", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(true, phase, "None", 0.02)); - REQUIRE_NOTHROW( - mpm_scheme->compute_particle_kinematics(false, phase, "None", 0.02)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + true, phase, "Cundall", 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + false, phase, "Cundall", 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics(true, phase, "None", + 0.02, step)); + REQUIRE_NOTHROW(mpm_scheme->compute_particle_kinematics( + false, phase, "None", 0.02, step)); // Update Stress Last REQUIRE_NOTHROW(mpm_scheme->postcompute_stress_strain(phase, true));