SHARE
TWEET

Untitled

a guest Jun 12th, 2019 75 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include "task2.h"
  2.  
  3. int flattenIndex(int x, int y, int size_x) {
  4.   // This is a utility function that maps a 2D index into a 1D index
  5.   // for our flat vector structures (normals, positions, masses...)
  6.   return x + y * size_x;
  7. }
  8.  
  9. void calcNormals(std::vector<float3>& normals,
  10.                  const std::vector<float3>& positions, int grid_x, int grid_y) {
  11.   // TODO: calculate the normal for each particle vertex
  12.   //
  13.   // To compute the particle normal consider the particles connected via
  14.   // structural springs. Use flattenIndex(int x, int y, int size_x) to convert
  15.   // from 2D coordinates to a 1D index for vectors normals and positions.
  16.   //
  17.   // Parameters:
  18.   // normals: resulting normals
  19.   // positions: particle positions
  20.   // grid_x, grid_y: particle-grid dimensions
  21.  
  22.  
  23.     for(int y = 0; y < grid_y; y++)
  24.     {
  25.         for(int x = 0; x < grid_x; x++)
  26.         {
  27.             float3 curr_index = positions.at(flattenIndex(x,y,grid_x));
  28.             float3 right_down = {0.0f,0.0f,0.0f};
  29.             float3 left_down = {0.0f,0.0f,0.0f};
  30.             float3 right_up = {0.0f,0.0f,0.0f};
  31.             float3 left_up = {0.0f,0.0f,0.0f};
  32.            // float3 temp_feder = {0.0f, 0.0f, 0.0f};
  33.             //float3 temp_feder2 = {0.0f, 0.0f, 0.0f};
  34.  
  35.  
  36.             if(x - 1 >= 0 && y + 1 < grid_y)
  37.             {
  38.                 float3 temp_feder = positions.at(flattenIndex(x,y+1,grid_x)) - curr_index;
  39.                 float3 temp_feder2 = positions.at(flattenIndex(x-1,y,grid_x)) - curr_index;
  40.                 left_down = normalize(cross(temp_feder,temp_feder2));
  41.  
  42.             }
  43.             if(x - 1 >= 0 && y - 1  >= 0)
  44.             {
  45.                 float3 temp_feder = positions.at(flattenIndex(x-1,y,grid_x)) - curr_index;
  46.                 float3 temp_feder2 = positions.at(flattenIndex(x,y - 1,grid_x)) - curr_index;
  47.                 left_up = normalize(cross(temp_feder,temp_feder2));
  48.  
  49.             }
  50.             if(x + 1 < grid_x && y -1 >= 0)
  51.             {
  52.                 float3 temp_feder = positions.at(flattenIndex(x,y-1,grid_x)) - curr_index;
  53.                 float3 temp_feder2 = positions.at(flattenIndex(x+1,y,grid_x)) - curr_index;
  54.                 right_up = normalize(cross(temp_feder,temp_feder2));
  55.  
  56.             }
  57.             if(x + 1 < grid_x && y + 1 < grid_y)
  58.             {
  59.                 float3 temp_feder = positions.at(flattenIndex(x + 1,y,grid_x)) - curr_index;
  60.                 float3 temp_feder2 = positions.at(flattenIndex(x,y + 1,grid_x)) - curr_index;
  61.                 right_down = normalize(cross(temp_feder,temp_feder2));
  62.  
  63.             }
  64.  
  65.            
  66.             int index = flattenIndex(x, y, grid_x);
  67.             normals.at(index) = normalize(left_down + left_up + right_up + right_down);
  68.  
  69.  
  70.  
  71.  
  72.  
  73.         }
  74.     }
  75.  
  76. }
  77.  
  78. float3 calcWindAcc(const float3& normal, const float3& wind_force,
  79.                    float inv_mass) {
  80.   // TODO: calculate the particle-dependent acceleration due to wind
  81.   //
  82.   // The wind strength should depend on the cosine (dot product) between the
  83.   // particle normal and the wind direction (note that the wind can hit the
  84.   // particle from the front or back). The resulting particle acceleration
  85.   // caused by the wind is the above scaling times the inverse mass of the
  86.   // particle times the wind force.
  87.   //
  88.   // Parameters:
  89.   // normal: particle normal
  90.   // wind_force: wind force in XYZ direction
  91.   // inv_mass: inverse mass of the particle
  92.     if(wind_force.x != 0.0f  && wind_force.y != 0.0f && wind_force.z != 0.0f)
  93.     {
  94.       float3 wind_force_normal = normalize(wind_force);
  95.       float save = std::fabs(dot(normal, wind_force_normal));
  96.       float3 acc =  inv_mass * wind_force * save;
  97.       return acc;
  98.     }
  99.  
  100.  
  101.  
  102.  
  103.   return float3(0.f, 0.f, 0.f);
  104. }
  105.  
  106. void verletIntegration(std::vector<float3>& p_next,
  107.                        const std::vector<float3>& p_cur,
  108.                        const std::vector<float3>& p_prev,
  109.                        const std::vector<float3>& normals,
  110.                        const std::vector<float>& inv_masses, float damp,
  111.                        const float3& gravity, const float3& wind_force,
  112.                        float dt) {
  113.   // TODO: for each particle perform the numerical integration
  114.   //
  115.   // Use the formula given in the assignment sheet to perform the Verlet
  116.   // Integration. Do NOT add spring forces, as they are taken care of by the
  117.   // fixup-step! Make sure that fixed particles (inverse mass = 0) do not change
  118.   // their position! Bonus: The particle-dependent accelerations should also
  119.   // consider the wind direction and strength => Use your implementation of
  120.   // calcWindAcc()
  121.   //
  122.   // Parameters:
  123.   // p_next: updated particle positions
  124.   // p_cur: current paricle positions
  125.   // p_prev: previous particle positions
  126.   // normals: particle normals
  127.   // inv_masses: inverse particle masses
  128.   // damp: damping factor
  129.   // gravity: gravity vector
  130.   // wind_force: wind force
  131.   // dt: timestep
  132.   for(int itr = 0; itr < p_next.size(); itr++)
  133.     p_next[itr] = inv_masses[itr] == 0.0f ? p_cur[itr] : (1 + damp) * p_cur[itr] - damp * p_prev[itr] +
  134.                                                       (gravity + calcWindAcc(normals[itr],wind_force,inv_masses[itr])) * powf(dt,2.0f);  
  135.        
  136.    
  137.  
  138. }
  139.  
  140. void fixupStep(std::vector<float3>& p_out, const std::vector<float3>& p_in,
  141.                int grid_x, int grid_y, float cloth_x, float cloth_y,
  142.                const std::vector<float>& inv_masses, float fixup_percent,
  143.                const std::vector<Sphere>& obstacles,
  144.                bool apply_structural_fixup, bool apply_shear_fixup,
  145.                bool apply_flexion_fixup) {
  146.     // TODO: implement the fixup-step
  147.     //
  148.     // For each particle iterate over all particles connected via springs and
  149.     // compute the movement required to restore each spring into its resting
  150.     // position. Based on the mass ratio between the particle pairs compute the
  151.     // part of the movement that would have to be carried out by the current
  152.     // particle. Sum these movements up and update the particle position with a
  153.     // fraction of that movement(fixup_percent of the full movement). Iterate over
  154.     // all obstacles and resolve the collisions.
  155.     //
  156.     // *Only* apply the relevant fixup (structural/shear/flexion) if the
  157.     // corresponding flag is set (apply_*_fixup)!
  158.     //
  159.     // Parameters:
  160.     // p_out: new particle positions
  161.     // p_in: current particle positions
  162.     // grid_x, grid_y: particle-grid dimensions
  163.     // cloth_x, cloth_y: cloth size
  164.     // inv_masses: inverse particle masses
  165.     // fix_percent: fixup percent [0..1]
  166.     // obstacles: obstacle vector
  167.     // apply_structural_fixup: only structural fixup if this is true
  168.     // apply_shear_fixup: only apply shear fixup if this is true
  169.     // apply_flexion_fixup: only apply flexion fixup if this is true
  170.     float flex_x = (cloth_x/(grid_x - 1.0f)) * 2.0f;
  171.     float flex_y = (cloth_y/(grid_y - 1.0f)) * 2.0f;
  172.     float shear = std::sqrt(std::pow(cloth_x/(grid_x - 1.0f), 2) + std::pow(cloth_y/(grid_y - 1.0f), 2));
  173.  
  174.  
  175.  
  176.     for(int y = 0; y < grid_y; y++)
  177.     {
  178.         for(int x = 0; x < grid_x; x++)
  179.         {
  180.             p_out.at(flattenIndex(x,y,grid_x)) = p_in.at(flattenIndex(x,y,grid_x));
  181.             float offset = 0.0f;
  182.             float distance = 0.0f;
  183.             float mass = 0.0f;
  184.  
  185.  
  186.             if(inv_masses.at(flattenIndex(x,y,grid_x)) > 0.0f)
  187.             {
  188.               /////////apply shear fixup /////////////////////////////
  189.               ///////////////////////////////////////////////////////
  190.               ///////////////////////////////////////////////////////
  191.  
  192.  
  193.                 if(apply_shear_fixup)
  194.                 {
  195.                   ///reset variables
  196.                     float3 feder = {0.0f,0.0f,0.0f};
  197.                     float3 particle1 = {0.0f,0.0f,0.0f};
  198.                     float3 particle2 = {0.0f,0.0f,0.0f};
  199.                     float3 particle3 = {0.0f,0.0f,0.0f};
  200.                     float3 particle4 = {0.0f,0.0f,0.0f};
  201.                     float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  202.                     if(x > 0 && y > 0)
  203.                     {
  204.                         feder = p_in.at(flattenIndex(x - 1, y - 1, grid_x));
  205.                         mass = inv_masses.at(flattenIndex(x - 1, y - 1, grid_x));
  206.                         // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  207.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  208.                         offset = distance - shear;
  209.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  210.  
  211.  
  212.                         particle1 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  213.  
  214.                     }
  215.  //////////
  216.                     if(x < grid_x - 1 && y < grid_y - 1)
  217.                     {
  218.                         feder = p_in.at(flattenIndex(x + 1, y + 1, grid_x));
  219.                         mass = inv_masses.at(flattenIndex(x + 1, y + 1, grid_x));
  220.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  221.                         offset = distance - shear;
  222.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  223.                         particle2 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  224.  
  225.                     }
  226.                     if(x > 0 && y < grid_y - 1)
  227.                     {
  228.                         feder = p_in.at(flattenIndex(x - 1, y + 1, grid_x));
  229.                         mass = inv_masses.at(flattenIndex(x - 1, y + 1, grid_x));
  230.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  231.                         offset = distance - shear;
  232.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  233.                         particle3 =  offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  234.  
  235.                     }
  236.                     if(x < grid_x - 1 && y > 0)
  237.                     {
  238.                         feder = p_in.at(flattenIndex(x + 1, y - 1, grid_x));
  239.                         mass = inv_masses.at(flattenIndex(x + 1, y - 1, grid_x));
  240.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  241.                         offset = distance - shear;
  242.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  243.                         particle4 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  244.  
  245.                     }
  246.                     p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) +  ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
  247.  
  248.                 }
  249.  ////////////////////////////////////////////////////////////////////////
  250.  
  251.                 ///apply flexion fixup ////////////////////////////////
  252.                 //////////////////////////////////////////////////////////////////////
  253.  
  254.                 if(apply_flexion_fixup)
  255.                 {
  256.                   ////reset variables
  257.                     float3 feder = {0.0f,0.0f,0.0f};
  258.                     float3 particle1 = {0.0f,0.0f,0.0f};
  259.                     float3 particle2 = {0.0f,0.0f,0.0f};
  260.                     float3 particle3 = {0.0f,0.0f,0.0f};
  261.                     float3 particle4 = {0.0f,0.0f,0.0f};
  262.  
  263.                     // offset = distance - shear;
  264.                     float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  265.                     if( x - 1 > 0 && x > 0)
  266.                     {
  267.                         feder = p_in.at(flattenIndex(x - 2, y, grid_x));
  268.                         mass = inv_masses.at(flattenIndex(x - 2, y, grid_x));
  269.                         // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  270.                         distance =length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  271.                         offset = distance - flex_x;
  272.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  273.  
  274.  
  275.                         particle1 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  276.  
  277.                     }
  278.  
  279.                     if(x < grid_x - 1 && x + 1 < grid_x - 1)
  280.                     {
  281.                         feder = p_in.at(flattenIndex(x + 2, y, grid_x));
  282.                         mass = inv_masses.at(flattenIndex(x + 2, y, grid_x));
  283.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  284.                         offset = distance - flex_x;
  285.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  286.  
  287.                         particle2 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  288.  
  289.                     }
  290.                     if(y + 1 < grid_y - 1 && y < grid_y - 1)
  291.                     {
  292.                         feder = p_in.at(flattenIndex(x, y + 2, grid_x));
  293.                         mass = inv_masses.at(flattenIndex(x, y + 2, grid_x));
  294.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  295.                         offset = distance - flex_y;
  296.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  297.                         particle3 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  298.  
  299.                     }
  300.                     if(y - 1 > 0 && y > 0)
  301.                     {
  302.                         feder = p_in.at(flattenIndex(x, y - 2, grid_x));
  303.                         mass = inv_masses.at(flattenIndex(x , y - 2, grid_x));
  304.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  305.                         offset = distance - flex_y;
  306.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  307.                         particle4 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  308.  
  309.                     }
  310.                     p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) +  ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
  311.  
  312.                 }
  313.  
  314.                 //////////////////////////////////////////////////////////////////////////////
  315.                 //////////apply structural fixup /////////////////////////////////////////////
  316.                 /////////////////////////////////////////////////////////////////////////////
  317.  
  318.                 if(apply_structural_fixup)
  319.                 {
  320.                   ////reset variables
  321.                     float3 feder = {0.0f,0.0f,0.0f};
  322.                     float3 particle1 = {0.0f,0.0f,0.0f};
  323.                     float3 particle2 = {0.0f,0.0f,0.0f};
  324.                     float3 particle3 = {0.0f,0.0f,0.0f};
  325.                     float3 particle4 = {0.0f,0.0f,0.0f};
  326.                     // offset = distance - shear;
  327.                     float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  328.                     if(y < grid_y - 1)
  329.                     {
  330.                         feder = p_in.at(flattenIndex(x, y + 1, grid_x));
  331.                         mass = inv_masses.at(flattenIndex(x, y + 1, grid_x));
  332.                         // float beginning_mass = inv_masses.at(flattenIndex(x,y,grid_x));
  333.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  334.                         offset = distance - (cloth_y/(grid_y - 1.0f));
  335.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  336.  
  337.  
  338.                         particle1 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  339.  
  340.                     }
  341.  
  342.                     if(x < grid_x - 1)
  343.                     {
  344.                         feder = p_in.at(flattenIndex(x + 1, y, grid_x));
  345.                         mass = inv_masses.at(flattenIndex(x + 1, y, grid_x));
  346.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  347.                         offset = distance - (cloth_x/(grid_x - 1.0f));
  348.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  349.                         particle2 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  350.  
  351.                     }
  352.                     if(x > 0)
  353.                     {
  354.                         feder = p_in.at(flattenIndex(x - 1 ,y, grid_x));
  355.                         mass = inv_masses.at(flattenIndex(x - 1, y, grid_x));
  356.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  357.                         offset = distance - (cloth_x/(grid_x - 1.0f));
  358.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  359.                         particle3 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  360.  
  361.                     }
  362.                     if(y > 0)
  363.                     {
  364.                         feder = p_in.at(flattenIndex(x, y - 1, grid_x));
  365.                         mass = inv_masses.at(flattenIndex(x , y - 1, grid_x));
  366.                         distance = length((p_in.at(flattenIndex(x,y,grid_x)) - feder));
  367.                         offset = distance - (cloth_y/(grid_y - 1.0f));
  368.                         float mass_end = (beginning_mass / (mass + beginning_mass));
  369.                         particle4 = offset *  normalize(feder - p_in.at(flattenIndex(x,y,grid_x))) * mass_end;
  370.                     }
  371.                     p_out.at(flattenIndex(x,y,grid_x)) = p_out.at(flattenIndex(x,y,grid_x)) +  ((particle1 + particle2 + particle3 + particle4)*fixup_percent);
  372.  
  373.                 }
  374.  
  375.  
  376.  
  377.  
  378.             }
  379.             ////////////////////////////////////////////////
  380.             //////////////// fix collission ////////////////
  381.             ////////////////////////////////////////////////
  382.             for(unsigned int itr = 0; itr < obstacles.size(); ++itr)
  383.             {
  384.                 fixCollision(p_out.at(flattenIndex(x,y,grid_x)), obstacles.at(itr));
  385.             }
  386.         }
  387.     }
  388.  
  389. }
  390.  
  391. void fixCollision(float3& particle_position, const Sphere& sphere) {
  392.   // TODO: resolve particle/obstacle collision
  393.   //
  394.   // Check if the particle is inside the sphere.
  395.   // If so, move the particle out of the sphere (along the sphere normal).
  396.   //
  397.   // Parameters:
  398.   // particle_position: position of the particle
  399.   // sphere: spherical obstacle containg the position and radius
  400.   if(length(particle_position - sphere.position) <  sphere.radius)
  401.   {
  402.     float3 temp = normalize(particle_position - sphere.position);
  403.     particle_position = sphere.position + temp * sphere.radius;
  404.    
  405.   }
  406.  
  407. }
  408.  
  409. void simulationStep(std::array<std::vector<float3>, 3>& p_buffers, int& next,
  410.                     int& current, int& previous, int grid_x, int grid_y,
  411.                     float cloth_x, float cloth_y,
  412.                     const std::vector<float>& inv_masses,
  413.                     std::vector<float3>& normals, const float3& gravity,
  414.                     const float3& wind_force, float damping,
  415.                     float fixup_percent, int fixup_iterations,
  416.                     const std::vector<Sphere>& obstacles, float dt,
  417.                     bool apply_structural_fixup, bool apply_shear_fixup,
  418.                     bool apply_flexion_fixup) {
  419.   // TODO: implement the simulation step
  420.   //
  421.   // 1. Calculate particle normals  => calcNormals()
  422.   // 2. Perform Verlet Integration  => verletIntegration()
  423.   // 3. Perform fixup steps         => fixupStep()
  424.   //
  425.   // Verlet Integration:
  426.   // You will need three position buffers, since the verlet integration needs
  427.   // the current und previous positions. The integrated position shall be
  428.   // written to the 'next' position buffer. Make sure to update the next,
  429.   // current, previous indices! e.g. after the integration step: next ->
  430.   // current, current -> last and last -> next
  431.   //
  432.   // Fixup:
  433.   // Repeatedly execute the FixUp step after each integration
  434.   // ('fixup_iterations' times). Make sure to update your position buffer
  435.   // indices in the right way!
  436.   //
  437.   // Make sure that 'current' points to the right position buffer before you
  438.   // return!
  439.   //
  440.   // Parameters:
  441.   // p_buffers: 3 position buffers
  442.   // next: index of next position buffer
  443.   // current: index of current position buffer
  444.   // previous: index of previous position buffer
  445.   // grid_x, grid_y: particle-grid dimensions
  446.   // cloth_x, cloth_y: cloth size
  447.   // inv_masses: inverse particle masses
  448.   // normals: particle normals
  449.   // gravity: gravity vector
  450.   // wind_force: wind force
  451.   // damp: damping factor for verlet integration
  452.   // fixup_percent: fixup percent [0..1]
  453.   // fixup_iterations: how often the fixup step should be executed
  454.   // obstacles: obstacle vector
  455.   // dt: timestep
  456.   //
  457.   // apply_structural_fixup: apply structural fixup
  458.   // (needs to be passed unmodified to fixupStep for automatic testing)
  459.   //
  460.   // apply_shear_fixup: apply shear fixup
  461.   // (needs to be passed unmodified to fixupStep for automatic testing)
  462.   //
  463.   // apply_flexion_fixup: apply flexion fixup
  464.   // (needs to be passed unmodified to fixupStep for automatic testing)
  465.   //
  466.  
  467.   // Simply make cloth fall down. You will want to replace that ;)
  468.   calcNormals(normals, p_buffers[current], grid_x, grid_y);
  469.   verletIntegration(p_buffers[next], p_buffers[current], p_buffers[previous], normals, inv_masses, damping,
  470.                      gravity, wind_force, dt);
  471.  
  472.  int save = next;
  473.  next = previous;
  474.  previous = current;
  475.  current = save;
  476.  
  477.  
  478.   for(unsigned int iteration = 0; iteration < fixup_iterations; iteration++)
  479.   {
  480.     fixupStep(p_buffers[next], p_buffers[current], grid_x, grid_y, cloth_x, cloth_y, inv_masses, fixup_percent,
  481.               obstacles, apply_structural_fixup, apply_shear_fixup, apply_flexion_fixup);
  482.      int save = current;
  483.      current = next;
  484.      next = save;
  485.    
  486.   }
  487.   //for (auto& p : p_buffers[current]) p += gravity * dt;
  488. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top