diff --git a/src/shared/kernels/all_kernels.h b/src/shared/kernels/all_kernels.h index 5e8e980fc2..b618a4e64b 100644 --- a/src/shared/kernels/all_kernels.h +++ b/src/shared/kernels/all_kernels.h @@ -34,5 +34,7 @@ #include "kernel_laguerre_gauss.h" #include "kernel_tabulated.h" #include "kernel_wenland_c2.h" +#include "anisotropic_kernel.hpp" + #endif // ALL_KERNELS_H \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/anisotropic_kernel.h b/src/shared/kernels/anisotropic_kernel.h similarity index 96% rename from tests/user_examples/extra_src/shared/anisotropic_kernel.h rename to src/shared/kernels/anisotropic_kernel.h index 80801c61f8..43ea4112bc 100644 --- a/tests/user_examples/extra_src/shared/anisotropic_kernel.h +++ b/src/shared/kernels/anisotropic_kernel.h @@ -26,7 +26,7 @@ * implemented in derived classes. The kernel function define the relevance * between two neighboring particles. Basically, the further the two * particles, the less relevance they have. - * @author Zhentong Wang and Xiaojing Tang + * @author Xiaojing Tang and Zhentong Wang * @version 0.1 * @version 0.3.0 */ @@ -34,12 +34,10 @@ #ifndef ANISOTROPIC_KERNEL_H #define ANISOTROPIC_KERNEL_H -#include "base_kernel_includes_nonisotropic.h" -#include "kernel_wenland_c2_anisotropic.h" +#include "base_kernel.h" namespace SPH { - /** * @class AnisotropicKernel * @brief Abstract base class of a general anisotropic SPH kernel function which @@ -90,9 +88,11 @@ class AnisotropicKernel : public KernelType this->factor_d2W_3D_ = this->factor_W_3D_; }; + /** Calculates the transform tensor form anisotropic space to isotropic space **/ Mat2d getCoordinateTransformationTensorG(Vec2d kernel_vector, Vec2d transform_vector); Mat3d getCoordinateTransformationTensorG(Vec3d kernel_vector, Vec3d transform_vector); + /** Calculates the unit vector between a pair of particles **/ virtual Vec2d e(const Real &distance, const Vec2d &displacement) const override; virtual Vec3d e(const Real &distance, const Vec3d &displacement) const override; diff --git a/tests/user_examples/extra_src/shared/anisotropic_kernel.hpp b/src/shared/kernels/anisotropic_kernel.hpp similarity index 98% rename from tests/user_examples/extra_src/shared/anisotropic_kernel.hpp rename to src/shared/kernels/anisotropic_kernel.hpp index 88e3cf4716..e0b9464c51 100644 --- a/tests/user_examples/extra_src/shared/anisotropic_kernel.hpp +++ b/src/shared/kernels/anisotropic_kernel.hpp @@ -26,8 +26,7 @@ Mat3d AnisotropicKernel::getCoordinateTransformationTensorG(Vec3d ke { Mat3d G_kernel_coordinate = Mat3d({{ Real(1.0) / (this->h_ * kernel_vector[0]), 0.0, 0.0}, {0.0, Real(1.0) / (this->h_ * kernel_vector[1]), 0.0}, - {0.0, 0.0, Real(1.0) / (this->h_ * kernel_vector[2])}} - ); + {0.0, 0.0, Real(1.0) / (this->h_ * kernel_vector[2])}} ); return G_kernel_coordinate; } //=========================================================================================// diff --git a/src/shared/kernels/base_kernel.h b/src/shared/kernels/base_kernel.h index bb632a6eba..ab9eb01f8b 100644 --- a/src/shared/kernels/base_kernel.h +++ b/src/shared/kernels/base_kernel.h @@ -87,14 +87,42 @@ class Kernel Real FactorW1D() const { return factor_W_1D_; }; Real FactorW2D() const { return factor_W_2D_; }; Real FactorW3D() const { return factor_W_3D_; }; + + /** + * unit vector pointing from j to i or inter-particle surface direction + */ + + virtual Vec2d e(const Real &distance, const Vec2d &displacement) const { return displacement / (distance + TinyReal); }; + virtual Vec3d e(const Real &distance, const Vec3d &displacement) const { return displacement / (distance + TinyReal); }; + + /** + * check if particles are within cutoff radius + */ + virtual bool checkIfWithinCutOffRadius(Vec2d displacement) + { + Real distance_metric = displacement.squaredNorm(); + if (distance_metric < CutOffRadiusSqr()) + return true; + else + return false; + }; + virtual bool checkIfWithinCutOffRadius(Vec3d displacement) + { + Real distance_metric = displacement.squaredNorm(); + if (distance_metric < CutOffRadiusSqr()) + return true; + else + return false; + }; + /** * Calculates the kernel value for the given displacement of two particles * r_ij pointing from particle j to particle i */ - Real W(const Real &r_ij, const Real &displacement) const; - Real W(const Real &r_ij, const Vec2d &displacement) const; - Real W(const Real &r_ij, const Vec3d &displacement) const; + virtual Real W(const Real &r_ij, const Real &displacement) const; + virtual Real W(const Real &r_ij, const Vec2d &displacement) const; + virtual Real W(const Real &r_ij, const Vec3d &displacement) const; /** this value could be use to calculate the value of W * they are realized in specific kernel implementations @@ -104,16 +132,16 @@ class Kernel virtual Real W_3D(const Real q) const = 0; /** Calculates the kernel value at the origin **/ - Real W0(const Real &point_i) const { return factor_W_1D_; }; - Real W0(const Vec2d &point_i) const { return factor_W_2D_; }; - Real W0(const Vec3d &point_i) const { return factor_W_3D_; }; + virtual Real W0(const Real &point_i) const { return factor_W_1D_; }; + virtual Real W0(const Vec2d &point_i) const { return factor_W_2D_; }; + virtual Real W0(const Vec3d &point_i) const { return factor_W_3D_; }; /** Calculates the kernel derivation for * the given distance of two particles */ - Real dW(const Real &r_ij, const Real &displacement) const; - Real dW(const Real &r_ij, const Vec2d &displacement) const; - Real dW(const Real &r_ij, const Vec3d &displacement) const; + virtual Real dW(const Real &r_ij, const Real &displacement) const; + virtual Real dW(const Real &r_ij, const Vec2d &displacement) const; + virtual Real dW(const Real &r_ij, const Vec3d &displacement) const; /** this value could be use to calculate the value of dW * they are realized in specific kernel implementations @@ -125,9 +153,9 @@ class Kernel /** Calculates the kernel second order derivation for * the given distance of two particles */ - Real d2W(const Real &r_ij, const Real &displacement) const; - Real d2W(const Real &r_ij, const Vec2d &displacement) const; - Real d2W(const Real &r_ij, const Vec3d &displacement) const; + virtual Real d2W(const Real &r_ij, const Real &displacement) const; + virtual Real d2W(const Real &r_ij, const Vec2d &displacement) const; + virtual Real d2W(const Real &r_ij, const Vec3d &displacement) const; /** this value could be use to calculate the value of d2W * they are realized in specific kernel implementations @@ -188,4 +216,4 @@ class Kernel void reduceTwice(); /** reduce for linear structures or filaments */ }; } // namespace SPH -#endif // BASE_KERNELS_H \ No newline at end of file +#endif // BASE_KERNELS_H diff --git a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp index 763b1e4e59..bf8182952b 100644 --- a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp @@ -79,8 +79,8 @@ void Integration1stHalfPK2::initialization(size_t index_i, Real dt) F_[index_i] += dF_dt_[index_i] * dt * 0.5; rho_[index_i] = rho0_ / F_[index_i].determinant(); // obtain the first Piola-Kirchhoff stress from the second Piola-Kirchhoff stress - // it seems using reproducing correction here increases convergence rate near the free surface - stress_PK1_B_[index_i] = elastic_solid_.StressPK1(F_[index_i], index_i) * B_[index_i]; + // it seems using reproducing correction here increases convergence rate near the free surface, note that the correction matrix is in a form of transpose + stress_PK1_B_[index_i] = elastic_solid_.StressPK1(F_[index_i], index_i) * B_[index_i].transpose(); } //=================================================================================================// Integration1stHalfKirchhoff:: diff --git a/src/shared/particle_neighborhood/neighborhood.cpp b/src/shared/particle_neighborhood/neighborhood.cpp index 2b4790e6b6..96e7433fd2 100644 --- a/src/shared/particle_neighborhood/neighborhood.cpp +++ b/src/shared/particle_neighborhood/neighborhood.cpp @@ -29,7 +29,7 @@ void NeighborBuilder::createNeighbor(Neighborhood &neighborhood, const Real &dis neighborhood.W_ij_.push_back(kernel_->W(distance, displacement)); neighborhood.dW_ijV_j_.push_back(kernel_->dW(distance, displacement) * Vol_j); neighborhood.r_ij_.push_back(distance); - neighborhood.e_ij_.push_back(displacement / (distance + TinyReal)); + neighborhood.e_ij_.push_back(kernel_->e(distance, displacement)); neighborhood.allocated_size_++; } //=================================================================================================// @@ -41,7 +41,7 @@ void NeighborBuilder::initializeNeighbor(Neighborhood &neighborhood, const Real neighborhood.W_ij_[current_size] = kernel_->W(distance, displacement); neighborhood.dW_ijV_j_[current_size] = kernel_->dW(distance, displacement) * Vol_j; neighborhood.r_ij_[current_size] = distance; - neighborhood.e_ij_[current_size] = displacement / (distance + TinyReal); + neighborhood.e_ij_[current_size] = kernel_->e(distance, displacement); } //=================================================================================================// void NeighborBuilder::createNeighbor(Neighborhood &neighborhood, const Real &distance, @@ -82,7 +82,7 @@ void NeighborBuilderInner::operator()(Neighborhood &neighborhood, size_t index_j = std::get<0>(list_data_j); Vecd displacement = pos_i - std::get<1>(list_data_j); Real distance_metric = displacement.squaredNorm(); - if (distance_metric < kernel_->CutOffRadiusSqr() && index_i != index_j) + if (kernel_->checkIfWithinCutOffRadius(displacement) && index_i != index_j) { neighborhood.current_size_ >= neighborhood.allocated_size_ ? createNeighbor(neighborhood, std::sqrt(distance_metric), displacement, index_j, std::get<2>(list_data_j)) diff --git a/tests/2d_examples/test_2d_anisotropic_beam/CMakeLists.txt b/tests/2d_examples/test_2d_anisotropic_beam/CMakeLists.txt new file mode 100644 index 0000000000..4872429c5b --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/CMakeLists.txt @@ -0,0 +1,25 @@ +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir + +set(CMAKE_VERBOSE_MAKEFILE on) + +STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} ) +PROJECT("${CURRENT_FOLDER}") + + + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) + +aux_source_directory(. DIR_SRCS) +ADD_EXECUTABLE(${PROJECT_NAME} ${DIR_SRCS}) + +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + + +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} sphinxsys_2d) diff --git a/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position.xml b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position.xml new file mode 100644 index 0000000000..aa309f10dc --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position.xml @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml new file mode 100644 index 0000000000..d4710a4b24 --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml @@ -0,0 +1,207 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_result.xml b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_result.xml new file mode 100644 index 0000000000..cdd197ce56 --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_result.xml @@ -0,0 +1,819 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_runtimes.dat b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_runtimes.dat new file mode 100644 index 0000000000..fd93545d1b --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +8 +5 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/regression_test_tool.py b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..b72295d423 --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/regression_test_tool.py @@ -0,0 +1,37 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_2d_anisotropic_beam +""" + +case_name = "test_2d_anisotropic_beam" +body_name = "BeamObserver" +parameter_name = "Position" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if converged == "true": + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp new file mode 100644 index 0000000000..39a510c19b --- /dev/null +++ b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp @@ -0,0 +1,311 @@ +/* ---------------------------------------------------------------------------* + * SPHinXsys: 2D anisotropic beam example-one body version * + * ----------------------------------------------------------------------------* + * This is a test cases using ASPH method with anisotropic kerne for * + * simulating solid. Particle space is anisotropic in different directions in this beam. * + * @author Xiaojing Tang + * ----------------------------------------------------------------------------*/ +#include "sphinxsys.h" +using namespace SPH; +//------------------------------------------------------------------------------ +// global parameters for the case +//------------------------------------------------------------------------------ +Real PL = 0.2; // beam length +Real PH = 0.02; // beam width +Real SL = 0.02; // constrained length + +// particle spacings and particle numbers +int y_num = 10; // particle number in y direction +Real ratio_ = 4.0; // anisotropic ratio, also dp_x / dp_y +Real resolution_ref = PH / y_num; // particle spacing in y direction +Real resolution_ref_large = ratio_ * resolution_ref; // large particle spacing, also the particle spacing in x direction +Real Total_PL = PL + SL; // total length +int x_num = Total_PL / resolution_ref_large; // particle number in x direction +// anisotropic parameters +Vec2d scaling_vector = Vec2d(1.0, 1.0 / ratio_); // scaling_vector for defining the anisotropic kernel +Real scaling_factor = 1.0 / ratio_; // scaling factor to calculate the time step + +//---------------------------------------------------------------------- +// particle generation considering the anisotropic resolution +//---------------------------------------------------------------------- +class AnisotropicParticleGenerator : public ParticleGenerator +{ + public: + AnisotropicParticleGenerator(SPHBody &sph_body) : ParticleGenerator(sph_body){}; + + virtual void initializeGeometricVariables() override + { + // set particles directly + for (int i = 0; i < x_num; i++) + { + for (int j = 0; j < y_num; j++) + { + Real x = -SL + (i + 0.5) * resolution_ref_large; + Real y = -PH / 2 + (j + 0.5) * resolution_ref; + initializePositionAndVolumetricMeasure(Vecd(x, y), (resolution_ref * resolution_ref_large)); + } + } + } +}; + +Real BW = resolution_ref * 4; // boundary width, at least three particles +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), + Vec2d(PL + 3.0 * BW, PL / 2.0)); +//---------------------------------------------------------------------- +// Material properties of the solid. +//---------------------------------------------------------------------- +Real rho0_s = 1.0e3; // reference density +Real Youngs_modulus = 2.0e6; // reference Youngs modulus +Real poisson = 0.3975; // Poisson ratio +//---------------------------------------------------------------------- +// Parameters for initial condition on velocity +//---------------------------------------------------------------------- +Real kl = 1.875; +Real M = sin(kl) + sinh(kl); +Real N = cos(kl) + cosh(kl); +Real Q = 2.0 * (cos(kl) * sinh(kl) - sin(kl) * cosh(kl)); +Real vf = 0.05; +Real R = PL / (0.5 * Pi); +//---------------------------------------------------------------------- +// Geometric shapes used in the system. +//---------------------------------------------------------------------- +// a beam base shape +std::vector beam_base_shape{ + Vecd(-SL, -PH / 2), Vecd(-SL, PH / 2), Vecd(0.0, PH / 2), + Vecd(0.0, -PH / 2), Vecd(-SL, -PH / 2)}; +// a beam shape +std::vector beam_shape{ + Vecd(0.0, -PH / 2), Vecd(0.0, PH / 2), Vecd(PL, PH / 2), Vecd(PL, -PH / 2), Vecd(0.0, -PH / 2)}; +// Beam observer location +StdVec observation_location = {Vecd(PL, 0.0)}; +//---------------------------------------------------------------------- +// Define the beam body +//---------------------------------------------------------------------- +class Beam : public MultiPolygonShape +{ + public: + explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(beam_base_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); + } +}; +//---------------------------------------------------------------------- +// application dependent initial condition +//---------------------------------------------------------------------- +class BeamInitialCondition + : public solid_dynamics::ElasticDynamicsInitialCondition +{ + public: + explicit BeamInitialCondition(SPHBody &sph_body) + : solid_dynamics::ElasticDynamicsInitialCondition(sph_body){}; + + void update(size_t index_i, Real dt) + { + /** initial velocity profile */ + Real x = pos_[index_i][0] / PL; + if (x > 0.0) + { + vel_[index_i][1] = vf * particles_->elastic_solid_.ReferenceSoundSpeed() * + (M * (cos(kl * x) - cosh(kl * x)) - N * (sin(kl * x) - sinh(kl * x))) / Q; + } + }; +}; +//---------------------------------------------------------------------- +// define the beam base which will be constrained. +//---------------------------------------------------------------------- +MultiPolygon createBeamConstrainShape() +{ + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(beam_base_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(beam_shape, ShapeBooleanOps::sub); + return multi_polygon; +}; +//---------------------------------------------------------------------- +// calculate correction matrix B to keep the accuracy +//---------------------------------------------------------------------- + +class AnisotropicCorrectConfiguration : public LocalDynamics, public GeneralDataDelegateInner +{ + public: + AnisotropicCorrectConfiguration(BaseInnerRelation &inner_relation, int beta = 0, Real alpha = Real(0)) + : LocalDynamics(inner_relation.getSPHBody()), + GeneralDataDelegateInner(inner_relation), + beta_(beta), alpha_(alpha), + B_(*particles_->getVariableByName("CorrectionMatrix")), + pos_(particles_->pos_) + { + particles_->registerVariable(show_neighbor_, "ShowingNeighbor", Real(0.0)); + } + virtual ~AnisotropicCorrectConfiguration(){}; + + protected: + protected: + int beta_; + Real alpha_; + StdLargeVec &B_; + StdLargeVec &pos_; + StdLargeVec show_neighbor_; + + void interaction(size_t index_i, Real dt = 0.0) + { + Matd local_configuration = Eps * Matd::Identity(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n]; + Vecd e_ij = inner_neighborhood.e_ij_[n]; + if (index_i == 67) + { + show_neighbor_[index_j] = 1.0; + }; + + Vecd gradW_ijV_j = dW_ijV_j * e_ij; + Vecd r_ji = pos_[index_i] - pos_[index_j]; + local_configuration -= r_ji * gradW_ijV_j.transpose(); + } + B_[index_i] = local_configuration; + }; + + void update(size_t index_i, Real dt) + { + Real det_sqr = pow(B_[index_i].determinant(), beta_); + Matd inverse = B_[index_i].inverse(); + B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr); + } +}; + +//------------------------------------------------------------------------------ +// the main program +//------------------------------------------------------------------------------ +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem with global controls. + //---------------------------------------------------------------------- + SPHSystem system(system_domain_bounds, resolution_ref_large); +#ifdef BOOST_AVAILABLE + // handle command line arguments + system.handleCommandlineOptions(ac, av); +#endif //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + SolidBody beam_body(system, makeShared("BeamBody")); + beam_body.sph_adaptation_->resetKernel>(scaling_vector); + beam_body.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); + beam_body.generateParticles(); + + ObserverBody beam_observer(system, "BeamObserver"); + beam_observer.sph_adaptation_->resetKernel>(scaling_vector); + beam_observer.generateParticles(observation_location); + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + //---------------------------------------------------------------------- + InnerRelation beam_body_inner(beam_body); + ContactRelation beam_observer_contact(beam_observer, {&beam_body}); + //----------------------------------------------------------------------------- + // this section define all numerical methods will be used in this case + //----------------------------------------------------------------------------- + SimpleDynamics beam_initial_velocity(beam_body); + // corrected strong configuration + InteractionWithUpdate beam_corrected_configuration(beam_body_inner); + beam_body.addBodyStateForRecording("ShowingNeighbor"); + // time step size calculation + ReduceDynamics computing_time_step_size(beam_body); + // stress relaxation for the beam + Dynamics1Level stress_relaxation_first_half(beam_body_inner); + Dynamics1Level stress_relaxation_second_half(beam_body_inner); + // clamping a solid body part. + BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); + SimpleDynamics constraint_beam_base(beam_base); + //----------------------------------------------------------------------------- + // outputs + //----------------------------------------------------------------------------- + IOEnvironment io_environment(system); + BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); + RegressionTestEnsembleAverage> + write_beam_tip_displacement("Position", io_environment, beam_observer_contact); + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + beam_initial_velocity.exec(); + beam_corrected_configuration.exec(); + //---------------------------------------------------------------------- + // Setup computing time-step controls. + //---------------------------------------------------------------------- + int ite = 0; + Real T0 = 1.0; + Real end_time = T0; + // time step size for output file + Real output_interval = 0.01 * T0; + Real Dt = 0.1 * output_interval; /**< Time period for data observing */ + Real dt = 0.0; // default acoustic time step sizes + + // statistics for computing time + TickCount t1 = TickCount::now(); + TimeInterval interval; + //----------------------------------------------------------------------------- + // from here the time stepping begins + //----------------------------------------------------------------------------- + write_beam_states.writeToFile(0); + write_beam_tip_displacement.writeToFile(0); + + // computation loop starts + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + // integrate time (loop) until the next output time + while (integration_time < output_interval) + { + + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + stress_relaxation_first_half.exec(dt); + constraint_beam_base.exec(); + stress_relaxation_second_half.exec(dt); + + ite++; + dt = scaling_factor * computing_time_step_size.exec(); + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + + if (ite % 100 == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " dt: " + << dt << "\n"; + } + } + } + + write_beam_tip_displacement.writeToFile(ite); + TickCount t2 = TickCount::now(); + write_beam_states.writeToFile(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (system.GenerateRegressionData()) + { + write_beam_tip_displacement.generateDataBase(Vec2d(1.0e-2, 1.0e-2), Vec2d(1.0e-2, 1.0e-2)); + } + else + { + write_beam_tip_displacement.testResult(); + } + + return 0; +} diff --git a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp index c90ff2545f..653f454f64 100644 --- a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp @@ -21,7 +21,7 @@ Real BW = resolution_ref * 4; // boundary width, at least three particles BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), Vec2d(PL + 3.0 * BW, PL / 2.0)); //---------------------------------------------------------------------- -// Material properties of the fluid. +// Material properties of the solid. //---------------------------------------------------------------------- Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus diff --git a/tests/user_examples/test_nonisotropic_kernel/CMakeLists.txt b/tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/CMakeLists.txt similarity index 72% rename from tests/user_examples/test_nonisotropic_kernel/CMakeLists.txt rename to tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/CMakeLists.txt index cb7a61ac7d..e6e6da0840 100644 --- a/tests/user_examples/test_nonisotropic_kernel/CMakeLists.txt +++ b/tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/CMakeLists.txt @@ -1,17 +1,16 @@ -STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} ) PROJECT("${CURRENT_FOLDER}") SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") - + aux_source_directory(. DIR_SRCS) ADD_EXECUTABLE(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS}) -target_link_libraries(${PROJECT_NAME} sphinxsys_2d GTest::gtest GTest::gtest_main) -target_link_libraries(${PROJECT_NAME} extra_sources_2d) +target_link_libraries(${PROJECT_NAME} sphinxsys_2d GTest::gtest GTest::gtest_main) set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") -add_test(NAME ${PROJECT_NAME} - COMMAND ${PROJECT_NAME} +add_test(NAME ${PROJECT_NAME}_particle_relaxation + COMMAND ${PROJECT_NAME} --r=true WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) diff --git a/tests/user_examples/test_nonisotropic_kernel/test_nonisotropic_kernel.cpp b/tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/test_anisotropic_kernel.cpp similarity index 96% rename from tests/user_examples/test_nonisotropic_kernel/test_nonisotropic_kernel.cpp rename to tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/test_anisotropic_kernel.cpp index f01a755af0..c31e5ba333 100644 --- a/tests/user_examples/test_nonisotropic_kernel/test_nonisotropic_kernel.cpp +++ b/tests/unit_tests_src/shared/test_kernels/test_anisotropic_kernel/test_anisotropic_kernel.cpp @@ -1,6 +1,4 @@ #include "anisotropic_kernel.hpp" -#include "base_kernel_includes_nonisotropic.h" -#include "kernel_wenland_c2_anisotropic.h" #include "sphinxsys.h" #include @@ -32,7 +30,7 @@ TEST(test_anisotropic_kernel, test_Laplacian) int x_num = PL / resolution_x; // Particle number in x direction , the same as particle number in y direction Vecd scaling_vector(1.0, 1.0 / ratio); // in x and y directions - AnisotropicKernel + AnisotropicKernel wendland(1.15 * resolution_x, scaling_vector, Vecd(0.0, 0.0)); // no rotation introduced Mat2d transform_tensor = wendland.getCoordinateTransformationTensorG(scaling_vector, Vecd(0.0, 0.0)); // tensor diff --git a/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.cpp b/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.cpp deleted file mode 100644 index 110baf3148..0000000000 --- a/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.cpp +++ /dev/null @@ -1,204 +0,0 @@ -/** - * @file base_kernel_includes_nonisotropic.cpp - * @author Luhui Han, Chi Zhang and Xiangyu Hu - */ - -#include "base_kernel_includes_nonisotropic.h" - -namespace SPH -{ -namespace Anisotropic -{ -//=================================================================================================// -Kernel::Kernel(Real h, Real kernel_size, Real truncation, const std::string &kernel_name) - : kernel_name_(kernel_name), h_(h), inv_h_(1.0 / h), - kernel_size_(kernel_size), - truncation_(truncation), rc_ref_(truncation * h), - rc_ref_sqr_(rc_ref_ * rc_ref_), - h_factor_W_1D_(std::bind(&Kernel::factorW1D, this, _1)), - h_factor_W_2D_(std::bind(&Kernel::factorW2D, this, _1)), - h_factor_W_3D_(std::bind(&Kernel::factorW3D, this, _1)), - h_factor_dW_1D_(std::bind(&Kernel::factordW1D, this, _1)), - h_factor_dW_2D_(std::bind(&Kernel::factordW2D, this, _1)), - h_factor_dW_3D_(std::bind(&Kernel::factordW3D, this, _1)), - h_factor_d2W_1D_(std::bind(&Kernel::factord2W1D, this, _1)), - h_factor_d2W_2D_(std::bind(&Kernel::factord2W2D, this, _1)), - h_factor_d2W_3D_(std::bind(&Kernel::factord2W3D, this, _1)){}; -//=================================================================================================// -void Kernel::setDerivativeParameters() -{ - factor_dW_1D_ = inv_h_ * factor_W_1D_; - factor_dW_2D_ = inv_h_ * factor_W_2D_; - factor_dW_3D_ = inv_h_ * factor_W_3D_; - factor_d2W_1D_ = inv_h_ * factor_dW_1D_; - factor_d2W_2D_ = inv_h_ * factor_dW_2D_; - factor_d2W_3D_ = inv_h_ * factor_dW_3D_; -} -//=================================================================================================// -void Kernel::resetSmoothingLength(Real h) -{ - Real ratio = h_ / h; - h_ /= ratio; - inv_h_ *= ratio; - rc_ref_ /= ratio; - rc_ref_sqr_ /= ratio * ratio; - - factor_W_1D_ *= ratio; - factor_W_2D_ *= ratio * ratio; - factor_W_3D_ *= ratio * ratio * ratio; - - setDerivativeParameters(); -} -//=================================================================================================// -Real Kernel::W(const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_W_1D_ * W_1D(q); -} -//=================================================================================================// -Real Kernel::W(const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_W_2D_ * W_2D(q); -} -//=================================================================================================// -Real Kernel::W(const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_W_3D_ * W_3D(q); -} -//=================================================================================================// -Real Kernel::dW(const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_dW_1D_ * dW_1D(q); -} -//=================================================================================================// -Real Kernel::dW(const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_dW_2D_ * dW_2D(q); -} -//=================================================================================================// -Real Kernel::dW(const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_dW_3D_ * dW_3D(q); -} -//=================================================================================================// -Real Kernel::d2W(const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_d2W_1D_ * d2W_1D(q); -} -//=================================================================================================// -Real Kernel::d2W(const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_d2W_2D_ * d2W_2D(q); -} -//=================================================================================================// -Real Kernel::d2W(const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_; - return factor_d2W_3D_ * d2W_3D(q); -} -//=================================================================================================// -Real Kernel::W(const Real &h_ratio, const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_W_1D_ * W_1D(q) * h_factor_W_1D_(h_ratio); -} -//=================================================================================================// -Real Kernel::W(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_W_2D_ * W_2D(q) * h_factor_W_2D_(h_ratio); -} -//=================================================================================================// -Real Kernel::W(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_W_3D_ * W_3D(q) * h_factor_W_3D_(h_ratio); -} -//=================================================================================================// -Real Kernel::W0(const Real &h_ratio, const Real &point_i) const -{ - return factor_W_1D_ * h_factor_W_1D_(h_ratio); -}; -//=================================================================================================// -Real Kernel::W0(const Real &h_ratio, const Vec2d &point_i) const -{ - return factor_W_2D_ * h_factor_W_2D_(h_ratio); -}; -//=================================================================================================// -Real Kernel::W0(const Real &h_ratio, const Vec3d &point_i) const -{ - return factor_W_3D_ * h_factor_W_3D_(h_ratio); -}; -//=================================================================================================// -Real Kernel::dW(const Real &h_ratio, const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_dW_1D_ * dW_1D(q) * h_factor_dW_1D_(h_ratio); -} -//=================================================================================================// -Real Kernel::dW(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_dW_2D_ * dW_2D(q) * h_factor_dW_2D_(h_ratio); -} -//=================================================================================================// -Real Kernel::dW(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_dW_3D_ * dW_3D(q) * h_factor_dW_1D_(h_ratio); -} -//=================================================================================================// -Real Kernel::d2W(const Real &h_ratio, const Real &r_ij, const Real &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_d2W_1D_ * d2W_1D(q) * h_factor_d2W_1D_(h_ratio); -} -//=================================================================================================// -Real Kernel::d2W(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_d2W_2D_ * d2W_2D(q) * h_factor_d2W_2D_(h_ratio); -} -//=================================================================================================// -Real Kernel::d2W(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const -{ - Real q = r_ij * inv_h_ * h_ratio; - return factor_d2W_3D_ * d2W_3D(q) * h_factor_d2W_3D_(h_ratio); -} -//=================================================================================================// -void Kernel::reduceOnce() -{ - factor_W_3D_ = factor_W_2D_; - factor_W_2D_ = factor_W_1D_; - factor_W_1D_ = 0.0; - setDerivativeParameters(); - - h_factor_W_3D_ = std::bind(&Kernel::factorW2D, this, _1); - h_factor_W_2D_ = std::bind(&Kernel::factorW1D, this, _1); - h_factor_dW_3D_ = std::bind(&Kernel::factordW2D, this, _1); - h_factor_dW_2D_ = std::bind(&Kernel::factordW1D, this, _1); - h_factor_d2W_3D_ = std::bind(&Kernel::factord2W2D, this, _1); - h_factor_d2W_2D_ = std::bind(&Kernel::factord2W1D, this, _1); -} -//=================================================================================================// -void Kernel::reduceTwice() -{ - factor_W_3D_ = factor_W_1D_; - factor_W_2D_ = 0.0; - factor_W_1D_ = 0.0; - setDerivativeParameters(); - - h_factor_W_3D_ = std::bind(&Kernel::factorW1D, this, _1); - h_factor_dW_3D_ = std::bind(&Kernel::factordW1D, this, _1); - h_factor_d2W_3D_ = std::bind(&Kernel::factord2W1D, this, _1); -} - -} // namespace Anisotropic -} // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.h b/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.h deleted file mode 100644 index 5ffbdf434c..0000000000 --- a/tests/user_examples/extra_src/shared/base_kernel_includes_nonisotropic.h +++ /dev/null @@ -1,221 +0,0 @@ -/* -----------------------------------------------------------------------------* - * SPHinXsys * - * -----------------------------------------------------------------------------* - * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * - * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * - * physical accurate simulation and aims to model coupled industrial dynamic * - * systems including fluid, solid, multi-body dynamics and beyond with SPH * - * (smoothed particle hydrodynamics), a meshless computational method using * - * particle discretization. * - * * - * SPHinXsys is partially funded by German Research Foundation * - * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * - * HU1527/12-1 and HU1527/12-4. * - * * - * Portions copyright (c) 2017-2023 Technical University of Munich and * - * the authors' affiliations. * - * * - * Licensed under the Apache License, Version 2.0 (the "License"); you may * - * not use this file except in compliance with the License. You may obtain a * - * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * - * * - * -----------------------------------------------------------------------------*/ -/** - * @file base_kernel_includes_nonisotropic.h - * @brief This is the base classes of kernel functions. Implementation will be - * implemented in derived classes. The kernal function define the relevance - * between two neighboring particles. Basically, the further the two - * particles, the less relevance they have. - * @author Luhui Han, Chi Zhang and Xiangyu Hu - * @version 0.1 - * @version 0.3.0 - * Add the reduced kernel for reduced dynamics of linear structure. - * -- Chi ZHANG - */ - -#ifndef BASE_KERNELS_ANISOTROPIC_H -#define BASE_KERNELS_ANISOTROPIC_H - -#include "base_data_package.h" - -#include -#include - -using namespace std::placeholders; - -namespace SPH -{ -namespace Anisotropic -{ -/** - * @class Kernel - * @brief Abstract base class of a general SPH kernel function which - * is a smoothed Dirac delta function, - * a kernel function is radial symmetric, and has a scaling factor. - * Based on difference data type in 2d or 3d buildings, - * the kernel is defined for 2 and 3 dimensions. - * The kernel gives value one at the origin. - * The naming of kernel function follows the stand SPH literature. - * Currently, only constant smoothing length is applied. - * Basically, one can assign different kernel for different particle interactions. - */ -class Kernel -{ - protected: - std::string kernel_name_; - Real h_, inv_h_; /**< reference smoothing length and its inverse **/ - - Real kernel_size_; /** FactorFunctor; - FactorFunctor h_factor_W_1D_, h_factor_W_2D_, h_factor_W_3D_; - FactorFunctor h_factor_dW_1D_, h_factor_dW_2D_, h_factor_dW_3D_; - FactorFunctor h_factor_d2W_1D_, h_factor_d2W_2D_, h_factor_d2W_3D_; - - Real factorW1D(const Real &h_ratio) const { return h_ratio; }; - Real factorW2D(const Real &h_ratio) const { return h_ratio * h_ratio; }; - Real factorW3D(const Real &h_ratio) const { return h_ratio * h_ratio * h_ratio; }; - Real factordW1D(const Real &h_ratio) const { return factorW1D(h_ratio) * h_ratio; }; - Real factordW2D(const Real &h_ratio) const { return factorW2D(h_ratio) * h_ratio; }; - Real factordW3D(const Real &h_ratio) const { return factorW3D(h_ratio) * h_ratio; }; - Real factord2W1D(const Real &h_ratio) const { return factordW1D(h_ratio) * h_ratio; }; - Real factord2W2D(const Real &h_ratio) const { return factordW2D(h_ratio) * h_ratio; }; - Real factord2W3D(const Real &h_ratio) const { return factordW3D(h_ratio) * h_ratio; }; - - public: - Real CutOffRadius(Real h_ratio) const { return rc_ref_ / h_ratio; }; - Real CutOffRadiusSqr(Real h_ratio) const { return rc_ref_sqr_ / (h_ratio * h_ratio); }; - - Real W(const Real &h_ratio, const Real &r_ij, const Real &displacement) const; - Real W(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const; - Real W(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const; - - /** Calculates the kernel value at the origin **/ - Real W0(const Real &h_ratio, const Real &point_i) const; - Real W0(const Real &h_ratio, const Vec2d &point_i) const; - Real W0(const Real &h_ratio, const Vec3d &point_i) const; - - /** Calculates the kernel derivation for the given distance of two particles **/ - Real dW(const Real &h_ratio, const Real &r_ij, const Real &displacement) const; - Real dW(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const; - Real dW(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const; - - /** Calculates the kernel second order derivation for the given distance of two particles **/ - Real d2W(const Real &h_ratio, const Real &r_ij, const Real &displacement) const; - Real d2W(const Real &h_ratio, const Real &r_ij, const Vec2d &displacement) const; - Real d2W(const Real &h_ratio, const Real &r_ij, const Vec3d &displacement) const; - //---------------------------------------------------------------------- - // Below are for reduced kernels. - //---------------------------------------------------------------------- - public: - void reduceOnce(); /** reduce for thin structures or films */ - void reduceTwice(); /** reduce for linear structures or filaments */ -}; -} // namespace Anisotropic -} // namespace SPH - -#endif // BASE_KERNELS_ANISOTROPIC_H \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.cpp b/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.cpp deleted file mode 100644 index f3035f3b31..0000000000 --- a/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.cpp +++ /dev/null @@ -1,70 +0,0 @@ -/** - * @file kernel_wenland.cpp - * @author Luhui Han, Chi Zhang, Yongchuan Yu and Xiangyu Hu - */ -#include "kernel_wenland_c2_anisotropic.h" - -#include - -namespace SPH -{ - -namespace Anisotropic -{ -//=================================================================================================// -KernelWendlandC2::KernelWendlandC2(Real h) - : Anisotropic::Kernel(h, 2.0, 2.0, "Wendland2CKernel") -{ - factor_W_1D_ = inv_h_ * 3.0 / 4.0; - factor_W_2D_ = inv_h_ * inv_h_ * 7.0 / (4.0 * Pi); - factor_W_3D_ = inv_h_ * inv_h_ * inv_h_ * 21.0 / (16.0 * Pi); - setDerivativeParameters(); -} -//=================================================================================================// -Real KernelWendlandC2::W_1D(const Real q) const -{ - return pow(1.0 - 0.5 * q, 4) * (1.0 + 2.0 * q); -} -//=================================================================================================// -Real KernelWendlandC2::W_2D(const Real q) const -{ - return W_1D(q); -} -//=================================================================================================// -Real KernelWendlandC2::W_3D(const Real q) const -{ - return W_2D(q); -} -//=================================================================================================// -Real KernelWendlandC2::dW_1D(const Real q) const -{ - return 0.625 * pow(q - 2.0, 3) * q; -} -//=================================================================================================// -Real KernelWendlandC2::dW_2D(const Real q) const -{ - return dW_1D(q); -} -//=================================================================================================// -Real KernelWendlandC2::dW_3D(const Real q) const -{ - return dW_2D(q); -} -//=================================================================================================// -Real KernelWendlandC2::d2W_1D(const Real q) const -{ - return 1.25 * pow(q - 2.0, 2) * (2.0 * q - 1.0); -} -//=================================================================================================// -Real KernelWendlandC2::d2W_2D(const Real q) const -{ - return d2W_1D(q); -} -//=================================================================================================// -Real KernelWendlandC2::d2W_3D(const Real q) const -{ - return d2W_2D(q); -} -//=================================================================================================// -} // namespace Anisotropic -} // namespace SPH diff --git a/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.h b/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.h deleted file mode 100644 index 8472d04d81..0000000000 --- a/tests/user_examples/extra_src/shared/kernel_wenland_c2_anisotropic.h +++ /dev/null @@ -1,66 +0,0 @@ -/* -------------------------------------------------------------------------* - * SPHinXsys * - * --------------------------------------------------------------------------* - * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * - * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * - * physical accurate simulation and aims to model coupled industrial dynamic * - * systems including fluid, solid, multi-body dynamics and beyond with SPH * - * (smoothed particle hydrodynamics), a meshless computational method using * - * particle discretization. * - * * - * SPHinXsys is partially funded by German Research Foundation * - * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1 * - * and HU1527/12-1. * - * * - * Portions copyright (c) 2017-2020 Technical University of Munich and * - * the authors' affiliations. * - * * - * Licensed under the Apache License, Version 2.0 (the "License"); you may * - * not use this file except in compliance with the License. You may obtain a * - * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * - * * - * --------------------------------------------------------------------------*/ -/** - * @file kernel_wenland_c2_anisotropic.h - * @brief This is the class for Wenland kernel. - * @details NThis kernel has compact support of 2h. - * The smoothing length h can be variable when variable h functions are applied. - * @author Luhui Han, Chi Zhang and Xiangyu Hu - */ - -#ifndef KERNEL_WENLAND_C2_Non_H -#define KERNEL_WENLAND_C2_Non_H - -#include "base_kernel_includes_nonisotropic.h" - -namespace SPH -{ -namespace Anisotropic -{ -/** - * @class KernelWendlandC2 - * @brief Kernel WendlandC2 - */ -class KernelWendlandC2 : public Anisotropic::Kernel -{ - public: - explicit KernelWendlandC2(Real h); - - /** Calculates the kernel value for - the given distance of two particles */ - virtual Real W_1D(const Real q) const override; - virtual Real W_2D(const Real q) const override; - virtual Real W_3D(const Real q) const override; - - virtual Real dW_1D(const Real q) const override; - virtual Real dW_2D(const Real q) const override; - virtual Real dW_3D(const Real q) const override; - - virtual Real d2W_1D(const Real q) const override; - virtual Real d2W_2D(const Real q) const override; - virtual Real d2W_3D(const Real q) const override; -}; -} // namespace Anisotropic - -} // namespace SPH -#endif // KERNEL_WENLAND_C2_H