Skip to content

Commit

Permalink
Merge branch 'thermal_optimization_clean' of https://github.com/Bo-Zh…
Browse files Browse the repository at this point in the history
…ang1995/SPHinXsys-Bo into thermal_optimization_clean
  • Loading branch information
Bo-Zhang1995 committed Sep 7, 2023
2 parents ee37623 + f107bcb commit 8feb1e3
Show file tree
Hide file tree
Showing 50 changed files with 2,430 additions and 648 deletions.
2 changes: 2 additions & 0 deletions src/shared/kernels/all_kernels.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
//=========================================================================================//
Expand Down
54 changes: 41 additions & 13 deletions src/shared/kernels/base_kernel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -188,4 +216,4 @@ class Kernel
void reduceTwice(); /** reduce for linear structures or filaments */
};
} // namespace SPH
#endif // BASE_KERNELS_H
#endif // BASE_KERNELS_H
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,10 @@ namespace SPH
{
//=================================================================================================//
CorrectedConfigurationInner::
CorrectedConfigurationInner(BaseInnerRelation &inner_relation, int beta, Real alpha)
CorrectedConfigurationInner(BaseInnerRelation &inner_relation, Real alpha)
: LocalDynamics(inner_relation.getSPHBody()),
GeneralDataDelegateInner(inner_relation),
beta_(beta), alpha_(alpha),
alpha_(alpha),
B_(*particles_->getVariableByName<Matd>("CorrectionMatrix")) {}
//=================================================================================================//
void CorrectedConfigurationInner::interaction(size_t index_i, Real dt)
Expand All @@ -22,20 +22,21 @@ void CorrectedConfigurationInner::interaction(size_t index_i, Real dt)
Vecd r_ji = inner_neighborhood.r_ij_[n] * inner_neighborhood.e_ij_[n];
local_configuration -= r_ji * gradW_ij.transpose();
}

B_[index_i] = local_configuration;
}
//=================================================================================================//
void CorrectedConfigurationInner::update(size_t index_i, Real dt)
{
Real det_sqr = pow(B_[index_i].determinant(), beta_);
Real det_sqr = alpha_;
Matd inverse = B_[index_i].inverse();
B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr);
Real weight1_ = B_[index_i].determinant() / (B_[index_i].determinant() + det_sqr);
Real weight2_ = det_sqr / (B_[index_i].determinant() + det_sqr);
B_[index_i] = weight1_ * inverse + weight2_ * Matd::Identity();
}
//=================================================================================================//
CorrectedConfigurationComplex::
CorrectedConfigurationComplex(ComplexRelation &complex_relation, int beta, Real alpha)
: CorrectedConfigurationInner(complex_relation.getInnerRelation(), beta, alpha),
CorrectedConfigurationComplex(ComplexRelation &complex_relation,Real alpha)
: CorrectedConfigurationInner(complex_relation.getInnerRelation(),alpha),
GeneralDataDelegateContactOnly(complex_relation.getContactRelation())
{
for (size_t k = 0; k != contact_particles_.size(); ++k)
Expand Down Expand Up @@ -64,4 +65,4 @@ void CorrectedConfigurationComplex::interaction(size_t index_i, Real dt)
}
//=================================================================================================//
} // namespace SPH
//=================================================================================================//
//=================================================================================================//
Original file line number Diff line number Diff line change
Expand Up @@ -36,11 +36,10 @@ namespace SPH
class CorrectedConfigurationInner : public LocalDynamics, public GeneralDataDelegateInner
{
public:
CorrectedConfigurationInner(BaseInnerRelation &inner_relation, int beta = 0, Real alpha = Real(0));
CorrectedConfigurationInner(BaseInnerRelation &inner_relation, Real alpha = Real(0));
virtual ~CorrectedConfigurationInner(){};

protected:
int beta_;
Real alpha_;
StdLargeVec<Matd> &B_;

Expand All @@ -51,7 +50,7 @@ class CorrectedConfigurationInner : public LocalDynamics, public GeneralDataDele
class CorrectedConfigurationComplex : public CorrectedConfigurationInner, public GeneralDataDelegateContactOnly
{
public:
CorrectedConfigurationComplex(ComplexRelation &complex_relation, int beta = 0, Real alpha = Real(0));
CorrectedConfigurationComplex(ComplexRelation &complex_relation,Real alpha = Real(0));
virtual ~CorrectedConfigurationComplex(){};

protected:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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::
Expand Down
6 changes: 3 additions & 3 deletions src/shared/particle_neighborhood/neighborhood.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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_++;
}
//=================================================================================================//
Expand All @@ -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,
Expand Down Expand Up @@ -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))
Expand Down
25 changes: 25 additions & 0 deletions tests/2d_examples/test_2d_anisotropic_beam/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
Original file line number Diff line number Diff line change
@@ -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>
Loading

0 comments on commit 8feb1e3

Please sign in to comment.