add nonisotropic kernel documents and cases

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

tests/user_examples/extra_src/shared/anisotropic_kernel.h → src/shared/kernels/anisotropic_kernel.h

@@ -26,20 +26,18 @@
  *			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
  */
 
 #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;
 

tests/user_examples/extra_src/shared/anisotropic_kernel.hpp → src/shared/kernels/anisotropic_kernel.hpp

@@ -26,8 +26,7 @@ Mat3d AnisotropicKernel<KernelType>::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;
 }
 //=========================================================================================//

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
+#endif // BASE_KERNELS_H

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::

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))

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)

tests/2d_examples/test_2d_anisotropic_beam/regression_test_tool/BeamObserver_Position.xml

@@ -0,0 +1,103 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<Position>
+    <Snapshot_0 Position_0="~[0.186832,-8.93823e-18]" />
+    <Snapshot_440 Position_0="~[0.185054,0.0256427]" />
+    <Snapshot_880 Position_0="~[0.179374,0.0489783]" />
+    <Snapshot_1311 Position_0="~[0.171672,0.0682099]" />
+    <Snapshot_1741 Position_0="~[0.163864,0.0827785]" />
+    <Snapshot_2168 Position_0="~[0.157772,0.0919152]" />
+    <Snapshot_2588 Position_0="~[0.15485,0.0958623]" />
+    <Snapshot_3008 Position_0="~[0.155972,0.0941709]" />
+    <Snapshot_3431 Position_0="~[0.160793,0.0876134]" />
+    <Snapshot_3861 Position_0="~[0.168152,0.0752882]" />
+    <Snapshot_4291 Position_0="~[0.176084,0.0583281]" />
+    <Snapshot_4728 Position_0="~[0.182717,0.0363509]" />
+    <Snapshot_5168 Position_0="~[0.186671,0.0114759]" />
+    <Snapshot_5608 Position_0="~[0.186436,-0.0142909]" />
+    <Snapshot_6048 Position_0="~[0.182134,-0.0388313]" />
+    <Snapshot_6484 Position_0="~[0.175303,-0.0603063]" />
+    <Snapshot_6914 Position_0="~[0.167296,-0.0768292]" />
+    <Snapshot_7344 Position_0="~[0.160114,-0.088551]" />
+    <Snapshot_7766 Position_0="~[0.155761,-0.0944869]" />
+    <Snapshot_8186 Position_0="~[0.155131,-0.0955362]" />
+    <Snapshot_8606 Position_0="~[0.158366,-0.0910661]" />
+    <Snapshot_9035 Position_0="~[0.164702,-0.081337]" />
+    <Snapshot_9465 Position_0="~[0.172624,-0.0661965]" />
+    <Snapshot_9897 Position_0="~[0.18019,-0.0464129]" />
+    <Snapshot_10337 Position_0="~[0.185459,-0.0227805]" />
+    <Snapshot_10777 Position_0="~[0.186875,0.00286214]" />
+    <Snapshot_11217 Position_0="~[0.184558,0.0283795]" />
+    <Snapshot_11656 Position_0="~[0.178566,0.0513461]" />
+    <Snapshot_12086 Position_0="~[0.170834,0.0699838]" />
+    <Snapshot_12516 Position_0="~[0.163201,0.0838515]" />
+    <Snapshot_12943 Position_0="~[0.157426,0.092405]" />
+    <Snapshot_13363 Position_0="~[0.154916,0.095717]" />
+    <Snapshot_13783 Position_0="~[0.156471,0.0935144]" />
+    <Snapshot_14208 Position_0="~[0.161708,0.0862054]" />
+    <Snapshot_14638 Position_0="~[0.169221,0.0733175]" />
+    <Snapshot_15068 Position_0="~[0.177024,0.0558139]" />
+    <Snapshot_15506 Position_0="~[0.18338,0.0335277]" />
+    <Snapshot_15946 Position_0="~[0.186837,0.00843373]" />
+    <Snapshot_16386 Position_0="~[0.186109,-0.0173243]" />
+    <Snapshot_16826 Position_0="~[0.181477,-0.0415059]" />
+    <Snapshot_17259 Position_0="~[0.174495,-0.0622211]" />
+    <Snapshot_17689 Position_0="~[0.16652,-0.0782383]" />
+    <Snapshot_18119 Position_0="~[0.159549,-0.0893291]" />
+    <Snapshot_18540 Position_0="~[0.155609,-0.0946944]" />
+    <Snapshot_18960 Position_0="~[0.155499,-0.0950447]" />
+    <Snapshot_19380 Position_0="~[0.159132,-0.0900016]" />
+    <Snapshot_19810 Position_0="~[0.165675,-0.0796811]" />
+    <Snapshot_20240 Position_0="~[0.173582,-0.0640518]" />
+    <Snapshot_20673 Position_0="~[0.180962,-0.0437496]" />
+    <Snapshot_21113 Position_0="~[0.18584,-0.0197787]" />
+    <Snapshot_21553 Position_0="~[0.186834,0.00581744]" />
+    <Snapshot_21993 Position_0="~[0.184042,0.0310311]" />
+    <Snapshot_22431 Position_0="~[0.177788,0.0535727]" />
+    <Snapshot_22861 Position_0="~[0.16997,0.0717138]" />
+    <Snapshot_23291 Position_0="~[0.162501,0.084926]" />
+    <Snapshot_23717 Position_0="~[0.157139,0.0928228]" />
+    <Snapshot_24137 Position_0="~[0.155095,0.0954915]" />
+    <Snapshot_24557 Position_0="~[0.157029,0.0927695]" />
+    <Snapshot_24983 Position_0="~[0.162575,0.0848013]" />
+    <Snapshot_25413 Position_0="~[0.170239,0.0713395]" />
+    <Snapshot_25843 Position_0="~[0.177943,0.0532723]" />
+    <Snapshot_26281 Position_0="~[0.183975,0.0307388]" />
+    <Snapshot_26721 Position_0="~[0.186947,0.00555204]" />
+    <Snapshot_27161 Position_0="~[0.18576,-0.0201512]" />
+    <Snapshot_27601 Position_0="~[0.180791,-0.0440362]" />
+    <Snapshot_28033 Position_0="~[0.173655,-0.0641199]" />
+    <Snapshot_28463 Position_0="~[0.16579,-0.0795469]" />
+    <Snapshot_28893 Position_0="~[0.159081,-0.0900057]" />
+    <Snapshot_29313 Position_0="~[0.155501,-0.0948223]" />
+    <Snapshot_29733 Position_0="~[0.155856,-0.0945236]" />
+    <Snapshot_30155 Position_0="~[0.159944,-0.088842]" />
+    <Snapshot_30585 Position_0="~[0.166713,-0.0779041]" />
+    <Snapshot_31015 Position_0="~[0.17455,-0.0618028]" />
+    <Snapshot_31449 Position_0="~[0.181685,-0.0410525]" />
+    <Snapshot_31889 Position_0="~[0.186171,-0.0167611]" />
+    <Snapshot_32329 Position_0="~[0.186719,0.00884826]" />
+    <Snapshot_32769 Position_0="~[0.183472,0.0337216]" />
+    <Snapshot_33206 Position_0="~[0.17701,0.055739]" />
+    <Snapshot_33636 Position_0="~[0.169131,0.0733605]" />
+    <Snapshot_34066 Position_0="~[0.161798,0.0859649]" />
+    <Snapshot_34491 Position_0="~[0.156849,0.0931979]" />
+    <Snapshot_34911 Position_0="~[0.155342,0.0951951]" />
+    <Snapshot_35331 Position_0="~[0.157688,0.0919067]" />
+    <Snapshot_35758 Position_0="~[0.163477,0.0833192]" />
+    <Snapshot_36188 Position_0="~[0.171225,0.0693064]" />
+    <Snapshot_36618 Position_0="~[0.178834,0.0506689]" />
+    <Snapshot_37057 Position_0="~[0.18455,0.0277726]" />
+    <Snapshot_37497 Position_0="~[0.187012,0.0025474]" />
+    <Snapshot_37937 Position_0="~[0.185369,-0.0229957]" />
+    <Snapshot_38377 Position_0="~[0.180053,-0.0465718]" />
+    <Snapshot_38808 Position_0="~[0.172753,-0.0660685]" />
+    <Snapshot_39238 Position_0="~[0.165026,-0.0808563]" />
+    <Snapshot_39667 Position_0="~[0.158681,-0.0906106]" />
+    <Snapshot_40087 Position_0="~[0.155483,-0.0948597]" />
+    <Snapshot_40507 Position_0="~[0.156252,-0.0939345]" />
+    <Snapshot_40930 Position_0="~[0.160756,-0.08762]" />
+    <Snapshot_41360 Position_0="~[0.167775,-0.0760275]" />
+    <Snapshot_41790 Position_0="~[0.175546,-0.0594198]" />
+    <Snapshot_42225 Position_0="~[0.182381,-0.038272]" />
+    <Snapshot_42665 Position_0="~[0.186445,-0.0137439]" />
+</Position>