diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 53aecc5..8150988 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -26,7 +26,7 @@ jobs: run: | sudo apt update sudo apt install --fix-missing -y doxygen graphviz clang-format clang-tidy cppcheck lcov - sudo apt install --fix-missing -y gcc g++ libspdlog-dev libcgal-dev freeglut3-dev libboost-all-dev libvtk9-dev qtbase5-dev xorg-dev libglu1-mesa-dev libglm-dev libglfw3-dev + sudo apt install --fix-missing -y gcc g++ libspdlog-dev libcgal-dev freeglut3-dev libboost-all-dev libvtk9-dev qtbase5-dev xorg-dev libglu1-mesa-dev libglm-dev libglfw3-dev libglew-dev - name: Build run: | @@ -48,7 +48,8 @@ jobs: run: | sudo apt update sudo apt install --fix-missing -y doxygen graphviz clang-format clang-tidy cppcheck lcov - sudo apt install --fix-missing -y gcc g++ libspdlog-dev libcgal-dev freeglut3-dev libboost-all-dev libvtk9-dev qtbase5-dev xorg-dev libglu1-mesa-dev libglm-dev libglfw3-dev + sudo apt install --fix-missing -y gcc g++ libspdlog-dev libcgal-dev freeglut3-dev libboost-all-dev libvtk9-dev qtbase5-dev xorg-dev libglu1-mesa-dev libglm-dev libglfw3-dev libglew-dev + sudo apt-get install --no-install-recommends -y xvfb - name: Build run: | @@ -56,6 +57,24 @@ jobs: cmake --build build --target all cmake --build build --target coverage + - name: Start Xvfb + run: | + nohup Xvfb :99 -screen 0 1024x768x24 & + echo "DISPLAY=:99" >> $GITHUB_ENV + + - name: Run demo2d + run: | + for i in {1..7}; do + ./build/bin/demo2d --test --instruction $i + done + env: + DISPLAY: :99 + + - name: Run MassSpring3D + run: ./build/bin/MassSpring3D --test --instruction 1 + env: + DISPLAY: :99 + - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v3 with: @@ -69,13 +88,41 @@ jobs: - name: Install dependencies run: | - brew install doxygen graphviz llvm cppcheck lcov spdlog boost freeglut glfw vtk cgal glm + brew install doxygen graphviz llvm cppcheck lcov spdlog boost freeglut glfw glew vtk cgal glm ln -s "$(brew --prefix llvm)/bin/clang-format" "/usr/local/bin/clang-format" ln -s "$(brew --prefix llvm)/bin/clang-tidy" "/usr/local/bin/clang-tidy" + + - name: Install XQuartz + run: | + brew install --cask xquartz + - name: Start Xvfb + run: | + sudo Xvfb :99 -screen 0 1024x768x24 & + echo "DISPLAY=:99" >> $GITHUB_ENV + - name: Build run: | cmake --preset=build cmake --build build --target all + + - name: Run demo2d + run: | + for i in {1..7}; do + ./build/bin/demo2d --test --instruction $i + done + env: + DISPLAY: :99 + + - name: Run MassSpring3D + run: ./build/bin/MassSpring3D --test --instruction 1 + env: + DISPLAY: :99 + + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v3 + with: + files: build/coverage/coverage.info + verbose: true # @bug 在 osx 下无法正常执行 lcov # cmake --build build --target coverage diff --git a/cmake/3rd.cmake b/cmake/3rd.cmake index 0b27f0d..19ca1c9 100644 --- a/cmake/3rd.cmake +++ b/cmake/3rd.cmake @@ -219,4 +219,4 @@ find_package(glfw3 REQUIRED) if (NOT glfw3_FOUND) message(FATAL_ERROR "glfw3 not found.\n" "Following https://www.glfw.org to install.") -endif () +endif () \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f97b571..bfbb4ca 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -81,7 +81,7 @@ packageProject( # 生成文件目录 BINARY_DIR ${PROJECT_BINARY_DIR} # 头文件路径 - INCLUDE_DIR ${PROJECT_SOURCE_DIR}/src/include ${PROJECT_SOURCE_DIR}/src/math/include ${PROJECT_SOURCE_DIR}/src/physics/include + INCLUDE_DIR ${PROJECT_SOURCE_DIR}/src/include ${PROJECT_SOURCE_DIR}/src/collision_detect/include ${PROJECT_SOURCE_DIR}/src/math/include ${PROJECT_SOURCE_DIR}/src/physics/include # 与 target 的 INSTALL_INTERFACE 一致 INCLUDE_DESTINATION include/${PROJECT_NAME}-${PROJECT_VERSION} # 头文件过滤 @@ -120,7 +120,7 @@ packageProject( # 生成文件目录 BINARY_DIR ${PROJECT_BINARY_DIR} # 头文件路径 - INCLUDE_DIR ${PROJECT_SOURCE_DIR}/src/include ${PROJECT_SOURCE_DIR}/src/math/include ${PROJECT_SOURCE_DIR}/src/physics/include + INCLUDE_DIR ${PROJECT_SOURCE_DIR}/src/include ${PROJECT_SOURCE_DIR}/src/collision_detect/include ${PROJECT_SOURCE_DIR}/src/math/include ${PROJECT_SOURCE_DIR}/src/physics/include # 与 target 的 INSTALL_INTERFACE 一致 INCLUDE_DESTINATION include/${PROJECT_NAME}-${PROJECT_VERSION} # 头文件过滤 diff --git a/src/dtkPhysMassSpring.cpp b/src/dtkPhysMassSpring.cpp index 487ecd9..1715492 100644 --- a/src/dtkPhysMassSpring.cpp +++ b/src/dtkPhysMassSpring.cpp @@ -32,883 +32,928 @@ using namespace std; namespace dtk { -dtkPhysMassSpring::dtkPhysMassSpring(double defaultMass, double defaultStiff, - double defaultDamp, - double defaultPointDamp, - double defaultPointResistence, - dtkDouble3 defaultGravityAccel) { - mDefaultMass = defaultMass; - mDefaultStiff = defaultStiff; - mDefaultDamp = defaultDamp; - mDefaultPointDamp = defaultPointDamp; - mDefaultPointResistence = defaultPointResistence; - mDefaultGravityAccel = defaultGravityAccel; + dtkPhysMassSpring::dtkPhysMassSpring(double defaultMass, double defaultStiff, + double defaultDamp, + double defaultPointDamp, + double defaultPointResistence, + dtkDouble3 defaultGravityAccel) { + mDefaultMass = defaultMass; + mDefaultStiff = defaultStiff; + mDefaultDamp = defaultDamp; + mDefaultPointDamp = defaultPointDamp; + mDefaultPointResistence = defaultPointResistence; + mDefaultGravityAccel = defaultGravityAccel; #ifdef DTK_CL - mUseMultiThread = false; - source_ = ""; - mMPPosData = NULL; - mMPVelData = NULL; - mMPAccelData = NULL; - mMPForceAccumData = NULL; - mMPGravityData = NULL; - mMPMassData = NULL; - mMPPosBufferData = NULL; - mMPVelBufferData = NULL; - mMPAccelBufferData = NULL; - mMPForceDecoratorData = NULL; - mMPActiveData = NULL; - - mSpringVertexIDData = NULL; - mSpringOriLengthData = NULL; - mSpringStiffData = NULL; - mSpringDampData = NULL; + mUseMultiThread = false; + source_ = ""; + mMPPosData = NULL; + mMPVelData = NULL; + mMPAccelData = NULL; + mMPForceAccumData = NULL; + mMPGravityData = NULL; + mMPMassData = NULL; + mMPPosBufferData = NULL; + mMPVelBufferData = NULL; + mMPAccelBufferData = NULL; + mMPForceDecoratorData = NULL; + mMPActiveData = NULL; + + mSpringVertexIDData = NULL; + mSpringOriLengthData = NULL; + mSpringStiffData = NULL; + mSpringDampData = NULL; #endif - mUnderControl = false; -} + mUnderControl = false; + } -dtkPhysMassSpring::~dtkPhysMassSpring() { - for (dtkID i = 0; i < mMassPoints.size(); i++) { - delete mMassPoints[i]; + dtkPhysMassSpring::dtkPhysMassSpring(double defaultMass, double defaultStiff, + double defaultDamp, + double defaultPointDamp, + double defaultPointResistence, + double defaultTimeStep, + dtkDouble3 defaultGravityAccel) { + mDefaultMass = defaultMass; + mDefaultStiff = defaultStiff; + mDefaultDamp = defaultDamp; + mDefaultPointDamp = defaultPointDamp; + mDefaultPointResistence = defaultPointResistence; + mDefaultTimeStep = defaultTimeStep; + mDefaultGravityAccel = defaultGravityAccel; +#ifdef DTK_CL + mUseMultiThread = false; + source_ = ""; + mMPPosData = NULL; + mMPVelData = NULL; + mMPAccelData = NULL; + mMPForceAccumData = NULL; + mMPGravityData = NULL; + mMPMassData = NULL; + mMPPosBufferData = NULL; + mMPVelBufferData = NULL; + mMPAccelBufferData = NULL; + mMPForceDecoratorData = NULL; + mMPActiveData = NULL; + + mSpringVertexIDData = NULL; + mSpringOriLengthData = NULL; + mSpringStiffData = NULL; + mSpringDampData = NULL; +#endif + mUnderControl = false; } - for (dtkID i = 0; i < mSprings.size(); i++) { - delete mSprings[i]; + dtkPhysMassSpring::~dtkPhysMassSpring() { + for (dtkID i = 0; i < mMassPoints.size(); i++) { + delete mMassPoints[i]; + } + + for (dtkID i = 0; i < mSprings.size(); i++) { + delete mSprings[i]; + } } -} -void dtkPhysMassSpring::SetPoints(dtkPoints::Ptr points) { - dtkAssert(points.get() != NULL, NULL_POINTER); - mPts = points; -} + void dtkPhysMassSpring::SetPoints(dtkPoints::Ptr points) { + dtkAssert(points.get() != NULL, NULL_POINTER); + mPts = points; + } -dtkPoints::Ptr dtkPhysMassSpring::GetPoints() { return mPts; } + dtkPoints::Ptr dtkPhysMassSpring::GetPoints() { return mPts; } -const GK::Point3 &dtkPhysMassSpring::GetPoint(dtkID id) const { - dtkAssert(mPts.get() != NULL, NULL_POINTER); - return mPts->GetPoint(id); -} + const GK::Point3& dtkPhysMassSpring::GetPoint(dtkID id) const { + dtkAssert(mPts.get() != NULL, NULL_POINTER); + return mPts->GetPoint(id); + } -void dtkPhysMassSpring::SetPoint(dtkID id, const GK::Point3 &coord) { - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - dtkAssert(mPts->SetPoint(id, coord), OUT_OF_RANGE); -} + void dtkPhysMassSpring::SetPoint(dtkID id, const GK::Point3& coord) { + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + dtkAssert(mPts->SetPoint(id, coord), OUT_OF_RANGE); + } -dtkID dtkPhysMassSpring::AddMassPoint(dtkID id, const double &mass, - const dtkT3 &vel, - double pointDamp, double pointResistence, - dtkDouble3 defaultGravityAccel) { - mMassPoints.push_back(new dtkPhysMassPoint( + dtkID dtkPhysMassSpring::AddMassPoint(dtkID id, const double& mass, + const dtkT3& vel, + double pointDamp, double pointResistence, + dtkDouble3 defaultGravityAccel) { + mMassPoints.push_back(new dtkPhysMassPoint( id, mPts, mass, vel, pointDamp, pointResistence, defaultGravityAccel)); #ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG - mAltitudeSpringForces.push_back(dtkT3(0, 0, 0)); + mAltitudeSpringForces.push_back(dtkT3(0, 0, 0)); #endif - return mMassPoints.size() - 1; -} + return mMassPoints.size() - 1; + } -dtkID dtkPhysMassSpring::AddSpring(dtkID p1, dtkID p2, const double &stiff, - const double &damp) { - if (p1 > p2) - swap(p1, p2); - if (mEdgeMap.find(dtkID2(p1, p2)) == mEdgeMap.end()) { - dtkPhysSpring *newSpring = + dtkID dtkPhysMassSpring::AddSpring(dtkID p1, dtkID p2, const double& stiff, + const double& damp, const double& rest_length_factor) { + if (p1 > p2) + swap(p1, p2); + if (mEdgeMap.find(dtkID2(p1, p2)) == mEdgeMap.end()) { + dtkPhysSpring* newSpring = new dtkPhysSpring(mMassPoints[p1], mMassPoints[p2], stiff, damp); - mSprings.push_back(newSpring); - mEdgeMap[dtkID2(p1, p2)] = newSpring; + newSpring->SetRestLength(newSpring->GetOriLength() * rest_length_factor); + mSprings.push_back(newSpring); + mEdgeMap[dtkID2(p1, p2)] = newSpring; + } + return mSprings.size() - 1; } - return mSprings.size() - 1; -} -bool dtkPhysMassSpring::Update_s(double timeslice, ItrMethod method, - bool limitDeformation) { - if (timeslice == 0) - return true; - switch (method) { - case Euler: - PreUpdate(timeslice, method); - UpdateStrings(timeslice, method, 0, limitDeformation); - UpdateMassPoints(timeslice, method, 0); - PostUpdate(method); - break; - case Mid: - for (dtkID i = 0; i < 2; i++) { - PreUpdate(timeslice, method, i); - UpdateStrings(timeslice, method, i, limitDeformation); - UpdateMassPoints(timeslice, method, i); - PostUpdate(method, i); - } - break; - case RK4: - for (dtkID i = 0; i < 4; i++) { - PreUpdate(timeslice, method, i); - UpdateStrings(timeslice, method, i, limitDeformation); - UpdateMassPoints(timeslice, method, i); - PostUpdate(method, i); - } - break; - case Heun: - for (dtkID i = 0; i < 2; i++) { - PreUpdate(timeslice, method, i); + bool dtkPhysMassSpring::Update_s(double timeslice, ItrMethod method, + bool limitDeformation) { + if (timeslice == 0) + return true; + switch (method) { + case Euler: + PreUpdate(timeslice, method); + UpdateStrings(timeslice, method, 0, limitDeformation); + UpdateMassPoints(timeslice, method, 0); + PostUpdate(method); + break; + case Mid: + for (dtkID i = 0; i < 2; i++) { + PreUpdate(timeslice, method, i); + UpdateStrings(timeslice, method, i, limitDeformation); + UpdateMassPoints(timeslice, method, i); + PostUpdate(method, i); + } + break; + case RK4: + for (dtkID i = 0; i < 4; i++) { + PreUpdate(timeslice, method, i); + UpdateStrings(timeslice, method, i, limitDeformation); + UpdateMassPoints(timeslice, method, i); + PostUpdate(method, i); + } + break; + case Heun: + for (dtkID i = 0; i < 2; i++) { + PreUpdate(timeslice, method, i); #ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG - if (i == 0) { - mAltitudeSpringForces.clear(); - for (dtkID j = 0; j < mMassPoints.size(); j++) - mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); + if (i == 0) { + mAltitudeSpringForces.clear(); + for (dtkID j = 0; j < mMassPoints.size(); j++) + mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); + } +#endif + UpdateStrings(timeslice, method, i, limitDeformation); + UpdateMassPoints(timeslice, method, i); + PostUpdate(method, i); } + case Collision: + for (dtkID i = 0; i < 2; i++) { + PreUpdate(timeslice, method, i); +#ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG + if (i == 0) { + mAltitudeSpringForces.clear(); + for (dtkID j = 0; j < mMassPoints.size(); j++) + mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); + } #endif - UpdateStrings(timeslice, method, i, limitDeformation); - UpdateMassPoints(timeslice, method, i); - PostUpdate(method, i); + UpdateStrings(timeslice, method, i, limitDeformation); + UpdateMassPoints(timeslice, method, i); + PostUpdate(method, i); + } + default: { + // not handle } - case Collision: - for (dtkID i = 0; i < 2; i++) { - PreUpdate(timeslice, method, i); + } + + return true; + } + + bool dtkPhysMassSpring::Update_iteration(double timeslice, ItrMethod method, + dtkID iteration, + bool limitDeformation) { + if (timeslice == 0) + return true; + PreUpdate(timeslice, method, iteration); #ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG - if (i == 0) { - mAltitudeSpringForces.clear(); - for (dtkID j = 0; j < mMassPoints.size(); j++) - mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); - } -#endif - UpdateStrings(timeslice, method, i, limitDeformation); - UpdateMassPoints(timeslice, method, i); - PostUpdate(method, i); + if (iteration == 0) { + mAltitudeSpringForces.clear(); + for (dtkID j = 0; j < mMassPoints.size(); j++) + mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); } - default: { - // not handle +#endif + UpdateStrings(timeslice, method, iteration, limitDeformation); + UpdateMassPoints(timeslice, method, iteration); + PostUpdate(method, iteration); + return true; } + + bool dtkPhysMassSpring::Update(double timeslice, ItrMethod method, + bool limitDeformation) { +#ifdef DTK_CL + if (mUseMultiThread) + return Update_mt(timeslice, method, limitDeformation); + else +#endif + return Update_s(timeslice, method, limitDeformation); } - return true; -} + bool dtkPhysMassSpring::UpdateStrings(double timeslice, ItrMethod method, + dtkID iteration, bool limitDeformation) { + for (dtkID i = 0; i < mSprings.size(); i++) + mSprings[i]->Update(timeslice, method, iteration, limitDeformation); -bool dtkPhysMassSpring::Update_iteration(double timeslice, ItrMethod method, - dtkID iteration, - bool limitDeformation) { - if (timeslice == 0) return true; - PreUpdate(timeslice, method, iteration); -#ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG - if (iteration == 0) { - mAltitudeSpringForces.clear(); - for (dtkID j = 0; j < mMassPoints.size(); j++) - mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); } -#endif - UpdateStrings(timeslice, method, iteration, limitDeformation); - UpdateMassPoints(timeslice, method, iteration); - PostUpdate(method, iteration); - return true; -} - -bool dtkPhysMassSpring::Update(double timeslice, ItrMethod method, - bool limitDeformation) { -#ifdef DTK_CL - if (mUseMultiThread) - return Update_mt(timeslice, method, limitDeformation); - else -#endif - return Update_s(timeslice, method, limitDeformation); -} -bool dtkPhysMassSpring::UpdateStrings(double timeslice, ItrMethod method, - dtkID iteration, bool limitDeformation) { - for (dtkID i = 0; i < mSprings.size(); i++) - mSprings[i]->Update(timeslice, method, iteration, limitDeformation); + bool dtkPhysMassSpring::UpdateMassPoints(double timeslice, ItrMethod method, + dtkID iteration) { + for (dtkID i = 0; i < mMassPoints.size(); i++) + mMassPoints[i]->Update(timeslice, method, iteration); - return true; -} + return true; + } -bool dtkPhysMassSpring::UpdateMassPoints(double timeslice, ItrMethod method, - dtkID iteration) { - for (dtkID i = 0; i < mMassPoints.size(); i++) - mMassPoints[i]->Update(timeslice, method, iteration); + bool dtkPhysMassSpring::ApplyImpulse(double timeslice) { + for (dtkID i = 0; i < mMassPoints.size(); i++) + mMassPoints[i]->ApplyImpulse(); - return true; -} + return true; + } -bool dtkPhysMassSpring::ApplyImpulse(double timeslice) { - for (dtkID i = 0; i < mMassPoints.size(); i++) - mMassPoints[i]->ApplyImpulse(); + bool dtkPhysMassSpring::PreUpdate(double timeslice, ItrMethod method, + dtkID iteration) { + return false; + } - return true; -} + bool dtkPhysMassSpring::PostUpdate(ItrMethod method, dtkID iteration) { + return false; + } -bool dtkPhysMassSpring::PreUpdate(double timeslice, ItrMethod method, - dtkID iteration) { - return false; -} + void dtkPhysMassSpring::SetSpringStiffness(int id, double newStiff) { + if (id >= ((int)mSprings.size()) || id < 0) + return; + mSprings[id]->SetStiffness(newStiff); + } -bool dtkPhysMassSpring::PostUpdate(ItrMethod method, dtkID iteration) { - return false; -} + void dtkPhysMassSpring::SetSpringDamp(int id, double newDamp) { + if (id >= ((int)mSprings.size()) || id < 0) + return; + mSprings[id]->SetDamp(newDamp); + } -void dtkPhysMassSpring::SetSpringStiffness(int id, double newStiff) { - if (id >= ((int)mSprings.size()) || id < 0) - return; - mSprings[id]->SetStiffness(newStiff); -} + void dtkPhysMassSpring::SetPointMass(int id, double newMass) { + if (id >= mMassPoints.size() || id < 0) + return; + mMassPoints[id]->SetMass(newMass); + } -void dtkPhysMassSpring::SetSpringDamp(int id, double newDamp) { - if (id >= ((int)mSprings.size()) || id < 0) - return; - mSprings[id]->SetDamp(newDamp); -} + void dtkPhysMassSpring::SetTimeStep(double timeStep) { + mDefaultTimeStep = timeStep; + } -void dtkPhysMassSpring::SetPointMass(int id, double newMass) { - if (id >= mMassPoints.size() || id < 0) - return; - mMassPoints[id]->SetMass(newMass); -} + double dtkPhysMassSpring::GetTimeStep() { + return mDefaultTimeStep; + } #ifdef DTK_CL -bool dtkPhysMassSpring::Update_mt(double timeslice, ItrMethod method, - bool limitDeformation) { - if (timeslice == 0) - return true; - for (dtkID i = 0; i < 2; i++) { - PreUpdate(timeslice, method, i); + bool dtkPhysMassSpring::Update_mt(double timeslice, ItrMethod method, + bool limitDeformation) { + if (timeslice == 0) + return true; + for (dtkID i = 0; i < 2; i++) { + PreUpdate(timeslice, method, i); #ifdef DTK_PHYSMASSSPRINGIMPL_DEBUG - if (i == 0) { - mAltitudeSpringForces.clear(); - for (dtkID j = 0; j < mMassPoints.size(); j++) - mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); - } + if (i == 0) { + mAltitudeSpringForces.clear(); + for (dtkID j = 0; j < mMassPoints.size(); j++) + mAltitudeSpringForces.push_back(mMassPoints[j]->GetForceAccum()); + } #endif - UpdateCopyToCLArray(); - RunCLKernels(timeslice, i); - UpdateCopyBack(); + UpdateCopyToCLArray(); + RunCLKernels(timeslice, i); + UpdateCopyBack(); - PostUpdate(method, i); + PostUpdate(method, i); + } + return true; } - return true; -} -template int dtkPhysMassSpring::SetupCL() { - cl_int status; - clGetPlatformIDs(1, &platform, NULL); - clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL); - context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL); - queue = clCreateCommandQueue(context, device, 0, &status); - if (status != CL_SUCCESS) { - cout << "create command error!" << endl; - } - const char *source = source_.c_str(); - size_t sourceSize[] = {strlen(source_.c_str())}; - program = clCreateProgramWithSource(context, 1, &source, sourceSize, &status); - if (status != CL_SUCCESS) { - cout << "create program with source error!" << endl; - } - const char *option = "-g -O0"; - status = clBuildProgram(program, 1, &device, option, NULL, NULL); - if (status != CL_SUCCESS) { - cout << "build program error!" << endl; - if (status == CL_INVALID_PROGRAM) - cout << "invalid program" << endl; - if (status == CL_INVALID_VALUE) - cout << "invalid value" << endl; - if (status == CL_INVALID_DEVICE) - cout << "invalid device" << endl; - if (status == CL_INVALID_BUILD_OPTIONS) - cout << "invalid build options" << endl; - if (status == CL_INVALID_OPERATION) - cout << "invalid operation" << endl; - if (status == CL_COMPILER_NOT_AVAILABLE) - cout << "compiler not available" << endl; - if (status == CL_BUILD_PROGRAM_FAILURE) - cout << "build program failure" << endl; - } - kernelSpring = clCreateKernel(program, "dtkPhysSpringUpdate", NULL); - kernelMP = clCreateKernel(program, "dtkPhysPointUpdate", NULL); - - mSpringOriLengthBuffer = clCreateBuffer( + template int dtkPhysMassSpring::SetupCL() { + cl_int status; + clGetPlatformIDs(1, &platform, NULL); + clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL); + context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL); + queue = clCreateCommandQueue(context, device, 0, &status); + if (status != CL_SUCCESS) { + cout << "create command error!" << endl; + } + const char* source = source_.c_str(); + size_t sourceSize[] = { strlen(source_.c_str()) }; + program = clCreateProgramWithSource(context, 1, &source, sourceSize, &status); + if (status != CL_SUCCESS) { + cout << "create program with source error!" << endl; + } + const char* option = "-g -O0"; + status = clBuildProgram(program, 1, &device, option, NULL, NULL); + if (status != CL_SUCCESS) { + cout << "build program error!" << endl; + if (status == CL_INVALID_PROGRAM) + cout << "invalid program" << endl; + if (status == CL_INVALID_VALUE) + cout << "invalid value" << endl; + if (status == CL_INVALID_DEVICE) + cout << "invalid device" << endl; + if (status == CL_INVALID_BUILD_OPTIONS) + cout << "invalid build options" << endl; + if (status == CL_INVALID_OPERATION) + cout << "invalid operation" << endl; + if (status == CL_COMPILER_NOT_AVAILABLE) + cout << "compiler not available" << endl; + if (status == CL_BUILD_PROGRAM_FAILURE) + cout << "build program failure" << endl; + } + kernelSpring = clCreateKernel(program, "dtkPhysSpringUpdate", NULL); + kernelMP = clCreateKernel(program, "dtkPhysPointUpdate", NULL); + + mSpringOriLengthBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_float), mSpringOriLengthData, 0); - mSpringStiffBuffer = clCreateBuffer( + mSpringStiffBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_float), mSpringStiffData, 0); - mSpringDampBuffer = clCreateBuffer( + mSpringDampBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_float), mSpringDampData, 0); - mSpringVertexIDBuffer = + mSpringVertexIDBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * 2.0 * sizeof(cl_int), mSpringVertexIDData, 0); - mMPPosBuffer = clCreateBuffer( + mMPPosBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3.0 * sizeof(cl_float), mMPPosData, 0); - mMPVelBuffer = clCreateBuffer( + mMPVelBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), mMPVelData, 0); - mMPAccelBuffer = clCreateBuffer( + mMPAccelBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), mMPAccelData, 0); - mMPForceAccumBuffer = + mMPForceAccumBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), mMPForceAccumData, 0); - mMPGravityBuffer = clCreateBuffer( + mMPGravityBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), mMPGravityData, 0); - mMPMassBuffer = clCreateBuffer( + mMPMassBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * sizeof(cl_float), mMPMassData, &status); - this->CheckCLStatus(status, "massbuffer create buffer error!"); - mMPPosBufferBuffer = + this->CheckCLStatus(status, "massbuffer create buffer error!"); + mMPPosBufferBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), mMPPosBufferData, 0); - mMPVelBufferBuffer = + mMPVelBufferBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), mMPVelBufferData, 0); - mMPAccelBufferBuffer = + mMPAccelBufferBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), mMPAccelBufferData, 0); - mMPForceDecoratorBuffer = + mMPForceDecoratorBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), mMPForceDecoratorData, 0); - mMPActiveBuffer = clCreateBuffer( + mMPActiveBuffer = clCreateBuffer( context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * sizeof(cl_bool), mMPActiveData, &status); - mMPAdjacentVertexIDsBuffer = + mMPAdjacentVertexIDsBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_int) * 2, mMPAdjacentVertexIDsData, &status); - mMPNumberOfAdjacentVertexBuffer = + mMPNumberOfAdjacentVertexBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfMassPoints() * sizeof(cl_int), mMPNumberOfAdjacentVertexData, &status); - this->CheckCLStatus(status, "numofadjvertexbuffer create buffer error!"); - mMPStartAdjacentVertexIDBuffer = + this->CheckCLStatus(status, "numofadjvertexbuffer create buffer error!"); + mMPStartAdjacentVertexIDBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_int), mMPStartAdjacentVertexIDData, &status); - mMPAdjacentSpringIDsBuffer = + mMPAdjacentSpringIDsBuffer = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, this->GetNumberOfSprings() * sizeof(cl_int) * 2, mMPAdjacentSpringIDsData, &status); - global_work_size = NWITEMS; - mUseMultiThread = true; -} - -template -int dtkPhysMassSpring::RunCLKernels(double timeslice, dtkID iteration) { - mTimesliceData = static_cast(timeslice); - mIterationData = iteration; - cl_int status; + global_work_size = NWITEMS; + mUseMultiThread = true; + } - status = - clSetKernelArg(kernelMP, 0, sizeof(cl_int), (void *)&mNumberOfMassPoints); - // this->CheckCLStatus(status, "set kernel arg error! MP 0"); - status = clSetKernelArg(kernelMP, 1, sizeof(cl_mem), (void *)&mMPPosBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 1"); - status = clSetKernelArg(kernelMP, 2, sizeof(cl_mem), (void *)&mMPVelBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 2"); - status = clSetKernelArg(kernelMP, 3, sizeof(cl_mem), (void *)&mMPAccelBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 3"); - status = - clSetKernelArg(kernelMP, 4, sizeof(cl_mem), (void *)&mMPForceAccumBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 4"); - status = - clSetKernelArg(kernelMP, 5, sizeof(cl_mem), (void *)&mMPGravityBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 5"); - status = clSetKernelArg(kernelMP, 6, sizeof(cl_mem), (void *)&mMPMassBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 6"); - status = - clSetKernelArg(kernelMP, 7, sizeof(cl_mem), (void *)&mMPPosBufferBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 7"); - status = - clSetKernelArg(kernelMP, 8, sizeof(cl_mem), (void *)&mMPVelBufferBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 8"); - status = clSetKernelArg(kernelMP, 9, sizeof(cl_mem), - (void *)&mMPAccelBufferBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 9"); - status = clSetKernelArg(kernelMP, 10, sizeof(cl_mem), - (void *)&mMPForceDecoratorBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 10"); - status = - clSetKernelArg(kernelMP, 11, sizeof(cl_mem), (void *)&mMPActiveBuffer); - // this->CheckCLStatus(status, "set kernel arg error! MP 11"); - status = clSetKernelArg(kernelMP, 12, sizeof(cl_int), &mIterationData); - // this->CheckCLStatus(status, "set kernel arg error! MP 12"); - status = clSetKernelArg(kernelMP, 13, sizeof(cl_float), &mTimesliceData); - // this->CheckCLStatus(status, "set kernel arg error! MP 13"); - status = clSetKernelArg(kernelMP, 14, sizeof(cl_mem), - (void *)&mMPAdjacentVertexIDsBuffer); - status = clSetKernelArg(kernelMP, 15, sizeof(cl_mem), - (void *)&mMPNumberOfAdjacentVertexBuffer); - status = clSetKernelArg(kernelMP, 16, sizeof(cl_mem), - (void *)&mMPStartAdjacentVertexIDBuffer); - status = clSetKernelArg(kernelMP, 17, sizeof(cl_mem), - (void *)&mMPAdjacentSpringIDsBuffer); - status = clSetKernelArg(kernelMP, 18, sizeof(cl_mem), - (void *)&mSpringOriLengthBuffer); - status = - clSetKernelArg(kernelMP, 19, sizeof(cl_mem), (void *)&mSpringStiffBuffer); - status = - clSetKernelArg(kernelMP, 20, sizeof(cl_mem), (void *)&mSpringDampBuffer); - - work_size = this->GetNumberOfMassPoints(); - size_t local_work_size = 256; - work_size = (work_size / local_work_size + 1) * local_work_size; - clEnqueueNDRangeKernel(queue, kernelMP, 1, NULL, &work_size, &local_work_size, - 0, NULL, NULL); - // clFinish( queue ); - - clEnqueueReadBuffer(queue, mMPPosBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), - mMPPosData, 0, NULL, NULL); // - clEnqueueReadBuffer(queue, mMPVelBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), - mMPVelData, 0, NULL, NULL); // - clEnqueueReadBuffer(queue, mMPPosBufferBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), - mMPPosBufferData, 0, NULL, NULL); // - clEnqueueReadBuffer(queue, mMPVelBufferBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), - mMPVelBufferData, 0, NULL, NULL); // - clFinish(queue); - - /* - clEnqueueReadBuffer(queue, mMPAccelBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 3.0 * sizeof(cl_float), mMPAccelData, - 0, NULL, NULL ); clEnqueueReadBuffer(queue, mMPForceAccumBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 3.0 * sizeof(cl_float), - mMPForceAccumData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, - mMPGravityBuffer, CL_TRUE, 0, this->GetNumberOfMassPoints() * 3.0 * - sizeof(cl_float), mMPGravityData, 0, NULL, NULL ); // - clEnqueueReadBuffer(queue, mMPMassBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * sizeof(cl_float), mMPMassData, 0, - NULL, NULL ); // clEnqueueReadBuffer(queue, mMPAccelBufferBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * 9.0 * sizeof(cl_float), - mMPAccelBufferData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, - mMPForceDecoratorBuffer, CL_TRUE, 0, this->GetNumberOfMassPoints() * 3.0 * - sizeof(cl_float), mMPForceDecoratorData, 0, NULL, NULL ); // - clEnqueueReadBuffer(queue, mMPActiveBuffer, CL_TRUE, 0, - this->GetNumberOfMassPoints() * sizeof(cl_bool), mMPActiveData, 0, - NULL, NULL ); // clEnqueueReadBuffer(queue, mSpringVertexIDBuffer, CL_TRUE, 0, - this->GetNumberOfSprings() * 2.0 * sizeof(cl_int), - mSpringVertexIDData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, - mSpringOriLengthBuffer, CL_TRUE, 0, this->GetNumberOfSprings() * - sizeof(cl_float), mSpringOriLengthData, 0, NULL, NULL ); // - clEnqueueReadBuffer(queue, mSpringStiffBuffer, CL_TRUE, 0, - this->GetNumberOfSprings() * sizeof(cl_float), mSpringStiffData, 0, - NULL, NULL ); // clEnqueueReadBuffer(queue, mSpringDampBuffer, CL_TRUE, 0, - this->GetNumberOfSprings() * sizeof(cl_float), mSpringDampData, 0, - NULL, NULL ); // - */ - /* - cout << "RunKernel Ouptput:" << endl; - cout << mSpringOriLengthData[0] << endl; - cout << mSpringStiffData[0] << endl; - cout << mSpringDampData[0] << endl; - cout << mSpringVertexIDData[0] << endl; - cout << mMPPosData[0] << endl; - cout << mMPVelData[0] << endl; - cout << mMPForceAccumData[0] << endl; - cout << mMPGravityData[0] << endl; - cout << mMPMassData[0] << endl; - cout << mMPPosBufferData[0] << endl; - cout << mMPVelBufferData[0] << endl; - cout << mMPAccelBufferData[0] << endl; - cout << mMPForceDecoratorData[0] << endl; - cout << mMPActiveData[0] << endl; - cout << mMPAccelData[0] << endl; - cout << "/RunKernel Ouptput" << endl; - */ + template + int dtkPhysMassSpring::RunCLKernels(double timeslice, dtkID iteration) { + mTimesliceData = static_cast(timeslice); + mIterationData = iteration; + cl_int status; + + status = + clSetKernelArg(kernelMP, 0, sizeof(cl_int), (void*)&mNumberOfMassPoints); + // this->CheckCLStatus(status, "set kernel arg error! MP 0"); + status = clSetKernelArg(kernelMP, 1, sizeof(cl_mem), (void*)&mMPPosBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 1"); + status = clSetKernelArg(kernelMP, 2, sizeof(cl_mem), (void*)&mMPVelBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 2"); + status = clSetKernelArg(kernelMP, 3, sizeof(cl_mem), (void*)&mMPAccelBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 3"); + status = + clSetKernelArg(kernelMP, 4, sizeof(cl_mem), (void*)&mMPForceAccumBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 4"); + status = + clSetKernelArg(kernelMP, 5, sizeof(cl_mem), (void*)&mMPGravityBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 5"); + status = clSetKernelArg(kernelMP, 6, sizeof(cl_mem), (void*)&mMPMassBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 6"); + status = + clSetKernelArg(kernelMP, 7, sizeof(cl_mem), (void*)&mMPPosBufferBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 7"); + status = + clSetKernelArg(kernelMP, 8, sizeof(cl_mem), (void*)&mMPVelBufferBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 8"); + status = clSetKernelArg(kernelMP, 9, sizeof(cl_mem), + (void*)&mMPAccelBufferBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 9"); + status = clSetKernelArg(kernelMP, 10, sizeof(cl_mem), + (void*)&mMPForceDecoratorBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 10"); + status = + clSetKernelArg(kernelMP, 11, sizeof(cl_mem), (void*)&mMPActiveBuffer); + // this->CheckCLStatus(status, "set kernel arg error! MP 11"); + status = clSetKernelArg(kernelMP, 12, sizeof(cl_int), &mIterationData); + // this->CheckCLStatus(status, "set kernel arg error! MP 12"); + status = clSetKernelArg(kernelMP, 13, sizeof(cl_float), &mTimesliceData); + // this->CheckCLStatus(status, "set kernel arg error! MP 13"); + status = clSetKernelArg(kernelMP, 14, sizeof(cl_mem), + (void*)&mMPAdjacentVertexIDsBuffer); + status = clSetKernelArg(kernelMP, 15, sizeof(cl_mem), + (void*)&mMPNumberOfAdjacentVertexBuffer); + status = clSetKernelArg(kernelMP, 16, sizeof(cl_mem), + (void*)&mMPStartAdjacentVertexIDBuffer); + status = clSetKernelArg(kernelMP, 17, sizeof(cl_mem), + (void*)&mMPAdjacentSpringIDsBuffer); + status = clSetKernelArg(kernelMP, 18, sizeof(cl_mem), + (void*)&mSpringOriLengthBuffer); + status = + clSetKernelArg(kernelMP, 19, sizeof(cl_mem), (void*)&mSpringStiffBuffer); + status = + clSetKernelArg(kernelMP, 20, sizeof(cl_mem), (void*)&mSpringDampBuffer); + + work_size = this->GetNumberOfMassPoints(); + size_t local_work_size = 256; + work_size = (work_size / local_work_size + 1) * local_work_size; + clEnqueueNDRangeKernel(queue, kernelMP, 1, NULL, &work_size, &local_work_size, + 0, NULL, NULL); + // clFinish( queue ); + + clEnqueueReadBuffer(queue, mMPPosBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), + mMPPosData, 0, NULL, NULL); // + clEnqueueReadBuffer(queue, mMPVelBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 3 * sizeof(cl_float), + mMPVelData, 0, NULL, NULL); // + clEnqueueReadBuffer(queue, mMPPosBufferBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), + mMPPosBufferData, 0, NULL, NULL); // + clEnqueueReadBuffer(queue, mMPVelBufferBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 9 * sizeof(cl_float), + mMPVelBufferData, 0, NULL, NULL); // + clFinish(queue); + + /* + clEnqueueReadBuffer(queue, mMPAccelBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 3.0 * sizeof(cl_float), mMPAccelData, + 0, NULL, NULL ); clEnqueueReadBuffer(queue, mMPForceAccumBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 3.0 * sizeof(cl_float), + mMPForceAccumData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, + mMPGravityBuffer, CL_TRUE, 0, this->GetNumberOfMassPoints() * 3.0 * + sizeof(cl_float), mMPGravityData, 0, NULL, NULL ); // + clEnqueueReadBuffer(queue, mMPMassBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * sizeof(cl_float), mMPMassData, 0, + NULL, NULL ); // clEnqueueReadBuffer(queue, mMPAccelBufferBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * 9.0 * sizeof(cl_float), + mMPAccelBufferData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, + mMPForceDecoratorBuffer, CL_TRUE, 0, this->GetNumberOfMassPoints() * 3.0 * + sizeof(cl_float), mMPForceDecoratorData, 0, NULL, NULL ); // + clEnqueueReadBuffer(queue, mMPActiveBuffer, CL_TRUE, 0, + this->GetNumberOfMassPoints() * sizeof(cl_bool), mMPActiveData, 0, + NULL, NULL ); // clEnqueueReadBuffer(queue, mSpringVertexIDBuffer, CL_TRUE, 0, + this->GetNumberOfSprings() * 2.0 * sizeof(cl_int), + mSpringVertexIDData, 0, NULL, NULL ); // clEnqueueReadBuffer(queue, + mSpringOriLengthBuffer, CL_TRUE, 0, this->GetNumberOfSprings() * + sizeof(cl_float), mSpringOriLengthData, 0, NULL, NULL ); // + clEnqueueReadBuffer(queue, mSpringStiffBuffer, CL_TRUE, 0, + this->GetNumberOfSprings() * sizeof(cl_float), mSpringStiffData, 0, + NULL, NULL ); // clEnqueueReadBuffer(queue, mSpringDampBuffer, CL_TRUE, 0, + this->GetNumberOfSprings() * sizeof(cl_float), mSpringDampData, 0, + NULL, NULL ); // + */ + /* + cout << "RunKernel Ouptput:" << endl; + cout << mSpringOriLengthData[0] << endl; + cout << mSpringStiffData[0] << endl; + cout << mSpringDampData[0] << endl; + cout << mSpringVertexIDData[0] << endl; + cout << mMPPosData[0] << endl; + cout << mMPVelData[0] << endl; + cout << mMPForceAccumData[0] << endl; + cout << mMPGravityData[0] << endl; + cout << mMPMassData[0] << endl; + cout << mMPPosBufferData[0] << endl; + cout << mMPVelBufferData[0] << endl; + cout << mMPAccelBufferData[0] << endl; + cout << mMPForceDecoratorData[0] << endl; + cout << mMPActiveData[0] << endl; + cout << mMPAccelData[0] << endl; + cout << "/RunKernel Ouptput" << endl; + */ - return 0; -} + return 0; + } -bool dtkPhysMassSpring::Open(const char *fileName) { - size_t size; - char *str; + bool dtkPhysMassSpring::Open(const char* fileName) { + size_t size; + char* str; - fstream f(fileName); //, (fstream::in | fstream::binary)); - if (!f) - return false; - else // if(f.is_open()) - { - size_t sizeFile; - f.seekg(0, fstream::end); - size = sizeFile = (size_t)f.tellg(); - f.seekg(0, fstream::beg); - - str = new char[size + 1]; - if (!str) { - f.close(); + fstream f(fileName); //, (fstream::in | fstream::binary)); + if (!f) return false; - } + else // if(f.is_open()) + { + size_t sizeFile; + f.seekg(0, fstream::end); + size = sizeFile = (size_t)f.tellg(); + f.seekg(0, fstream::beg); + + str = new char[size + 1]; + if (!str) { + f.close(); + return false; + } - f.read(str, sizeFile); - f.close(); - str[size] = '\0'; + f.read(str, sizeFile); + f.close(); + str[size] = '\0'; - source_ = str; - delete[] str; - return true; + source_ = str; + delete[] str; + return true; + } } -} - -void dtkPhysMassSpring::InitialCLArray() { - size_t sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_float); - mMPPosData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPVelData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPAccelData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPForceAccumData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPForceDecoratorData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPGravityData = (cl_float *)malloc(sizeInBytes * 3.0); - mMPMassData = (cl_float *)malloc(sizeInBytes); - mMPPosBufferData = (cl_float *)malloc(sizeInBytes * 9.0); - mMPVelBufferData = (cl_float *)malloc(sizeInBytes * 9.0); - mMPAccelBufferData = (cl_float *)malloc(sizeInBytes * 9.0); - - sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_bool); - mMPActiveData = (cl_bool *)malloc(sizeInBytes); - - sizeInBytes = this->GetNumberOfSprings() * sizeof(cl_int); - mSpringVertexIDData = (cl_int *)malloc(sizeInBytes * 2.0); - - sizeInBytes = this->GetNumberOfSprings() * sizeof(cl_float); - mSpringOriLengthData = (cl_float *)malloc(sizeInBytes); - mSpringStiffData = (cl_float *)malloc(sizeInBytes); - mSpringDampData = (cl_float *)malloc(sizeInBytes); - - sizeInBytes = this->GetNumberOfSprings() * 2.0 * sizeof(cl_int); - mMPAdjacentVertexIDsData = (cl_int *)malloc(sizeInBytes); - mMPAdjacentSpringIDsData = (cl_int *)malloc(sizeInBytes); - sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_int); - mMPNumberOfAdjacentVertexData = (cl_int *)malloc(sizeInBytes); - mMPStartAdjacentVertexIDData = (cl_int *)malloc(sizeInBytes); - - vector> tmpAdjacentVertexIDs; - vector> tmpAdjacentSpringIDs; - tmpAdjacentVertexIDs.resize(this->GetNumberOfMassPoints()); - tmpAdjacentSpringIDs.resize(this->GetNumberOfMassPoints()); - - for (dtkID i = 0; i < mSprings.size(); i++) { - int id1, id2; - id1 = mSpringVertexIDData[i * 2] = + + void dtkPhysMassSpring::InitialCLArray() { + size_t sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_float); + mMPPosData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPVelData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPAccelData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPForceAccumData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPForceDecoratorData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPGravityData = (cl_float*)malloc(sizeInBytes * 3.0); + mMPMassData = (cl_float*)malloc(sizeInBytes); + mMPPosBufferData = (cl_float*)malloc(sizeInBytes * 9.0); + mMPVelBufferData = (cl_float*)malloc(sizeInBytes * 9.0); + mMPAccelBufferData = (cl_float*)malloc(sizeInBytes * 9.0); + + sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_bool); + mMPActiveData = (cl_bool*)malloc(sizeInBytes); + + sizeInBytes = this->GetNumberOfSprings() * sizeof(cl_int); + mSpringVertexIDData = (cl_int*)malloc(sizeInBytes * 2.0); + + sizeInBytes = this->GetNumberOfSprings() * sizeof(cl_float); + mSpringOriLengthData = (cl_float*)malloc(sizeInBytes); + mSpringStiffData = (cl_float*)malloc(sizeInBytes); + mSpringDampData = (cl_float*)malloc(sizeInBytes); + + sizeInBytes = this->GetNumberOfSprings() * 2.0 * sizeof(cl_int); + mMPAdjacentVertexIDsData = (cl_int*)malloc(sizeInBytes); + mMPAdjacentSpringIDsData = (cl_int*)malloc(sizeInBytes); + sizeInBytes = this->GetNumberOfMassPoints() * sizeof(cl_int); + mMPNumberOfAdjacentVertexData = (cl_int*)malloc(sizeInBytes); + mMPStartAdjacentVertexIDData = (cl_int*)malloc(sizeInBytes); + + vector> tmpAdjacentVertexIDs; + vector> tmpAdjacentSpringIDs; + tmpAdjacentVertexIDs.resize(this->GetNumberOfMassPoints()); + tmpAdjacentSpringIDs.resize(this->GetNumberOfMassPoints()); + + for (dtkID i = 0; i < mSprings.size(); i++) { + int id1, id2; + id1 = mSpringVertexIDData[i * 2] = mSprings[i]->GetFirstVertex()->GetPointID(); - id2 = mSpringVertexIDData[i * 2 + 1] = + id2 = mSpringVertexIDData[i * 2 + 1] = mSprings[i]->GetSecondVertex()->GetPointID(); - tmpAdjacentVertexIDs[id1].push_back( + tmpAdjacentVertexIDs[id1].push_back( mSprings[i]->GetSecondVertex()->GetPointID()); - tmpAdjacentVertexIDs[id2].push_back( + tmpAdjacentVertexIDs[id2].push_back( mSprings[i]->GetFirstVertex()->GetPointID()); - tmpAdjacentSpringIDs[id1].push_back(i); - tmpAdjacentSpringIDs[id2].push_back(i); - - mSpringOriLengthData[i] = mSprings[i]->GetOriLength(); - mSpringStiffData[i] = mSprings[i]->GetStiffness(); - mSpringDampData[i] = mSprings[i]->GetDamp(); - } - - int startADJ = 0; - int startSADJ = 0; - for (dtkID i = 0; i < mMassPoints.size(); i++) { - dtkT3 pos = mMassPoints[i]->GetPosition(); - dtkT3 vel = mMassPoints[i]->GetVel(); - dtkT3 accel = mMassPoints[i]->GetAccel(); - dtkT3 forceAccum = mMassPoints[i]->GetForceAccum(); - dtkT3 gravity = mMassPoints[i]->GetGravity(); - dtkT3 forceDecorator = mMassPoints[i]->GetForceDecorator(); - std::vector> posBuffer; - posBuffer.resize(3); - posBuffer = mMassPoints[i]->GetPosBuffer(); - std::vector> velBuffer; - velBuffer.resize(3); - velBuffer = mMassPoints[i]->GetVelBuffer(); - std::vector> accelBuffer; - accelBuffer.resize(3); - accelBuffer = mMassPoints[i]->GetAccelBuffer(); - bool activity = mMassPoints[i]->GetActive(); - mMPMassData[i] = mMassPoints[i]->GetMass(); - mMPActiveData[i] = activity; - for (int j = 0; j < 3; j++) { - mMPPosData[i * 3 + j] = pos[j]; - mMPVelData[i * 3 + j] = vel[j]; - mMPAccelData[i * 3 + j] = accel[j]; - mMPForceAccumData[i * 3 + j] = forceAccum[j]; - mMPGravityData[i * 3 + j] = gravity[j]; - mMPForceDecoratorData[i * 3 + j] = forceDecorator[j]; + tmpAdjacentSpringIDs[id1].push_back(i); + tmpAdjacentSpringIDs[id2].push_back(i); + + mSpringOriLengthData[i] = mSprings[i]->GetOriLength(); + mSpringStiffData[i] = mSprings[i]->GetStiffness(); + mSpringDampData[i] = mSprings[i]->GetDamp(); } - for (int j = 0; j < 3; j++) { - for (int m = 0; m < 3; m++) { - mMPPosBufferData[i * 9 + j * 3 + m] = posBuffer[j][m]; - mMPVelBufferData[i * 9 + j * 3 + m] = velBuffer[j][m]; - mMPAccelBufferData[i * 9 + j * 3 + m] = accelBuffer[j][m]; + + int startADJ = 0; + int startSADJ = 0; + for (dtkID i = 0; i < mMassPoints.size(); i++) { + dtkT3 pos = mMassPoints[i]->GetPosition(); + dtkT3 vel = mMassPoints[i]->GetVel(); + dtkT3 accel = mMassPoints[i]->GetAccel(); + dtkT3 forceAccum = mMassPoints[i]->GetForceAccum(); + dtkT3 gravity = mMassPoints[i]->GetGravity(); + dtkT3 forceDecorator = mMassPoints[i]->GetForceDecorator(); + std::vector> posBuffer; + posBuffer.resize(3); + posBuffer = mMassPoints[i]->GetPosBuffer(); + std::vector> velBuffer; + velBuffer.resize(3); + velBuffer = mMassPoints[i]->GetVelBuffer(); + std::vector> accelBuffer; + accelBuffer.resize(3); + accelBuffer = mMassPoints[i]->GetAccelBuffer(); + bool activity = mMassPoints[i]->GetActive(); + mMPMassData[i] = mMassPoints[i]->GetMass(); + mMPActiveData[i] = activity; + for (int j = 0; j < 3; j++) { + mMPPosData[i * 3 + j] = pos[j]; + mMPVelData[i * 3 + j] = vel[j]; + mMPAccelData[i * 3 + j] = accel[j]; + mMPForceAccumData[i * 3 + j] = forceAccum[j]; + mMPGravityData[i * 3 + j] = gravity[j]; + mMPForceDecoratorData[i * 3 + j] = forceDecorator[j]; + } + for (int j = 0; j < 3; j++) { + for (int m = 0; m < 3; m++) { + mMPPosBufferData[i * 9 + j * 3 + m] = posBuffer[j][m]; + mMPVelBufferData[i * 9 + j * 3 + m] = velBuffer[j][m]; + mMPAccelBufferData[i * 9 + j * 3 + m] = accelBuffer[j][m]; + } } + + mMPStartAdjacentVertexIDData[i] = startADJ; + mMPNumberOfAdjacentVertexData[i] = tmpAdjacentVertexIDs[i].size(); + for (int j = 0; j < tmpAdjacentVertexIDs[i].size(); j++) + mMPAdjacentVertexIDsData[startADJ++] = tmpAdjacentVertexIDs[i][j]; + for (int j = 0; j < tmpAdjacentSpringIDs[i].size(); j++) + mMPAdjacentSpringIDsData[startSADJ++] = tmpAdjacentSpringIDs[i][j]; } - mMPStartAdjacentVertexIDData[i] = startADJ; - mMPNumberOfAdjacentVertexData[i] = tmpAdjacentVertexIDs[i].size(); - for (int j = 0; j < tmpAdjacentVertexIDs[i].size(); j++) - mMPAdjacentVertexIDsData[startADJ++] = tmpAdjacentVertexIDs[i][j]; - for (int j = 0; j < tmpAdjacentSpringIDs[i].size(); j++) - mMPAdjacentSpringIDsData[startSADJ++] = tmpAdjacentSpringIDs[i][j]; - } - - assert(startADJ == this->GetNumberOfSprings() * 2); - assert(startSADJ == this->GetNumberOfSprings() * 2); - - mNumberOfMassPoints = this->GetNumberOfMassPoints(); - /* - mSpringOriLengthData[0] = 1; - mSpringStiffData[0] = 2; - mSpringDampData[0] = 3; - mSpringVertexIDData[0] = 4; - mMPPosData[0] = 5; - mMPVelData[0] = 6; - mMPForceAccumData[0] = 7; - mMPGravityData[0] = 8; - mMPMassData[0] = 9; - mMPPosBufferData[0] = 10; - mMPVelBufferData[0] = 11; - mMPAccelBufferData[0] = 12; - mMPForceDecoratorData[0] = 13; - mMPActiveData[0] = true; - mMPAccelData[0] = 15; - */ -} - -void dtkPhysMassSpring::UpdateCopyToCLArray() { - for (dtkID i = 0; i < mMassPoints.size(); i++) { - // dtkT3 vel = mMassPoints[i]->GetVel(); - // dtkT3 accel = mMassPoints[i]->GetAccel(); - // dtkT3 pos = mMassPoints[i]->GetPosition(); - dtkT3 forceAccum = mMassPoints[i]->GetForceAccum(); + assert(startADJ == this->GetNumberOfSprings() * 2); + assert(startSADJ == this->GetNumberOfSprings() * 2); + + mNumberOfMassPoints = this->GetNumberOfMassPoints(); /* - dtkT3 forceDecorator = mMassPoints[i]->GetForceDecorator(); - vector< dtkT3 > posBuffer = mMassPoints[i]->GetPosBuffer(); - vector< dtkT3 > velBuffer = mMassPoints[i]->GetVelBuffer(); - vector< dtkT3 > accelBuffer = mMassPoints[i]->GetAccelBuffer(); + mSpringOriLengthData[0] = 1; + mSpringStiffData[0] = 2; + mSpringDampData[0] = 3; + mSpringVertexIDData[0] = 4; + mMPPosData[0] = 5; + mMPVelData[0] = 6; + mMPForceAccumData[0] = 7; + mMPGravityData[0] = 8; + mMPMassData[0] = 9; + mMPPosBufferData[0] = 10; + mMPVelBufferData[0] = 11; + mMPAccelBufferData[0] = 12; + mMPForceDecoratorData[0] = 13; + mMPActiveData[0] = true; + mMPAccelData[0] = 15; */ - bool activity = mMassPoints[i]->IsActive(); - mMPActiveData[i] = activity; - for (int j = 0; j < 3; j++) { - // mMPVelData[i * 3 + j] = vel[j]; - // mMPAccelData[i * 3 + j] = accel[j]; - // mMPPosData[i * 3 + j] = pos[j]; - mMPForceAccumData[i * 3 + j] = static_cast(forceAccum[j]); + } + + void dtkPhysMassSpring::UpdateCopyToCLArray() { + for (dtkID i = 0; i < mMassPoints.size(); i++) { + // dtkT3 vel = mMassPoints[i]->GetVel(); + // dtkT3 accel = mMassPoints[i]->GetAccel(); + // dtkT3 pos = mMassPoints[i]->GetPosition(); + dtkT3 forceAccum = mMassPoints[i]->GetForceAccum(); + /* + dtkT3 forceDecorator = mMassPoints[i]->GetForceDecorator(); + vector< dtkT3 > posBuffer = mMassPoints[i]->GetPosBuffer(); + vector< dtkT3 > velBuffer = mMassPoints[i]->GetVelBuffer(); + vector< dtkT3 > accelBuffer = mMassPoints[i]->GetAccelBuffer(); + */ + bool activity = mMassPoints[i]->IsActive(); + mMPActiveData[i] = activity; + for (int j = 0; j < 3; j++) { + // mMPVelData[i * 3 + j] = vel[j]; + // mMPAccelData[i * 3 + j] = accel[j]; + // mMPPosData[i * 3 + j] = pos[j]; + mMPForceAccumData[i * 3 + j] = static_cast(forceAccum[j]); + } + /* + for(int j = 0;j < 4;j++) + { + for(int m = 0;m < 3;m++) + { + mMPPosBufferData[i * 12 + j * 3 + m] = posBuffer[j][m]; + mMPVelBufferData[i * 12 + j * 3 + m] = velBuffer[j][m]; + mMPAccelBufferData[i * 12 + j * 3 + m] = accelBuffer[j][m]; + } + } + */ } /* - for(int j = 0;j < 4;j++) - { - for(int m = 0;m < 3;m++) - { - mMPPosBufferData[i * 12 + j * 3 + m] = posBuffer[j][m]; - mMPVelBufferData[i * 12 + j * 3 + m] = velBuffer[j][m]; - mMPAccelBufferData[i * 12 + j * 3 + m] = accelBuffer[j][m]; - } + cl_int status; + status = clEnqueueWriteBuffer(queue, + mMPPosBuffer, + CL_TRUE, + 0, + sizeof(cl_float) * this->GetNumberOfMassPoints() * 3.0, + mMPPosData, + 0, NULL, NULL); + + status = clEnqueueWriteBuffer(queue, + mMPActiveBuffer, + CL_TRUE, + 0, + sizeof(cl_bool) * this->GetNumberOfMassPoints(), + mMPActiveData, + 0, NULL, &writeEvt); + this->CheckCLStatus( status, "write to spring vertex buffer error!"); + status = WaitForEventAndRelease(&writeEvt); + this->CheckCLStatus( status, "wait for evt and release error!"); + */ + clEnqueueWriteBuffer(queue, mMPForceAccumBuffer, CL_TRUE, 0, + sizeof(cl_float) * this->GetNumberOfMassPoints() * 3, + mMPForceAccumData, 0, NULL, NULL); + // clFinish(queue); + } + + void dtkPhysMassSpring::UpdateCopyBack() { + for (dtkID i = 0; i < mMassPoints.size(); i++) { + dtkT3 vel(mMPVelData[i * 3], mMPVelData[i * 3 + 1], + mMPVelData[i * 3 + 2]); + dtkT3 pos(mMPPosData[i * 3], mMPPosData[i * 3 + 1], + mMPPosData[i * 3 + 2]); + mMassPoints[i]->SetPosition(pos); + mMassPoints[i]->SetVel(vel); + dtkT3 velBuffer, posBuffer; + for (int m = 0; m < 3; m++) { + velBuffer[m] = mMPVelBufferData[i * 9 + m]; + posBuffer[m] = mMPPosBufferData[i * 9 + m]; + } + mMassPoints[i]->SetVel(velBuffer, 0); + mMassPoints[i]->SetPosition(posBuffer, 0); + mMassPoints[i]->SetForceAccum(dtkT3(0, 0, 0)); } - */ } - /* -cl_int status; -status = clEnqueueWriteBuffer(queue, - mMPPosBuffer, - CL_TRUE, - 0, - sizeof(cl_float) * this->GetNumberOfMassPoints() * 3.0, - mMPPosData, - 0, NULL, NULL); - -status = clEnqueueWriteBuffer(queue, - mMPActiveBuffer, - CL_TRUE, - 0, - sizeof(cl_bool) * this->GetNumberOfMassPoints(), - mMPActiveData, - 0, NULL, &writeEvt); -this->CheckCLStatus( status, "write to spring vertex buffer error!"); -status = WaitForEventAndRelease(&writeEvt); -this->CheckCLStatus( status, "wait for evt and release error!"); -*/ - clEnqueueWriteBuffer(queue, mMPForceAccumBuffer, CL_TRUE, 0, - sizeof(cl_float) * this->GetNumberOfMassPoints() * 3, - mMPForceAccumData, 0, NULL, NULL); - // clFinish(queue); -} - -void dtkPhysMassSpring::UpdateCopyBack() { - for (dtkID i = 0; i < mMassPoints.size(); i++) { - dtkT3 vel(mMPVelData[i * 3], mMPVelData[i * 3 + 1], - mMPVelData[i * 3 + 2]); - dtkT3 pos(mMPPosData[i * 3], mMPPosData[i * 3 + 1], - mMPPosData[i * 3 + 2]); - mMassPoints[i]->SetPosition(pos); - mMassPoints[i]->SetVel(vel); - dtkT3 velBuffer, posBuffer; - for (int m = 0; m < 3; m++) { - velBuffer[m] = mMPVelBufferData[i * 9 + m]; - posBuffer[m] = mMPPosBufferData[i * 9 + m]; + + void dtkPhysMassSpring::CheckCLStatus(cl_int status, string message) { + if (status != CL_SUCCESS) + cout << message << endl; + return; + } + + int dtkPhysMassSpring::WaitForEventAndRelease(cl_event* event) { + cl_int status = CL_SUCCESS; + cl_int eventStatus = CL_QUEUED; + while (eventStatus != CL_COMPLETE) { + status = clGetEventInfo(*event, CL_EVENT_COMMAND_EXECUTION_STATUS, + sizeof(cl_int), &eventStatus, NULL); + CheckCLStatus(status, "clGetEventInfo error!"); } - mMassPoints[i]->SetVel(velBuffer, 0); - mMassPoints[i]->SetPosition(posBuffer, 0); - mMassPoints[i]->SetForceAccum(dtkT3(0, 0, 0)); - } -} - -void dtkPhysMassSpring::CheckCLStatus(cl_int status, string message) { - if (status != CL_SUCCESS) - cout << message << endl; - return; -} - -int dtkPhysMassSpring::WaitForEventAndRelease(cl_event *event) { - cl_int status = CL_SUCCESS; - cl_int eventStatus = CL_QUEUED; - while (eventStatus != CL_COMPLETE) { - status = clGetEventInfo(*event, CL_EVENT_COMMAND_EXECUTION_STATUS, - sizeof(cl_int), &eventStatus, NULL); - CheckCLStatus(status, "clGetEventInfo error!"); - } - status = clReleaseEvent(*event); - CheckCLStatus(status, "clReleaseEvent error!"); - return CL_SUCCESS; -} - -void dtkPhysMassSpring::UseMultiThread(const char *fileName) { - this->Open(fileName); - this->InitialCLArray(); - this->SetupCL(); -} + status = clReleaseEvent(*event); + CheckCLStatus(status, "clReleaseEvent error!"); + return CL_SUCCESS; + } + + void dtkPhysMassSpring::UseMultiThread(const char* fileName) { + this->Open(fileName); + this->InitialCLArray(); + this->SetupCL(); + } #endif -size_t dtkPhysMassSpring::FindTwins(Ptr ms, double distance) { - size_t count = 0; + size_t dtkPhysMassSpring::FindTwins(Ptr ms, double distance) { + size_t count = 0; - for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { - dtkPhysMassPoint *point1 = this->GetMassPoint(i); - if (point1->HasTwin()) - continue; + for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { + dtkPhysMassPoint* point1 = this->GetMassPoint(i); + if (point1->HasTwin()) + continue; - double minDistance = dtkDoubleMax; - dtkT3 pos1 = point1->GetPosition(); + double minDistance = dtkDoubleMax; + dtkT3 pos1 = point1->GetPosition(); - dtkPhysMassPoint *point2; - dtkT3 pos2; - dtkID minDistanceID = 0; + dtkPhysMassPoint* point2; + dtkT3 pos2; + dtkID minDistanceID = 0; - for (dtkID j = 0; j < ms->GetNumberOfMassPoints(); j++) { - point2 = ms->GetMassPoint(j); - if (point2->HasTwin()) - continue; - double tempDistance; - pos2 = point2->GetPosition(); - tempDistance = length(pos1 - pos2); + for (dtkID j = 0; j < ms->GetNumberOfMassPoints(); j++) { + point2 = ms->GetMassPoint(j); + if (point2->HasTwin()) + continue; + double tempDistance; + pos2 = point2->GetPosition(); + tempDistance = length(pos1 - pos2); - if (tempDistance < minDistance) { - minDistance = tempDistance; - minDistanceID = j; + if (tempDistance < minDistance) { + minDistance = tempDistance; + minDistanceID = j; + } + } + if (minDistance < distance) { + point1->AddTwin(ms->GetMassPoint(minDistanceID)); + ms->GetMassPoint(minDistanceID)->AddTwin(point1); + count++; } } - if (minDistance < distance) { - point1->AddTwin(ms->GetMassPoint(minDistanceID)); - ms->GetMassPoint(minDistanceID)->AddTwin(point1); - count++; - } - } - return count; -} - -void dtkPhysMassSpring::ConvertImpulseToForce(double timeslice) { - mImpulseForce = dtkT3(0, 0, 0); - for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { - dtkPhysMassPoint *point = this->GetMassPoint(i); - mImpulseForce = mImpulseForce + point->GetAndClearImpulse(); + return count; } - mImpulseForce = mImpulseForce * (1.0 / timeslice); -} -void dtkPhysMassSpring::RegisterLabel(dtkID label) { - mLabels.push_back(label); - mTransportForces[label] = dtkT3(0, 0, 0); -} + void dtkPhysMassSpring::ConvertImpulseToForce(double timeslice) { + mImpulseForce = dtkT3(0, 0, 0); + for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { + dtkPhysMassPoint* point = this->GetMassPoint(i); + mImpulseForce = mImpulseForce + point->GetAndClearImpulse(); + } + mImpulseForce = mImpulseForce * (1.0 / timeslice); + } -void dtkPhysMassSpring::TransportForce(double timeslice) { - for (std::map>::iterator itr = mTransportForces.begin(); - itr != mTransportForces.end(); itr++) { - mTransportForces[itr->first] = dtkT3(0, 0, 0); + void dtkPhysMassSpring::RegisterLabel(dtkID label) { + mLabels.push_back(label); + mTransportForces[label] = dtkT3(0, 0, 0); } - for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { - dtkPhysMassPoint *point = this->GetMassPoint(i); - if (!point->IsActive()) { - dtkID label = point->GetLabel(); - if (label != 0) { - mTransportForces[label] = + + void dtkPhysMassSpring::TransportForce(double timeslice) { + for (std::map>::iterator itr = mTransportForces.begin(); + itr != mTransportForces.end(); itr++) { + mTransportForces[itr->first] = dtkT3(0, 0, 0); + } + for (dtkID i = 0; i < this->GetNumberOfMassPoints(); i++) { + dtkPhysMassPoint* point = this->GetMassPoint(i); + if (!point->IsActive()) { + dtkID label = point->GetLabel(); + if (label != 0) { + mTransportForces[label] = mTransportForces[label] + point->GetForceAccum(); + } } } } -} -void dtkPhysMassSpring::SetTriangleMesh( + void dtkPhysMassSpring::SetTriangleMesh( dtkStaticTriangleMesh::Ptr newTriangleMesh) { - this->SetPoints(newTriangleMesh->GetPoints()); - for (dtkID i = 0; i < newTriangleMesh->GetNumberOfPoints(); i++) { - this->AddMassPoint(i, mDefaultMass, dtkT3(0, 0, 0), - mDefaultPointDamp, mDefaultPointResistence, - mDefaultGravityAccel); - } - const std::vector &ec = newTriangleMesh->GetECTable(); - set edges; - for (dtkID i = 0; i < ec.size(); i++) { - if (ec[i][0] > ec[i][1]) - edges.insert(dtkID2(ec[i][1], ec[i][0])); - else - edges.insert(dtkID2(ec[i][0], ec[i][1])); - - if (ec[i][0] > ec[i][2]) - edges.insert(dtkID2(ec[i][2], ec[i][0])); - else - edges.insert(dtkID2(ec[i][0], ec[i][2])); + this->SetPoints(newTriangleMesh->GetPoints()); + for (dtkID i = 0; i < newTriangleMesh->GetNumberOfPoints(); i++) { + this->AddMassPoint(i, mDefaultMass, dtkT3(0, 0, 0), + mDefaultPointDamp, mDefaultPointResistence, + mDefaultGravityAccel); + } + const std::vector& ec = newTriangleMesh->GetECTable(); + set edges; + for (dtkID i = 0; i < ec.size(); i++) { + if (ec[i][0] > ec[i][1]) + edges.insert(dtkID2(ec[i][1], ec[i][0])); + else + edges.insert(dtkID2(ec[i][0], ec[i][1])); + + if (ec[i][0] > ec[i][2]) + edges.insert(dtkID2(ec[i][2], ec[i][0])); + else + edges.insert(dtkID2(ec[i][0], ec[i][2])); + + if (ec[i][2] > ec[i][1]) + edges.insert(dtkID2(ec[i][1], ec[i][2])); + else + edges.insert(dtkID2(ec[i][2], ec[i][1])); + } - if (ec[i][2] > ec[i][1]) - edges.insert(dtkID2(ec[i][1], ec[i][2])); - else - edges.insert(dtkID2(ec[i][2], ec[i][1])); + for (set::iterator itr = edges.begin(); itr != edges.end(); itr++) { + this->AddSpring(itr->a, itr->b, mDefaultStiff, mDefaultDamp); + } } - for (set::iterator itr = edges.begin(); itr != edges.end(); itr++) { - this->AddSpring(itr->a, itr->b, mDefaultStiff, mDefaultDamp); + dtkPhysSpring* dtkPhysMassSpring::GetSpringByPoints(dtkID2 edge) { + map::iterator it; + it = mEdgeMap.find(edge); + if (it != mEdgeMap.end()) + return it->second; + else + return 0; } -} -dtkPhysSpring *dtkPhysMassSpring::GetSpringByPoints(dtkID2 edge) { - map::iterator it; - it = mEdgeMap.find(edge); - if (it != mEdgeMap.end()) - return it->second; - else - return 0; -} - -void dtkPhysMassSpring::DeleteSpring(dtkID id1, dtkID id2) { - if (id1 > id2) - swap(id1, id2); - dtkPhysSpring *spring = mEdgeMap[dtkID2(id1, id2)]; - for (dtkID i = 0; i < mSprings.size(); i++) { - if (mSprings[i] == spring) { - mSprings.erase(mSprings.begin() + i); - break; + void dtkPhysMassSpring::DeleteSpring(dtkID id1, dtkID id2) { + if (id1 > id2) + swap(id1, id2); + dtkPhysSpring* spring = mEdgeMap[dtkID2(id1, id2)]; + for (dtkID i = 0; i < mSprings.size(); i++) { + if (mSprings[i] == spring) { + mSprings.erase(mSprings.begin() + i); + break; + } } + delete spring; } - delete spring; -} -void dtkPhysMassSpring::AbandonTwins() { - for (dtkID i = 0; i < GetNumberOfMassPoints(); i++) { - mMassPoints[i]->AbandonTwins(); + void dtkPhysMassSpring::AbandonTwins() { + for (dtkID i = 0; i < GetNumberOfMassPoints(); i++) { + mMassPoints[i]->AbandonTwins(); + } } -} } // namespace dtk diff --git a/src/dtkStaticTriangleMesh.cpp b/src/dtkStaticTriangleMesh.cpp index 6a7c985..eedc431 100644 --- a/src/dtkStaticTriangleMesh.cpp +++ b/src/dtkStaticTriangleMesh.cpp @@ -18,368 +18,397 @@ #include #include -// Debug + // Debug #include using namespace std; #include "dtkStaticTriangleMesh.h" namespace dtk { -dtkStaticTriangleMesh::dtkStaticTriangleMesh() : mModified(false) { - // nothing. -} - -dtkStaticTriangleMesh::~dtkStaticTriangleMesh() { - // nothing. -} - -void dtkStaticTriangleMesh::SetPoints(dtkPoints::Ptr points) { - dtkAssert(points.get() != NULL, NULL_POINTER); - mPts = points; -} - -dtkPoints::Ptr dtkStaticTriangleMesh::GetPoints() { return mPts; } - -const GK::Point3 &dtkStaticTriangleMesh::GetPoint(dtkID id) const { - dtkAssert(mPts.get() != NULL, NULL_POINTER); - return mPts->GetPoint(id); -} - -void dtkStaticTriangleMesh::SetPoint(dtkID id, const GK::Point3 &coord) { - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - dtkAssert(mPts->SetPoint(id, coord), OUT_OF_RANGE); -} - -void dtkStaticTriangleMesh::Clear() { - mEC.clear(); - mV2E.clear(); - if (mPts.get() != NULL) - mV2E.resize(mPts->GetMaxID() + 1, dtkErrorID); - mE2E.clear(); - mB2E.clear(); -} - -dtkID dtkStaticTriangleMesh::InsertTriangle(dtkID v0, dtkID v1, dtkID v2) { - // You must call SetPoints before excuting this function. - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - - const GK::Point3 &coord0 = mPts->GetPoint(v0); - const GK::Point3 &coord1 = mPts->GetPoint(v1); - const GK::Point3 &coord2 = mPts->GetPoint(v2); - - dtkSign orient = GK::Orient3D(coord0, coord1, coord2); - // dtkAssert(orient != dtkSign::ZERO, DEGENERACY_TRIANGLE); - if (orient == dtkSign::ZERO) - return dtkErrorID; - // Disable this line to ensure clockwise, but harm other functionality - // if (orient == dtkSign::NEGATIVE) std::swap(v1, v2); - - // Insert into EC Table. - dtkID3 entry = dtkID3(v0, v1, v2); - mEC.push_back(entry); - - mModified = true; - - return (dtkID)(mEC.size() - 1); -} - -dtkID dtkStaticTriangleMesh::InsertTriangle(dtkID v[3]) { - return InsertTriangle(v[0], v[1], v[2]); -} - -dtkID dtkStaticTriangleMesh::InsertTriangle(const dtkID3 &tri) { - return InsertTriangle(tri.a, tri.b, tri.c); -} - -dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(dtkID v0, dtkID v1, - dtkID v2) { - // You must call SetPoints before excuting this function. - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - - if (IsRepeat(v0, v1, v2) == true) { + dtkStaticTriangleMesh::dtkStaticTriangleMesh() : mModified(false) { + // nothing. + } + + dtkStaticTriangleMesh::~dtkStaticTriangleMesh() { + // nothing. + } + + void dtkStaticTriangleMesh::SetPoints(dtkPoints::Ptr points) { + dtkAssert(points.get() != NULL, NULL_POINTER); + mPts = points; + } + + dtkPoints::Ptr dtkStaticTriangleMesh::GetPoints() { return mPts; } + + const GK::Point3& dtkStaticTriangleMesh::GetPoint(dtkID id) const { + dtkAssert(mPts.get() != NULL, NULL_POINTER); + return mPts->GetPoint(id); + } + + void dtkStaticTriangleMesh::SetPoint(dtkID id, const GK::Point3& coord) { + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + dtkAssert(mPts->SetPoint(id, coord), OUT_OF_RANGE); + } + + void dtkStaticTriangleMesh::Clear() { + mEC.clear(); + mV2E.clear(); + if (mPts.get() != NULL) + mV2E.resize(mPts->GetMaxID() + 1, dtkErrorID); + mE2E.clear(); + mB2E.clear(); + } + + dtkID dtkStaticTriangleMesh::InsertTriangle(dtkID v0, dtkID v1, dtkID v2) { + // You must call SetPoints before excuting this function. + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + + const GK::Point3& coord0 = mPts->GetPoint(v0); + const GK::Point3& coord1 = mPts->GetPoint(v1); + const GK::Point3& coord2 = mPts->GetPoint(v2); + + dtkSign orient = GK::Orient3D(coord0, coord1, coord2); + // dtkAssert(orient != dtkSign::ZERO, DEGENERACY_TRIANGLE); + if (orient == dtkSign::ZERO) + return dtkErrorID; + // Disable this line to ensure clockwise, but harm other functionality + // if (orient == dtkSign::NEGATIVE) std::swap(v1, v2); + + // Insert into EC Table. + dtkID3 entry = dtkID3(v0, v1, v2); + mEC.push_back(entry); + + mModified = true; + return (dtkID)(mEC.size() - 1); } - const GK::Point3 &coord0 = mPts->GetPoint(v0); - const GK::Point3 &coord1 = mPts->GetPoint(v1); - const GK::Point3 &coord2 = mPts->GetPoint(v2); - - dtkSign orient = GK::Orient3D(coord0, coord1, coord2); - // dtkAssert(orient != dtkSign::ZERO, DEGENERACY_TRIANGLE); - if (orient == dtkSign::ZERO) - return dtkErrorID; - // Disable this line to ensure clockwise, but harm other functionality - if (orient == dtkSign::NEGATIVE) - std::swap(v1, v2); - - // Insert into EC Table. - dtkID3 entry = dtkID3(v0, v1, v2); - mEC.push_back(entry); - - mModified = true; - - return (dtkID)(mEC.size() - 1); -} - -dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(dtkID v[3]) { - return InsertTriangleNotRepeat(v[0], v[1], v[2]); -} - -dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(const dtkID3 &tri) { - return InsertTriangleNotRepeat(tri.a, tri.b, tri.c); -} - -bool dtkStaticTriangleMesh::IsRepeat(dtkID v0, dtkID v1, dtkID v2) { - for (unsigned int i = 0; i < mEC.size(); i++) { - if ((mEC[i].a == v0 && mEC[i].b == v1 && mEC[i].c == v2) || - (mEC[i].a == v0 && mEC[i].b == v2 && mEC[i].c == v1) || - (mEC[i].a == v1 && mEC[i].b == v0 && mEC[i].c == v2) || - (mEC[i].a == v1 && mEC[i].b == v2 && mEC[i].c == v0) || - (mEC[i].a == v2 && mEC[i].b == v0 && mEC[i].c == v1) || - (mEC[i].a == v2 && mEC[i].b == v1 && mEC[i].c == v0)) { - return true; - } + dtkID dtkStaticTriangleMesh::InsertTriangle(dtkID v[3]) { + return InsertTriangle(v[0], v[1], v[2]); } - return false; -} - -dtkID dtkStaticTriangleMesh::RemoveTriangle(dtkID v0, dtkID v1, dtkID v2) { - // You must call SetPoints before excuting this function. - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - - dtkID3 remove = dtkID3(v0, v1, v2); - std::vector::iterator iter; - iter = find(mEC.begin(), mEC.end(), remove); - - mEC.erase(iter); - mModified = true; - - return ((dtkID)mEC.size() - 1); -} - -dtkID dtkStaticTriangleMesh::RemoveTriangle(dtkID v[3]) { - return RemoveTriangle(v[0], v[1], v[2]); -} - -dtkID dtkStaticTriangleMesh::RemoveTriangle(const dtkID3 &tri) { - return RemoveTriangle(tri.a, tri.b, tri.c); -} - -dtkID dtkStaticTriangleMesh::RemoveTriangleDisordered(dtkID v0, dtkID v1, - dtkID v2) { - for (unsigned int i = 0; i < mEC.size(); i++) { - if (mEC[i].a == v0 && mEC[i].b == v1 && mEC[i].c == v2) { - return RemoveTriangle(v0, v1, v2); - } else if (mEC[i].a == v0 && mEC[i].b == v2 && mEC[i].c == v1) { - return RemoveTriangle(v0, v2, v1); - } else if (mEC[i].a == v1 && mEC[i].b == v0 && mEC[i].c == v2) { - return RemoveTriangle(v1, v0, v2); - } else if (mEC[i].a == v1 && mEC[i].b == v2 && mEC[i].c == v0) { - return RemoveTriangle(v1, v2, v0); - } else if (mEC[i].a == v2 && mEC[i].b == v0 && mEC[i].c == v1) { - return RemoveTriangle(v2, v0, v1); - } else if (mEC[i].a == v2 && mEC[i].b == v1 && mEC[i].c == v0) { - return RemoveTriangle(v2, v1, v0); + + dtkID dtkStaticTriangleMesh::InsertTriangle(const dtkID3& tri) { + return InsertTriangle(tri.a, tri.b, tri.c); + } + + dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(dtkID v0, dtkID v1, + dtkID v2) { + // You must call SetPoints before excuting this function. + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + + if (IsRepeat(v0, v1, v2) == true) { + return (dtkID)(mEC.size() - 1); } + + const GK::Point3& coord0 = mPts->GetPoint(v0); + const GK::Point3& coord1 = mPts->GetPoint(v1); + const GK::Point3& coord2 = mPts->GetPoint(v2); + + dtkSign orient = GK::Orient3D(coord0, coord1, coord2); + // dtkAssert(orient != dtkSign::ZERO, DEGENERACY_TRIANGLE); + if (orient == dtkSign::ZERO) + return dtkErrorID; + // Disable this line to ensure clockwise, but harm other functionality + if (orient == dtkSign::NEGATIVE) + std::swap(v1, v2); + + // Insert into EC Table. + dtkID3 entry = dtkID3(v0, v1, v2); + mEC.push_back(entry); + + mModified = true; + + return (dtkID)(mEC.size() - 1); + } + + dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(dtkID v[3]) { + return InsertTriangleNotRepeat(v[0], v[1], v[2]); + } + + dtkID dtkStaticTriangleMesh::InsertTriangleNotRepeat(const dtkID3& tri) { + return InsertTriangleNotRepeat(tri.a, tri.b, tri.c); } - return 0; -} - -dtkID dtkStaticTriangleMesh::modifyElement(dtkID srcIndex, int srcOrder, - dtkID dest) { - assert(srcIndex < mEC.size()); - assert(srcOrder <= 2 && srcOrder >= 0); - assert(dest < mPts->GetNumberOfPoints()); - - mEC[srcIndex][srcOrder] = dest; - - return srcIndex; -} - -dtkID dtkStaticTriangleMesh::modifyElement(dtkID srcIndex, dtkID3 dest) { - dtkAssert(srcIndex < mEC.size()); - dtkAssert(dest.a < mPts->GetNumberOfPoints() && - dest.b < mPts->GetNumberOfPoints() && - dest.c < mPts->GetNumberOfPoints()); - - mEC[srcIndex] = dest; - - return srcIndex; -} - -void dtkStaticTriangleMesh::Modified() { mModified = true; } - -bool dtkStaticTriangleMesh::IsModified() const { return mModified; } - -void dtkStaticTriangleMesh::Rebuild() { - const static dtkID2 edges[3] = {dtkID2(1, 2), dtkID2(2, 3), dtkID2(3, 1)}; - const static dtkFloat3 upz = dtkFloat3(0.0f, 0.0f, 1.0f); - - if (!mModified) - return; - dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); - - std::vector(mEC).swap(mEC); - - // clear mV2E, mE2E, mB2E, rebulid mV2E, mE2E, mB2E according the new mEC. - dtkID maxID = mPts->GetMaxID(); - mV2E.clear(); - mV2E.resize(maxID + 1, dtkErrorID); - mE2E.clear(); - mE2E.resize(mEC.size(), dtkID3(dtkErrorID, dtkErrorID, dtkErrorID)); - mB2E.clear(); - - // Rebuild. dtkDWORD is 64-bits data type, hes represent half-edges. - std::map hes; - // once a half-edge is finished, it will add into finishedEdge. - std::set finishedEdge; - - // Find the adjacent triangle of each triangle. - // And find the top-most vertex and its incident triangle. - dtkID firstTri; - GK::Float maxZ = -dtkFloatMax; - for (size_t i = 0; i < mEC.size(); ++i) { - const dtkID3 &tri = mEC[i]; - - for (int j = 0; j < 3; ++j) { - // get three half-edges from a triangle with three vertexs. - dtkID2 edge = sort(mEC[i][edges[j].a - 1], mEC[i][edges[j].b - 1]); - // combine the from point and to point into key. - dtkDWORD key = Merge(edge.a, edge.b); - - // One edge must incident to no more than 2 triangles. - dtkAssert(finishedEdge.find(key) == finishedEdge.end(), - MUST_BE_MANIFOLD_MESH); - - std::map::iterator pos; - pos = hes.find(key); - if (pos != - hes.end()) // a pair edge found, because edge with a, b are sorted. - { - // return the mapped value. - dtkID he = pos->second; - dtkID d, s; - // decode the code of half-edge - // d present triangle index, s present half-edge index in local - // triangle. - dtkStaticTriangleMesh::DecodeHE(he, d, s); - - dtkAssert(mE2E[d][s - 1] == dtkErrorID, MUST_BE_MANIFOLD_MESH); - dtkAssert(mE2E[i][j] == dtkErrorID, MUST_BE_MANIFOLD_MESH); - - //------------------------------------------------------------------------- - // there should be: - // dthID he_ij = dtkStaticTriangleMesh::EncodeHE((dtkID)i, (dtkID)(j + - // 1)); mE2E[d][s-1] = he_ij; mE2E[i][j] = he; - //------------------------------------------------------------------------- - mE2E[d][s - 1] = (dtkID)i; - mE2E[i][j] = d; - - hes.erase(pos); - finishedEdge.insert(key); - } else { - // add the new half-edge into hes. - dtkID he = dtkStaticTriangleMesh::EncodeHE((dtkID)i, (dtkID)(j + 1)); - hes.insert(std::make_pair(key, he)); - } - //------------------------------------------------ - // for what ??? - const GK::Point3 &coord = mPts->GetPoint(tri[j]); - if (coord.z() > maxZ) { - maxZ = coord.z(); - firstTri = (dtkID)i; + bool dtkStaticTriangleMesh::IsRepeat(dtkID v0, dtkID v1, dtkID v2) { + for (unsigned int i = 0; i < mEC.size(); i++) { + if ((mEC[i].a == v0 && mEC[i].b == v1 && mEC[i].c == v2) || + (mEC[i].a == v0 && mEC[i].b == v2 && mEC[i].c == v1) || + (mEC[i].a == v1 && mEC[i].b == v0 && mEC[i].c == v2) || + (mEC[i].a == v1 && mEC[i].b == v2 && mEC[i].c == v0) || + (mEC[i].a == v2 && mEC[i].b == v0 && mEC[i].c == v1) || + (mEC[i].a == v2 && mEC[i].b == v1 && mEC[i].c == v0)) { + return true; } - //------------------------------------------------- } + return false; + } + + dtkID dtkStaticTriangleMesh::RemoveTriangle(dtkID v0, dtkID v1, dtkID v2) { + // You must call SetPoints before excuting this function. + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + + dtkID3 remove = dtkID3(v0, v1, v2); + std::vector::iterator iter; + iter = find(mEC.begin(), mEC.end(), remove); + + mEC.erase(iter); + mModified = true; + + return ((dtkID)mEC.size() - 1); } - // Set the border HEs - // the way to justify board edge: no pair half-edge, which stay in hes. - for (std::map::iterator iter = hes.begin(); - iter != hes.end(); ++iter) { - dtkDWORD e = iter->first; - dtkID he = iter->second; + dtkID dtkStaticTriangleMesh::RemoveTriangle(dtkID v[3]) { + return RemoveTriangle(v[0], v[1], v[2]); + } - dtkID c, i; - dtkStaticTriangleMesh::DecodeHE(he, c, i); + dtkID dtkStaticTriangleMesh::RemoveTriangle(const dtkID3& tri) { + return RemoveTriangle(tri.a, tri.b, tri.c); + } - mB2E.push_back(he); - mE2E[c][i - 1] = dtkStaticTriangleMesh::EncodeHE((dtkID)mB2E.size() - 1, 0); - // i % 3 ??? - mV2E[mEC[c][i % 3]] = he; + dtkID dtkStaticTriangleMesh::RemoveTriangleDisordered(dtkID v0, dtkID v1, + dtkID v2) { + for (unsigned int i = 0; i < mEC.size(); i++) { + if (mEC[i].a == v0 && mEC[i].b == v1 && mEC[i].c == v2) { + return RemoveTriangle(v0, v1, v2); + } + else if (mEC[i].a == v0 && mEC[i].b == v2 && mEC[i].c == v1) { + return RemoveTriangle(v0, v2, v1); + } + else if (mEC[i].a == v1 && mEC[i].b == v0 && mEC[i].c == v2) { + return RemoveTriangle(v1, v0, v2); + } + else if (mEC[i].a == v1 && mEC[i].b == v2 && mEC[i].c == v0) { + return RemoveTriangle(v1, v2, v0); + } + else if (mEC[i].a == v2 && mEC[i].b == v0 && mEC[i].c == v1) { + return RemoveTriangle(v2, v0, v1); + } + else if (mEC[i].a == v2 && mEC[i].b == v1 && mEC[i].c == v0) { + return RemoveTriangle(v2, v1, v0); + } + } + return 0; } - std::vector setted; - std::queue detected; - setted.resize(mEC.size(), false); + dtkID dtkStaticTriangleMesh::modifyElement(dtkID srcIndex, int srcOrder, + dtkID dest) { + assert(srcIndex < mEC.size()); + assert(srcOrder <= 2 && srcOrder >= 0); + assert(dest < mPts->GetNumberOfPoints()); - const GK::Point3 &coord0 = mPts->GetPoint(mEC[firstTri][0]); - const GK::Point3 &coord1 = mPts->GetPoint(mEC[firstTri][1]); - const GK::Point3 &coord2 = mPts->GetPoint(mEC[firstTri][2]); + mEC[srcIndex][srcOrder] = dest; - // if the triangle's direction is negetive, modify the mE2E and mEC.but why - // not modify the mB2E? - if (GK::Orient3D(coord0, coord1, coord2) == dtkSign::NEGATIVE) { - std::swap(mEC[firstTri][1], mEC[firstTri][2]); - std::swap(mE2E[firstTri][0], mE2E[firstTri][2]); + return srcIndex; } - setted[firstTri] = true; - detected.push(firstTri); - while (!detected.empty()) { - dtkID curTri = detected.front(); - detected.pop(); + dtkID dtkStaticTriangleMesh::modifyElement(dtkID srcIndex, dtkID3 dest) { + dtkAssert(srcIndex < mEC.size()); + dtkAssert(dest.a < mPts->GetNumberOfPoints() && + dest.b < mPts->GetNumberOfPoints() && + dest.c < mPts->GetNumberOfPoints()); - for (int i = 0; i < 3; ++i) { - dtkID neighbour = mE2E[curTri][i]; - dtkID s, t, u; + mEC[srcIndex] = dest; - s = mEC[curTri][(i + 0) % 3]; - t = mEC[curTri][(i + 1) % 3]; - u = mEC[curTri][(i + 2) % 3]; - dtkDWORD e = Merge(s, t); + return srcIndex; + } + + void dtkStaticTriangleMesh::update_normals() { + mVertexNormals.resize(mPts->GetNumberOfPoints(), dtkDouble3(0.0f, 0.0f, 0.0f)); + for (size_t i = 0; i < mEC.size(); ++i) { + const dtkID3& tri = mEC[i]; + const GK::Point3& p0 = GetPoint(tri[0]); + const GK::Point3& p1 = GetPoint(tri[1]); + const GK::Point3& p2 = GetPoint(tri[2]); - dtkID2 sorted = sort(s, t); + GK::Vector3 v0 = p1 - p0; + GK::Vector3 v1 = p2 - p0; + GK::Vector3 faceNormal = GK::Normalize(GK::CrossProduct(v0, v1)); - if (hes.find(Merge(sorted.a, sorted.b)) != hes.end()) - continue; // Is boundary edge, has been dealt before. + dtkDouble3 fn = dtkDouble3(faceNormal.x(), faceNormal.y(), faceNormal.z()); - if (!setted[neighbour]) { - dtkDWORD e0, e1, e2; - e0 = Merge(mEC[neighbour][1], mEC[neighbour][0]); - e1 = Merge(mEC[neighbour][2], mEC[neighbour][1]); - e2 = Merge(mEC[neighbour][0], mEC[neighbour][2]); + mVertexNormals[tri[0]] += fn; + mVertexNormals[tri[1]] += fn; + mVertexNormals[tri[2]] += fn; + } + for (size_t i = 0; i < mVertexNormals.size(); ++i) { + mVertexNormals[i] = normalize(mVertexNormals[i]); + } + } - // Need to change the orientation. - if (e != e0 && e != e1 && e != e2) { - std::swap(mEC[neighbour][1], mEC[neighbour][2]); - std::swap(mE2E[neighbour][0], mE2E[neighbour][2]); + void dtkStaticTriangleMesh::Modified() { mModified = true; } + + bool dtkStaticTriangleMesh::IsModified() const { return mModified; } + + void dtkStaticTriangleMesh::Rebuild() { + const static dtkID2 edges[3] = { dtkID2(1, 2), dtkID2(2, 3), dtkID2(3, 1) }; + const static dtkFloat3 upz = dtkFloat3(0.0f, 0.0f, 1.0f); + + if (!mModified) + return; + dtkAssert(mPts.get() != NULL, ILLEGAL_STATE); + + std::vector(mEC).swap(mEC); + + // clear mV2E, mE2E, mB2E, rebulid mV2E, mE2E, mB2E according the new mEC. + dtkID maxID = mPts->GetMaxID(); + mV2E.clear(); + mV2E.resize(maxID + 1, dtkErrorID); + mE2E.clear(); + mE2E.resize(mEC.size(), dtkID3(dtkErrorID, dtkErrorID, dtkErrorID)); + mB2E.clear(); + + // Rebuild. dtkDWORD is 64-bits data type, hes represent half-edges. + std::map hes; + // once a half-edge is finished, it will add into finishedEdge. + std::set finishedEdge; + + // Find the adjacent triangle of each triangle. + // And find the top-most vertex and its incident triangle. + dtkID firstTri; + GK::Float maxZ = -dtkFloatMax; + for (size_t i = 0; i < mEC.size(); ++i) { + const dtkID3& tri = mEC[i]; + + for (int j = 0; j < 3; ++j) { + // get three half-edges from a triangle with three vertexs. + dtkID2 edge = sort(mEC[i][edges[j].a - 1], mEC[i][edges[j].b - 1]); + // combine the from point and to point into key. + dtkDWORD key = Merge(edge.a, edge.b); + + // One edge must incident to no more than 2 triangles. + dtkAssert(finishedEdge.find(key) == finishedEdge.end(), + MUST_BE_MANIFOLD_MESH); + + std::map::iterator pos; + pos = hes.find(key); + if (pos != + hes.end()) // a pair edge found, because edge with a, b are sorted. + { + // return the mapped value. + dtkID he = pos->second; + dtkID d, s; + // decode the code of half-edge + // d present triangle index, s present half-edge index in local + // triangle. + dtkStaticTriangleMesh::DecodeHE(he, d, s); + + dtkAssert(mE2E[d][s - 1] == dtkErrorID, MUST_BE_MANIFOLD_MESH); + dtkAssert(mE2E[i][j] == dtkErrorID, MUST_BE_MANIFOLD_MESH); + + //------------------------------------------------------------------------- + // there should be: + // dthID he_ij = dtkStaticTriangleMesh::EncodeHE((dtkID)i, (dtkID)(j + + // 1)); mE2E[d][s-1] = he_ij; mE2E[i][j] = he; + //------------------------------------------------------------------------- + mE2E[d][s - 1] = (dtkID)i; + mE2E[i][j] = d; + + hes.erase(pos); + finishedEdge.insert(key); + } + else { + // add the new half-edge into hes. + dtkID he = dtkStaticTriangleMesh::EncodeHE((dtkID)i, (dtkID)(j + 1)); + hes.insert(std::make_pair(key, he)); } - setted[neighbour] = true; - detected.push(neighbour); + //------------------------------------------------ + // for what ??? + const GK::Point3& coord = mPts->GetPoint(tri[j]); + if (coord.z() > maxZ) { + maxZ = coord.z(); + firstTri = (dtkID)i; + } + //------------------------------------------------- } + } + + // Set the border HEs + // the way to justify board edge: no pair half-edge, which stay in hes. + for (std::map::iterator iter = hes.begin(); + iter != hes.end(); ++iter) { + dtkDWORD e = iter->first; + dtkID he = iter->second; - dtkID idx; - if (t == mEC[neighbour][0]) - idx = 1; - if (t == mEC[neighbour][1]) - idx = 2; - if (t == mEC[neighbour][2]) - idx = 3; + dtkID c, i; + dtkStaticTriangleMesh::DecodeHE(he, c, i); - mE2E[curTri][i] = dtkStaticTriangleMesh::EncodeHE(neighbour, idx); + mB2E.push_back(he); + mE2E[c][i - 1] = dtkStaticTriangleMesh::EncodeHE((dtkID)mB2E.size() - 1, 0); + // i % 3 ??? + mV2E[mEC[c][i % 3]] = he; } - for (int i = 0; i < 3; ++i) { - dtkID he = mV2E[mEC[curTri][i]]; + std::vector setted; + std::queue detected; + setted.resize(mEC.size(), false); - if (!IsBoundaryHE(he)) - mV2E[mEC[curTri][i]] = dtkStaticTriangleMesh::EncodeHE(curTri, i + 1); + const GK::Point3& coord0 = mPts->GetPoint(mEC[firstTri][0]); + const GK::Point3& coord1 = mPts->GetPoint(mEC[firstTri][1]); + const GK::Point3& coord2 = mPts->GetPoint(mEC[firstTri][2]); + + // if the triangle's direction is negetive, modify the mE2E and mEC.but why + // not modify the mB2E? + if (GK::Orient3D(coord0, coord1, coord2) == dtkSign::NEGATIVE) { + std::swap(mEC[firstTri][1], mEC[firstTri][2]); + std::swap(mE2E[firstTri][0], mE2E[firstTri][2]); } - } + setted[firstTri] = true; + + detected.push(firstTri); + while (!detected.empty()) { + dtkID curTri = detected.front(); + detected.pop(); + + for (int i = 0; i < 3; ++i) { + dtkID neighbour = mE2E[curTri][i]; + dtkID s, t, u; + + s = mEC[curTri][(i + 0) % 3]; + t = mEC[curTri][(i + 1) % 3]; + u = mEC[curTri][(i + 2) % 3]; + dtkDWORD e = Merge(s, t); - mModified = false; -} + dtkID2 sorted = sort(s, t); + + if (hes.find(Merge(sorted.a, sorted.b)) != hes.end()) + continue; // Is boundary edge, has been dealt before. + + if (!setted[neighbour]) { + dtkDWORD e0, e1, e2; + e0 = Merge(mEC[neighbour][1], mEC[neighbour][0]); + e1 = Merge(mEC[neighbour][2], mEC[neighbour][1]); + e2 = Merge(mEC[neighbour][0], mEC[neighbour][2]); + + // Need to change the orientation. + if (e != e0 && e != e1 && e != e2) { + std::swap(mEC[neighbour][1], mEC[neighbour][2]); + std::swap(mE2E[neighbour][0], mE2E[neighbour][2]); + } + setted[neighbour] = true; + + detected.push(neighbour); + } + + dtkID idx; + if (t == mEC[neighbour][0]) + idx = 1; + if (t == mEC[neighbour][1]) + idx = 2; + if (t == mEC[neighbour][2]) + idx = 3; + + mE2E[curTri][i] = dtkStaticTriangleMesh::EncodeHE(neighbour, idx); + } + + for (int i = 0; i < 3; ++i) { + dtkID he = mV2E[mEC[curTri][i]]; + + if (!IsBoundaryHE(he)) + mV2E[mEC[curTri][i]] = dtkStaticTriangleMesh::EncodeHE(curTri, i + 1); + } + } + + mModified = false; + } } // namespace dtk diff --git a/src/include/dtk.h b/src/include/dtk.h index 092f0c2..575c9a3 100644 --- a/src/include/dtk.h +++ b/src/include/dtk.h @@ -175,7 +175,7 @@ #include "dtkUtility.h" // #include "dtkVolume.h" -#include "dtkArray_will_del.h" +// #include "dtkArray_will_del.h" // Common #ifdef DTK_CUDA diff --git a/src/include/dtkPhysMassSpring.h b/src/include/dtkPhysMassSpring.h index 8ba47f5..8955a6c 100644 --- a/src/include/dtkPhysMassSpring.h +++ b/src/include/dtkPhysMassSpring.h @@ -35,256 +35,279 @@ using namespace std; #define NWITEMS 256 namespace dtk { -class dtkPhysMassSpring : public boost::noncopyable { -public: - typedef std::shared_ptr Ptr; - - static Ptr - New(double defaultMass = 2, double defaultK = 2000, double defaultDamp = 5880, - double defaultPointDamp = 1.0, double defaultPointResistence = 2.5, - dtkDouble3 defaultGravityAccel = dtkDouble3( - 0, 0, 0)) // K = k0 * cross section area, damp is additional damp - { - return Ptr(new dtkPhysMassSpring(defaultMass, defaultK, defaultDamp, - defaultPointDamp, defaultPointResistence, - defaultGravityAccel)); - } - -public: - virtual ~dtkPhysMassSpring(); - dtkPhysMassSpring(double defaultMass, double defaultK, double defaultDamp, - double defaultPointDamp, - double defaultPointResistence = 2.5, - dtkDouble3 defaultGravityAccel = dtkDouble3(0, 0, 0)); - void SetPoints(dtkPoints::Ptr points); - dtkPoints::Ptr GetPoints(); - - const GK::Point3 &GetPoint(dtkID id) const; - void SetPoint(dtkID id, const GK::Point3 &coord); - - dtkPhysMassPoint *GetMassPoint(dtkID id) { return mMassPoints[id]; }; - dtkID GetNumberOfMassPoints() { return mMassPoints.size(); } - dtkID AddMassPoint(dtkID id, const double &mass = 1.0, - const dtkT3 &vel = dtkT3(0, 0, 0), - double pointDamp = 1.0, double pointResistence = 2.5, - dtkDouble3 defaultGravityAccel = dtkDouble3(0, 0, 0)); - - dtkPhysSpring *GetSpring(dtkID id) { return mSprings[id]; }; - dtkID GetNumberOfSprings() { return (dtkID)mSprings.size(); } - dtkID AddSpring(dtkID p1, dtkID p2, const double &stiff = 0, - const double &damp = 0); - void DeleteSpring(dtkID p1, dtkID p2); - - void SetSpringStiffness(int id, double newK); - void SetSpringDamp(int id, double newDamp); - void SetPointMass(int id, double newMass); - - /** - * @brief 更新弹簧及质点状态 - * @param[in] timeslice : 更新时间间隔 - * @param[in] method : 迭代更新算法 - * @param[in] limitDeformation : 弹簧弹性限度 - * @note 更新弹簧及质点状态 - * @return - * true update successfully \n - * false update failure \n - */ - bool Update(double timeslice, ItrMethod method = Euler, - bool limitDeformation = false); - - /** - * @brief 单线程更新弹簧及质点状态 - * @param[in] timeslice : 更新时间间隔 - * @param[in] method : 迭代更新算法 - * @param[in] limitDeformation : 弹簧弹性限度 - * @note 单线程更新弹簧及质点状态 - * @return - * true update successfully \n - * false update failure \n - */ - bool Update_s(double timeslice, ItrMethod method = Euler, + class dtkPhysMassSpring : public boost::noncopyable { + public: + typedef std::shared_ptr Ptr; + + static Ptr + New(double defaultMass = 2, double defaultK = 2000, double defaultDamp = 5880, + double defaultPointDamp = 1.0, double defaultPointResistence = 2.5, + dtkDouble3 defaultGravityAccel = dtkDouble3( + 0, 0, 0)) // K = k0 * cross section area, damp is additional damp + { + return Ptr(new dtkPhysMassSpring(defaultMass, defaultK, defaultDamp, + defaultPointDamp, defaultPointResistence, + defaultGravityAccel)); + } + + static Ptr + New(double defaultMass = 2, double defaultK = 2000, double defaultDamp = 5880, + double defaultPointDamp = 1.0, double defaultPointResistence = 2.5, double defaultTimeStep = 0.01, + dtkDouble3 defaultGravityAccel = dtkDouble3( + 0, 0, 0)) // K = k0 * cross section area, damp is additional damp + { + return Ptr(new dtkPhysMassSpring(defaultMass, defaultK, defaultDamp, + defaultPointDamp, defaultPointResistence, defaultTimeStep, + defaultGravityAccel)); + } + + public: + virtual ~dtkPhysMassSpring(); + dtkPhysMassSpring(double defaultMass, double defaultK, double defaultDamp, + double defaultPointDamp, + double defaultPointResistence = 2.5, + dtkDouble3 defaultGravityAccel = dtkDouble3(0, 0, 0)); + dtkPhysMassSpring(double defaultMass, double defaultK, double defaultDamp, + double defaultPointDamp, + double defaultPointResistence = 2.5, + double defaultTimeStep = 0.01, + dtkDouble3 defaultGravityAccel = dtkDouble3(0, 0, 0)); + void SetPoints(dtkPoints::Ptr points); + dtkPoints::Ptr GetPoints(); + + const GK::Point3& GetPoint(dtkID id) const; + void SetPoint(dtkID id, const GK::Point3& coord); + + dtkPhysMassPoint* GetMassPoint(dtkID id) { return mMassPoints[id]; }; + dtkID GetNumberOfMassPoints() { return mMassPoints.size(); } + dtkID AddMassPoint(dtkID id, const double& mass = 1.0, + const dtkT3& vel = dtkT3(0, 0, 0), + double pointDamp = 1.0, double pointResistence = 2.5, + dtkDouble3 defaultGravityAccel = dtkDouble3(0, 0, 0)); + + dtkPhysSpring* GetSpring(dtkID id) { return mSprings[id]; }; + dtkID GetNumberOfSprings() { return (dtkID)mSprings.size(); } + dtkID AddSpring(dtkID p1, dtkID p2, const double& stiff = 0, + const double& damp = 0, const double& rest_length_factor = 1.0); + void DeleteSpring(dtkID p1, dtkID p2); + + void SetSpringStiffness(int id, double newK); + void SetSpringDamp(int id, double newDamp); + void SetPointMass(int id, double newMass); + void SetTimeStep(double timeStep); + double GetTimeStep(); + double GetDefaultStiffness() const { return mDefaultStiff; } + double GetDefaultDamp() const { return mDefaultDamp; } + double GetDefaultPointDamp() const { return mDefaultPointDamp; } + dtk::dtkDouble3 GetDefaultGravityAccel() const { return mDefaultGravityAccel; } + + /** + * @brief 更新弹簧及质点状态 + * @param[in] timeslice : 更新时间间隔 + * @param[in] method : 迭代更新算法 + * @param[in] limitDeformation : 弹簧弹性限度 + * @note 更新弹簧及质点状态 + * @return + * true update successfully \n + * false update failure \n + */ + bool Update(double timeslice, ItrMethod method = Euler, bool limitDeformation = false); - virtual bool PreUpdate(double timeslice, ItrMethod method = Euler, - dtkID iteration = 0); - virtual bool PostUpdate(ItrMethod method = Euler, dtkID iteration = 0); - // 更新迭代 - bool Update_iteration(double timeslice, ItrMethod method = Euler, - dtkID iteration = 0, bool limitDeformation = false); + /** + * @brief 单线程更新弹簧及质点状态 + * @param[in] timeslice : 更新时间间隔 + * @param[in] method : 迭代更新算法 + * @param[in] limitDeformation : 弹簧弹性限度 + * @note 单线程更新弹簧及质点状态 + * @return + * true update successfully \n + * false update failure \n + */ + bool Update_s(double timeslice, ItrMethod method = Euler, + bool limitDeformation = false); + virtual bool PreUpdate(double timeslice, ItrMethod method = Euler, + dtkID iteration = 0); + virtual bool PostUpdate(ItrMethod method = Euler, dtkID iteration = 0); - // 应用冲量 - bool ApplyImpulse(double timeslice); + // 更新迭代 + bool Update_iteration(double timeslice, ItrMethod method = Euler, + dtkID iteration = 0, bool limitDeformation = false); - /** - * @brief 更新弹簧状态，点位置，速度等 - * @param[in] timeslice : 更新时间间隔 - * @param[in] method : 迭代更新算法 - * @param[in] limitDeformation : 弹簧弹性限度 - * @note 更新弹簧状态，点位置，速度等 - * @return - * true update successfully \n - * false update failure \n - */ + // 应用冲量 + bool ApplyImpulse(double timeslice); - bool UpdateStrings(double timeslice, ItrMethod method = Euler, - dtkID iteration = 0, bool limitDeformation = false); + /** + * @brief 更新弹簧状态，点位置，速度等 + * @param[in] timeslice : 更新时间间隔 + * @param[in] method : 迭代更新算法 + * @param[in] limitDeformation : 弹簧弹性限度 + * @note 更新弹簧状态，点位置，速度等 + * @return + * true update successfully \n + * false update failure \n + */ - /** - * @brief 更新质点状态，点位置，速度等 - * @param[in] timeslice : 更新时间间隔 - * @param[in] method : 迭代更新算法 - * @param[in] limitDeformation : 弹簧弹性限度 - * @note 更新质点状态，点位置，速度等 - * @return - * true update successfully \n - * false update failure \n - */ + bool UpdateStrings(double timeslice, ItrMethod method = Euler, + dtkID iteration = 0, bool limitDeformation = false); - bool UpdateMassPoints(double timeslice, ItrMethod method = Euler, - dtkID iteration = 0); + /** + * @brief 更新质点状态，点位置，速度等 + * @param[in] timeslice : 更新时间间隔 + * @param[in] method : 迭代更新算法 + * @param[in] limitDeformation : 弹簧弹性限度 + * @note 更新质点状态，点位置，速度等 + * @return + * true update successfully \n + * false update failure \n + */ - // 从三角网格添加质量弹簧 + bool UpdateMassPoints(double timeslice, ItrMethod method = Euler, + dtkID iteration = 0); - void SetTriangleMesh(dtkStaticTriangleMesh::Ptr newTriangleMesh); + // 从三角网格添加质量弹簧 - // 寻找周围距离小于distance的twins点. + void SetTriangleMesh(dtkStaticTriangleMesh::Ptr newTriangleMesh); - size_t FindTwins(Ptr, double distance); - void AbandonTwins(); + // 寻找周围距离小于distance的twins点. - void RegisterLabel(dtkID label); + size_t FindTwins(Ptr, double distance); + void AbandonTwins(); - void TransportForce(double timeslice); + void RegisterLabel(dtkID label); - const dtkT3 &GetTransportForce(dtkID label) { - return mTransportForces[label]; - } + void TransportForce(double timeslice); - void ConvertImpulseToForce(double timeslice); + const dtkT3& GetTransportForce(dtkID label) { + return mTransportForces[label]; + } - const dtkT3 &GetImpulseForce() { return mImpulseForce; } + void ConvertImpulseToForce(double timeslice); - bool IsUnderControl() { return mUnderControl; } + const dtkT3& GetImpulseForce() { return mImpulseForce; } - void SetUnderControl(bool underControl) { mUnderControl = underControl; } + bool IsUnderControl() { return mUnderControl; } - dtkPhysSpring *GetSpringByPoints(dtkID2); + void SetUnderControl(bool underControl) { mUnderControl = underControl; } + + dtkPhysSpring* GetSpringByPoints(dtkID2); #ifdef DTK_CL - /** - * OpenCL related initialisations. - * Set up Context, Device list, Command Queue, Memory buffers - * Build CL kernel program executable - * @return 1 on success and 0 on failure - */ - bool Update_mt(double timeslice, ItrMethod method = Euler, - bool limitDeformation = false); - int SetupCL(); - int RunCLKernels(double timeslice, dtkID iteration); - bool Open(const char *fileName); - void InitialCLArray(); - void UpdateCopyToCLArray(); - void UpdateCopyBack(); - void CheckCLStatus(cl_int, std::string); - int WaitForEventAndRelease(cl_event *); - void UseMultiThread(const char *fileName); + /** + * OpenCL related initialisations. + * Set up Context, Device list, Command Queue, Memory buffers + * Build CL kernel program executable + * @return 1 on success and 0 on failure + */ + bool Update_mt(double timeslice, ItrMethod method = Euler, + bool limitDeformation = false); + int SetupCL(); + int RunCLKernels(double timeslice, dtkID iteration); + bool Open(const char* fileName); + void InitialCLArray(); + void UpdateCopyToCLArray(); + void UpdateCopyBack(); + void CheckCLStatus(cl_int, std::string); + int WaitForEventAndRelease(cl_event*); + void UseMultiThread(const char* fileName); #endif #ifdef DTK_DEBUG -public: - std::vector> mAltitudeSpringForces; - std::vector> mTotalSpringForces; + public: + std::vector> mAltitudeSpringForces; + std::vector> mTotalSpringForces; #endif -protected: - dtkPoints::Ptr mPts; /**< 点集 */ - std::vector mMassPoints; /**< 质点集 */ - std::vector mSprings; /**< 弹簧 */ + protected: + dtkPoints::Ptr mPts; /**< 点集 */ + std::vector mMassPoints; /**< 质点集 */ + std::vector mSprings; /**< 弹簧 */ - // 边集 - std::map + // 边集 + std::map mEdgeMap; /**< map from spring specified by dtkID2 to spring id */ - double mDefaultMass; /**< 质量 */ - double mDefaultStiff; /**< 弹簧刚性系数，弹性系数 */ - double mDefaultDamp; /**< 弹簧阻尼 */ - double mDefaultPointDamp; /**< 点阻尼 */ - double mDefaultPointResistence; /**< 阻力 */ - dtkDouble3 mDefaultGravityAccel; /**< 重力加速度 */ + double mDefaultMass; /**< 质量 */ + double mDefaultStiff; /**< 弹簧刚性系数，弹性系数 */ + double mDefaultDamp; /**< 弹簧阻尼 */ + double mDefaultPointDamp; /**< 点阻尼 */ + double mDefaultPointResistence; /**< 阻力 */ + double mDefaultTimeStep; /**< 时间步长 */ + dtkDouble3 mDefaultGravityAccel; /**< 重力加速度 */ - bool mUnderControl; /**< 弹簧受控 */ + bool mUnderControl; /**< 弹簧受控 */ - std::vector mLabels; // 标记点集, 可给予Transport力 + std::vector mLabels; // 标记点集, 可给予Transport力 - std::map> mTransportForces; // - dtkT3 mImpulseForce; // 瞬时力 + std::map> mTransportForces; // + dtkT3 mImpulseForce; // 瞬时力 #ifdef DTK_CL - bool mUseMultiThread; - // OpenCL - cl_float mTimesliceData; - cl_int mIterationData; - cl_int mNumberOfMassPoints; - - cl_int *mSpringVertexIDData; - cl_mem mSpringVertexIDBuffer; - cl_float *mSpringOriLengthData; - cl_mem mSpringOriLengthBuffer; - cl_float *mSpringStiffData; - cl_mem mSpringStiffBuffer; - cl_float *mSpringDampData; - cl_mem mSpringDampBuffer; - - cl_float *mMPPosData; - cl_mem mMPPosBuffer; - cl_float *mMPVelData; - cl_mem mMPVelBuffer; - cl_float *mMPAccelData; - cl_mem mMPAccelBuffer; - cl_float *mMPForceAccumData; - cl_mem mMPForceAccumBuffer; - cl_float *mMPGravityData; - cl_mem mMPGravityBuffer; - cl_float *mMPMassData; - cl_mem mMPMassBuffer; - cl_float *mMPPosBufferData; - cl_mem mMPPosBufferBuffer; - cl_float *mMPVelBufferData; - cl_mem mMPVelBufferBuffer; - cl_float *mMPAccelBufferData; - cl_mem mMPAccelBufferBuffer; - cl_float *mMPForceDecoratorData; - cl_mem mMPForceDecoratorBuffer; - cl_bool *mMPActiveData; - cl_mem mMPActiveBuffer; - - // new spring calculate - cl_int *mMPAdjacentVertexIDsData; - cl_mem mMPAdjacentVertexIDsBuffer; - cl_int *mMPNumberOfAdjacentVertexData; - cl_mem mMPNumberOfAdjacentVertexBuffer; - cl_int *mMPStartAdjacentVertexIDData; - cl_mem mMPStartAdjacentVertexIDBuffer; - - cl_int *mMPAdjacentSpringIDsData; - cl_mem mMPAdjacentSpringIDsBuffer; - - size_t maxWorkGroupSize; - size_t *maxWorkItemSizes; - - cl_context context; - cl_platform_id platform; - cl_device_id device; - - cl_command_queue queue; - cl_program program; - cl_kernel kernelSpring; - cl_kernel kernelMP; - size_t global_work_size; - size_t work_size; - std::string source_; + bool mUseMultiThread; + // OpenCL + cl_float mTimesliceData; + cl_int mIterationData; + cl_int mNumberOfMassPoints; + + cl_int* mSpringVertexIDData; + cl_mem mSpringVertexIDBuffer; + cl_float* mSpringOriLengthData; + cl_mem mSpringOriLengthBuffer; + cl_float* mSpringStiffData; + cl_mem mSpringStiffBuffer; + cl_float* mSpringDampData; + cl_mem mSpringDampBuffer; + + cl_float* mMPPosData; + cl_mem mMPPosBuffer; + cl_float* mMPVelData; + cl_mem mMPVelBuffer; + cl_float* mMPAccelData; + cl_mem mMPAccelBuffer; + cl_float* mMPForceAccumData; + cl_mem mMPForceAccumBuffer; + cl_float* mMPGravityData; + cl_mem mMPGravityBuffer; + cl_float* mMPMassData; + cl_mem mMPMassBuffer; + cl_float* mMPPosBufferData; + cl_mem mMPPosBufferBuffer; + cl_float* mMPVelBufferData; + cl_mem mMPVelBufferBuffer; + cl_float* mMPAccelBufferData; + cl_mem mMPAccelBufferBuffer; + cl_float* mMPForceDecoratorData; + cl_mem mMPForceDecoratorBuffer; + cl_bool* mMPActiveData; + cl_mem mMPActiveBuffer; + + // new spring calculate + cl_int* mMPAdjacentVertexIDsData; + cl_mem mMPAdjacentVertexIDsBuffer; + cl_int* mMPNumberOfAdjacentVertexData; + cl_mem mMPNumberOfAdjacentVertexBuffer; + cl_int* mMPStartAdjacentVertexIDData; + cl_mem mMPStartAdjacentVertexIDBuffer; + + cl_int* mMPAdjacentSpringIDsData; + cl_mem mMPAdjacentSpringIDsBuffer; + + size_t maxWorkGroupSize; + size_t* maxWorkItemSizes; + + cl_context context; + cl_platform_id platform; + cl_device_id device; + + cl_command_queue queue; + cl_program program; + cl_kernel kernelSpring; + cl_kernel kernelMP; + size_t global_work_size; + size_t work_size; + std::string source_; #endif -}; + }; }; // namespace dtk #endif /* SIMPLEPHYSICSENGINE_DTKPHYSMASSSPRING_H */ diff --git a/src/include/dtkPhysSpring.h b/src/include/dtkPhysSpring.h index 0d03e54..d1735e0 100644 --- a/src/include/dtkPhysSpring.h +++ b/src/include/dtkPhysSpring.h @@ -21,33 +21,36 @@ #include "dtkPhysMassPoint.h" namespace dtk { -class dtkPhysSpring { -public: - dtkPhysSpring(dtkPhysMassPoint *p1, dtkPhysMassPoint *p2, - const double &stiff = 5, const double &damp = 0); - - dtkPhysMassPoint *GetFirstVertex() { return mVerteces[0]; } - dtkPhysMassPoint *GetSecondVertex() { return mVerteces[1]; } - - // 更新力到质点,迭代更新 - bool Update(double timeslice, ItrMethod method = Euler, dtkID iteration = 0, - bool limitDeformation = false); - void SetStiffness(double newStiffness) { mStiffness = newStiffness; } - void SetDamp(double newDamp) { mDamp = newDamp; } - double GetOriLength() { return mOriLength; } - double GetStiffness() { return mStiffness; } - double GetDamp() { return mDamp; } - -public: - virtual ~dtkPhysSpring(); - -private: - dtkPhysMassPoint *mVerteces[2]; /**< 弹簧两个质点 */ - - double mOriLength; /**< 长度 */ - double mStiffness; /**< 刚度 */ - double mDamp; /**< 阻尼 */ -}; + class dtkPhysSpring { + public: + dtkPhysSpring(dtkPhysMassPoint* p1, dtkPhysMassPoint* p2, + const double& stiff = 5, const double& damp = 0); + + dtkPhysMassPoint* GetFirstVertex() { return mVerteces[0]; } + dtkPhysMassPoint* GetSecondVertex() { return mVerteces[1]; } + + // 更新力到质点,迭代更新 + bool Update(double timeslice, ItrMethod method = Euler, dtkID iteration = 0, + bool limitDeformation = false); + void SetStiffness(double newStiffness) { mStiffness = newStiffness; } + void SetDamp(double newDamp) { mDamp = newDamp; } + double GetOriLength() { return mOriLength; } + void SetRestLength(const double& newRestLength) { mRestLength = newRestLength; } + double GetRestLength() { return mRestLength; } + double GetStiffness() { return mStiffness; } + double GetDamp() { return mDamp; } + + public: + virtual ~dtkPhysSpring(); + + private: + dtkPhysMassPoint* mVerteces[2]; /**< 弹簧两个质点 */ + + double mOriLength; /**< 原始长度 */ + double mRestLength; /**< 静息长度 */ + double mStiffness; /**< 刚度 */ + double mDamp; /**< 阻尼 */ + }; } // namespace dtk #endif /* SIMPLEPHYSICSENGINE_DTKPHYSSPRING_H */ diff --git a/src/include/dtkStaticTriangleMesh.h b/src/include/dtkStaticTriangleMesh.h index f3ec7d1..d3efae1 100644 --- a/src/include/dtkStaticTriangleMesh.h +++ b/src/include/dtkStaticTriangleMesh.h @@ -26,148 +26,156 @@ #include "dtkConfig.h" #include "dtkIDTypes.h" #include "dtkPoints.h" +#include "dtkMatrix.h" +#include "dtkTx.h" namespace dtk { -// This class is implemented following "Compact Array-Based Mesh Data -// Structures" -/** - * @class - * @brief 三角网格 - * @author <> - * @note This class is implemented following "Compact Array-Based Mesh Data - * Structures" - * - */ -class dtkStaticTriangleMesh : public boost::noncopyable { -public: - typedef std::shared_ptr Ptr; - - static dtkStaticTriangleMesh::Ptr New() { - return dtkStaticTriangleMesh::Ptr(new dtkStaticTriangleMesh()); - } - -public: - virtual ~dtkStaticTriangleMesh(); - - void SetPoints(dtkPoints::Ptr points); - dtkPoints::Ptr GetPoints(); - - const GK::Point3 &GetPoint(dtkID id) const; - size_t GetNumberOfPoints() { return mPts->GetNumberOfPoints(); } - void SetPoint(dtkID id, const GK::Point3 &coord); - - void Clear(); - dtkID InsertTriangle(dtkID v0, dtkID v1, dtkID v2); - dtkID InsertTriangle(dtkID v[3]); - dtkID InsertTriangle(const dtkID3 &tri); - - dtkID InsertTriangleNotRepeat(dtkID v0, dtkID v1, dtkID v2); - dtkID InsertTriangleNotRepeat(dtkID v[3]); - dtkID InsertTriangleNotRepeat(const dtkID3 &tri); - - bool IsRepeat(dtkID v0, dtkID v1, dtkID v2); - - dtkID RemoveTriangle(dtkID v0, dtkID v1, dtkID v2); - dtkID RemoveTriangle(dtkID v[3]); - dtkID RemoveTriangle(const dtkID3 &tri); - - dtkID RemoveTriangleDisordered(dtkID v0, dtkID v1, dtkID v2); - - dtkID modifyElement(dtkID srcIndex, int srcOrder, dtkID dest); - dtkID modifyElement(dtkID srcIndex, dtkID3 dest); - - void Modified(); - bool IsModified() const; - void Rebuild(); - - bool IsValidVertex(dtkID &vertexID) const { - if (vertexID < 0 || vertexID > mV2E.size()) - return false; - return (mV2E[vertexID] != dtkErrorID); - } - - const std::vector &GetECTable() const { return mEC; } - const std::vector &GetV2ETable() const { return mV2E; } - const std::vector &GetE2ETable() const { return mE2E; } - const std::vector &GetB2ETable() const { return mB2E; } - -private: - dtkStaticTriangleMesh(); - -public: - // Generate the HF. - static dtkID EncodeHE(dtkID tri, dtkID idx) { return (tri << 2) | idx; } - - // Get the components of a HE. - // Return true if it is on the boundary - static bool DecodeHE(dtkID he, dtkID &tri, dtkID &idx) { - static const dtkWORD mask = 0xFFFFFFFC; - tri = (he & mask) >> 2; - idx = (he & ~mask); - - return (idx == 0); - } - - static bool IsBoundaryHE(const dtkID &he) { - static const dtkWORD mask = 0x00000003; - return ((he & mask) == 0); - } - - static dtkID2 GetHEVertex(const dtkID &he, const std::vector &ec, - const std::vector &b2e) { - dtkID2 retVal; - dtkID c, i; - dtkStaticTriangleMesh::DecodeHE(he, c, i); - - if (i == 0) { - dtkID twin = b2e[c]; - retVal = GetHEVertex(twin, ec, b2e); - std::swap(retVal.a, retVal.b); - } else { - retVal.a = ec[c][i - 1]; - retVal.b = ec[c][i % 3]; + // This class is implemented following "Compact Array-Based Mesh Data + // Structures" + /** + * @class + * @brief 三角网格 + * @author <> + * @note This class is implemented following "Compact Array-Based Mesh Data + * Structures" + * + */ + class dtkStaticTriangleMesh : public boost::noncopyable { + public: + typedef std::shared_ptr Ptr; + + static dtkStaticTriangleMesh::Ptr New() { + return dtkStaticTriangleMesh::Ptr(new dtkStaticTriangleMesh()); } - return retVal; - } + public: + virtual ~dtkStaticTriangleMesh(); + + void SetPoints(dtkPoints::Ptr points); + dtkPoints::Ptr GetPoints(); + + const GK::Point3& GetPoint(dtkID id) const; + size_t GetNumberOfPoints() { return mPts->GetNumberOfPoints(); } + size_t GetNumberOfTriangles() { return mEC.size(); } + void SetPoint(dtkID id, const GK::Point3& coord); + + void Clear(); + dtkID InsertTriangle(dtkID v0, dtkID v1, dtkID v2); + dtkID InsertTriangle(dtkID v[3]); + dtkID InsertTriangle(const dtkID3& tri); + + dtkID InsertTriangleNotRepeat(dtkID v0, dtkID v1, dtkID v2); + dtkID InsertTriangleNotRepeat(dtkID v[3]); + dtkID InsertTriangleNotRepeat(const dtkID3& tri); + + bool IsRepeat(dtkID v0, dtkID v1, dtkID v2); - static dtkID3 SortTriVertex(dtkID v0, dtkID v1, dtkID v2) { - if (v2 < v1) - std::swap(v1, v2); - if (v1 < v0) - std::swap(v0, v1); - if (v2 < v1) - std::swap(v1, v2); + dtkID RemoveTriangle(dtkID v0, dtkID v1, dtkID v2); + dtkID RemoveTriangle(dtkID v[3]); + dtkID RemoveTriangle(const dtkID3& tri); - return dtkID3(v0, v1, v2); - } + dtkID RemoveTriangleDisordered(dtkID v0, dtkID v1, dtkID v2); -private: - static dtkDWORD Merge(const dtkID &from, const dtkID &to) { - dtkDWORD key(0); + dtkID modifyElement(dtkID srcIndex, int srcOrder, dtkID dest); + dtkID modifyElement(dtkID srcIndex, dtkID3 dest); - key = ((dtkDWORD)from) << 32; - key |= to; + void update_normals(); - return key; - } + void Modified(); + bool IsModified() const; + void Rebuild(); - static void Split(const dtkDWORD &key, dtkID &from, dtkID &to) { - static const dtkDWORD mask = 0x00000000FFFFFFFF; + bool IsValidVertex(dtkID& vertexID) const { + if (vertexID < 0 || vertexID > mV2E.size()) + return false; + return (mV2E[vertexID] != dtkErrorID); + } + + const std::vector& GetECTable() const { return mEC; } + const std::vector& GetV2ETable() const { return mV2E; } + const std::vector& GetE2ETable() const { return mE2E; } + const std::vector& GetB2ETable() const { return mB2E; } + const std::vector& GetVertexNormals() const { return mVertexNormals; } + + private: + dtkStaticTriangleMesh(); + + public: + // Generate the HF. + static dtkID EncodeHE(dtkID tri, dtkID idx) { return (tri << 2) | idx; } + + // Get the components of a HE. + // Return true if it is on the boundary + static bool DecodeHE(dtkID he, dtkID& tri, dtkID& idx) { + static const dtkWORD mask = 0xFFFFFFFC; + tri = (he & mask) >> 2; + idx = (he & ~mask); + + return (idx == 0); + } - from = static_cast((key & ~mask) >> 32); - to = static_cast(key & mask); - } + static bool IsBoundaryHE(const dtkID& he) { + static const dtkWORD mask = 0x00000003; + return ((he & mask) == 0); + } + + static dtkID2 GetHEVertex(const dtkID& he, const std::vector& ec, + const std::vector& b2e) { + dtkID2 retVal; + dtkID c, i; + dtkStaticTriangleMesh::DecodeHE(he, c, i); + + if (i == 0) { + dtkID twin = b2e[c]; + retVal = GetHEVertex(twin, ec, b2e); + std::swap(retVal.a, retVal.b); + } + else { + retVal.a = ec[c][i - 1]; + retVal.b = ec[c][i % 3]; + } + + return retVal; + } + + static dtkID3 SortTriVertex(dtkID v0, dtkID v1, dtkID v2) { + if (v2 < v1) + std::swap(v1, v2); + if (v1 < v0) + std::swap(v0, v1); + if (v2 < v1) + std::swap(v1, v2); + + return dtkID3(v0, v1, v2); + } + + private: + static dtkDWORD Merge(const dtkID& from, const dtkID& to) { + dtkDWORD key(0); + + key = ((dtkDWORD)from) << 32; + key |= to; + + return key; + } + + static void Split(const dtkDWORD& key, dtkID& from, dtkID& to) { + static const dtkDWORD mask = 0x00000000FFFFFFFF; + + from = static_cast((key & ~mask) >> 32); + to = static_cast(key & mask); + } -private: - bool mModified; - dtkPoints::Ptr mPts; + private: + bool mModified; + dtkPoints::Ptr mPts; - std::vector mEC; // present a triangle - std::vector mV2E; - std::vector mE2E; - std::vector mB2E; -}; + std::vector mVertexNormals; + std::vector mEC; // present a triangle + std::vector mV2E; + std::vector mE2E; + std::vector mB2E; + }; } // namespace dtk #endif /* SIMPLEPHYSICSENGINE_DTKSTATICTRIANGLEMESH_H */ diff --git a/src/math/include/dtkMatrix.h b/src/math/include/dtkMatrix.h index f2f3597..97e6490 100644 --- a/src/math/include/dtkMatrix.h +++ b/src/math/include/dtkMatrix.h @@ -27,288 +27,288 @@ #include "dtkConfig.h" namespace dtk { -/** - * @class - * @brief 矩阵 - * @author <> - * @note - * 实现任意维的矩阵 - */ -template > struct dtkMatrix { - // STL-friendly typedefs - - typedef typename MatrixT::iterator iterator; - typedef typename MatrixT::const_iterator const_iterator; - typedef typename MatrixT::size_type size_type; - typedef long difference_type; - typedef T &reference; - typedef const T &const_reference; - typedef T value_type; - typedef T *pointer; - typedef const T *const_pointer; - typedef typename MatrixT::reverse_iterator reverse_iterator; - typedef typename MatrixT::const_reverse_iterator const_reverse_iterator; - - // the actual representation - - unsigned int ni, nj; - MatrixT a; - - // the interface - - dtkMatrix(void) : ni(0), nj(0) {} - - dtkMatrix(int ni_, int nj_) : ni(ni_), nj(nj_), a(ni_ * nj_) { - assert(ni_ >= 0 && nj >= 0); - } - - dtkMatrix(int ni_, int nj_, MatrixT &a_) : ni(ni_), nj(nj_), a(a_) { - assert(ni_ >= 0 && nj >= 0); - } - - dtkMatrix(int ni_, int nj_, const T &value) - : ni(ni_), nj(nj_), a(ni_ * nj_, value) { - assert(ni_ >= 0 && nj >= 0); - } - - dtkMatrix(int ni_, int nj_, const T &value, size_type max_n_) - : ni(ni_), nj(nj_), a(ni_ * nj_, value, max_n_) { - assert(ni_ >= 0 && nj >= 0); - } - - void init(const T *values) { - for (unsigned int i = 0; i < ni * nj; i++) - a[i] = values[i]; - } - - // dtkMatrix(int ni_, int nj_, T* data_) - // : ni(ni_), nj(nj_), a(ni_*nj_, data_) - //{ assert(ni_>=0 && nj>=0); } - - // dtkMatrix(int ni_, int nj_, T* data_, size_type max_n_) - // : ni(ni_), nj(nj_), a(ni_*nj_, data_, max_n_) - //{ assert(ni_>=0 && nj>=0); } - - template - dtkMatrix(dtkMatrix &other) + /** + * @class + * @brief 矩阵 + * @author <> + * @note + * 实现任意维的矩阵 + */ + template > struct dtkMatrix { + // STL-friendly typedefs + + typedef typename MatrixT::iterator iterator; + typedef typename MatrixT::const_iterator const_iterator; + typedef typename MatrixT::size_type size_type; + typedef long difference_type; + typedef T& reference; + typedef const T& const_reference; + typedef T value_type; + typedef T* pointer; + typedef const T* const_pointer; + typedef typename MatrixT::reverse_iterator reverse_iterator; + typedef typename MatrixT::const_reverse_iterator const_reverse_iterator; + + // the actual representation + + unsigned int ni, nj; + MatrixT a; + + // the interface + + dtkMatrix(void) : ni(0), nj(0) {} + + dtkMatrix(int ni_, int nj_) : ni(ni_), nj(nj_), a(ni_* nj_) { + assert(ni_ >= 0 && nj >= 0); + } + + dtkMatrix(int ni_, int nj_, MatrixT& a_) : ni(ni_), nj(nj_), a(a_) { + assert(ni_ >= 0 && nj >= 0); + } + + dtkMatrix(int ni_, int nj_, const T& value) + : ni(ni_), nj(nj_), a(ni_* nj_, value) { + assert(ni_ >= 0 && nj >= 0); + } + + dtkMatrix(int ni_, int nj_, const T& value, size_type max_n_) + : ni(ni_), nj(nj_), a(ni_* nj_, value, max_n_) { + assert(ni_ >= 0 && nj >= 0); + } + + void init(const T* values) { + for (unsigned int i = 0; i < ni * nj; i++) + a[i] = values[i]; + } + + // dtkMatrix(int ni_, int nj_, T* data_) + // : ni(ni_), nj(nj_), a(ni_*nj_, data_) + //{ assert(ni_>=0 && nj>=0); } + + // dtkMatrix(int ni_, int nj_, T* data_, size_type max_n_) + // : ni(ni_), nj(nj_), a(ni_*nj_, data_, max_n_) + //{ assert(ni_>=0 && nj>=0); } + + template + dtkMatrix(dtkMatrix& other) : ni(other.ni), nj(other.nj), a(other.a) {} - ~dtkMatrix(void) { + ~dtkMatrix(void) { #ifndef NDEBUG - ni = nj = 0; + ni = nj = 0; #endif - } - - const T &operator()(unsigned int i, unsigned int j) const { - assert(i >= 0 && i < ni && j >= 0 && j < nj); - return a[i + ni * j]; - } - - T &operator()(unsigned int i, unsigned int j) { - assert(i >= 0 && i < ni && j >= 0 && j < nj); - return a[i + ni * j]; - } - - const T &operator[](unsigned int i) const { - assert(i >= 0 && i < ni * nj); - return a[i]; - } - - T &operator[](unsigned int i) { - assert(i >= 0 && i < ni * nj); - return a[i]; - } - - bool operator==(const dtkMatrix &x) const { - return ni == x.ni && nj == x.nj && a == x.a; - } - - bool operator!=(const dtkMatrix &x) const { - return ni != x.ni || nj != x.nj || a != x.a; - } - - bool operator<(const dtkMatrix &x) const { - if (ni < x.ni) - return true; - else if (ni > x.ni) - return false; - if (nj < x.nj) - return true; - else if (nj > x.nj) - return false; - return a < x.a; - } - - bool operator>(const dtkMatrix &x) const { - if (ni > x.ni) - return true; - else if (ni < x.ni) - return false; - if (nj > x.nj) - return true; - else if (nj < x.nj) - return false; - return a > x.a; - } - - bool operator<=(const dtkMatrix &x) const { - if (ni < x.ni) - return true; - else if (ni > x.ni) - return false; - if (nj < x.nj) - return true; - else if (nj > x.nj) - return false; - return a <= x.a; - } - - bool operator>=(const dtkMatrix &x) const { - if (ni > x.ni) - return true; - else if (ni < x.ni) - return false; - if (nj > x.nj) - return true; - else if (nj < x.nj) - return false; - return a >= x.a; - } - - void assign(const T &value) { a.assign(value); } - - void assign(int ni_, int nj_, const T &value) { - a.assign(ni_ * nj_, value); - ni = ni_; - nj = nj_; - } - - void assign(int ni_, int nj_, const T *copydata) { - a.assign(ni_ * nj_, copydata); - ni = ni_; - nj = nj_; - } - - const T &at(int i, int j) const { - assert(i >= 0 && i < ni && j >= 0 && j < nj); - return a[i + ni * j]; - } - - T &at(int i, int j) { - assert(i >= 0 && i < ni && j >= 0 && j < nj); - return a[i + ni * j]; - } - - const T &back(void) const { - assert(a.size()); - return a.back(); - } - - T &back(void) { - assert(a.size()); - return a.back(); - } - - const_iterator begin(void) const { return a.begin(); } - - iterator begin(void) { return a.begin(); } - - size_type capacity(void) const { return a.capacity(); } - - void clear(void) { - a.clear(); - ni = nj = 0; - } - - bool empty(void) const { return a.empty(); } - - const_iterator end(void) const { return a.end(); } - - iterator end(void) { return a.end(); } - - void fill(int ni_, int nj_, const T &value) { - a.fill(ni_ * nj_, value); - ni = ni_; - nj = nj_; - } - - const T &front(void) const { - assert(a.size()); - return a.front(); - } - - T &front(void) { - assert(a.size()); - return a.front(); - } - - size_type max_size(void) const { return a.max_size(); } - - reverse_iterator rbegin(void) { return reverse_iterator(end()); } - - const_reverse_iterator rbegin(void) const { - return const_reverse_iterator(end()); - } - - reverse_iterator rend(void) { return reverse_iterator(begin()); } - - const_reverse_iterator rend(void) const { - return const_reverse_iterator(begin()); - } - - void reserve(int reserve_ni, int reserve_nj) { - a.reserve(reserve_ni * reserve_nj); - } - - void resize(int ni_, int nj_) { - assert(ni_ >= 0 && nj_ >= 0); - a.resize(ni_ * nj_); - ni = ni_; - nj = nj_; - } - - void resize(int ni_, int nj_, const T &value) { - assert(ni_ >= 0 && nj_ >= 0); - a.resize(ni_ * nj_, value); - ni = ni_; - nj = nj_; - } + } + + const T& operator()(unsigned int i, unsigned int j) const { + assert(i >= 0 && i < ni && j >= 0 && j < nj); + return a[i + ni * j]; + } + + T& operator()(unsigned int i, unsigned int j) { + assert(i >= 0 && i < ni && j >= 0 && j < nj); + return a[i + ni * j]; + } + + const T& operator[](unsigned int i) const { + assert(i >= 0 && i < ni * nj); + return a[i]; + } + + T& operator[](unsigned int i) { + assert(i >= 0 && i < ni * nj); + return a[i]; + } + + bool operator==(const dtkMatrix& x) const { + return ni == x.ni && nj == x.nj && a == x.a; + } + + bool operator!=(const dtkMatrix& x) const { + return ni != x.ni || nj != x.nj || a != x.a; + } + + bool operator<(const dtkMatrix& x) const { + if (ni < x.ni) + return true; + else if (ni > x.ni) + return false; + if (nj < x.nj) + return true; + else if (nj > x.nj) + return false; + return a < x.a; + } + + bool operator>(const dtkMatrix& x) const { + if (ni > x.ni) + return true; + else if (ni < x.ni) + return false; + if (nj > x.nj) + return true; + else if (nj < x.nj) + return false; + return a > x.a; + } + + bool operator<=(const dtkMatrix& x) const { + if (ni < x.ni) + return true; + else if (ni > x.ni) + return false; + if (nj < x.nj) + return true; + else if (nj > x.nj) + return false; + return a <= x.a; + } + + bool operator>=(const dtkMatrix& x) const { + if (ni > x.ni) + return true; + else if (ni < x.ni) + return false; + if (nj > x.nj) + return true; + else if (nj < x.nj) + return false; + return a >= x.a; + } + + void assign(const T& value) { a.assign(value); } + + void assign(int ni_, int nj_, const T& value) { + a.assign(ni_ * nj_, value); + ni = ni_; + nj = nj_; + } + + void assign(int ni_, int nj_, const T* copydata) { + a.assign(ni_ * nj_, copydata); + ni = ni_; + nj = nj_; + } + + const T& at(int i, int j) const { + assert(i >= 0 && i < ni && j >= 0 && j < nj); + return a[i + ni * j]; + } + + T& at(int i, int j) { + assert(i >= 0 && i < ni && j >= 0 && j < nj); + return a[i + ni * j]; + } + + const T& back(void) const { + assert(a.size()); + return a.back(); + } + + T& back(void) { + assert(a.size()); + return a.back(); + } + + const_iterator begin(void) const { return a.begin(); } + + iterator begin(void) { return a.begin(); } + + size_type capacity(void) const { return a.capacity(); } + + void clear(void) { + a.clear(); + ni = nj = 0; + } + + bool empty(void) const { return a.empty(); } + + const_iterator end(void) const { return a.end(); } + + iterator end(void) { return a.end(); } + + void fill(int ni_, int nj_, const T& value) { + a.fill(ni_ * nj_, value); + ni = ni_; + nj = nj_; + } + + const T& front(void) const { + assert(a.size()); + return a.front(); + } + + T& front(void) { + assert(a.size()); + return a.front(); + } + + size_type max_size(void) const { return a.max_size(); } + + reverse_iterator rbegin(void) { return reverse_iterator(end()); } + + const_reverse_iterator rbegin(void) const { + return const_reverse_iterator(end()); + } + + reverse_iterator rend(void) { return reverse_iterator(begin()); } + + const_reverse_iterator rend(void) const { + return const_reverse_iterator(begin()); + } + + void reserve(int reserve_ni, int reserve_nj) { + a.reserve(reserve_ni * reserve_nj); + } + + void resize(int ni_, int nj_) { + assert(ni_ >= 0 && nj_ >= 0); + a.resize(ni_ * nj_); + ni = ni_; + nj = nj_; + } + + void resize(int ni_, int nj_, const T& value) { + assert(ni_ >= 0 && nj_ >= 0); + a.resize(ni_ * nj_, value); + ni = ni_; + nj = nj_; + } - void set_zero(void) { a.set_zero(); } + void set_zero(void) { a.set_zero(); } - size_type size(void) const { return a.size(); } + size_type size(void) const { return a.size(); } - void swap(dtkMatrix &x) { - std::swap(ni, x.ni); - std::swap(nj, x.nj); - a.swap(x.a); - } + void swap(dtkMatrix& x) { + std::swap(ni, x.ni); + std::swap(nj, x.nj); + a.swap(x.a); + } - void trim(void) { a.trim(); } -}; + void trim(void) { a.trim(); } + }; -// some common arrays -typedef dtkMatrix> dtkMatrixDouble; -typedef dtkMatrix> dtkMatrixFloat; -typedef dtkMatrix> dtkMatrixLLong; -typedef dtkMatrix> + // some common arrays + typedef dtkMatrix> dtkMatrixDouble; + typedef dtkMatrix> dtkMatrixFloat; + typedef dtkMatrix> dtkMatrixLLong; + typedef dtkMatrix> dtkMatrixULLong; -typedef dtkMatrix> dtkMatrixInt; -typedef dtkMatrix> dtkMatrixUInt; -typedef dtkMatrix> dtkMatrixShort; -typedef dtkMatrix> dtkMatrixUShort; -typedef dtkMatrix> dtkMatrixChar; -typedef dtkMatrix> dtkMatrixUChar; - -typedef glm::mat2x2 dtkMatrix22; -typedef glm::mat2x3 dtkMatrix23; -typedef glm::mat3x2 dtkMatrix32; -typedef glm::mat3x3 dtkMatrix33; -typedef glm::mat2x4 dtkMatrix24; -typedef glm::mat4x2 dtkMatrix42; -typedef glm::mat3x4 dtkMatrix34; -typedef glm::mat4x3 dtkMatrix43; -typedef glm::mat4x4 dtkMatrix44; + typedef dtkMatrix> dtkMatrixInt; + typedef dtkMatrix> dtkMatrixUInt; + typedef dtkMatrix> dtkMatrixShort; + typedef dtkMatrix> dtkMatrixUShort; + typedef dtkMatrix> dtkMatrixChar; + typedef dtkMatrix> dtkMatrixUChar; + + typedef glm::mat2x2 dtkMatrix22; + typedef glm::mat2x3 dtkMatrix23; + typedef glm::mat3x2 dtkMatrix32; + typedef glm::mat3x3 dtkMatrix33; + typedef glm::mat2x4 dtkMatrix24; + typedef glm::mat4x2 dtkMatrix42; + typedef glm::mat3x4 dtkMatrix34; + typedef glm::mat4x3 dtkMatrix43; + typedef glm::mat4x4 dtkMatrix44; } // namespace dtk diff --git a/src/math/include/dtkSparseMatrix.h b/src/math/include/dtkSparseMatrix.h new file mode 100644 index 0000000..644aaa1 --- /dev/null +++ b/src/math/include/dtkSparseMatrix.h @@ -0,0 +1,129 @@ +/** + * @file dtkSparseMatrix.h + * @brief dtkSparseMatrix 头文件 + * @author Chance.H (chancewithchance@outlook.com) + * @version 1.0 + * @date 2024-7-30 + * @copyright MIT LICENSE + * https://github.com/Simple-XX/SimplePhysicsEngine + */ +#ifndef SIMPLEPHYSICSENGINE_DTKSPARSEMATRIX_H +#define SIMPLEPHYSICSENGINE_DTKSPARSEMATRIX_H +#include +#include +#include + +namespace dtk { + + /** + * @class + * @brief 稀疏矩阵 + * @author <> + * @note + * 实现稀疏矩阵 + */ + template + class dtkSparseMatrix { + public: + typedef unsigned int Index; + typedef T ValueType; + + struct Triplet { + Index row, col; + ValueType value; + + Triplet(Index r, Index c, const ValueType& v) : row(r), col(c), value(v) {} + }; + + dtkSparseMatrix() : rows_(0), cols_(0) {} + + dtkSparseMatrix(Index rows, Index cols) : rows_(rows), cols_(cols) {} + + void resize(Index rows, Index cols) { + rows_ = rows; + cols_ = cols; + triplets_.clear(); + } + + void insert(Index row, Index col, const ValueType& value) { + assert(row < rows_ && col < cols_); + triplets_.emplace_back(row, col, value); + } + + void setFromTriplets(const std::vector& triplets) { + triplets_ = triplets; + } + + void reserve(Index reserveSize) { + triplets_.reserve(reserveSize); + } + + void clear() { + triplets_.clear(); + } + + Index rows() const { return rows_; } + Index cols() const { return cols_; } + Index nonZeros() const { return triplets_.size(); } + + const std::vector& triplets() const { return triplets_; } + + dtkSparseMatrix transpose() const { + dtkSparseMatrix result(cols_, rows_); + for (const auto& triplet : triplets_) { + result.insert(triplet.col, triplet.row, triplet.value); + } + return result; + } + + template + std::vector operator*(const VectorT& rhs) const { + assert(rhs.size() == cols_); + std::vector result(rows_, T(0)); + for (const auto& triplet : triplets_) { + result[triplet.row] += triplet.value * rhs[triplet.col]; + } + return result; + } + + dtkSparseMatrix operator+(const dtkSparseMatrix& rhs) const { + assert(rows_ == rhs.rows_ && cols_ == rhs.cols_); + dtkSparseMatrix result(rows_, cols_); + result.triplets_ = triplets_; + for (const auto& triplet : rhs.triplets_) { + auto it = std::find_if(result.triplets_.begin(), result.triplets_.end(), + [&](const Triplet& t) { return t.row == triplet.row && t.col == triplet.col; }); + if (it != result.triplets_.end()) { + it->value += triplet.value; + } + else { + result.triplets_.emplace_back(triplet.row, triplet.col, triplet.value); + } + } + return result; + } + + dtkSparseMatrix operator-(const dtkSparseMatrix& rhs) const { + assert(rows_ == rhs.rows_ && cols_ == rhs.cols_); + dtkSparseMatrix result(rows_, cols_); + result.triplets_ = triplets_; + for (const auto& triplet : rhs.triplets_) { + auto it = std::find_if(result.triplets_.begin(), result.triplets_.end(), + [&](const Triplet& t) { return t.row == triplet.row && t.col == triplet.col; }); + if (it != result.triplets_.end()) { + it->value -= triplet.value; + } + else { + result.triplets_.emplace_back(triplet.row, triplet.col, -triplet.value); + } + } + return result; + } + + private: + Index rows_, cols_; + std::vector triplets_; + }; + +} // namespace dtk +#endif /* SIMPLEPHYSICSENGINE_DTKSPARSEMATRIX_H */ \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index af2f399..30e9c63 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ enable_language(CXX) include_directories( ${SimplePhysicsEngine_SOURCE_DIR}/src/include + ${SimplePhysicsEngine_SOURCE_DIR}/src/collision_detect/include ${SimplePhysicsEngine_SOURCE_DIR}/src/math/include ${SimplePhysicsEngine_SOURCE_DIR}/src/physics/include ) diff --git a/test/system_test/CMakeLists.txt b/test/system_test/CMakeLists.txt index 68b5296..e1e5a56 100644 --- a/test/system_test/CMakeLists.txt +++ b/test/system_test/CMakeLists.txt @@ -12,4 +12,5 @@ project( add_subdirectory(demo2d) add_subdirectory(demo3d) +add_subdirectory(MassSpring3D) # add_subdirectory(demo_old) diff --git a/test/system_test/MassSpring3D/CMakeLists.txt b/test/system_test/MassSpring3D/CMakeLists.txt new file mode 100644 index 0000000..bd253d3 --- /dev/null +++ b/test/system_test/MassSpring3D/CMakeLists.txt @@ -0,0 +1,26 @@ +# This file is a part of Simple-XX/SimplePhysicsEngine +# (https://github.com/Simple-XX/SimplePhysicsEngine). +# +# CMakeLists.txt for Simple-XX/SimplePhysicsEngine. + +project(MassSpring3D) + +add_executable(${PROJECT_NAME} + main.cpp + ClothSimulation.cpp + Shader.cpp + Renderer.cpp + dtkPhysMassSpringSolver.cpp +) + +target_include_directories(${PROJECT_NAME} PRIVATE + include +) + +target_link_libraries(${PROJECT_NAME} GLEW::GLEW) + +if (APPLE) + target_compile_definitions(${PROJECT_NAME} PRIVATE + GL_SILENCE_DEPRECATION + ) +endif () diff --git a/test/system_test/MassSpring3D/ClothSimulation.cpp b/test/system_test/MassSpring3D/ClothSimulation.cpp new file mode 100644 index 0000000..0eb417e --- /dev/null +++ b/test/system_test/MassSpring3D/ClothSimulation.cpp @@ -0,0 +1,325 @@ +/** + * @file ClothSimulation.cpp + * @brief ClothSimulation 实现 + * @author Chance + * @version 1.0 + * @date 2024-07-12 + * @copyright MIT LICENSE + * https://github.com/Simple-XX/SimplePhysicsEngine + * @par change log: + * + *
DateAuthorDescription + *
2024-07-12Chance创建文件 + *
+ */ +#include "ClothSimulation.h" + +ClothSimulation::ClothSimulation(const unsigned int& windowWidth, const unsigned int& windowHeight, const dtk::dtkDouble3& gravity) : Scene(windowWidth, windowHeight), _gravity(gravity) { +}; + +const ClothSimulation::ClothMesh ClothSimulation::GetClothMesh() const { + return _cloth_mesh; +}; + +void ClothSimulation::move(const dtk::dtkDouble2& v) { + +}; + +void ClothSimulation::CleanUp() { + delete g_phongShader; + delete g_pickShader; + delete g_render_target; +}; + +void ClothSimulation::Init() { + InitShader(); + InitCloth(); + InitScene(); +}; + +void ClothSimulation::Update(float dt) { + if (IsVisible() == false || IsPause() == true) return; + + _solver->solve(_iter_num); + _solver->solve(_iter_num); + + _solver->satisfy(); + + // mCollisionDetectResponse->Update(dt, ); + // _cloth_mesh->SetPoints(_solver->getCurrentState()); + float* pointData = _solver->getCurrentState().data(); + for (dtk::dtkID id = 0; id < _cloth_mesh->GetNumberOfPoints(); id++) { + _cloth_mesh->SetPoint(id, dtk::GK::Point3(pointData[3 * id + 0], pointData[3 * id + 1], pointData[3 * id + 2])); + } + _cloth_mesh->update_normals(); + + UpdateRenderTarget(); +}; + +void ClothSimulation::Render() { + if (!IsVisible()) return; + Renderer renderer; + renderer.setProgram(g_phongShader); + renderer.setModelview(g_ModelViewMatrix); + renderer.setProjection(g_ProjectionMatrix); + g_phongShader->setAlbedo(g_albedo); + g_phongShader->setAmbient(g_ambient); + g_phongShader->setLight(g_light); + renderer.setProgramInput(g_render_target); + renderer.setElementCount(_cloth_mesh->GetNumberOfTriangles() * 3); + renderer.draw(); +}; + +void ClothSimulation::InitShader() { + GLShader basic_vert(GL_VERTEX_SHADER); + GLShader phong_frag(GL_FRAGMENT_SHADER); + GLShader pick_frag(GL_FRAGMENT_SHADER); + + auto ibasic = std::ifstream("./test/system_test/MassSpring3D/shaders/basic.vshader"); + auto iphong = std::ifstream("./test/system_test/MassSpring3D/shaders/phong.fshader"); + auto ifrag = std::ifstream("./test/system_test/MassSpring3D/shaders/pick.fshader"); + + basic_vert.compile(ibasic); + phong_frag.compile(iphong); + pick_frag.compile(ifrag); + + g_phongShader = new PhongShader; + g_pickShader = new PickShader; + g_phongShader->link(basic_vert, phong_frag); + g_pickShader->link(basic_vert, pick_frag); +} + +void ClothSimulation::InitCloth() { + _cloth_mesh = dtkFactory::CreateClothMesh(SystemParam::w, SystemParam::n); + + ClothDrop(); + + g_render_target = new ProgramInput; + UpdateRenderTarget(); // set position data + const std::vector& mEc = _cloth_mesh->GetECTable(); + unsigned int* indexBuffer = (unsigned int*)&mEc[0]; + unsigned int indexBufferSize = mEc.size() * 3; + g_render_target->setIndexData(indexBuffer, indexBufferSize); +} + +void ClothSimulation::InitScene() { + g_ModelViewMatrix = glm::lookAt( + glm::vec3(0.618, -0.786, 0.3f) * g_camera_distance, + glm::vec3(0.0f, 0.0f, -1.0f), + glm::vec3(0.0f, 0.0f, 1.0f) + ) * glm::translate(glm::mat4(1), glm::vec3(0.0f, 0.0f, SystemParam::w / 4)); + g_ProjectionMatrix = glm::perspective(PI / 4.0f, g_windowWidth * 1.0f / g_windowHeight, 0.01f, 1000.0f); +}; + +void ClothSimulation::SetParameters() { + +}; + +void ClothSimulation::ClothDrop() { + _system = dtkFactory::CreateClothMassSpringSystem(_cloth_mesh); + // _sphere_mesh = dtkFactory::CreateSphereMesh(dtk::dtkDouble3(0, 0, -1), 0.64, 20); + + _solver = dtkFactory::CreateClothMassSpringSolver(_system); + + mCollisionDetectResponse = dtk::dtkPhysMassSpringCollisionResponse::New(); + mCollisionDetectResponse->SetMassSpring(0, _system); +}; + +void ClothSimulation::UpdateClothMesh() { + const VectorXf& current_state = _solver->getCurrentState(); + for (dtk::dtkID id = 0; id < _cloth_mesh->GetNumberOfPoints(); id++) { + dtk::dtkDouble3 p(current_state[3 * id + 0], current_state[3 * id + 1], current_state[3 * id + 2]); + _cloth_mesh->SetPoint(id, dtk::GK::Point3(p.x, p.y, p.z)); + } + _cloth_mesh->update_normals(); +} + +void ClothSimulation::UpdateRenderTarget() { + dtk::dtkPoints::Ptr mPts = _cloth_mesh->GetPoints(); + + unsigned int vertexBufferSize = mPts->GetNumberOfPoints() * 3; + float* vertexBuffer = new float[vertexBufferSize]; + const std::vector& normalData = _cloth_mesh->GetVertexNormals(); + float* normalBuffer = new float[vertexBufferSize]; + for (dtk::dtkID i = 0;i < mPts->GetNumberOfPoints();i++) { + // vertexBuffer[3 * i + 0] = (float)mPts->GetPoint(i)[0]; + // vertexBuffer[3 * i + 1] = (float)mPts->GetPoint(i)[1]; + // vertexBuffer[3 * i + 2] = (float)mPts->GetPoint(i)[2]; + normalBuffer[3 * i + 0] = (float)normalData[i].x; + normalBuffer[3 * i + 1] = (float)normalData[i].y; + normalBuffer[3 * i + 2] = (float)normalData[i].z; + } + + // std::cout << "vertexBuffer:" << std::endl; + // for (int i = 0;i < vertexBufferSize;i++) { + // std::cout << vertexBuffer[i] << " "; + // } + // std::cout << std::endl; + + g_render_target->setPositionData(_solver->getCurrentState().data(), vertexBufferSize); + g_render_target->setNormalData(normalBuffer, vertexBufferSize); +}; + +dtk::dtkStaticTriangleMesh::Ptr dtkFactory::CreateClothMesh(float w, int n) { + dtk::dtkStaticTriangleMesh::Ptr result = dtk::dtkStaticTriangleMesh::New(); + + unsigned int idx = 0; // vertex index + const float d = w / (n - 1); // step distance + const float ud = 1.0f / (n - 1); // unit step distance + const glm::vec3 o = glm::vec3(-w / 2.0f, w / 2.0f, 0.0f); // origin + const glm::vec3 ux = glm::vec3(1.0f, 0.0f, 0.0f); // unit x direction + const glm::vec3 uy = glm::vec3(0.0f, -1.0f, 0.0f); // unit y direction + + dtk::dtkPointsVector::Ptr vertices = dtk::dtkPointsVector::New(); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + glm::vec3 p = o + d * j * ux + d * i * uy; // vertex position + dtk::GK::Point3 p3(p.x, p.y, p.z); + dtk::dtkID pid = j + i * n; + vertices->InsertPoint(pid, p3); // add vertex + } + } + result->SetPoints(vertices); + + for (int i = 0;i < n;i++) { + for (int j = 0;j < n;j++) { + //add connectivity + if (i > 0 && j < n - 1) { + result->InsertTriangleNotRepeat( + j + i * n, + j + 1 + (i - 1) * n, + j + (i - 1) * n + ); + } + + if (j > 0 && i > 0) { + result->InsertTriangleNotRepeat( + j + i * n, + j + (i - 1) * n, + j - 1 + i * n + ); + } + } + } + + result->update_normals(); + return result; +} + +dtk::dtkStaticTriangleMesh::Ptr dtkFactory::CreateSphereMesh(dtk::dtkDouble3 center, float radius, int n) { + dtk::dtkStaticTriangleMesh::Ptr result = dtk::dtkStaticTriangleMesh::New(); + dtk::dtkPointsVector::Ptr vertices = dtk::dtkPointsVector::New(); + + // Generate vertices + for (int i = 0; i <= n; ++i) { + float theta = i * glm::pi() / n; + for (int j = 0; j <= n; ++j) { + float phi = j * 2 * glm::pi() / n; + float x = center.x + radius * sin(theta) * cos(phi); + float y = center.y + radius * sin(theta) * sin(phi); + float z = center.z + radius * cos(theta); + dtk::GK::Point3 p3(x, y, z); + dtk::dtkID pid = j + i * (n + 1); + vertices->InsertPoint(pid, p3); + } + } + result->SetPoints(vertices); + + // Generate triangles + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + dtk::dtkID p1 = j + i * (n + 1); + dtk::dtkID p2 = j + (i + 1) * (n + 1); + dtk::dtkID p3 = (j + 1) + i * (n + 1); + dtk::dtkID p4 = (j + 1) + (i + 1) * (n + 1); + + if (i != 0) { + result->InsertTriangleNotRepeat(p1, p2, p3); + } + if (i != (n - 1)) { + result->InsertTriangleNotRepeat(p3, p2, p4); + } + } + } + + return result; +} + +dtk::dtkPhysMassSpring::Ptr dtkFactory::CreateClothMassSpringSystem(const dtk::dtkStaticTriangleMesh::Ptr& mesh) { + dtk::dtkDouble3 gravity(0, 0, -SystemParam::g); + dtk::dtkPhysMassSpring::Ptr system = dtk::dtkPhysMassSpring::New(SystemParam::m, SystemParam::k, SystemParam::b, SystemParam::a, SystemParam::r, SystemParam::h, gravity); + // system->SetTriangleMesh(mesh); + + // n must be odd + assert(SystemParam::n % 2 == 1); + + // compute n_points and n_springs + unsigned int n_points = SystemParam::n * SystemParam::n; + unsigned int n_springs = (SystemParam::n - 1) * (5 * SystemParam::n - 2); + system->SetPoints(mesh->GetPoints()); + for (unsigned int id = 0; id < n_points; id++) { + system->AddMassPoint(id, SystemParam::m, dtk::dtkT3(0, 0, 0), SystemParam::a, SystemParam::c, gravity); + } + + unsigned int n = SystemParam::n; + double rest_length_factor = 1.05; + for (unsigned int i = 0; i < SystemParam::n; i++) { + for (unsigned int j = 0; j < SystemParam::n; j++) { + // bottom right corner + if (i == n - 1 && j == n - 1) { + continue; + } + + if (i == n - 1) { + // structural spring + system->AddSpring(n * i + j, n * i + j + 1, SystemParam::k, SystemParam::a, rest_length_factor); + + // bending spring + if (j % 2 == 0) { + system->AddSpring(n * i + j, n * i + j + 2, SystemParam::k, SystemParam::a, rest_length_factor); + } + continue; + } + + // right edge + if (j == n - 1) { + // structural spring + system->AddSpring(n * i + j, n * (i + 1) + j, SystemParam::k, SystemParam::a, rest_length_factor); + + // bending spring + if (i % 2 == 0) { + system->AddSpring(n * i + j, n * (i + 2) + j, SystemParam::k, SystemParam::a, rest_length_factor); + } + continue; + } + + // structural springs + system->AddSpring(n * i + j, n * i + j + 1, SystemParam::k, SystemParam::a, rest_length_factor); + + system->AddSpring(n * i + j, n * (i + 1) + j, SystemParam::k, SystemParam::a, rest_length_factor); + + // shearing springs + system->AddSpring(n * i + j, n * (i + 1) + j + 1, SystemParam::k, SystemParam::a, rest_length_factor); + + system->AddSpring(n * (i + 1) + j, n * i + j + 1, SystemParam::k, SystemParam::a, rest_length_factor); + + // bending springs + if (j % 2 == 0) { + system->AddSpring(n * i + j, n * i + j + 2, SystemParam::k, SystemParam::a, rest_length_factor); + } + if (i % 2 == 0) { + system->AddSpring(n * i + j, n * (i + 2) + j, SystemParam::k, SystemParam::a, rest_length_factor); + } + } + } + + // TODO: setting gravity as an external force + + return system; +} + +dtk::dtkPhysMassSpringSolver::Ptr dtkFactory::CreateClothMassSpringSolver(const dtk::dtkPhysMassSpring::Ptr& system) { + dtk::dtkPhysMassSpringSolver::Ptr solver = dtk::dtkPhysMassSpringSolver::New(system); + return solver; +} + diff --git a/test/system_test/MassSpring3D/Renderer.cpp b/test/system_test/MassSpring3D/Renderer.cpp new file mode 100644 index 0000000..e8f0f46 --- /dev/null +++ b/test/system_test/MassSpring3D/Renderer.cpp @@ -0,0 +1,49 @@ +#include "Renderer.h" +#include + +Renderer::Renderer() {} + +void Renderer::setProgram(GLProgram* program) { + assert(program != nullptr); + assert((*program) > 0); + this->program = program; +} + +void Renderer::setProgramInput(ProgramInput* input) { + this->input = input; +} + +void Renderer::setModelview(const glm::mat4& mv) { + assert(program != nullptr); + assert((*program) > 0); + glUseProgram(*program); + program->setModelView(mv); + glUseProgram(0); +} + +void Renderer::setProjection(const glm::mat4& p) { + assert(program != nullptr); + assert((*program) > 0); + glUseProgram(*program); + program->setProjection(p); + glUseProgram(0); +} + +void Renderer::setElementCount(unsigned int n_elements) { this->n_elements = n_elements; } + +void Renderer::draw() { + assert(program != nullptr); + assert((*program) > 0); + glUseProgram(*program); + glBindVertexArray(*input); + glEnableVertexAttribArray(0); + glEnableVertexAttribArray(1); + glEnableVertexAttribArray(2); + glDrawElements(GL_TRIANGLES, n_elements, GL_UNSIGNED_INT, 0); + glDisableVertexAttribArray(0); + glDisableVertexAttribArray(1); + glDisableVertexAttribArray(2); + glBindVertexArray(0); + glUseProgram(0); +} + diff --git a/test/system_test/MassSpring3D/Scene.cpp b/test/system_test/MassSpring3D/Scene.cpp new file mode 100644 index 0000000..e69de29 diff --git a/test/system_test/MassSpring3D/Shader.cpp b/test/system_test/MassSpring3D/Shader.cpp new file mode 100644 index 0000000..234fc65 --- /dev/null +++ b/test/system_test/MassSpring3D/Shader.cpp @@ -0,0 +1,195 @@ +#include "Shader.h" +#include +#include +#include + +// GLSHADER /////////////////////////////////////////////////////////////////////////////////// +GLShader::GLShader(GLenum shaderType) { + // if (!glewIsSupported("GL_VERSION_2_0")) { + // printf("OpenGL 2.0 not supported\n"); + // } + handle = glCreateShader(shaderType); +}; + +GLShader::operator GLuint() const { + return handle; +} + +GLShader::~GLShader() { + glDeleteShader(handle); +} + +void GLShader::compile(const char* source) { + GLint compiled = 0; // Compiled flag + const char* ptrs[] = { source }; + const GLint lens[] = { static_cast(std::strlen(source)) }; + glShaderSource(handle, 1, ptrs, lens); + glCompileShader(handle); + glGetShaderiv(handle, GL_COMPILE_STATUS, &compiled); + if (!compiled) { + GLint logSize = 0; + glGetShaderiv(handle, GL_INFO_LOG_LENGTH, &logSize); + std::vector errorLog(logSize); + glGetShaderInfoLog(handle, logSize, &logSize, &errorLog[0]); + std::cerr << &errorLog[0] << std::endl; + throw std::runtime_error("Failed to compile shader."); + } +} + +void GLShader::compile(std::ifstream& source) { + std::vector text; + source.seekg(0, std::ios_base::end); + std::streampos fileSize = source.tellg(); + text.resize(fileSize); + + source.seekg(0, std::ios_base::beg); + source.read(&text[0], fileSize); + compile(&text[0]); +} + +// GLPROGRAM ////////////////////////////////////////////////////////////////////////////////// +GLProgram::GLProgram() : handle(glCreateProgram()) {} + +void GLProgram::link(const GLShader& vshader, const GLShader& fshader) { + GLint linked = 0; // Linked flag + glAttachShader(handle, vshader); + glAttachShader(handle, fshader); + glLinkProgram(handle); + glDetachShader(handle, vshader); + glDetachShader(handle, fshader); + glGetProgramiv(handle, GL_LINK_STATUS, &linked); + if (!linked) + throw std::runtime_error("Failed to link shaders."); + + // get camera uniforms + uModelViewMatrix = glGetUniformLocation(handle, "uModelViewMatrix"); + uProjectionMatrix = glGetUniformLocation(handle, "uProjectionMatrix"); + + // post link + postLink(); + +} + +void GLProgram::postLink() {} + +void GLProgram::setUniformMat4(GLuint unif, glm::mat4 m) { + glUseProgram(handle); + glUniformMatrix4fv(unif, + 1, GL_FALSE, glm::value_ptr(m[0])); + glUseProgram(0); +} + + +void GLProgram::setModelView(glm::mat4 m) { + setUniformMat4(uModelViewMatrix, m); +} + +void GLProgram::setProjection(glm::mat4 m) { + setUniformMat4(uProjectionMatrix, m); +} + +GLProgram::operator GLuint() const { return handle; } + +GLProgram::~GLProgram() { glDeleteProgram(handle); } + +// PROGRAM INPUT ////////////////////////////////////////////////////////////////////////////// + +ProgramInput::ProgramInput() { + // generate buffers + glGenBuffers(4, &vbo[0]); + + // generate vertex array object + glGenVertexArrays(1, &handle); + glBindVertexArray(handle); + + // positions + glBindBuffer(GL_ARRAY_BUFFER, vbo[0]); + glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, 0); + + // normals + glBindBuffer(GL_ARRAY_BUFFER, vbo[1]); + glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, 0); + + // texture coordinates + glBindBuffer(GL_ARRAY_BUFFER, vbo[2]); + glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, 0, 0); + + // indices + glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbo[3]); + + glBindVertexArray(0); + +} + +void ProgramInput::bufferData(unsigned int index, void* buff, size_t size) { + glBindBuffer(GL_ARRAY_BUFFER, vbo[index]); + glBufferData(GL_ARRAY_BUFFER, size, + buff, GL_STATIC_DRAW); +} +void ProgramInput::setPositionData(float* buff, unsigned int len) { + bufferData(0, buff, sizeof(float) * len); +} + +void ProgramInput::setNormalData(float* buff, unsigned int len) { + bufferData(1, buff, sizeof(float) * len); +} + +void ProgramInput::setTextureData(float* buff, unsigned int len) { + bufferData(2, buff, sizeof(float) * len); +} + +void ProgramInput::setIndexData(unsigned int* buff, unsigned int len) { + bufferData(3, buff, sizeof(unsigned int) * len); +} + +ProgramInput::operator GLuint() const { + return handle; +} + + +ProgramInput::~ProgramInput() { + glDeleteBuffers(4, vbo); + glDeleteVertexArrays(1, &handle); +} + +// SHADER PROGRAMS //////////////////////////////////////////////////////////////////////////// +PhongShader::PhongShader() : GLProgram() {} +void PhongShader::postLink() { + // get uniforms + uAlbedo = glGetUniformLocation(handle, "uAlbedo"); + uAmbient = glGetUniformLocation(handle, "uAmbient"); + uLight = glGetUniformLocation(handle, "uLight"); +} +void PhongShader::setAlbedo(const glm::vec3& albedo) { + assert(uAlbedo >= 0); + glUseProgram(*this); + glUniform3f(uAlbedo, albedo[0], albedo[1], albedo[2]); + glUseProgram(0); +} +void PhongShader::setAmbient(const glm::vec3& ambient) { + assert(uAmbient >= 0); + glUseProgram(*this); + glUniform3f(uAmbient, ambient[0], ambient[1], ambient[2]); + glUseProgram(0); +} +void PhongShader::setLight(const glm::vec3& light) { + assert(uLight >= 0); + glUseProgram(*this); + glUniform3f(uLight, light[0], light[1], light[2]); + glUseProgram(0); +} + + +PickShader::PickShader() : GLProgram() {} + +void PickShader::postLink() { + // get uniforms + uTessFact = glGetUniformLocation(handle, "uTessFact"); +} + +void PickShader::setTessFact(unsigned int n) { + assert(uTessFact >= 0); + glUseProgram(*this); + glUniform1i(uTessFact, n); + glUseProgram(0); +} \ No newline at end of file diff --git a/test/system_test/MassSpring3D/dtkPhysMassSpringSolver.cpp b/test/system_test/MassSpring3D/dtkPhysMassSpringSolver.cpp new file mode 100644 index 0000000..6b488d6 --- /dev/null +++ b/test/system_test/MassSpring3D/dtkPhysMassSpringSolver.cpp @@ -0,0 +1,179 @@ +#include "dtkPhysMassSpringSolver.h" + +dtk::dtkPhysMassSpringSolver::dtkPhysMassSpringSolver() { +} + +dtk::dtkPhysMassSpringSolver::dtkPhysMassSpringSolver(const dtk::dtkPhysMassSpring::Ptr& massSpringSystem) { + _system = massSpringSystem; + _time_step = _system->GetTimeStep(); + + // current state + _current_state.resize(3 * _system->GetNumberOfMassPoints()); + for (dtk::dtkID i = 0; i < _system->GetNumberOfMassPoints(); i++) { + dtk::dtkPhysMassPoint* massPoint = _system->GetMassPoint(i); + _current_state[3 * i] = massPoint->GetPosition()[0]; + _current_state[3 * i + 1] = massPoint->GetPosition()[1]; + _current_state[3 * i + 2] = massPoint->GetPosition()[2]; + } + _prev_state = _current_state; + _initial_state = _current_state; + _spring_directions.resize(3 * _system->GetNumberOfSprings()); + + // M + TripletList MTriplets; + _M.resize(3 * _system->GetNumberOfMassPoints(), 3 * _system->GetNumberOfMassPoints()); + for (dtk::dtkID i = 0; i < _system->GetNumberOfMassPoints(); i++) { + for (dtk::dtkID j = 0; j < 3; j++) { + MTriplets.push_back(Triplet(3 * i + j, 3 * i + j, _system->GetMassPoint(i)->GetMass())); + } + } + _M.setFromTriplets(MTriplets.begin(), MTriplets.end()); + + // L + TripletList LTriplets; + _L.resize(3 * _system->GetNumberOfMassPoints(), 3 * _system->GetNumberOfMassPoints()); + for (dtk::dtkID i = 0;i < _system->GetNumberOfSprings();i++) { + double stiffness = _system->GetSpring(i)->GetStiffness(); + dtk::dtkPhysSpring* spring = _system->GetSpring(i); + dtk::dtkID id1 = spring->GetFirstVertex()->GetPointID();// TODO GetFirstVertex修改成GetFirstPoint(); + dtk::dtkID id2 = spring->GetSecondVertex()->GetPointID(); + for (dtk::dtkID j = 0; j < 3; j++) { + LTriplets.push_back(Triplet(3 * id1 + j, 3 * id1 + j, 1 * stiffness)); + LTriplets.push_back(Triplet(3 * id1 + j, 3 * id2 + j, -1 * stiffness)); + LTriplets.push_back(Triplet(3 * id2 + j, 3 * id1 + j, -1 * stiffness)); + LTriplets.push_back(Triplet(3 * id2 + j, 3 * id2 + j, 1 * stiffness)); + } + } + _L.setFromTriplets(LTriplets.begin(), LTriplets.end()); + + // J + TripletList JTriplets; + _J.resize(3 * _system->GetNumberOfMassPoints(), 3 * _system->GetNumberOfSprings()); + for (dtk::dtkID i = 0;i < _system->GetNumberOfSprings();i++) { + double stiffness = _system->GetSpring(i)->GetStiffness(); + dtk::dtkPhysSpring* spring = _system->GetSpring(i); + dtk::dtkID id1 = spring->GetFirstVertex()->GetPointID(); + dtk::dtkID id2 = spring->GetSecondVertex()->GetPointID(); + for (unsigned int j = 0; j < 3; j++) { + JTriplets.push_back( + Triplet(3 * id1 + j, 3 * i + j, 1 * stiffness)); + JTriplets.push_back( + Triplet(3 * id2 + j, 3 * i + j, -1 * stiffness)); + } + } + _J.setFromTriplets(JTriplets.begin(), JTriplets.end()); + + // pre-factor + double h2 = _system->GetTimeStep() * _system->GetTimeStep(); + SparseMatrix A = _M + h2 * _L; + _system_matrix.compute(A); + + + // std::cout << "M: " << std::endl; + // printSparseMatrix(_M); + // std::cout << "L: " << std::endl; + // printSparseMatrix(_L); + // std::cout << "J: " << std::endl; + // printSparseMatrix(_J); + // std::cout << "A: " << std::endl; + // printSparseMatrix(A); +} + +void dtk::dtkPhysMassSpringSolver::solve(unsigned int iter_num) { + float damping_factor = _system->GetDefaultPointDamp(); + + // update inertial term + _inertial_term = _M * ((damping_factor + 1) * (_current_state)-damping_factor * _prev_state); + _prev_state = _current_state; + + // std::cout << "damping_factor: " << damping_factor << std::endl; + // std::cout << "M: " << std::endl; + // printSparseMatrix(_M); + // std::cout << "_current_state: " << std::endl << _current_state << std::endl; + // std::cout << "_prev_state: " << std::endl << _prev_state << std::endl; + // std::cout << "_internal_term: " << std::endl << _inertial_term << std::endl; + + // perform steps + for (unsigned int i = 0; i < iter_num; i++) { + localStep(); + globalStep(); + } + + // std::cout << "_current_state: " << std::endl << _current_state << std::endl; +} + +void dtk::dtkPhysMassSpringSolver::localStep() { + for (dtk::dtkID id = 0;id < _system->GetNumberOfSprings();id++) { + dtk::dtkPhysSpring* spring = _system->GetSpring(id); + dtk::dtkID id1 = spring->GetFirstVertex()->GetPointID(); + dtk::dtkID id2 = spring->GetSecondVertex()->GetPointID(); + double rest_length = spring->GetRestLength(); + Vector3f p12( + _current_state[3 * id1 + 0] - _current_state[3 * id2 + 0], + _current_state[3 * id1 + 1] - _current_state[3 * id2 + 1], + _current_state[3 * id1 + 2] - _current_state[3 * id2 + 2] + ); + + p12.normalize(); + _spring_directions[3 * id + 0] = rest_length * p12[0]; + _spring_directions[3 * id + 1] = rest_length * p12[1]; + _spring_directions[3 * id + 2] = rest_length * p12[2]; + } + // std::cout << "_spring_directions: " << std::endl << _spring_directions << std::endl; +} + +void dtk::dtkPhysMassSpringSolver::globalStep() { + float h2 = _system->GetTimeStep() * _system->GetTimeStep(); // shorthand + + // compute right hand side + // dtk::dtkDouble3 fext(0, 0, 0); + dtk::dtkDouble3 fext = _system->GetDefaultGravityAccel(); // TODO: change to _system->GetImpulseForce() + VectorXf fext_force = VectorXf(Vector3f(fext.x, fext.y, fext.z).replicate(_system->GetNumberOfMassPoints(), 1)); + + // // 生成随机外力 + // std::random_device rd; + // std::mt19937 gen(rd()); + // std::uniform_real_distribution<> dis(-1.0, 1.0); + // dtk::dtkDouble3 random_force(dis(gen), dis(gen), dis(gen)); + // // 随机选择一个质点 + // std::uniform_int_distribution<> point_dis(0, _system->GetNumberOfMassPoints() - 1); + // int random_point = point_dis(gen); + // fext_force.segment<3>(random_point * 3) += Vector3f(random_force.x, random_force.y, random_force.z); + + VectorXf b = _inertial_term + h2 * _J * _spring_directions + h2 * fext_force; + // std::cout << "b: " << std::endl << b << std::endl; + + // solve system and update state + _current_state = _system_matrix.solve(b); + // std::cout << "_current_state: " << std::endl << _current_state << std::endl; +} + +void dtk::dtkPhysMassSpringSolver::satisfy(ClothDropType type) { + if (type == Sphere) { + const float radius = 0.64f; + const Eigen::Vector3f center(0, 0, -1); + + for (int i = 0; i < _system->GetNumberOfMassPoints(); i++) { + Vector3f p( + _current_state[3 * i + 0] - center[0], + _current_state[3 * i + 1] - center[1], + _current_state[3 * i + 2] - center[2] + ); + + if (p.norm() < radius) { + p.normalize(); + p = radius * p; + } + else continue; + + for (int j = 0; j < 3; j++) { + _current_state[3 * i + j] = p[j] + center[j]; + } + } + } + else { + int fix_idx = 0; + for (int i = 0; i < 3; i++) + _current_state[fix_idx + i] = _initial_state[fix_idx + i]; + } +} \ No newline at end of file diff --git a/test/system_test/MassSpring3D/include/ClothSimulation.h b/test/system_test/MassSpring3D/include/ClothSimulation.h new file mode 100644 index 0000000..17b2c1d --- /dev/null +++ b/test/system_test/MassSpring3D/include/ClothSimulation.h @@ -0,0 +1,123 @@ +/** + * @file ClothSimulation.h + * @brief ClothSimulation 头文件 + * @author Chance + * @version 1.0 + * @date 2024-07-12 + * @copyright MIT LICENSE + * https://github.com/Simple-XX/SimplePhysicsEngine + * @par change log: + * + *
DateAuthorDescription + *
2024-07-12Chance创建文件 + *
+ */ +#ifndef SIMPLEPHYSICSENGINE_CLOTHSIMULATION_H +#define SIMPLEPHYSICSENGINE_CLOTHSIMULATION_H + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "dtkStaticTriangleMesh.h" +#include "dtkPhysMassSpring.h" +#include "dtkMatrix.h" + +#include "Scene.h" +#include "Shader.h" +#include "Renderer.h" +#include "dtkCollisionDetectHierarchyKDOPS.h" +#include "dtkPhysCore.h" +#include "dtkPhysMassSpringSolver.h" +#include "dtkPhysMassSpringCollisionResponse.h" + +namespace SystemParam { + static const int n = 33; // must be odd, n * n = n_vertices | 61 + static const float w = 2.0f; // width | 2.0f + static const float h = 0.008f; // time step, smaller for better results | 0.008f = 0.016f/2 + static const float r = w / (n - 1) * 1.05f; // spring rest legnth + static const float k = 1.0f; // spring stiffness | 1.0f; + static const float m = 0.05f / (n * n); // point mass | 0.25f + static const float a = 0.993f; // point damping, close to 1.0 | 0.993f + static const float b = 5880.0f; // damping | 5880.0f + static const float c = 2.5f; // point resistence | 2.5f + static const float g = 9.8f * m; // gravitational force | 9.8f +} + +class ClothSimulation : public Scene { + typedef Eigen::VectorXf VectorXf; +public: + ClothSimulation(const unsigned int& windowWidth, const unsigned int& windowHeight, const dtk::dtkDouble3& gravity); + ~ClothSimulation() = default; + + using ClothMesh = dtk::dtkStaticTriangleMesh::Ptr; + using SphereMesh = dtk::dtkStaticTriangleMesh::Ptr; + using ClothMassSpring = dtk::dtkPhysMassSpring::Ptr; + using ClothMassSpringSolver = dtk::dtkPhysMassSpringSolver::Ptr; + + const ClothMesh GetClothMesh() const; + + void move(const dtk::dtkDouble2& v); + + void CleanUp(); + + void Init(); + void Update(float dt); + void Render(); + void UpdateRenderTarget(); + +private: + void InitShader(); + void InitCloth(); + void InitScene(); + void SetParameters(); + void ClothDrop(); + + void UpdateClothMesh(); + + // Shader + PhongShader* g_phongShader; // linked phong shader + PickShader* g_pickShader; // linked pick shader + + // Camera + dtk::dtkMatrix44 g_ModelViewMatrix; + dtk::dtkMatrix44 g_ProjectionMatrix; + const float g_camera_distance = 4.2f; + const float PI = glm::pi(); + + ProgramInput* g_render_target; // vertex, index + + // @todo: change to dtk vector + const glm::vec3 g_albedo = glm::vec3(0.0f, 0.2f, 0.9f); + const glm::vec3 g_ambient = glm::vec3(0.01f, 0.01f, 0.01f); + const glm::vec3 g_light = glm::vec3(1.0f, 1.0f, -1.0f); + + // 重力 + dtk::dtkDouble3 _gravity; + + ClothMesh _cloth_mesh; + SphereMesh _sphere_mesh; + ClothMassSpring _system; + ClothMassSpringSolver _solver; + + dtk::dtkPhysMassSpringCollisionResponse::Ptr mCollisionDetectResponse; + + const int _iter_num = 5; +}; + +class dtkFactory { +public: + static dtk::dtkStaticTriangleMesh::Ptr CreateClothMesh(float w, int n); + static dtk::dtkPhysMassSpring::Ptr CreateClothMassSpringSystem(const dtk::dtkStaticTriangleMesh::Ptr& mesh); + static dtk::dtkPhysMassSpringSolver::Ptr CreateClothMassSpringSolver(const dtk::dtkPhysMassSpring::Ptr& system); + static dtk::dtkStaticTriangleMesh::Ptr CreateSphereMesh(dtk::dtkDouble3 center, float radius, int n); +}; + +#endif /* SIMPLEPHYSICSENGINE_CLOTHSIMULATION_H */ diff --git a/test/system_test/MassSpring3D/include/Renderer.h b/test/system_test/MassSpring3D/include/Renderer.h new file mode 100644 index 0000000..525184f --- /dev/null +++ b/test/system_test/MassSpring3D/include/Renderer.h @@ -0,0 +1,22 @@ +#pragma once +#include +#include "Shader.h" + + +class Renderer { +protected: + GLProgram* program; + ProgramInput* input; + unsigned int n_elements; + +public: + Renderer(); + + void setProgram(GLProgram* program); + void setProgramInput(ProgramInput* input); + void setModelview(const glm::mat4& mv); + void setProjection(const glm::mat4& p); + void setElementCount(unsigned int n_elements); + + void draw(); +}; \ No newline at end of file diff --git a/test/system_test/MassSpring3D/include/Scene.h b/test/system_test/MassSpring3D/include/Scene.h new file mode 100644 index 0000000..dbc7a14 --- /dev/null +++ b/test/system_test/MassSpring3D/include/Scene.h @@ -0,0 +1,45 @@ +#include +#include + +enum SceneState { SCENE_ACTIVE, SCENE_PAUSE }; +enum SceneVisibility { SCENE_VISIBLE, SCENE_HIDDEN }; + +class Scene { + SceneState State; + SceneVisibility Visibility; +public: + Scene() : State(SCENE_ACTIVE), Visibility(SCENE_VISIBLE) {}; + Scene(const unsigned int& width, const unsigned int& height) : State(SCENE_PAUSE), Visibility(SCENE_HIDDEN), g_windowWidth(width), g_windowHeight(height) {}; + virtual void Init() = 0; + virtual void Update(float dt) = 0; + virtual void Render() = 0; + + bool IsPause() const { + return State == SCENE_PAUSE; + } + + void SetPause(bool pause) { + if (pause) + State = SCENE_PAUSE; + else + State = SCENE_ACTIVE; + } + + bool IsVisible() const { + return Visibility == SCENE_VISIBLE; + } + + void SetVisible(bool visible) { + if (visible) { + Visibility = SCENE_VISIBLE; + State = SCENE_ACTIVE; + } + else { + Visibility = SCENE_HIDDEN; + State = SCENE_PAUSE; + } + } + + unsigned int g_windowWidth = 800; + unsigned int g_windowHeight = 600; +}; \ No newline at end of file diff --git a/test/system_test/MassSpring3D/include/Shader.h b/test/system_test/MassSpring3D/include/Shader.h new file mode 100644 index 0000000..652eb7f --- /dev/null +++ b/test/system_test/MassSpring3D/include/Shader.h @@ -0,0 +1,87 @@ +#pragma once +#include +#include +#include + +class NonCopyable { +private: + NonCopyable(const NonCopyable& other) = delete; + NonCopyable& operator=(const NonCopyable& other) = delete; + +public: + NonCopyable() {} +}; + +class GLShader : public NonCopyable { +private: + GLuint handle; // Shader handle + +public: + GLShader(GLenum shaderType); + /*GLShader(GLenum shaderType, const char* source); + GLShader(GLenum shaderType, std::ifstream& source);*/ + void compile(const char* source); + void compile(std::ifstream& source); + operator GLuint() const; // cast to GLuint + ~GLShader(); +}; + +class GLProgram : public NonCopyable { +protected: + GLuint handle; + GLuint uModelViewMatrix, uProjectionMatrix; + void setUniformMat4(GLuint unif, glm::mat4 m); + +public: + GLProgram(); + virtual void link(const GLShader& vshader, const GLShader& fshader); + virtual void postLink(); + operator GLuint() const; // cast to GLuint + void setModelView(glm::mat4 m); + void setProjection(glm::mat4 m); + + ~GLProgram(); +}; + +class ProgramInput : public NonCopyable { +private: + GLuint handle; // vertex array object handle + GLuint vbo[4]; // vertex buffer object handles | position, normal, texture, index + void bufferData(unsigned int index, void* buff, size_t size); + +public: + ProgramInput(); + + void setPositionData(float* buff, unsigned int len); + void setNormalData(float* buff, unsigned int len); + void setTextureData(float* buff, unsigned int len); + void setIndexData(unsigned int* buff, unsigned int len); + + operator GLuint() const; // cast to GLuint + + ~ProgramInput(); +}; + +class PhongShader : public GLProgram { +private: + // Albedo | Ambient Light | Light Direction + GLuint uAlbedo, uAmbient, uLight; + +public: + PhongShader(); + virtual void postLink(); + void setAlbedo(const glm::vec3& albedo); + void setAmbient(const glm::vec3& ambient); + void setLight(const glm::vec3& light); +}; + + +class PickShader : public GLProgram { +private: + GLuint uTessFact; + +public: + PickShader(); + virtual void postLink(); + void setTessFact(unsigned int n); +}; \ No newline at end of file diff --git a/test/system_test/MassSpring3D/include/dtkPhysMassSpringSolver.h b/test/system_test/MassSpring3D/include/dtkPhysMassSpringSolver.h new file mode 100644 index 0000000..1e30e73 --- /dev/null +++ b/test/system_test/MassSpring3D/include/dtkPhysMassSpringSolver.h @@ -0,0 +1,72 @@ +#ifndef SIMPLEPHYSICSENGINE_DTKPHYSMASSSPRINGSOLVER_H +#define SIMPLEPHYSICSENGINE_DTKPHYSMASSSPRINGSOLVER_H +#include "dtkPhysMassPoint.h" +#include "dtkPhysMassSpring.h" + +#include +#include +#include + +namespace dtk { + class dtkPhysMassSpringSolver { + private: + typedef Eigen::Vector3f Vector3f; + typedef Eigen::VectorXf VectorXf; + typedef Eigen::SparseMatrix SparseMatrix; + typedef Eigen::SimplicialLLT > Cholesky; + typedef Eigen::Map Map; + typedef std::pair Edge; + typedef Eigen::Triplet Triplet; + typedef std::vector TripletList; + public: + typedef std::shared_ptr Ptr; + + static dtkPhysMassSpringSolver::Ptr New() { + return Ptr(new dtkPhysMassSpringSolver()); + } + static dtkPhysMassSpringSolver::Ptr New(const dtkPhysMassSpring::Ptr& massSpring) { + return Ptr(new dtkPhysMassSpringSolver(massSpring)); + } + + void solve(unsigned int iter_num); + enum ClothDropType { + Sphere, + Hang + }; + + + void satisfy(ClothDropType type = Sphere); + VectorXf getCurrentState() const { return _current_state; }; + + static void printSparseMatrix(const SparseMatrix& matrix) { + for (int k = 0; k < matrix.outerSize(); ++k) { + for (SparseMatrix::InnerIterator it(matrix, k); it; ++it) { + std::cout << "Element at (" << it.row() << ", " << it.col() << ") = " << it.value() << std::endl; + } + } + } + private: + dtkPhysMassSpringSolver(); + dtkPhysMassSpringSolver(const dtkPhysMassSpring::Ptr& massSpring); + + void localStep(); + void globalStep(); + + Cholesky _system_matrix; + dtkPhysMassSpring::Ptr _system; + + // M, L, J matrices + SparseMatrix _M;/**< 质量稀疏矩阵 */ + SparseMatrix _L;/**< Laplace稀疏矩阵，表示弹簧质点系统中质点与质点之间的连接关系，以及连接弹簧的刚度 */ + SparseMatrix _J;/**< Jacobian稀疏矩阵，表示弹簧质点系统中每个质点与弹簧的连接关系，以及连接弹簧的刚度 */ + + VectorXf _initial_state; // q(0) + VectorXf _current_state; // q(n) + VectorXf _prev_state; // q(n-1) + VectorXf _spring_directions; // d, spring directions + VectorXf _inertial_term; /**< 惯性项 = M * y, y = (a + 1) * q(n) - a * q(n - 1), a = damp_factor */ + + float _time_step; + }; +} +#endif // SIMPLEPHYSICSENGINE_DTKPHYSMASSSPRINGSOLVER_H \ No newline at end of file diff --git a/test/system_test/MassSpring3D/main.cpp b/test/system_test/MassSpring3D/main.cpp new file mode 100644 index 0000000..d9ff20f --- /dev/null +++ b/test/system_test/MassSpring3D/main.cpp @@ -0,0 +1,250 @@ +/** + * @file ClothSimulation.cpp + * @brief ClothSimulation main 实现 + * @author Chance + * @version 1.0 + * @date 2024-07-12 + * @copyright MIT LICENSE + * https://github.com/Simple-XX/SimplePhysicsEngine + * @par change log: + * + *
DateAuthorDescription + *
2024-07-12Chance创建文件 + *
+ */ + +#include +#include +#include + +#include +#include +#include +#include + // #include + +#include "ClothSimulation.h" +#include "dtkStaticTriangleMesh.h" +#include "Shader.h" +#include "Renderer.h" + +static auto last_clock = std::chrono::high_resolution_clock::now(); +static auto init_clock = std::chrono::high_resolution_clock::now(); + +// The Width of the screen +const unsigned int WINDOW_WIDTH = 800; +// The height of the screen +const unsigned int WINDOW_HEIGHT = 600; +static ClothSimulation scene(WINDOW_WIDTH, WINDOW_HEIGHT, { 0, -9.8 }); +// static ProgramInput* g_render_target; // vertex, index + +bool TEST = false; +char INSTRUCTION = '0'; + +static void checkGlErrors(); + +static void draw_text(int x, int y, const char* format, ...) { + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + int w = glutGet(GLUT_WINDOW_WIDTH); + int h = glutGet(GLUT_WINDOW_HEIGHT); + gluOrtho2D(0, w, h, 0); + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + + glColor3f(0.9f, 0.9f, 0.9f); + glRasterPos2i(x, y); + + char buffer[256]; + va_list args; + va_start(args, format); + int len = vsnprintf(buffer, 256, format, args); + va_end(args); + for (int i = 0; i < len; ++i) { + glutBitmapCharacter(GLUT_BITMAP_9_BY_15, buffer[i]); + } + + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); +} + +void display() { + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + glTranslatef(0.0f, -8.0f, -25.0f); + + auto now = std::chrono::high_resolution_clock::now(); + auto dt = std::chrono::duration_cast>( + now - last_clock) + .count(); + if (TEST) { + auto total_time = std::chrono::duration_cast>( + now - init_clock) + .count(); + if (total_time > 5.0) { + std::cout << "TEST " << INSTRUCTION << "." << std::endl; + exit(1); + } + } + last_clock = now; + + int h = glutGet(GLUT_WINDOW_HEIGHT); + int w = glutGet(GLUT_WINDOW_WIDTH); + + draw_text(5, 20, "dtk @SoftBody simulation"); + + draw_text(5, 40, "Push [1-1] to switch scene"); + draw_text(w - 150, h - 20, "refer: apollonia"); + + if (scene.IsPause()) + draw_text(5, h - 20, "dt: %.2f ms PAUSED", dt * 1000); + else + draw_text(5, h - 20, "dt: %.2f ms", dt * 1000); + + scene.Update(std::min(dt, 0.08)); + scene.Render(); + + GLenum err; + while ((err = glGetError()) != GL_NO_ERROR) { + std::cerr << "OpenGL error: " << err << std::endl; + } + + glutSwapBuffers(); + checkGlErrors(); +} + +void reshape(int width, int height) { + glViewport(0, 0, width, height); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluPerspective(45.0, width / (float)height, 0.1, 100.0); +} + +void mouse(int button, int state, int x, int y) {} + +void move_pos(const dtk::dtkDouble2& v) { scene.move(v); } + +void keyboard(unsigned char key, int x, int y) { + switch (key) { + case '1': + scene.SetVisible(!scene.IsVisible()); + break; + case 'w': + move_pos(dtk::dtkDouble2(0, 1)); + break; + case 'a': + move_pos(dtk::dtkDouble2(-1, 0)); + break; + case 's': + move_pos(dtk::dtkDouble2(0, -1)); + break; + case 'd': + move_pos(dtk::dtkDouble2(1, 0)); + break; + case ' ': + scene.SetPause(!scene.IsPause()); + break; + case 27: + exit(0); + break; + default: + break; + } +} + +void motion(int x, int y) {} + +void special(int key, int x, int y) {} + +void idle() { + display(); +} + +static void initGlutState(int argc, char** argv, const char* window_title = "", const unsigned int window_width = 800, const unsigned int window_height = 600) { + glutInit(&argc, argv); + glutInitWindowSize(window_width, window_height); + glutInitWindowPosition(50, 50); + glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH); /// TODO + glutCreateWindow(window_title); + glutDisplayFunc(&display); + glutReshapeFunc(&reshape); + glutMouseFunc(&mouse); + glutMotionFunc(&motion); + glutSpecialFunc(&special); + glutKeyboardFunc(&keyboard); + glutIdleFunc(&idle); +} + +static void initGlewState() { + GLenum err = glewInit(); + // if (!glewIsSupported("GL_VERSION_2_0")) { + // printf("OpenGL 2.0 not supported\n"); + // exit(1); + // } + if (err != GLEW_OK) { + std::cerr << "Error initializing GLEW: " << glewGetErrorString(err) << std::endl; + exit(0); + } +} + +static void checkGlErrors() { + const GLenum errCode = glGetError(); + + if (errCode != GL_NO_ERROR) { + std::string error("GL Error: "); + error += reinterpret_cast(gluErrorString(errCode)); + std::cerr << error << std::endl; + throw std::runtime_error(error); + } +} + +static void initGLState() { + glClearColor(0.25f, 0.25f, 0.25f, 0); + glClearDepth(1.); + glPixelStorei(GL_UNPACK_ALIGNMENT, 1); + glPixelStorei(GL_PACK_ALIGNMENT, 1); + glEnable(GL_DEPTH_TEST); + glDepthFunc(GL_LESS); + glReadBuffer(GL_BACK); + glEnable(GL_FRAMEBUFFER_SRGB); + + checkGlErrors(); +} + +int main(int argc, char* argv[]) { + // Parse command line arguments + for (int i = 1; i < argc; ++i) { + if (std::strcmp(argv[i], "--test") == 0) { + TEST = true; + } + else if (std::strcmp(argv[i], "--instruction") == 0 && i + 1 < argc) { + INSTRUCTION = argv[++i][0]; + } + } + + try { + const char* window_title = "SimplePhysicsEngine-ST-MassSpring3D"; + initGlutState(argc, argv, window_title, WINDOW_WIDTH, WINDOW_HEIGHT); + initGlewState(); + initGLState(); + + scene.Init(); + checkGlErrors(); + + if (INSTRUCTION != '0') { + keyboard(INSTRUCTION, 0, 0); + } + glutMainLoop(); + + scene.CleanUp(); + return 0; + } + catch (const std::runtime_error& e) { + std::cout << "Exception caught: " << e.what() << std::endl; + return -1; + } +} \ No newline at end of file diff --git a/test/system_test/MassSpring3D/shaders/basic.vshader b/test/system_test/MassSpring3D/shaders/basic.vshader new file mode 100644 index 0000000..b6128bf --- /dev/null +++ b/test/system_test/MassSpring3D/shaders/basic.vshader @@ -0,0 +1,19 @@ +#version 430 + +uniform mat4 uModelViewMatrix; +// uniform mat4 uNormalMatrix; +uniform mat4 uProjectionMatrix; + +layout(location=0) in vec3 aPosition; +layout(location=1) in vec3 aNormal; +layout(location=2) in vec2 aTexCoord; + +out vec3 vNormal; +out vec2 vTexCoord; + +void main(){ + vNormal = aNormal; + vTexCoord = aTexCoord; + vec4 position = vec4(aPosition, 1.0); + gl_Position = uProjectionMatrix * uModelViewMatrix * position; +} \ No newline at end of file diff --git a/test/system_test/MassSpring3D/shaders/phong.fshader b/test/system_test/MassSpring3D/shaders/phong.fshader new file mode 100644 index 0000000..c473e87 --- /dev/null +++ b/test/system_test/MassSpring3D/shaders/phong.fshader @@ -0,0 +1,22 @@ +#version 430 + +in vec3 vNormal; +out vec4 fragColor; + +uniform vec3 uAlbedo; +uniform vec3 uAmbient; +uniform vec3 uLight; +void main(){ + vec3 normal = normalize(vNormal); + vec3 albedo = uAlbedo; + if(!gl_FrontFacing) { + normal = -normal; + albedo = 1.0 - albedo; + } + + vec3 toLight = normalize(-uLight); + float diffuse = max(0, dot(toLight, normal)); + + vec3 color = diffuse * albedo + uAmbient * albedo; + fragColor = vec4(color, 1.0); +} \ No newline at end of file diff --git a/test/system_test/MassSpring3D/shaders/pick.fshader b/test/system_test/MassSpring3D/shaders/pick.fshader new file mode 100644 index 0000000..9e21c5b --- /dev/null +++ b/test/system_test/MassSpring3D/shaders/pick.fshader @@ -0,0 +1,12 @@ +#version 430 + +in vec2 vTexCoord; +out vec4 fragColor; + +uniform int uTessFact; + +void main(){ + float vx = round(vTexCoord.x * (uTessFact - 1)); + float vy = round(vTexCoord.y * (uTessFact - 1)); + fragColor = vec4(vx / (uTessFact - 1), vy / (uTessFact - 1), 0.2, 1.0); +} \ No newline at end of file diff --git a/test/system_test/demo2d/main.cpp b/test/system_test/demo2d/main.cpp index 485c9e3..17137ed 100644 --- a/test/system_test/demo2d/main.cpp +++ b/test/system_test/demo2d/main.cpp @@ -24,14 +24,18 @@ #include "FemSimulation.h" static auto last_clock = std::chrono::high_resolution_clock::now(); +static auto init_clock = std::chrono::high_resolution_clock::now(); // The Width of the screen const unsigned int SCREEN_WIDTH = 800; // The height of the screen const unsigned int SCREEN_HEIGHT = 600; -static FemSimulation world({0, -9.8}); +static FemSimulation world({ 0, -9.8 }); -static void draw_text(int x, int y, const char *format, ...) { +bool TEST = false; +char INSTRUCTION = '0'; + +static void draw_text(int x, int y, const char* format, ...) { glMatrixMode(GL_PROJECTION); glPushMatrix(); glLoadIdentity(); @@ -59,7 +63,7 @@ static void draw_text(int x, int y, const char *format, ...) { glPopMatrix(); } -static void draw_divided_polys(dtk::dtkPolygonRigidBody &body) { +static void draw_divided_polys(dtk::dtkPolygonRigidBody& body) { int len1 = body.mPolygonList.size(); for (int i = 0; i < len1; i++) { glBegin(GL_LINE_LOOP); @@ -72,10 +76,11 @@ static void draw_divided_polys(dtk::dtkPolygonRigidBody &body) { } } -static void draw_body(const dtk::dtkPolygonRigidBody &body) { +static void draw_body(const dtk::dtkPolygonRigidBody& body) { if (std::isinf(body.get_mass())) { glColor3f(1.0f, 1.0f, 1.0f); - } else { + } + else { glColor3f(0.8f, 0.8f, 0.0f); } glBegin(GL_LINE_LOOP); @@ -85,7 +90,7 @@ static void draw_body(const dtk::dtkPolygonRigidBody &body) { } glEnd(); if (!std::isinf(body.get_mass())) { - auto &pos = body.get_position(); + auto& pos = body.get_position(); auto v = body.local_to_world(body.get_velocity() * 0.2); glBegin(GL_LINES); glColor3f(0.0f, 1.0f, 0.0f); @@ -99,7 +104,7 @@ static void draw_body(const dtk::dtkPolygonRigidBody &body) { } } -static void draw_mesh(const Mesh &mesh) { +static void draw_mesh(const Mesh& mesh) { glColor3f(0.8f, 0.8f, 0.0f); glBegin(GL_LINES); for (int i = 0; i < mesh.n_fem_element_; ++i) { @@ -116,7 +121,7 @@ static void draw_mesh(const Mesh &mesh) { glEnd(); } -static void draw_mesh_shell(const Mesh &mesh) { +static void draw_mesh_shell(const Mesh& mesh) { glColor3f(0.8f, 0.0f, 0.0f); glBegin(GL_LINE_LOOP); @@ -128,12 +133,12 @@ static void draw_mesh_shell(const Mesh &mesh) { // draw_divided_polys(*(mesh.shell)); } -static void draw_joint(const dtk::dtkRevoluteJoint &joint) { +static void draw_joint(const dtk::dtkRevoluteJoint& joint) { auto centroid_a = - joint.get_a()->local_to_world(joint.get_a()->get_centroid()); + joint.get_a()->local_to_world(joint.get_a()->get_centroid()); auto anchor_a = joint.world_anchor_a(); auto centroid_b = - joint.get_b()->local_to_world(joint.get_b()->get_centroid()); + joint.get_b()->local_to_world(joint.get_b()->get_centroid()); auto anchor_b = joint.world_anchor_b(); glColor3f(0.6f, 0.6f, 0.6f); @@ -149,9 +154,9 @@ static void draw_joint(const dtk::dtkRevoluteJoint &joint) { glEnd(); } -static void draw_arbiter(const CollisionPair::ptr &pair) { - auto &contacts = pair->get_contacts(); - for (auto &contact : contacts) { +static void draw_arbiter(const CollisionPair::ptr& pair) { + auto& contacts = pair->get_contacts(); + for (auto& contact : contacts) { auto pos = contact.position; auto ra = pos + dtk::normalize(contact.ra) * 0.2; auto rb = pos + dtk::normalize(contact.rb) * 0.2; @@ -179,12 +184,12 @@ void drawParticle(float x, float y, float s) { glEnd(); } -void draw_sph(SPHSolver &sph) { +void draw_sph(SPHSolver& sph) { glMatrixMode(GL_MODELVIEW); glLoadIdentity(); glTranslatef(-1.5f, -1.0f, -5.0f); // for each (auto & particle in sph.particles) - for (auto &particle : sph.particles) { + for (auto& particle : sph.particles) { glColor3f(0.5, 1.0, 1.0); drawParticle(particle.position[0], particle.position[1], 0.01); } @@ -196,19 +201,19 @@ double random(double low, double high) { static void test_polygon() { dtkFactory::make_fence(world); - world.add(dtkFactory::make_polygon(200, {{-4, 0}, {4, 0}, {4, 1}, {-4, 1}}, - {0, 0})); - world.add(dtkFactory::make_polygon(200, {{-1, -2}, {1, -2}, {1, 2}, {-1, 2}}, - {0, 4})); - world.add(dtkFactory::make_polygon(200, {{2, -1}, {0, 1}, {-2, -1}, {0, 0}}, - {0, 10})); + world.add(dtkFactory::make_polygon(200, { {-4, 0}, {4, 0}, {4, 1}, {-4, 1} }, + { 0, 0 })); + world.add(dtkFactory::make_polygon(200, { {-1, -2}, {1, -2}, {1, 2}, {-1, 2} }, + { 0, 4 })); + world.add(dtkFactory::make_polygon(200, { {2, -1}, {0, 1}, {-2, -1}, {0, 0} }, + { 0, 10 })); } static void test_stack() { dtkFactory::make_fence(world); for (int i = 0; i < 10; ++i) { double x = random(-0.1 * i, 0.1 * i); - auto body = dtkFactory::make_box(1, 1, 1, {x, 0.51f + 1.05f * i}); + auto body = dtkFactory::make_box(1, 1, 1, { x, 0.51f + 1.05f * i }); body->set_friction(0.2); world.add(body); } @@ -232,24 +237,24 @@ static void test_pyramid() { } static void test_joint() { - auto ground = dtkFactory::make_box(dtk::inf, 100, 20, {0, -10}); + auto ground = dtkFactory::make_box(dtk::inf, 100, 20, { 0, -10 }); world.add(ground); - auto box1 = dtkFactory::make_box(500, 1, 1, {13.5, 11}); + auto box1 = dtkFactory::make_box(500, 1, 1, { 13.5, 11 }); world.add(box1); - auto joint1 = dtkFactory::make_revolute_joint(ground, box1, {4.5, 11}); + auto joint1 = dtkFactory::make_revolute_joint(ground, box1, { 4.5, 11 }); world.add(joint1); for (size_t i = 0; i < 5; ++i) { - auto box2 = dtkFactory::make_box(100, 1, 1, {3.5f - i, 2}); + auto box2 = dtkFactory::make_box(100, 1, 1, { 3.5f - i, 2 }); world.add(box2); - auto joint2 = dtkFactory::make_revolute_joint(ground, box2, {3.5f - i, 11}); + auto joint2 = dtkFactory::make_revolute_joint(ground, box2, { 3.5f - i, 11 }); world.add(joint2); } } static void test_chain() { - auto ground = dtkFactory::make_box(dtk::inf, 100, 20, {0, -10}); + auto ground = dtkFactory::make_box(dtk::inf, 100, 20, { 0, -10 }); ground->set_friction(0.4); world.add(ground); @@ -257,11 +262,11 @@ static void test_chain() { const double y = 12.0f; auto last = ground; for (int i = 0; i < 15; ++i) { - auto box = dtkFactory::make_box(mass, 0.75, 0.25, {0.5f + i, y}); + auto box = dtkFactory::make_box(mass, 0.75, 0.25, { 0.5f + i, y }); box->set_friction(0.4); world.add(box); auto joint = - dtkFactory::make_revolute_joint(last, box, dtk::dtkDouble2(i, y)); + dtkFactory::make_revolute_joint(last, box, dtk::dtkDouble2(i, y)); world.add(joint); last = box; } @@ -271,8 +276,8 @@ static void test_mesh() { dtkFactory::make_fence(world); auto mesh = dtkFactory::make_mesh(15, 3, dtk::dtkDouble2(0.0, 6.0)); world.add(mesh); - auto box = dtkFactory::make_polygon(200, {{-1, 0}, {1, 0}, {1, 2}, {-1, 2}}, - {0.2, 0}); + auto box = dtkFactory::make_polygon(200, { {-1, 0}, {1, 0}, {1, 2}, {-1, 2} }, + { 0.2, 0 }); world.add(box); } @@ -289,8 +294,18 @@ void display() { auto now = std::chrono::high_resolution_clock::now(); auto dt = std::chrono::duration_cast>( - now - last_clock) - .count(); + now - last_clock) + .count(); + if (TEST) { + auto total_time = std::chrono::duration_cast>( + now - init_clock) + .count(); + if (total_time > 5.0) { + std::cout << "TEST " << INSTRUCTION << "." << std::endl; + exit(0); + } + } + last_clock = now; int h = glutGet(GLUT_WINDOW_HEIGHT); @@ -308,24 +323,24 @@ void display() { world.step(std::min(dt, 0.01)); - for (auto &body : world.get_bodies()) { + for (auto& body : world.get_bodies()) { draw_body(*std::dynamic_pointer_cast(body).get()); } - for (auto &joint : world.get_joints()) { + for (auto& joint : world.get_joints()) { draw_joint(*std::dynamic_pointer_cast(joint).get()); } - for (auto &arbiter : world.get_arbiters()) { + for (auto& arbiter : world.get_arbiters()) { draw_arbiter(arbiter.second); } - for (auto &mesh : world.get_meshes()) { + for (auto& mesh : world.get_meshes()) { draw_mesh(*std::dynamic_pointer_cast(mesh).get()); draw_mesh_shell(*std::dynamic_pointer_cast(mesh).get()); } - for (auto &sph : world.get_sphs()) { + for (auto& sph : world.get_sphs()) { draw_sph(*std::dynamic_pointer_cast(sph).get()); } @@ -341,7 +356,7 @@ void reshape(int width, int height) { void mouse(int button, int state, int x, int y) {} -void move_pos(const dtk::dtkDouble2 &v) { world.move(v); } +void move_pos(const dtk::dtkDouble2& v) { world.move(v); } void keyboard(unsigned char key, int x, int y) { switch (key) { @@ -402,7 +417,17 @@ void special(int key, int x, int y) {} void idle() { display(); } -int main(int argc, char *argv[]) { +int main(int argc, char* argv[]) { + // Parse command line arguments + for (int i = 1; i < argc; ++i) { + if (std::strcmp(argv[i], "--test") == 0) { + TEST = true; + } + else if (std::strcmp(argv[i], "--instruction") == 0 && i + 1 < argc) { + INSTRUCTION = argv[++i][0]; + } + } + glutInit(&argc, argv); glutInitWindowSize(SCREEN_WIDTH, SCREEN_HEIGHT); glutInitWindowPosition(50, 50); @@ -418,6 +443,9 @@ int main(int argc, char *argv[]) { glutIdleFunc(&idle); /// @todo not found // glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION); + if (INSTRUCTION != '0') { + keyboard(INSTRUCTION, 0, 0); + } glutMainLoop(); return 0; }