Snippets
Created by
aokomoriuta
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 | #define PARTICLE_SIMULATOR_TWO_DIMENSION
#pragma warning(push, 0)
#include <particle_simulator.hpp>
#pragma warning(pop)
struct ParticleNumberDensity
{
double val;
ParticleNumberDensity()
: val(0)
{}
ParticleNumberDensity(const double v)
: val(v)
{}
void clear()
{
val = 0;
}
};
// 粒子
class Particle final
{
private:
// 位置ベクトル
PS::F64vec x;
// 粒子数密度
ParticleNumberDensity n;
public:
static constexpr double r_e = 2.4;
Particle()
{}
Particle(const Particle& src)
: x(src.x),
n(src.n)
{}
Particle(Particle&& src)
: x(std::move(src.x)),
n(src.n)
{}
Particle& operator=(const Particle& src)
{
Copy(src);
return *this;
}
Particle& operator=(Particle&& src)
{
this->x = std::move(src.x);
this->n = src.n;
return *this;
}
// 距離から重み関数を計算する
// @param r 距離
static double W(const double r)
{
// 影響半径内ならr_e/r-1を返す(ただし距離0の場合は0)
return ((0 < r) && (r < r_e)) ? (r_e / r - 1) : 0;
}
void Copy(const Particle& src)
{
this->x = src.x;
this->n = src.n;
}
void copyFromFP(const Particle& src)
{
Copy(src);
}
////////////////
// プロパティ //
////////////////
// 位置ベクトル
const auto& X() const
{
return x;
}
auto& X()
{
return x;
}
auto getPos() const
{
return X();
}
void setPos(const decltype(x)& src)
{
x = src;
}
// 粒子数密度
const auto& N() const
{
return n;
}
auto& N()
{
return n;
}
void copyFromForce(const ParticleNumberDensity& src)
{
n.val = src.val;
}
auto getRSearch() const
{
return r_e;
}
};
// エントリポイント
int main(int argc, char *argv[])
{
const double l_0 = 1.0;
ParticleSimulator::Initialize(argc, argv);
ParticleSimulator::ParticleSystem<Particle> particles;
particles.initialize();
const int L = 9;
const int H = 9;
const auto n = L*H;
particles.setNumberOfParticleLocal(n);
// 水を追加
for(int i = 0; i < L; i++)
{
for(int j = 0; j < H; j++)
{
const auto idx = i*H + j;
particles[idx].X().x = i*l_0;
particles[idx].X().y = j*l_0;
}
}
ParticleSimulator::DomainInfo domain;
domain.initialize();
domain.setBoundaryCondition(ParticleSimulator::BOUNDARY_CONDITION_OPEN);
domain.setPosRootDomain(ParticleSimulator::F64vec(0, 0), ParticleSimulator::F64vec(10, 10)); // 開放条件の場合呼ぶ必要はないが後のことも考えて一応
ParticleSimulator::TreeForForceShort<ParticleNumberDensity, Particle, Particle>::Gather particleNumberDensity;
particleNumberDensity.initialize(n);
domain.decomposeDomainAll(particles);
particles.exchangeParticle(domain);
particleNumberDensity.calcForceAllAndWriteBack([](
const Particle particle_i[], const std::int32_t n_i,
const Particle particle_j[], const std::int32_t n_j,
ParticleNumberDensity pnd[])
{
for(auto i = decltype(n_i)(0); i < n_i; i++)
{
double sum = 0;
for(auto j = decltype(n_j)(0); j < n_j; j++)
{
const auto dr = (particle_j[j].X() - particle_i[i].X());
const auto r = std::sqrt(dr*dr);
sum += Particle::W(r);
}
pnd[i] = sum;
}
}, particles, domain);
for(auto i = decltype(n)(0); i < n; i++)
{
std::cout << i << ", " <<
particles[i].X().x << ", " << particles[i].X().y << ", " <<
particles[i].N().val << std::endl;
}
PS::Finalize();
// 終了
std::cout << "finished" << std::endl;
return 0;
}
|
Comments (0)
You can clone a snippet to your computer for local editing. Learn more.