原始人のプログラミング日記

println(++C++); // println(C + 1);, C += 2;

C++でスライムを作ってみた

結論から言うと

結果


f:id:kawazuini:20170717205416g:plain
わー!すごーい!

ソースコード

PhysicalSlime.h

/**
 * @file   PhysicalSlime.h
 * @brief  PhysicalSlime
 * @author Maeda Takumi
 */
#ifndef PHYSICALSLIME_H
#define PHYSICALSLIME_H

#include "main.h"

class PhysicalPoint;

/**
 * @brief  \~english  physical slime
 * @brief  \~japanese 物理的スライム
 * @author \~ Maeda Tkumi
 */
class PhysicalSlime : public KDrawer, private KUpdater {
private:
    /* 中心座標     */ KVector mPosition;
    /* 半径         */ float mRadius;
    /* 水平分割数   */ int mStack;
    /* 垂直分割数   */ int mSlice;
    /* 有効頂点数   */ int mVertexSize;

    /* 質量         */ float mMass;
    /* 頂点毎の質量 */ float mPointMass;

    /* 頂点配列     */ Vector<Vector<KVector>> mVertex;
    /* 法線配列     */ Vector<Vector<KVector>> mNormal;
    /* 頂点毎の質点 */ Vector<Vector<PhysicalPoint*>> mPoints;

    /* 理想体積     */ float mIdealVolume;
    /* 粘度         */ float mViscosity;

    // 体積を計算します。
    float calcVolume() const;
    // 頂点配列のみを移動します(重心の移動に使用します)。
    void translateVertex(const KVector& aPosition);
    // 有効な頂点かどうかを返します。
    bool isValidIndex(const int& aStack, const int& aSlice) const;
    // 有効な垂直分割のインデックスを取得します。
    int getValidSlice(const int& aStack, const int& aSlice) const;
public:
    /**
     * \~english
     * @param aPosition  inital coordinate
     * @param aRadius    inital radius
     * @param aStack     Number of horizonal divisions
     * @param aSlice     Number of vertical divisions
     * @param aMass      mass
     * @param aVisvosity viscosity
     * \~japanese
     * @param aPosition  初期座標
     * @param aRadius    初期半径
     * @param aStack     水平分割数
     * @param aSlice     垂直分割数
     * @param aMass      質量
     * @param aVisvosity 粘度
     */
    PhysicalSlime(
            const KVector& aPosition,
            const float& aRadius,
            const int& aStack,
            const int& aSlice,
            const float& aMass,
            const float& aVisvosity
            );
    virtual ~PhysicalSlime();

    void draw() const override;
    void update() override;

    /**
     * \~english
     * @brief translate center coordinate into argument.
     * @param aPosition new center coordinate
     * \~japanese
     * @brief 中心座標を指定座標に遷移させます。
     * @param aPosition 設定する中心座標
     */
    virtual void translate(const KVector& aPosition);

    /// @brief \~english  move each vertex randomly.
    /// @brief \~japanese 各頂点をランダムに移動します。
    void spike();
    /**
     * \~english
     * @brief apply force.
     * @param aPos   point of application of force
     * @param aForce force vector
     * \~japanese
     * @brief 力を加えます。
     * @param aPos   力の作用点
     * @param aForce 力ベクトル
     */
    void applyForce(
            const KVector& aPos,
            const KVector& aForce
            );
};

#endif /* PHYSICALSLIME_H */

PhysicalSlime.cpp

/**
 * @file   PhysicalSlime.cpp
 * @brief  PhysicalSlime
 * @author Maeda Takumi
 */
#include "PhysicalSlime.h"

#include "PhysicalPoint.h"

PhysicalSlime::PhysicalSlime(
        const KVector& aPosition,
        const float& aRadius,
        const int& aStack,
        const int& aSlice,
        const float& aMass,
        const float& aVisvosity
        ) :
mRadius(aRadius),
mStack(aStack),
mSlice(aSlice),
mVertexSize((mStack - 1) * mSlice + 2),
mMass(aMass),
mPointMass(mMass / mVertexSize),
mVertex(mStack + 1, Vector<KVector>(mSlice)),
mNormal(mStack + 1, Vector<KVector>(mSlice)),
mPoints(mStack + 1, Vector<PhysicalPoint*>(mSlice, NULL)),
mIdealVolume(4.0f / 3.0f * Math::PI * powf(mRadius, 3.0f)),
mViscosity(aVisvosity) {
    const float thetaStack(Math::PI / mStack);
    const float thetaSlice((Math::PI * 2) / mStack);
    // 頂点と法線と質点の割り当て
    for (int i = 0; i <= mStack; ++i) {
        float theta(thetaStack * i);
        float h(mRadius * cos(theta)); // 輪切り高さ
        float r(mRadius * sin(theta)); // 輪切り半径
        for (int j = 0; j < mSlice; ++j) {
            if (!(i == 0 || i == mStack) || j == 0) { // 重なる頂点は作らない
                float phi(thetaSlice * j);
                KVector vertex(r * sin(phi), h, r * cos(phi));
                mVertex[i][j] = vertex;
                mNormal[i][j] = vertex.normalization();
                mPoints[i][j] = new PhysicalPoint(mPointMass, vertex);
            }
        }
    }
    translate(aPosition);

    spike();
}

PhysicalSlime::~PhysicalSlime() {
    for (Vector<PhysicalPoint*>& i : mPoints) {
        for (PhysicalPoint*& j : i) {
            if (j) delete j;
        }
    }
}

void PhysicalSlime::draw() const {
    glColor4ub(0xfc, 0xe2, 0xc4, 0x11);

    glBegin(GL_TRIANGLE_STRIP);
    for (int i = 0; i < mStack; ++i) {
        for (int j = 0; j <= mSlice; ++j) {
            int slice(getValidSlice(i, j));
            KVector n(mNormal[i][slice]);
            KVector v(mPoints[i][slice]->position());
            glNormal3f(DEPLOY_VEC(n));
            glVertex3f(DEPLOY_VEC(v));

            slice = getValidSlice(i + 1, j);
            n = mNormal[i + 1][slice];
            v = mPoints[i + 1][slice]->position();
            glNormal3f(DEPLOY_VEC(n));
            glVertex3f(DEPLOY_VEC(v));
        }
    }
    glEnd();
}

void PhysicalSlime::update() {
    const float cVol(calcVolume());

    float correction((mIdealVolume - cVol) / mIdealVolume);
    correction = Math::min(Math::max(-1.0f, correction), 1.0f) / mVertexSize;
    println(correction * mVertexSize);

    KVector centroid;
    for (int i = 0; i <= mStack; ++i) {
        for (int j = 0; j < mSlice; ++j) {
            if (isValidIndex(i, j)) {
                KVector position(mPoints[i][j]->position());

                mNormal[i][j] = (position - mPosition).normalization();
                mPoints[i][j]->applyForce(mNormal[i][j] * mRadius * correction);

                KVector resilience(mVertex[i][j] - position);
                mPoints[i][j]->applyForce(resilience / sqrtf(mVertexSize));

                // 速度を0にする力(柔らかさを変更できるようにする)
                mPoints[i][j]->applyForce(-mPoints[i][j]->velocity() * mPointMass * 1.0_s * mViscosity);

                centroid += position;
            }
        }
    }

    // 重心の更新
    translateVertex(centroid / mVertexSize);
}

float PhysicalSlime::calcVolume() const {
    const KVector y_axis(mVertex[mStack][0] - mVertex[0][0]);

    float volume(0.0f);

    // 上面・底面(輪切りにして計算(積分のイメージ))
    float tArea(0.0f), bArea(0.0f);
    KVector tCenter = mVertex[mStack][0], bCenter;
    for (int i = 1; i <= mStack; ++i) {
        bArea = 0.0f;
        bCenter = KVector();
        if (i != mStack) {
            for (int j = 0; j < mSlice; ++j) bCenter += mPoints[i][getValidSlice(i, j)]->position();
            bCenter /= mSlice;

            KVector point, next = mPoints[i][0]->position();
            for (int j = 0; j < mSlice; ++j) {
                point = next;
                next = mPoints[i][getValidSlice(i, j)]->position();
                KVector a(point - bCenter), b(next - bCenter); // 中心から各頂点までのベクトル
                bArea += a.length() * b.length() * a.angle(b);
            }
        }
        // (半径^2 * θ) / 2(扇形の面積)
        bArea /= 2; // まとめて割る

        KVector heightV((bCenter - tCenter).extractParallel(y_axis));
        float aveArea((tArea + bArea) / (tArea && bArea ? 2 : 3));
        volume += Math::sign(heightV.dot(y_axis)) * aveArea * heightV.length();

        tArea = bArea;
        tCenter = bCenter;
    }
    return volume;
}

void PhysicalSlime::translateVertex(const KVector& aPosition) {
    const KVector move(aPosition - mPosition);
    mPosition += move;
    for (Vector<KVector>& i : mVertex) {
        for (KVector& j : i) {
            j += move;
        }
    }
}

bool PhysicalSlime::isValidIndex(const int& aStack, const int& aSlice) const {
    return mPoints[aStack][aSlice];
}

int PhysicalSlime::getValidSlice(const int& aStack, const int& aSlice) const {
    if (aSlice == mSlice || !isValidIndex(aStack, aSlice)) return 0;
    return aSlice;
}

void PhysicalSlime::translate(const KVector & aPosition) {
    const KVector move(aPosition - mPosition);
    mPosition += move;
    for (int i = 0; i <= mStack; ++i) {
        for (int j = 0; j < mSlice; ++j) {
            if (isValidIndex(i, j)) {
                mVertex[i][j] += move;
                mPoints[i][j]->translate(mPoints[i][j]->position() + move);
            }
        }
    }
}

void PhysicalSlime::spike() {
    for (int i = 0; i <= mStack; ++i) {
        for (int j = 0; j < mSlice; ++j) {
            if (isValidIndex(i, j)) {
                KVector randVec(random(300) - 150, random(300) - 150, random(300) - 150);
                mPoints[i][j]->translate(mPoints[i][j]->position() + randVec.normalization() * mRadius * 0.5);
            }
        }
    }
}

void PhysicalSlime::applyForce(const KVector& aPos, const KVector& aForce) {
    // 力の作用店に最も近い頂点を探す。
    PhysicalPoint * nearPoint(NULL);
    float nearLen(0xffffff);
    for (int i = 0; i <= mStack; ++i) {
        for (int j = 0; j < mSlice; ++j) {
            if (isValidIndex(i, j)) {
                float len((mPoints[i][j]->position() - aPos).length());
                if (len < nearLen) {
                    nearLen = len;
                    nearPoint = mPoints[i][j];
                }
            }
        }
    }
    if (nearPoint) { // 力の届く範囲内
        // 最近頂点からの距離により伝わる力を調整する。
        // 粘度が強いほど、力を伝わりやすくする。
        KVector nearPos(nearPoint->position());
        float diameter(mRadius * 2.0f);
        for (int i = 0; i <= mStack; ++i) {
            for (int j = 0; j < mSlice; ++j) {
                if (isValidIndex(i, j)) {
                    float len((mPoints[i][j]->position() - nearPos).length());
                    float correction((Math::max((diameter - len) / diameter, 0.0f) + mViscosity) / mVertexSize);
                    mPoints[i][j]->applyForce(-(nearPos - aPos).normalization() * aForce.length() * correction);
                }
            }
        }
    }
}

解説

注意

このプログラムはしっかりとした物理法則にしたがうものではなく、学術的な理論とは程遠いによって、それなりにスライムらしい挙動を実装しただけです。

概要

このプログラムの要点は、毎フレーム各頂点に3つの力を加えていることです(実質それだけです)。

  1. 体積を一定にしようとする力
  2. 球体になろうとする力
  3. 慣性を打ち消す力

1つめの力は、理想体積と計算で求める体積との比率により大きさが変動する力です。
2つめの力は、現在の頂点位置から理想的な球形への頂点位置と戻ろうとする力です。
3つめの力は、球形を安定させる力です。

1と2の力は物理的に理解できるのですが、3の力が物理法則に従っているかどうかが分かりません。
かといって3の力がない場合、球形が安定せず体積が無限に大きく(小さく)なってしまいます。

終わりに

試行錯誤が多く完成したときはテンションが異常に上がったのですが、いざ解説を書こうとすると大変でテンション下がりました。
解説はおざなりでしたが、パラメータや加える力の式を変えることで大きく挙動が変わると思うので今後とも作りこんでいこうと思います。