Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include "task2.h"
- int flattenIndex(int x, int y, int size_x) {
- // This is a utility function that maps a 2D index into a 1D index
- // for our flat vector structures (normals, positions, masses...)
- return x + y * size_x;
- }
- void calcNormals(std::vector<float3>& normals,
- const std::vector<float3>& positions, int grid_x, int grid_y) {
- // TODO: calculate the normal for each particle vertex
- //
- // To compute the particle normal consider the particles connected via
- // structural springs. Use flattenIndex(int x, int y, int size_x) to convert
- // from 2D coordinates to a 1D index for vectors normals and positions.
- //
- // Parameters:
- // normals: resulting normals
- // positions: particle positions
- // grid_x, grid_y: particle-grid dimensions
- for(int y = 0; y < grid_y; y++)
- {
- for(int x = 0; x < grid_x; x++)
- {
- float3 curr_index = positions.at(flattenIndex(x,y,grid_x));
- float3 right_down = {0.0f,0.0f,0.0f};
- float3 left_down = {0.0f,0.0f,0.0f};
- float3 right_up = {0.0f,0.0f,0.0f};
- float3 left_up = {0.0f,0.0f,0.0f};
- // float3 temp_feder = {0.0f, 0.0f, 0.0f};
- //float3 temp_feder2 = {0.0f, 0.0f, 0.0f};
- if(x - 1 >= 0 && y + 1 < grid_y)
- {
- float3 temp_feder = positions.at(flattenIndex(x,y+1,grid_x)) - curr_index;
- float3 temp_feder2 = positions.at(flattenIndex(x-1,y,grid_x)) - curr_index;
- left_down = normalize(cross(temp_feder,temp_feder2));
- }
- if(x - 1 >= 0 && y - 1 >= 0)
- {
- float3 temp_feder = positions.at(flattenIndex(x-1,y,grid_x)) - curr_index;
- float3 temp_feder2 = positions.at(flattenIndex(x,y - 1,grid_x)) - curr_index;
- left_up = normalize(cross(temp_feder,temp_feder2));
- }
- if(x + 1 < grid_x && y -1 >= 0)
- {
- float3 temp_feder = positions.at(flattenIndex(x,y-1,grid_x)) - curr_index;
- float3 temp_feder2 = positions.at(flattenIndex(x+1,y,grid_x)) - curr_index;
- right_up = normalize(cross(temp_feder,temp_feder2));
- }
- if(x + 1 < grid_x && y + 1 < grid_y)
- {
- float3 temp_feder = positions.at(flattenIndex(x + 1,y,grid_x)) - curr_index;
- float3 temp_feder2 = positions.at(flattenIndex(x,y + 1,grid_x)) - curr_index;
- right_down = normalize(cross(temp_feder,temp_feder2));
- }
- int index = flattenIndex(x, y, grid_x);
- normals.at(index) = normalize(left_down + left_up + right_up + right_down);
- }
- }
- }
- float3 calcWindAcc(const float3& normal, const float3& wind_force,
- float inv_mass) {
- // TODO: calculate the particle-dependent acceleration due to wind
- //
- // The wind strength should depend on the cosine (dot product) between the
- // particle normal and the wind direction (note that the wind can hit the
- // particle from the front or back). The resulting particle acceleration
- // caused by the wind is the above scaling times the inverse mass of the
- // particle times the wind force.
- //
- // Parameters:
- // normal: particle normal
- // wind_force: wind force in XYZ direction
- // inv_mass: inverse mass of the particle
- if(wind_force.x != 0.0f && wind_force.y != 0.0f && wind_force.z != 0.0f)
- {
- float3 wind_force_normal = normalize(wind_force);
- float save = std::fabs(dot(normal, wind_force_normal));
- float3 acc = inv_mass * wind_force * save;
- return acc;
- }
- return float3(0.f, 0.f, 0.f);
- }
- void verletIntegration(std::vector<float3>& p_next,
- const std::vector<float3>& p_cur,
- const std::vector<float3>& p_prev,
- const std::vector<float3>& normals,
- const std::vector<float>& inv_masses, float damp,
- const float3& gravity, const float3& wind_force,
- float dt) {
- // TODO: for each particle perform the numerical integration
- //
- // Use the formula given in the assignment sheet to perform the Verlet
- // Integration. Do NOT add spring forces, as they are taken care of by the
- // fixup-step! Make sure that fixed particles (inverse mass = 0) do not change
- // their position! Bonus: The particle-dependent accelerations should also
- // consider the wind direction and strength => Use your implementation of
- // calcWindAcc()
- //
- // Parameters:
- // p_next: updated particle positions
- // p_cur: current paricle positions
- // p_prev: previous particle positions
- // normals: particle normals
- // inv_masses: inverse particle masses
- // damp: damping factor
- // gravity: gravity vector
- // wind_force: wind force
- // dt: timestep
- for(int itr = 0; itr < p_next.size(); itr++)
- p_next[itr] = inv_masses[itr] == 0.0f ? p_cur[itr] : (1 + damp) * p_cur[itr] - damp * p_prev[itr] +
- (gravity + calcWindAcc(normals[itr],wind_force,inv_masses[itr])) * powf(dt,2.0f);
- }
- void fixupStep(std::vector<float3>& p_out, const std::vector<float3>& p_in,
- int grid_x, int grid_y, float cloth_x, float cloth_y,
- const std::vector<float>& inv_masses, float fixup_percent,
- const std::vector<Sphere>& obstacles,
- bool apply_structural_fixup, bool apply_shear_fixup,
- bool apply_flexion_fixup) {
- // TODO: implement the fixup-step
- //
- // For each particle iterate over all particles connected via springs and
- // compute the movement required to restore each spring into its resting
- // position. Based on the mass ratio between the particle pairs compute the
- // part of the movement that would have to be carried out by the current
- // particle. Sum these movements up and update the particle position with a
- // fraction of that movement(fixup_percent of the full movement). Iterate over
- // all obstacles and resolve the collisions.
- //
- // *Only* apply the relevant fixup (structural/shear/flexion) if the
- // corresponding flag is set (apply_*_fixup)!
- //
- // Parameters:
- // p_out: new particle positions
- // p_in: current particle positions
- // grid_x, grid_y: particle-grid dimensions
- // cloth_x, cloth_y: cloth size
- // inv_masses: inverse particle masses
- // fix_percent: fixup percent [0..1]
- // obstacles: obstacle vector
- // apply_structural_fixup: only structural fixup if this is true
- // apply_shear_fixup: only apply shear fixup if this is true
- // apply_flexion_fixup: only apply flexion fixup if this is true
- float flex_x = (cloth_x/(grid_x - 1.0f)) * 2.0f;
- float flex_y = (cloth_y/(grid_y - 1.0f)) * 2.0f;
- float shear = std::sqrt(std::pow(cloth_x/(grid_x - 1.0f), 2) + std::pow(cloth_y/(grid_y - 1.0f), 2));
- for(int y = 0; y < grid_y; y++)
- {
- for(int x = 0; x < grid_x; x++)
- {
- p_out.at(flattenIndex(x,y,grid_x)) = p_in.at(flattenIndex(x,y,grid_x));
- float offset = 0.0f;
- float distance = 0.0f;
- float mass = 0.0f;
- if(inv_masses.at(flattenIndex(x,y,grid_x)) > 0.0f)
- {
- /////////apply shear fixup /////////////////////////////
- ///////////////////////////////////////////////////////
- ///////////////////////////////////////////////////////
- if(apply_shear_fixup)
- {
- ///reset variables
- float3 feder = {0.0f,0.0f,0.0f};
- float3 particle1 = {0.0f,0.0f,0.0f};
- float3 particle2 = {0.0f,0.0f,0.0f};
- float3 particle3 = {0.0f,0.0f,0.0f};
- float3 particle4 = {0.0f,0.0f,0.0f};
- float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- if(x > 0 && y > 0)
- {
- feder = p_in.at(flattenIndex(x - 1, y - 1, grid_x));
- mass = inv_masses.at(flattenIndex(x - 1, y - 1, grid_x));
- // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - shear;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle1 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- //////////
- if(x < grid_x - 1 && y < grid_y - 1)
- {
- feder = p_in.at(flattenIndex(x + 1, y + 1, grid_x));
- mass = inv_masses.at(flattenIndex(x + 1, y + 1, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - shear;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle2 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(x > 0 && y < grid_y - 1)
- {
- feder = p_in.at(flattenIndex(x - 1, y + 1, grid_x));
- mass = inv_masses.at(flattenIndex(x - 1, y + 1, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - shear;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle3 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(x < grid_x - 1 && y > 0)
- {
- feder = p_in.at(flattenIndex(x + 1, y - 1, grid_x));
- mass = inv_masses.at(flattenIndex(x + 1, y - 1, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - shear;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle4 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) + ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
- }
- ////////////////////////////////////////////////////////////////////////
- ///apply flexion fixup ////////////////////////////////
- //////////////////////////////////////////////////////////////////////
- if(apply_flexion_fixup)
- {
- ////reset variables
- float3 feder = {0.0f,0.0f,0.0f};
- float3 particle1 = {0.0f,0.0f,0.0f};
- float3 particle2 = {0.0f,0.0f,0.0f};
- float3 particle3 = {0.0f,0.0f,0.0f};
- float3 particle4 = {0.0f,0.0f,0.0f};
- // offset = distance - shear;
- float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- if( x - 1 > 0 && x > 0)
- {
- feder = p_in.at(flattenIndex(x - 2, y, grid_x));
- mass = inv_masses.at(flattenIndex(x - 2, y, grid_x));
- // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- distance =length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - flex_x;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle1 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(x < grid_x - 1 && x + 1 < grid_x - 1)
- {
- feder = p_in.at(flattenIndex(x + 2, y, grid_x));
- mass = inv_masses.at(flattenIndex(x + 2, y, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - flex_x;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle2 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(y + 1 < grid_y - 1 && y < grid_y - 1)
- {
- feder = p_in.at(flattenIndex(x, y + 2, grid_x));
- mass = inv_masses.at(flattenIndex(x, y + 2, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - flex_y;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle3 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(y - 1 > 0 && y > 0)
- {
- feder = p_in.at(flattenIndex(x, y - 2, grid_x));
- mass = inv_masses.at(flattenIndex(x , y - 2, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - flex_y;
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle4 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) + ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
- }
- //////////////////////////////////////////////////////////////////////////////
- //////////apply structural fixup /////////////////////////////////////////////
- /////////////////////////////////////////////////////////////////////////////
- if(apply_structural_fixup)
- {
- ////reset variables
- float3 feder = {0.0f,0.0f,0.0f};
- float3 particle1 = {0.0f,0.0f,0.0f};
- float3 particle2 = {0.0f,0.0f,0.0f};
- float3 particle3 = {0.0f,0.0f,0.0f};
- float3 particle4 = {0.0f,0.0f,0.0f};
- // offset = distance - shear;
- float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- if(y < grid_y - 1)
- {
- feder = p_in.at(flattenIndex(x, y + 1, grid_x));
- mass = inv_masses.at(flattenIndex(x, y + 1, grid_x));
- // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - (cloth_y/(grid_y - 1.0f));
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle1 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(x < grid_x - 1)
- {
- feder = p_in.at(flattenIndex(x + 1, y, grid_x));
- mass = inv_masses.at(flattenIndex(x + 1, y, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - (cloth_x/(grid_x - 1.0f));
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle2 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(x > 0)
- {
- feder = p_in.at(flattenIndex(x - 1 ,y, grid_x));
- mass = inv_masses.at(flattenIndex(x - 1, y, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - (cloth_x/(grid_x - 1.0f));
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle3 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- if(y > 0)
- {
- feder = p_in.at(flattenIndex(x, y - 1, grid_x));
- mass = inv_masses.at(flattenIndex(x , y - 1, grid_x));
- distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
- offset = distance - (cloth_y/(grid_y - 1.0f));
- float mass_end = (beginning_mass / (mass + beginning_mass));
- particle4 = offset * normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
- }
- p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) + ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
- }
- }
- ////////////////////////////////////////////////
- //////////////// fix collission ////////////////
- ////////////////////////////////////////////////
- for(unsigned int itr = 0; itr < obstacles.size(); ++itr)
- {
- fixCollision(p_out.at(flattenIndex(x,y,grid_x)), obstacles.at(itr));
- }
- }
- }
- }
- void fixCollision(float3& particle_position, const Sphere& sphere) {
- // TODO: resolve particle/obstacle collision
- //
- // Check if the particle is inside the sphere.
- // If so, move the particle out of the sphere (along the sphere normal).
- //
- // Parameters:
- // particle_position: position of the particle
- // sphere: spherical obstacle containg the position and radius
- if(length(particle_position - sphere.position) < sphere.radius)
- {
- float3 temp = normalize(particle_position - sphere.position);
- particle_position = sphere.position + temp * sphere.radius;
- }
- }
- void simulationStep(std::array<std::vector<float3>, 3>& p_buffers, int& next,
- int& current, int& previous, int grid_x, int grid_y,
- float cloth_x, float cloth_y,
- const std::vector<float>& inv_masses,
- std::vector<float3>& normals, const float3& gravity,
- const float3& wind_force, float damping,
- float fixup_percent, int fixup_iterations,
- const std::vector<Sphere>& obstacles, float dt,
- bool apply_structural_fixup, bool apply_shear_fixup,
- bool apply_flexion_fixup) {
- // TODO: implement the simulation step
- //
- // 1. Calculate particle normals => calcNormals()
- // 2. Perform Verlet Integration => verletIntegration()
- // 3. Perform fixup steps => fixupStep()
- //
- // Verlet Integration:
- // You will need three position buffers, since the verlet integration needs
- // the current und previous positions. The integrated position shall be
- // written to the 'next' position buffer. Make sure to update the next,
- // current, previous indices! e.g. after the integration step: next ->
- // current, current -> last and last -> next
- //
- // Fixup:
- // Repeatedly execute the FixUp step after each integration
- // ('fixup_iterations' times). Make sure to update your position buffer
- // indices in the right way!
- //
- // Make sure that 'current' points to the right position buffer before you
- // return!
- //
- // Parameters:
- // p_buffers: 3 position buffers
- // next: index of next position buffer
- // current: index of current position buffer
- // previous: index of previous position buffer
- // grid_x, grid_y: particle-grid dimensions
- // cloth_x, cloth_y: cloth size
- // inv_masses: inverse particle masses
- // normals: particle normals
- // gravity: gravity vector
- // wind_force: wind force
- // damp: damping factor for verlet integration
- // fixup_percent: fixup percent [0..1]
- // fixup_iterations: how often the fixup step should be executed
- // obstacles: obstacle vector
- // dt: timestep
- //
- // apply_structural_fixup: apply structural fixup
- // (needs to be passed unmodified to fixupStep for automatic testing)
- //
- // apply_shear_fixup: apply shear fixup
- // (needs to be passed unmodified to fixupStep for automatic testing)
- //
- // apply_flexion_fixup: apply flexion fixup
- // (needs to be passed unmodified to fixupStep for automatic testing)
- //
- // Simply make cloth fall down. You will want to replace that ;)
- calcNormals(normals, p_buffers[current], grid_x, grid_y);
- verletIntegration(p_buffers[next], p_buffers[current], p_buffers[previous], normals, inv_masses, damping,
- gravity, wind_force, dt);
- int save = next;
- next = previous;
- previous = current;
- current = save;
- for(unsigned int iteration = 0; iteration < fixup_iterations; iteration++)
- {
- fixupStep(p_buffers[next], p_buffers[current], grid_x, grid_y, cloth_x, cloth_y, inv_masses, fixup_percent,
- obstacles, apply_structural_fixup, apply_shear_fixup, apply_flexion_fixup);
- int save = current;
- current = next;
- next = save;
- }
- //for (auto& p : p_buffers[current]) p += gravity * dt;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement