Reimplemented the voxel texture type for displaying hair continuum grids.
authorLukas Tönne <lukas.toenne@gmail.com>
Fri, 31 Oct 2014 19:29:51 +0000 (20:29 +0100)
committerLukas Tönne <lukas.toenne@gmail.com>
Tue, 20 Jan 2015 08:30:05 +0000 (09:30 +0100)
Conflicts:
source/blender/physics/intern/BPH_mass_spring.cpp

source/blender/physics/BPH_mass_spring.h
source/blender/physics/intern/BPH_mass_spring.cpp
source/blender/physics/intern/hair_volume.c
source/blender/physics/intern/implicit.h
source/blender/render/CMakeLists.txt
source/blender/render/SConscript
source/blender/render/intern/source/voxeldata.c

index 166d0402d1c1ee1b1fc1ff8b93be792f0f193a9f..c189dea28a4987c974a10ea833d15afbb506f3a9 100644 (file)
@@ -33,7 +33,10 @@ extern "C" {
 #endif
 
 struct Implicit_Data;
+struct Object;
 struct ClothModifierData;
+struct ListBase;
+struct VoxelData;
 
 typedef enum eMassSpringSolverStatus {
        BPH_SOLVER_SUCCESS              = (1 << 0),
@@ -50,7 +53,7 @@ void BPH_cloth_solver_free(struct ClothModifierData *clmd);
 int BPH_cloth_solve(struct Object *ob, float frame, struct ClothModifierData *clmd, struct ListBase *effectors);
 void BKE_cloth_solver_set_positions(struct ClothModifierData *clmd);
 
-bool implicit_hair_volume_get_texture_data(struct Object *UNUSED(ob), struct ClothModifierData *clmd, struct ListBase *UNUSED(effectors), struct VoxelData *vd);
+bool BPH_cloth_solver_get_texture_data(struct Object *ob, struct ClothModifierData *clmd, struct VoxelData *vd);
 
 #ifdef __cplusplus
 }
index a894603b884ce0f2a095d07496289840d80de79e..7a2d85f140db20e289430b65619512395a753d45 100644 (file)
-/*
- * ***** BEGIN GPL LICENSE BLOCK *****
- *
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- *
- * The Original Code is Copyright (C) Blender Foundation
- * All rights reserved.
- *
- * The Original Code is: all of this file.
- *
- * Contributor(s): Lukas Toenne
- *
- * ***** END GPL LICENSE BLOCK *****
- */
-
-/** \file blender/blenkernel/intern/BPH_mass_spring.c
- *  \ingroup bph
- */
-
-extern "C" {
-#include "MEM_guardedalloc.h"
-
-#include "DNA_cloth_types.h"
-#include "DNA_scene_types.h"
-#include "DNA_object_force.h"
-#include "DNA_object_types.h"
-#include "DNA_meshdata_types.h"
-#include "DNA_modifier_types.h"
-
-#include "BLI_math.h"
-#include "BLI_linklist.h"
-#include "BLI_utildefines.h"
-
-#include "BKE_cloth.h"
-#include "BKE_collision.h"
-#include "BKE_effect.h"
-}
-
-#include "BPH_mass_spring.h"
-#include "implicit.h"
-
-/* Number of off-diagonal non-zero matrix blocks.
- * Basically there is one of these for each vertex-vertex interaction.
- */
-static int cloth_count_nondiag_blocks(Cloth *cloth)
-{
-       LinkNode *link;
-       int nondiag = 0;
-       
-       for (link = cloth->springs; link; link = link->next) {
-               ClothSpring *spring = (ClothSpring *)link->link;
-               switch (spring->type) {
-                       case CLOTH_SPRING_TYPE_BENDING_ANG:
-                               /* angular bending combines 3 vertices */
-                               nondiag += 3;
-                               break;
-                               
-                       default:
-                               /* all other springs depend on 2 vertices only */
-                               nondiag += 1;
-                               break;
-               }
-       }
-       
-       return nondiag;
-}
-
-int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd)
-{
-       Cloth *cloth = clmd->clothObject;
-       ClothVertex *verts = cloth->verts;
-       const float ZERO[3] = {0.0f, 0.0f, 0.0f};
-       Implicit_Data *id;
-       unsigned int i, nondiag;
-       
-       nondiag = cloth_count_nondiag_blocks(cloth);
-       cloth->implicit = id = BPH_mass_spring_solver_create(cloth->numverts, nondiag);
-       
-       for (i = 0; i < cloth->numverts; i++) {
-               BPH_mass_spring_set_vertex_mass(id, i, verts[i].mass);
-       }
-       
-       for (i = 0; i < cloth->numverts; i++) {
-               BPH_mass_spring_set_motion_state(id, i, verts[i].x, ZERO);
-       }
-       
-       return 1;
-}
-
-void BPH_cloth_solver_free(ClothModifierData *clmd)
-{
-       Cloth *cloth = clmd->clothObject;
-       
-       if (cloth->implicit) {
-               BPH_mass_spring_solver_free(cloth->implicit);
-               cloth->implicit = NULL;
-       }
-}
-
-void BKE_cloth_solver_set_positions(ClothModifierData *clmd)
-{
-       Cloth *cloth = clmd->clothObject;
-       ClothVertex *verts = cloth->verts;
-       unsigned int numverts = cloth->numverts, i;
-       ClothHairRoot *cloth_roots = clmd->roots;
-       Implicit_Data *id = cloth->implicit;
-       
-       for (i = 0; i < numverts; i++) {
-               ClothHairRoot *root = &cloth_roots[i];
-               
-               BPH_mass_spring_set_rest_transform(id, i, root->rot);
-               BPH_mass_spring_set_motion_state(id, i, verts[i].x, verts[i].v);
-       }
-}
-
-static bool collision_response(ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, float dt, float restitution, float r_impulse[3])
-{
-       Cloth *cloth = clmd->clothObject;
-       int index = collpair->ap1;
-       bool result = false;
-       
-       float v1[3], v2_old[3], v2_new[3], v_rel_old[3], v_rel_new[3];
-       float epsilon2 = BLI_bvhtree_getepsilon(collmd->bvhtree);
-
-       float margin_distance = collpair->distance - epsilon2;
-       float mag_v_rel;
-       
-       zero_v3(r_impulse);
-       
-       if (margin_distance > 0.0f)
-               return false; /* XXX tested before already? */
-       
-       /* only handle static collisions here */
-       if ( collpair->flag & COLLISION_IN_FUTURE )
-               return false;
-       
-       /* velocity */
-       copy_v3_v3(v1, cloth->verts[index].v);
-       collision_get_collider_velocity(v2_old, v2_new, collmd, collpair);
-       /* relative velocity = velocity of the cloth point relative to the collider */
-       sub_v3_v3v3(v_rel_old, v1, v2_old);
-       sub_v3_v3v3(v_rel_new, v1, v2_new);
-       /* normal component of the relative velocity */
-       mag_v_rel = dot_v3v3(v_rel_old, collpair->normal);
-       
-       /* only valid when moving toward the collider */
-       if (mag_v_rel < -ALMOST_ZERO) {
-               float v_nor_old, v_nor_new;
-               float v_tan_old[3], v_tan_new[3];
-               float bounce, repulse;
-               
-               /* Collision response based on
-                * "Simulating Complex Hair with Robust Collision Handling" (Choe, Choi, Ko, ACM SIGGRAPH 2005)
-                * http://graphics.snu.ac.kr/publications/2005-choe-HairSim/Choe_2005_SCA.pdf
-                */
-               
-               v_nor_old = mag_v_rel;
-               v_nor_new = dot_v3v3(v_rel_new, collpair->normal);
-               
-               madd_v3_v3v3fl(v_tan_old, v_rel_old, collpair->normal, -v_nor_old);
-               madd_v3_v3v3fl(v_tan_new, v_rel_new, collpair->normal, -v_nor_new);
-               
-               bounce = -v_nor_old * restitution;
-               
-               repulse = -margin_distance / dt; /* base repulsion velocity in normal direction */
-               /* XXX this clamping factor is quite arbitrary ...
-                * not sure if there is a more scientific approach, but seems to give good results
-                */
-               CLAMP(repulse, 0.0f, 4.0f * bounce);
-               
-               if (margin_distance < -epsilon2) {
-                       mul_v3_v3fl(r_impulse, collpair->normal, max_ff(repulse, bounce) - v_nor_new);
-               }
-               else {
-                       bounce = 0.0f;
-                       mul_v3_v3fl(r_impulse, collpair->normal, repulse - v_nor_new);
-               }
-               
-               result = true;
-       }
-       
-       return result;
-}
-
-/* Init constraint matrix
- * This is part of the modified CG method suggested by Baraff/Witkin in
- * "Large Steps in Cloth Simulation" (Siggraph 1998)
- */
-static void cloth_setup_constraints(ClothModifierData *clmd, ColliderContacts *contacts, int totcolliders, float dt)
-{
-       Cloth *cloth = clmd->clothObject;
-       Implicit_Data *data = cloth->implicit;
-       ClothVertex *verts = cloth->verts;
-       int numverts = cloth->numverts;
-       int i, j, v;
-       
-       const float ZERO[3] = {0.0f, 0.0f, 0.0f};
-       
-       for (v = 0; v < numverts; v++) {
-               if (verts[v].flags & CLOTH_VERT_FLAG_PINNED) {
-                       /* pinned vertex constraints */
-                       BPH_mass_spring_add_constraint_ndof0(data, v, ZERO); /* velocity is defined externally */
-               }
-               
-               verts[v].impulse_count = 0;
-       }
-
-       for (i = 0; i < totcolliders; ++i) {
-               ColliderContacts *ct = &contacts[i];
-               for (j = 0; j < ct->totcollisions; ++j) {
-                       CollPair *collpair = &ct->collisions[j];
-//                     float restitution = (1.0f - clmd->coll_parms->damping) * (1.0f - ct->ob->pd->pdef_sbdamp);
-                       float restitution = 0.0f;
-                       int v = collpair->face1;
-                       float impulse[3];
-                       
-                       /* pinned verts handled separately */
-                       if (verts[v].flags & CLOTH_VERT_FLAG_PINNED)
-                               continue;
-                       
-                       /* XXX cheap way of avoiding instability from multiple collisions in the same step
-                        * this should eventually be supported ...
-                        */
-                       if (verts[v].impulse_count > 0)
-                               continue;
-                       
-                       /* calculate collision response */
-                       if (!collision_response(clmd, ct->collmd, collpair, dt, restitution, impulse))
-                               continue;
-                       
-                       BPH_mass_spring_add_constraint_ndof2(data, v, collpair->normal, impulse);
-                       ++verts[v].impulse_count;
-                       
-                       BKE_sim_debug_data_add_dot(clmd->debug_data, collpair->pa, 0, 1, 0, "collision", hash_collpair(936, collpair));
-//                     BKE_sim_debug_data_add_dot(clmd->debug_data, collpair->pb, 1, 0, 0, "collision", hash_collpair(937, collpair));
-//                     BKE_sim_debug_data_add_line(clmd->debug_data, collpair->pa, collpair->pb, 0.7, 0.7, 0.7, "collision", hash_collpair(938, collpair));
-                       
-                       { /* DEBUG */
-                               float nor[3];
-                               mul_v3_v3fl(nor, collpair->normal, -collpair->distance);
-                               BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pa, nor, 1, 1, 0, "collision", hash_collpair(939, collpair));
-//                             BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, impulse, 1, 1, 0, "collision", hash_collpair(940, collpair));
-//                             BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, collpair->normal, 1, 1, 0, "collision", hash_collpair(941, collpair));
-                       }
-               }
-       }
-}
-
-/* computes where the cloth would be if it were subject to perfectly stiff edges
- * (edge distance constraints) in a lagrangian solver.  then add forces to help
- * guide the implicit solver to that state.  this function is called after
- * collisions*/
-static int UNUSED_FUNCTION(cloth_calc_helper_forces)(Object *UNUSED(ob), ClothModifierData *clmd, float (*initial_cos)[3], float UNUSED(step), float dt)
-{
-       Cloth *cloth= clmd->clothObject;
-       float (*cos)[3] = (float (*)[3])MEM_callocN(sizeof(float)*3*cloth->numverts, "cos cloth_calc_helper_forces");
-       float *masses = (float *)MEM_callocN(sizeof(float)*cloth->numverts, "cos cloth_calc_helper_forces");
-       LinkNode *node;
-       ClothSpring *spring;
-       ClothVertex *cv;
-       int i, steps;
-       
-       cv = cloth->verts;
-       for (i=0; i<cloth->numverts; i++, cv++) {
-               copy_v3_v3(cos[i], cv->tx);
-               
-               if (cv->goal == 1.0f || len_squared_v3v3(initial_cos[i], cv->tx) != 0.0f) {
-                       masses[i] = 1e+10;
-               }
-               else {
-                       masses[i] = cv->mass;
-               }
-       }
-       
-       steps = 55;
-       for (i=0; i<steps; i++) {
-               for (node=cloth->springs; node; node=node->next) {
-                       /* ClothVertex *cv1, *cv2; */ /* UNUSED */
-                       int v1, v2;
-                       float len, c, l, vec[3];
-                       
-                       spring = (ClothSpring *)node->link;
-                       if (spring->type != CLOTH_SPRING_TYPE_STRUCTURAL && spring->type != CLOTH_SPRING_TYPE_SHEAR) 
-                               continue;
-                       
-                       v1 = spring->ij; v2 = spring->kl;
-                       /* cv1 = cloth->verts + v1; */ /* UNUSED */
-                       /* cv2 = cloth->verts + v2; */ /* UNUSED */
-                       len = len_v3v3(cos[v1], cos[v2]);
-                       
-                       sub_v3_v3v3(vec, cos[v1], cos[v2]);
-                       normalize_v3(vec);
-                       
-                       c = (len - spring->restlen);
-                       if (c == 0.0f)
-                               continue;
-                       
-                       l = c / ((1.0f / masses[v1]) + (1.0f / masses[v2]));
-                       
-                       mul_v3_fl(vec, -(1.0f / masses[v1]) * l);
-                       add_v3_v3(cos[v1], vec);
-       
-                       sub_v3_v3v3(vec, cos[v2], cos[v1]);
-                       normalize_v3(vec);
-                       
-                       mul_v3_fl(vec, -(1.0f / masses[v2]) * l);
-                       add_v3_v3(cos[v2], vec);
-               }
-       }
-       
-       cv = cloth->verts;
-       for (i=0; i<cloth->numverts; i++, cv++) {
-               float vec[3];
-               
-               /*compute forces*/
-               sub_v3_v3v3(vec, cos[i], cv->tx);
-               mul_v3_fl(vec, cv->mass*dt*20.0f);
-               add_v3_v3(cv->tv, vec);
-               //copy_v3_v3(cv->tx, cos[i]);
-       }
-       
-       MEM_freeN(cos);
-       MEM_freeN(masses);
-       
-       return 1;
-}
-
-BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s, float time)
-{
-       Cloth *cloth = clmd->clothObject;
-       ClothSimSettings *parms = clmd->sim_parms;
-       Implicit_Data *data = cloth->implicit;
-       ClothVertex *verts = cloth->verts;
-       
-       bool no_compress = parms->flags & CLOTH_SIMSETTINGS_FLAG_NO_SPRING_COMPRESS;
-       
-       zero_v3(s->f);
-       zero_m3(s->dfdx);
-       zero_m3(s->dfdv);
-       
-       s->flags &= ~CLOTH_SPRING_FLAG_NEEDED;
-       
-       // calculate force of structural + shear springs
-       if ((s->type & CLOTH_SPRING_TYPE_STRUCTURAL) || (s->type & CLOTH_SPRING_TYPE_SHEAR) || (s->type & CLOTH_SPRING_TYPE_SEWING) ) {
-#ifdef CLOTH_FORCE_SPRING_STRUCTURAL
-               float k, scaling;
-               
-               s->flags |= CLOTH_SPRING_FLAG_NEEDED;
-               
-               scaling = parms->structural + s->stiffness * fabsf(parms->max_struct - parms->structural);
-               k = scaling / (parms->avg_spring_len + FLT_EPSILON);
-               
-               if (s->type & CLOTH_SPRING_TYPE_SEWING) {
-                       // TODO: verify, half verified (couldn't see error)
-                       // sewing springs usually have a large distance at first so clamp the force so we don't get tunnelling through colission objects
-                       BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, parms->max_sewing, s->f, s->dfdx, s->dfdv);
-               }
-               else {
-                       BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, 0.0f, s->f, s->dfdx, s->dfdv);
-               }
-#endif
-       }
-       else if (s->type & CLOTH_SPRING_TYPE_GOAL) {
-#ifdef CLOTH_FORCE_SPRING_GOAL
-               float goal_x[3], goal_v[3];
-               float k, scaling;
-               
-               s->flags |= CLOTH_SPRING_FLAG_NEEDED;
-               
-               // current_position = xold + t * (newposition - xold)
-               interp_v3_v3v3(goal_x, verts[s->ij].xold, verts[s->ij].xconst, time);
-               sub_v3_v3v3(goal_v, verts[s->ij].xconst, verts[s->ij].xold); // distance covered over dt==1
-               
-               scaling = parms->goalspring + s->stiffness * fabsf(parms->max_struct - parms->goalspring);
-               k = verts[s->ij].goal * scaling / (parms->avg_spring_len + FLT_EPSILON);
-               
-               BPH_mass_spring_force_spring_goal(data, s->ij, goal_x, goal_v, k, parms->goalfrict * 0.01f, s->f, s->dfdx, s->dfdv);
-#endif
-       }
-       else if (s->type & CLOTH_SPRING_TYPE_BENDING) {  /* calculate force of bending springs */
-#ifdef CLOTH_FORCE_SPRING_BEND
-               float kb, cb, scaling;
-               
-               s->flags |= CLOTH_SPRING_FLAG_NEEDED;
-               
-               scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);
-               kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
-               
-               scaling = parms->bending_damping;
-               cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
-               
-               BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);
-#endif
-       }
-       else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) {
-#ifdef CLOTH_FORCE_SPRING_BEND
-               float kb, cb, scaling;
-               
-               s->flags |= CLOTH_SPRING_FLAG_NEEDED;
-               
-               scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);
-               kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
-               
-               scaling = parms->bending_damping;
-               cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));
-               
-               /* XXX assuming same restlen for ij and jk segments here, this can be done correctly for hair later */
-               BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb);
-               
-               {
-                       float x_kl[3], x_mn[3], v[3], d[3];
-                       
-                       BPH_mass_spring_get_motion_state(data, s->kl, x_kl, v);
-                       BPH_mass_spring_get_motion_state(data, s->mn, x_mn, v);
-                       
-                       BKE_sim_debug_data_add_dot(clmd->debug_data, x_kl, 0.9, 0.9, 0.9, "target", hash_vertex(7980, s->kl));
-                       BKE_sim_debug_data_add_line(clmd->debug_data, x_kl, x_mn, 0.8, 0.8, 0.8, "target", hash_vertex(7981, s->kl));
-                       
-                       copy_v3_v3(d, s->target);
-                       BKE_sim_debug_data_add_vector(clmd->debug_data, x_kl, d, 0.8, 0.8, 0.2, "target", hash_vertex(7982, s->kl));
-                       
-//                     copy_v3_v3(d, s->target_ij);
-//                     BKE_sim_debug_data_add_vector(clmd->debug_data, x, d, 1, 0.4, 0.4, "target", hash_vertex(7983, s->kl));
-               }
-#endif
-       }
-}
-
-static void hair_get_boundbox(ClothModifierData *clmd, float gmin[3], float gmax[3])
-{
-       Cloth *cloth = clmd->clothObject;
-       Implicit_Data *data = cloth->implicit;
-       unsigned int numverts = cloth->numverts;
-       int i;
-       
-       INIT_MINMAX(gmin, gmax);
-       for (i = 0; i < numverts; i++) {
-               float x[3];
-               BPH_mass_spring_get_motion_state(data, i, x, NULL);
-               DO_MINMAX(x, gmin, gmax);
-       }
-}
-
-static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time)
-{
-       /* Collect forces and derivatives:  F, dFdX, dFdV */
-       Cloth *cloth = clmd->clothObject;
-       Implicit_Data *data = cloth->implicit;
-       unsigned int i  = 0;
-       float           drag    = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
-       float           gravity[3] = {0.0f, 0.0f, 0.0f};
-       MFace           *mfaces         = cloth->mfaces;
-       unsigned int numverts = cloth->numverts;
-       ClothVertex *vert;
-       
-#ifdef CLOTH_FORCE_GRAVITY
-       /* global acceleration (gravitation) */
-       if (clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
-               /* scale gravity force */
-               mul_v3_v3fl(gravity, clmd->scene->physics_settings.gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity);
-       }
-       vert = cloth->verts;
-       for (i = 0; i < cloth->numverts; i++, vert++) {
-               BPH_mass_spring_force_gravity(data, i, vert->mass, gravity);
-       }
-#endif
-
-       /* cloth_calc_volume_force(clmd); */
-
-#ifdef CLOTH_FORCE_DRAG
-       BPH_mass_spring_force_drag(data, drag);
-#endif
-       
-       /* handle external forces like wind */
-       if (effectors) {
-               /* cache per-vertex forces to avoid redundant calculation */
-               float (*winvec)[3] = (float (*)[3])MEM_callocN(sizeof(float) * 3 * numverts, "effector forces");
-               for (i = 0; i < cloth->numverts; i++) {
-                       float x[3], v[3];
-                       EffectedPoint epoint;
-                       
-                       BPH_mass_spring_get_motion_state(data, i, x, v);
-                       pd_point_from_loc(clmd->scene, x, v, i, &epoint);
-                       pdDoEffectors(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL);
-               }
-               
-               for (i = 0; i < cloth->numfaces; i++) {
-                       MFace *mf = &mfaces[i];
-                       BPH_mass_spring_force_face_wind(data, mf->v1, mf->v2, mf->v3, mf->v4, winvec);
-               }
-
-               /* Hair has only edges */
-               if (cloth->numfaces == 0) {
-                       for (LinkNode *link = cloth->springs; link; link = link->next) {
-                               ClothSpring *spring = (ClothSpring *)link->link;
-                               if (spring->type == CLOTH_SPRING_TYPE_STRUCTURAL)
-                                       BPH_mass_spring_force_edge_wind(data, spring->ij, spring->kl, winvec);
-                       }
-               }
-
-               MEM_freeN(winvec);
-       }
-       
-       // calculate spring forces
-       for (LinkNode *link = cloth->springs; link; link = link->next) {
-               ClothSpring *spring = (ClothSpring *)link->link;
-               // only handle active springs
-               if (!(spring->flags & CLOTH_SPRING_FLAG_DEACTIVATE))
-                       cloth_calc_spring_force(clmd, spring, time);
-       }
-}
-
-#if 0
-static void cloth_calc_volume_force(ClothModifierData *clmd)
-{
-       ClothSimSettings *parms = clmd->sim_parms;
-       Cloth *cloth = clmd->clothObject;
-       Implicit_Data *data = cloth->implicit;
-       int numverts = cloth->numverts;
-       ClothVertex *vert;
-       
-       /* 2.0f is an experimental value that seems to give good results */
-       float smoothfac = 2.0f * parms->velocity_smooth;
-       float collfac = 2.0f * parms->collider_friction;
-       float pressfac = parms->pressure;
-       float minpress = parms->pressure_threshold;
-       float gmin[3], gmax[3];
-       int i;
-       
-       hair_get_boundbox(clmd, gmin, gmax);
-       
-       /* gather velocities & density */
-       if (smoothfac > 0.0f || pressfac > 0.0f) {
-               HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);
-               
-               vert = cloth->verts;
-               for (i = 0; i < numverts; i++, vert++) {
-                       float x[3], v[3];
-                       
-                       if (vert->solver_index < 0) {
-                               copy_v3_v3(x, vert->x);
-                               copy_v3_v3(v, vert->v);
-                       }
-                       else {
-                               BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);
-                       }
-                       BPH_hair_volume_add_vertex(vertex_grid, x, v);
-               }
-               BPH_hair_volume_normalize_vertex_grid(vertex_grid);
-               
-#if 0
-               /* apply velocity filter */
-               BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);
-#endif
-               
-               vert = cloth->verts;
-               for (i = 0; i < numverts; i++, vert++) {
-                       float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3];
-                       
-                       if (vert->solver_index < 0)
-                               continue;
-                       
-                       /* calculate volumetric forces */
-                       BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);
-                       BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv);
-                       /* apply on hair data */
-                       BPH_mass_spring_force_extern(data, vert->solver_index, f, dfdx, dfdv);
-               }
-               
-               BPH_hair_volume_free_vertex_grid(vertex_grid);
-       }
-}
-#endif
-
-static void cloth_continuum_step(ClothModifierData *clmd)
-{
-       ClothSimSettings *parms = clmd->sim_parms;
-       Cloth *cloth = clmd->clothObject;
-       Implicit_Data *data = cloth->implicit;
-       int numverts = cloth->numverts;
-       ClothVertex *vert;
-       
-       const float fluid_factor = 0.95f; /* blend between PIC and FLIP methods */
-       float smoothfac = parms->velocity_smooth;
-       float pressfac = parms->pressure;
-       float minpress = parms->pressure_threshold;
-       float gmin[3], gmax[3];
-       int i;
-       
-       /* clear grid info */
-       zero_v3_int(clmd->hair_grid_res);
-       zero_v3(clmd->hair_grid_min);
-       zero_v3(clmd->hair_grid_max);
-       
-       hair_get_boundbox(clmd, gmin, gmax);
-       
-       /* gather velocities & density */
-       if (smoothfac > 0.0f || pressfac > 0.0f) {
-               HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);
-               
-               vert = cloth->verts;
-               for (i = 0; i < numverts; i++, vert++) {
-                       float x[3], v[3];
-                       
-                       BPH_mass_spring_get_motion_state(data, i, x, v);
-                               BPH_mass_spring_get_position(data, i, x);
-                               BPH_mass_spring_get_new_velocity(data, i, v);
-                       BPH_hair_volume_add_vertex(vertex_grid, x, v);
-               }
-               BPH_hair_volume_normalize_vertex_grid(vertex_grid);
-               
-#if 0
-               /* apply velocity filter */
-               BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);
-#endif
-               
-               vert = cloth->verts;
-               for (i = 0; i < numverts; i++, vert++) {
-                       float x[3], v[3], nv[3];
-                       
-                       /* calculate volumetric velocity influence */
-                       BPH_mass_spring_get_position(data, i, x);
-                       BPH_mass_spring_get_new_velocity(data, i, v);
-                       
-                       BPH_hair_volume_grid_velocity(vertex_grid, x, v, fluid_factor, nv);
-                       
-                       interp_v3_v3v3(nv, v, nv, smoothfac);
-                       
-                       /* apply on hair data */
-                       BPH_mass_spring_set_new_velocity(data, i, nv);
-               }
-               
-               /* store basic grid info in the modifier data */
-               BPH_hair_volume_grid_geometry(vertex_grid, NULL, clmd->hair_grid_res, clmd->hair_grid_min, clmd->hair_grid_max);
-               
-#if 0 /* DEBUG hair velocity vector field */
-               {
-                       const int size = 64;
-                       int i, j;
-                       float offset[3], a[3], b[3];
-                       const int axis = 0;
-                       const float shift = 0.45f;
-                       
-                       copy_v3_v3(offset, clmd->hair_grid_min);
-                       zero_v3(a);
-                       zero_v3(b);
-                       
-                       offset[axis] = interpf(clmd->hair_grid_max[axis], clmd->hair_grid_min[axis], shift);
-                       a[(axis+1) % 3] = clmd->hair_grid_max[(axis+1) % 3] - clmd->hair_grid_min[(axis+1) % 3];
-                       b[(axis+2) % 3] = clmd->hair_grid_max[(axis+2) % 3] - clmd->hair_grid_min[(axis+2) % 3];
-                       
-                       for (j = 0; j < size; ++j) {
-                               for (i = 0; i < size; ++i) {
-                                       float x[3], v[3], nv[3];
-                                       madd_v3_v3v3fl(x, offset, a, (float)i / (float)(size-1));
-                                       madd_v3_v3fl(x, b, (float)j / (float)(size-1));
-                                       zero_v3(v);
-                                       
-                                       BPH_hair_volume_grid_velocity(vertex_grid, x, v, 0.0f, nv);
-                                       
-                                       BKE_sim_debug_data_add_vector(clmd->debug_data, x, nv, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));
-                               }
-                       }
-               }
-#endif
-               
-               BPH_hair_volume_free_vertex_grid(vertex_grid);
-       }
-}
-
-static void cloth_clear_result(ClothModifierData *clmd)
-{
-       ClothSolverResult *sres = clmd->solver_result;
-       
-       sres->status = 0;
-       sres->max_error = sres->min_error = sres->avg_error = 0.0f;
-       sres->max_iterations = sres->min_iterations = 0;
-       sres->avg_iterations = 0.0f;
-}
-
-static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *result, int steps)
-{
-       ClothSolverResult *sres = clmd->solver_result;
-       
-       if (sres->status) { /* already initialized ? */
-               /* error only makes sense for successful iterations */
-               if (result->status == BPH_SOLVER_SUCCESS) {
-                       sres->min_error = min_ff(sres->min_error, result->error);
-                       sres->max_error = max_ff(sres->max_error, result->error);
-                       sres->avg_error += result->error / (float)steps;
-               }
-               
-               sres->min_iterations = min_ii(sres->min_iterations, result->iterations);
-               sres->max_iterations = max_ii(sres->max_iterations, result->iterations);
-               sres->avg_iterations += (float)result->iterations / (float)steps;
-       }
-       else {
-               /* error only makes sense for successful iterations */
-               if (result->status == BPH_SOLVER_SUCCESS) {
-                       sres->min_error = sres->max_error = result->error;
-                       sres->avg_error += result->error / (float)steps;
-               }
-               
-               sres->min_iterations = sres->max_iterations  = result->iterations;
-               sres->avg_iterations += (float)result->iterations / (float)steps;
-       }
-       
-       sres->status |= result->status;
-}
-
-int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors)
-{
-       unsigned int i=0;
-       float step=0.0f, tf=clmd->sim_parms->timescale;
-       Cloth *cloth = clmd->clothObject;
-       ClothVertex *verts = cloth->verts/*, *cv*/;
-       unsigned int numverts = cloth->numverts;
-       float dt = clmd->sim_parms->timescale / clmd->sim_parms->stepsPerFrame;
-       Implicit_Data *id = cloth->implicit;
-       ColliderContacts *contacts = NULL;
-       int totcolliders = 0;
-       
-       BPH_mass_spring_solver_debug_data(id, clmd->debug_data);
-       
-       BKE_sim_debug_data_clear_category(clmd->debug_data, "collision");
-       
-       if (!clmd->solver_result)
-               clmd->solver_result = (ClothSolverResult *)MEM_callocN(sizeof(ClothSolverResult), "cloth solver result");
-       cloth_clear_result(clmd);
-       
-       if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) { /* do goal stuff */
-               for (i = 0; i < numverts; i++) {
-                       // update velocities with constrained velocities from pinned verts
-                       if (verts [i].flags & CLOTH_VERT_FLAG_PINNED) {
-                               float v[3];
-                               sub_v3_v3v3(v, verts[i].xconst, verts[i].xold);
-                               // mul_v3_fl(v, clmd->sim_parms->stepsPerFrame);
-                               BPH_mass_spring_set_velocity(id, i, v);
-                       }
-               }
-       }
-       
-       if (clmd->debug_data) {
-               for (i = 0; i < numverts; i++) {
-//                     BKE_sim_debug_data_add_dot(clmd->debug_data, verts[i].x, 1.0f, 0.1f, 1.0f, "points", hash_vertex(583, i));
-               }
-       }
-       
-       while (step < tf) {
-               ImplicitSolverResult result;
-               
-               /* initialize forces to zero */
-               BPH_mass_spring_clear_forces(id);
-               BPH_mass_spring_clear_constraints(id);
-               
-               /* copy velocities for collision */
-               for (i = 0; i < numverts; i++) {
-                       BPH_mass_spring_get_motion_state(id, i, NULL, verts[i].tv);
-                       copy_v3_v3(verts[i].v, verts[i].tv);
-               }
-               
-               /* determine contact points */
-               if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED) {
-                       if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_POINTS) {
-                               cloth_find_point_contacts(ob, clmd, 0.0f, tf, &contacts, &totcolliders);
-                       }
-               }
-               
-               /* setup vertex constraints for pinned vertices and contacts */
-               cloth_setup_constraints(clmd, contacts, totcolliders, dt);
-               
-               // damping velocity for artistic reasons
-               // this is a bad way to do it, should be removed imo - lukas_t
-               if (clmd->sim_parms->vel_damping != 1.0f) {
-                       for (i = 0; i < numverts; i++) {
-                               float v[3];
-                               BPH_mass_spring_get_motion_state(id, i, NULL, v);
-                               mul_v3_fl(v, clmd->sim_parms->vel_damping);
-                               BPH_mass_spring_set_velocity(id, i, v);
-                       }
-               }
-               
-               // calculate forces
-               cloth_calc_force(clmd, frame, effectors, step);
-               
-               // calculate new velocity and position
-               BPH_mass_spring_solve_velocities(id, dt, &result);
-               cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);
-               
-               cloth_continuum_step(clmd);
-               
-               BPH_mass_spring_solve_positions(id, dt);
-               
-               BPH_mass_spring_apply_result(id);
-               
-               /* move pinned verts to correct position */
-               for (i = 0; i < numverts; i++) {
-                       if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) {
-                               if (verts[i].flags & CLOTH_VERT_FLAG_PINNED) {
-                                       float x[3];
-                                       interp_v3_v3v3(x, verts[i].xold, verts[i].xconst, step + dt);
-                                       BPH_mass_spring_set_position(id, i, x);
-                               }
-                       }
-                       
-                       BPH_mass_spring_get_motion_state(id, i, verts[i].txold, NULL);
-                       
-//                     if (!(verts[i].flags & CLOTH_VERT_FLAG_PINNED) && i > 0) {
-//                             BKE_sim_debug_data_add_line(clmd->debug_data, id->X[i], id->X[i-1], 0.6, 0.3, 0.3, "hair", hash_vertex(4892, i));
-//                             BKE_sim_debug_data_add_line(clmd->debug_data, id->Xnew[i], id->Xnew[i-1], 1, 0.5, 0.5, "hair", hash_vertex(4893, i));
-//                     }
-//                     BKE_sim_debug_data_add_vector(clmd->debug_data, id->X[i], id->V[i], 0, 0, 1, "velocity", hash_vertex(3158, i));
-               }
-               
-               /* free contact points */
-               if (contacts) {
-                       cloth_free_contacts(contacts, totcolliders);
-               }
-               
-               step += dt;
-       }
-       
-       /* copy results back to cloth data */
-       for (i = 0; i < numverts; i++) {
-               BPH_mass_spring_get_motion_state(id, i, verts[i].x, verts[i].v);
-               copy_v3_v3(verts[i].txold, verts[i].x);
-       }
-       
-       BPH_mass_spring_solver_debug_data(id, NULL);
-       
-       return 1;
-}
+/*\r
+ * ***** BEGIN GPL LICENSE BLOCK *****\r
+ *\r
+ * This program is free software; you can redistribute it and/or\r
+ * modify it under the terms of the GNU General Public License\r
+ * as published by the Free Software Foundation; either version 2\r
+ * of the License, or (at your option) any later version.\r
+ *\r
+ * This program is distributed in the hope that it will be useful,\r
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
+ * GNU General Public License for more details.\r
+ *\r
+ * You should have received a copy of the GNU General Public License\r
+ * along with this program; if not, write to the Free Software Foundation,\r
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.\r
+ *\r
+ * The Original Code is Copyright (C) Blender Foundation\r
+ * All rights reserved.\r
+ *\r
+ * The Original Code is: all of this file.\r
+ *\r
+ * Contributor(s): Lukas Toenne\r
+ *\r
+ * ***** END GPL LICENSE BLOCK *****\r
+ */\r
+\r
+/** \file blender/blenkernel/intern/BPH_mass_spring.c\r
+ *  \ingroup bph\r
+ */\r
+\r
+extern "C" {\r
+#include "MEM_guardedalloc.h"\r
+\r
+#include "DNA_cloth_types.h"\r
+#include "DNA_scene_types.h"\r
+#include "DNA_object_force.h"\r
+#include "DNA_object_types.h"\r
+#include "DNA_meshdata_types.h"\r
+#include "DNA_modifier_types.h"\r
+\r
+#include "BLI_math.h"\r
+#include "BLI_linklist.h"\r
+#include "BLI_utildefines.h"\r
+\r
+#include "BKE_cloth.h"\r
+#include "BKE_collision.h"\r
+#include "BKE_effect.h"\r
+}\r
+\r
+#include "BPH_mass_spring.h"\r
+#include "implicit.h"\r
+\r
+/* Number of off-diagonal non-zero matrix blocks.\r
+ * Basically there is one of these for each vertex-vertex interaction.\r
+ */\r
+static int cloth_count_nondiag_blocks(Cloth *cloth)\r
+{\r
+       LinkNode *link;\r
+       int nondiag = 0;\r
+       \r
+       for (link = cloth->springs; link; link = link->next) {\r
+               ClothSpring *spring = (ClothSpring *)link->link;\r
+               switch (spring->type) {\r
+                       case CLOTH_SPRING_TYPE_BENDING_ANG:\r
+                               /* angular bending combines 3 vertices */\r
+                               nondiag += 3;\r
+                               break;\r
+                               \r
+                       default:\r
+                               /* all other springs depend on 2 vertices only */\r
+                               nondiag += 1;\r
+                               break;\r
+               }\r
+       }\r
+       \r
+       return nondiag;\r
+}\r
+\r
+int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       ClothVertex *verts = cloth->verts;\r
+       const float ZERO[3] = {0.0f, 0.0f, 0.0f};\r
+       Implicit_Data *id;\r
+       unsigned int i, nondiag;\r
+       \r
+       nondiag = cloth_count_nondiag_blocks(cloth);\r
+       cloth->implicit = id = BPH_mass_spring_solver_create(cloth->numverts, nondiag);\r
+       \r
+       for (i = 0; i < cloth->numverts; i++) {\r
+               BPH_mass_spring_set_vertex_mass(id, i, verts[i].mass);\r
+       }\r
+       \r
+       for (i = 0; i < cloth->numverts; i++) {\r
+               BPH_mass_spring_set_motion_state(id, i, verts[i].x, ZERO);\r
+       }\r
+       \r
+       return 1;\r
+}\r
+\r
+void BPH_cloth_solver_free(ClothModifierData *clmd)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       \r
+       if (cloth->implicit) {\r
+               BPH_mass_spring_solver_free(cloth->implicit);\r
+               cloth->implicit = NULL;\r
+       }\r
+}\r
+\r
+void BKE_cloth_solver_set_positions(ClothModifierData *clmd)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       ClothVertex *verts = cloth->verts;\r
+       unsigned int numverts = cloth->numverts, i;\r
+       ClothHairRoot *cloth_roots = clmd->roots;\r
+       Implicit_Data *id = cloth->implicit;\r
+       \r
+       for (i = 0; i < numverts; i++) {\r
+               ClothHairRoot *root = &cloth_roots[i];\r
+               \r
+               BPH_mass_spring_set_rest_transform(id, i, root->rot);\r
+               BPH_mass_spring_set_motion_state(id, i, verts[i].x, verts[i].v);\r
+       }\r
+}\r
+\r
+static bool collision_response(ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, float dt, float restitution, float r_impulse[3])\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       int index = collpair->ap1;\r
+       bool result = false;\r
+       \r
+       float v1[3], v2_old[3], v2_new[3], v_rel_old[3], v_rel_new[3];\r
+       float epsilon2 = BLI_bvhtree_getepsilon(collmd->bvhtree);\r
+\r
+       float margin_distance = collpair->distance - epsilon2;\r
+       float mag_v_rel;\r
+       \r
+       zero_v3(r_impulse);\r
+       \r
+       if (margin_distance > 0.0f)\r
+               return false; /* XXX tested before already? */\r
+       \r
+       /* only handle static collisions here */\r
+       if ( collpair->flag & COLLISION_IN_FUTURE )\r
+               return false;\r
+       \r
+       /* velocity */\r
+       copy_v3_v3(v1, cloth->verts[index].v);\r
+       collision_get_collider_velocity(v2_old, v2_new, collmd, collpair);\r
+       /* relative velocity = velocity of the cloth point relative to the collider */\r
+       sub_v3_v3v3(v_rel_old, v1, v2_old);\r
+       sub_v3_v3v3(v_rel_new, v1, v2_new);\r
+       /* normal component of the relative velocity */\r
+       mag_v_rel = dot_v3v3(v_rel_old, collpair->normal);\r
+       \r
+       /* only valid when moving toward the collider */\r
+       if (mag_v_rel < -ALMOST_ZERO) {\r
+               float v_nor_old, v_nor_new;\r
+               float v_tan_old[3], v_tan_new[3];\r
+               float bounce, repulse;\r
+               \r
+               /* Collision response based on\r
+                * "Simulating Complex Hair with Robust Collision Handling" (Choe, Choi, Ko, ACM SIGGRAPH 2005)\r
+                * http://graphics.snu.ac.kr/publications/2005-choe-HairSim/Choe_2005_SCA.pdf\r
+                */\r
+               \r
+               v_nor_old = mag_v_rel;\r
+               v_nor_new = dot_v3v3(v_rel_new, collpair->normal);\r
+               \r
+               madd_v3_v3v3fl(v_tan_old, v_rel_old, collpair->normal, -v_nor_old);\r
+               madd_v3_v3v3fl(v_tan_new, v_rel_new, collpair->normal, -v_nor_new);\r
+               \r
+               bounce = -v_nor_old * restitution;\r
+               \r
+               repulse = -margin_distance / dt; /* base repulsion velocity in normal direction */\r
+               /* XXX this clamping factor is quite arbitrary ...\r
+                * not sure if there is a more scientific approach, but seems to give good results\r
+                */\r
+               CLAMP(repulse, 0.0f, 4.0f * bounce);\r
+               \r
+               if (margin_distance < -epsilon2) {\r
+                       mul_v3_v3fl(r_impulse, collpair->normal, max_ff(repulse, bounce) - v_nor_new);\r
+               }\r
+               else {\r
+                       bounce = 0.0f;\r
+                       mul_v3_v3fl(r_impulse, collpair->normal, repulse - v_nor_new);\r
+               }\r
+               \r
+               result = true;\r
+       }\r
+       \r
+       return result;\r
+}\r
+\r
+/* Init constraint matrix\r
+ * This is part of the modified CG method suggested by Baraff/Witkin in\r
+ * "Large Steps in Cloth Simulation" (Siggraph 1998)\r
+ */\r
+static void cloth_setup_constraints(ClothModifierData *clmd, ColliderContacts *contacts, int totcolliders, float dt)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       Implicit_Data *data = cloth->implicit;\r
+       ClothVertex *verts = cloth->verts;\r
+       int numverts = cloth->numverts;\r
+       int i, j, v;\r
+       \r
+       const float ZERO[3] = {0.0f, 0.0f, 0.0f};\r
+       \r
+       for (v = 0; v < numverts; v++) {\r
+               if (verts[v].flags & CLOTH_VERT_FLAG_PINNED) {\r
+                       /* pinned vertex constraints */\r
+                       BPH_mass_spring_add_constraint_ndof0(data, v, ZERO); /* velocity is defined externally */\r
+               }\r
+               \r
+               verts[v].impulse_count = 0;\r
+       }\r
+\r
+       for (i = 0; i < totcolliders; ++i) {\r
+               ColliderContacts *ct = &contacts[i];\r
+               for (j = 0; j < ct->totcollisions; ++j) {\r
+                       CollPair *collpair = &ct->collisions[j];\r
+//                     float restitution = (1.0f - clmd->coll_parms->damping) * (1.0f - ct->ob->pd->pdef_sbdamp);\r
+                       float restitution = 0.0f;\r
+                       int v = collpair->face1;\r
+                       float impulse[3];\r
+                       \r
+                       /* pinned verts handled separately */\r
+                       if (verts[v].flags & CLOTH_VERT_FLAG_PINNED)\r
+                               continue;\r
+                       \r
+                       /* XXX cheap way of avoiding instability from multiple collisions in the same step\r
+                        * this should eventually be supported ...\r
+                        */\r
+                       if (verts[v].impulse_count > 0)\r
+                               continue;\r
+                       \r
+                       /* calculate collision response */\r
+                       if (!collision_response(clmd, ct->collmd, collpair, dt, restitution, impulse))\r
+                               continue;\r
+                       \r
+                       BPH_mass_spring_add_constraint_ndof2(data, v, collpair->normal, impulse);\r
+                       ++verts[v].impulse_count;\r
+                       \r
+                       BKE_sim_debug_data_add_dot(clmd->debug_data, collpair->pa, 0, 1, 0, "collision", hash_collpair(936, collpair));\r
+//                     BKE_sim_debug_data_add_dot(clmd->debug_data, collpair->pb, 1, 0, 0, "collision", hash_collpair(937, collpair));\r
+//                     BKE_sim_debug_data_add_line(clmd->debug_data, collpair->pa, collpair->pb, 0.7, 0.7, 0.7, "collision", hash_collpair(938, collpair));\r
+                       \r
+                       { /* DEBUG */\r
+                               float nor[3];\r
+                               mul_v3_v3fl(nor, collpair->normal, -collpair->distance);\r
+                               BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pa, nor, 1, 1, 0, "collision", hash_collpair(939, collpair));\r
+//                             BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, impulse, 1, 1, 0, "collision", hash_collpair(940, collpair));\r
+//                             BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, collpair->normal, 1, 1, 0, "collision", hash_collpair(941, collpair));\r
+                       }\r
+               }\r
+       }\r
+}\r
+\r
+/* computes where the cloth would be if it were subject to perfectly stiff edges\r
+ * (edge distance constraints) in a lagrangian solver.  then add forces to help\r
+ * guide the implicit solver to that state.  this function is called after\r
+ * collisions*/\r
+static int UNUSED_FUNCTION(cloth_calc_helper_forces)(Object *UNUSED(ob), ClothModifierData *clmd, float (*initial_cos)[3], float UNUSED(step), float dt)\r
+{\r
+       Cloth *cloth= clmd->clothObject;\r
+       float (*cos)[3] = (float (*)[3])MEM_callocN(sizeof(float)*3*cloth->numverts, "cos cloth_calc_helper_forces");\r
+       float *masses = (float *)MEM_callocN(sizeof(float)*cloth->numverts, "cos cloth_calc_helper_forces");\r
+       LinkNode *node;\r
+       ClothSpring *spring;\r
+       ClothVertex *cv;\r
+       int i, steps;\r
+       \r
+       cv = cloth->verts;\r
+       for (i=0; i<cloth->numverts; i++, cv++) {\r
+               copy_v3_v3(cos[i], cv->tx);\r
+               \r
+               if (cv->goal == 1.0f || len_squared_v3v3(initial_cos[i], cv->tx) != 0.0f) {\r
+                       masses[i] = 1e+10;\r
+               }\r
+               else {\r
+                       masses[i] = cv->mass;\r
+               }\r
+       }\r
+       \r
+       steps = 55;\r
+       for (i=0; i<steps; i++) {\r
+               for (node=cloth->springs; node; node=node->next) {\r
+                       /* ClothVertex *cv1, *cv2; */ /* UNUSED */\r
+                       int v1, v2;\r
+                       float len, c, l, vec[3];\r
+                       \r
+                       spring = (ClothSpring *)node->link;\r
+                       if (spring->type != CLOTH_SPRING_TYPE_STRUCTURAL && spring->type != CLOTH_SPRING_TYPE_SHEAR) \r
+                               continue;\r
+                       \r
+                       v1 = spring->ij; v2 = spring->kl;\r
+                       /* cv1 = cloth->verts + v1; */ /* UNUSED */\r
+                       /* cv2 = cloth->verts + v2; */ /* UNUSED */\r
+                       len = len_v3v3(cos[v1], cos[v2]);\r
+                       \r
+                       sub_v3_v3v3(vec, cos[v1], cos[v2]);\r
+                       normalize_v3(vec);\r
+                       \r
+                       c = (len - spring->restlen);\r
+                       if (c == 0.0f)\r
+                               continue;\r
+                       \r
+                       l = c / ((1.0f / masses[v1]) + (1.0f / masses[v2]));\r
+                       \r
+                       mul_v3_fl(vec, -(1.0f / masses[v1]) * l);\r
+                       add_v3_v3(cos[v1], vec);\r
+       \r
+                       sub_v3_v3v3(vec, cos[v2], cos[v1]);\r
+                       normalize_v3(vec);\r
+                       \r
+                       mul_v3_fl(vec, -(1.0f / masses[v2]) * l);\r
+                       add_v3_v3(cos[v2], vec);\r
+               }\r
+       }\r
+       \r
+       cv = cloth->verts;\r
+       for (i=0; i<cloth->numverts; i++, cv++) {\r
+               float vec[3];\r
+               \r
+               /*compute forces*/\r
+               sub_v3_v3v3(vec, cos[i], cv->tx);\r
+               mul_v3_fl(vec, cv->mass*dt*20.0f);\r
+               add_v3_v3(cv->tv, vec);\r
+               //copy_v3_v3(cv->tx, cos[i]);\r
+       }\r
+       \r
+       MEM_freeN(cos);\r
+       MEM_freeN(masses);\r
+       \r
+       return 1;\r
+}\r
+\r
+BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s, float time)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       ClothSimSettings *parms = clmd->sim_parms;\r
+       Implicit_Data *data = cloth->implicit;\r
+       ClothVertex *verts = cloth->verts;\r
+       \r
+       bool no_compress = parms->flags & CLOTH_SIMSETTINGS_FLAG_NO_SPRING_COMPRESS;\r
+       \r
+       zero_v3(s->f);\r
+       zero_m3(s->dfdx);\r
+       zero_m3(s->dfdv);\r
+       \r
+       s->flags &= ~CLOTH_SPRING_FLAG_NEEDED;\r
+       \r
+       // calculate force of structural + shear springs\r
+       if ((s->type & CLOTH_SPRING_TYPE_STRUCTURAL) || (s->type & CLOTH_SPRING_TYPE_SHEAR) || (s->type & CLOTH_SPRING_TYPE_SEWING) ) {\r
+#ifdef CLOTH_FORCE_SPRING_STRUCTURAL\r
+               float k, scaling;\r
+               \r
+               s->flags |= CLOTH_SPRING_FLAG_NEEDED;\r
+               \r
+               scaling = parms->structural + s->stiffness * fabsf(parms->max_struct - parms->structural);\r
+               k = scaling / (parms->avg_spring_len + FLT_EPSILON);\r
+               \r
+               if (s->type & CLOTH_SPRING_TYPE_SEWING) {\r
+                       // TODO: verify, half verified (couldn't see error)\r
+                       // sewing springs usually have a large distance at first so clamp the force so we don't get tunnelling through colission objects\r
+                       BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, parms->max_sewing, s->f, s->dfdx, s->dfdv);\r
+               }\r
+               else {\r
+                       BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, 0.0f, s->f, s->dfdx, s->dfdv);\r
+               }\r
+#endif\r
+       }\r
+       else if (s->type & CLOTH_SPRING_TYPE_GOAL) {\r
+#ifdef CLOTH_FORCE_SPRING_GOAL\r
+               float goal_x[3], goal_v[3];\r
+               float k, scaling;\r
+               \r
+               s->flags |= CLOTH_SPRING_FLAG_NEEDED;\r
+               \r
+               // current_position = xold + t * (newposition - xold)\r
+               interp_v3_v3v3(goal_x, verts[s->ij].xold, verts[s->ij].xconst, time);\r
+               sub_v3_v3v3(goal_v, verts[s->ij].xconst, verts[s->ij].xold); // distance covered over dt==1\r
+               \r
+               scaling = parms->goalspring + s->stiffness * fabsf(parms->max_struct - parms->goalspring);\r
+               k = verts[s->ij].goal * scaling / (parms->avg_spring_len + FLT_EPSILON);\r
+               \r
+               BPH_mass_spring_force_spring_goal(data, s->ij, goal_x, goal_v, k, parms->goalfrict * 0.01f, s->f, s->dfdx, s->dfdv);\r
+#endif\r
+       }\r
+       else if (s->type & CLOTH_SPRING_TYPE_BENDING) {  /* calculate force of bending springs */\r
+#ifdef CLOTH_FORCE_SPRING_BEND\r
+               float kb, cb, scaling;\r
+               \r
+               s->flags |= CLOTH_SPRING_FLAG_NEEDED;\r
+               \r
+               scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);\r
+               kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));\r
+               \r
+               scaling = parms->bending_damping;\r
+               cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));\r
+               \r
+               BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->restlen, kb, cb, s->f, s->dfdx, s->dfdv);\r
+#endif\r
+       }\r
+       else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) {\r
+#ifdef CLOTH_FORCE_SPRING_BEND\r
+               float kb, cb, scaling;\r
+               \r
+               s->flags |= CLOTH_SPRING_FLAG_NEEDED;\r
+               \r
+               scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending);\r
+               kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));\r
+               \r
+               scaling = parms->bending_damping;\r
+               cb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON));\r
+               \r
+               /* XXX assuming same restlen for ij and jk segments here, this can be done correctly for hair later */\r
+               BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb);\r
+               \r
+               {\r
+                       float x_kl[3], x_mn[3], v[3], d[3];\r
+                       \r
+                       BPH_mass_spring_get_motion_state(data, s->kl, x_kl, v);\r
+                       BPH_mass_spring_get_motion_state(data, s->mn, x_mn, v);\r
+                       \r
+                       BKE_sim_debug_data_add_dot(clmd->debug_data, x_kl, 0.9, 0.9, 0.9, "target", hash_vertex(7980, s->kl));\r
+                       BKE_sim_debug_data_add_line(clmd->debug_data, x_kl, x_mn, 0.8, 0.8, 0.8, "target", hash_vertex(7981, s->kl));\r
+                       \r
+                       copy_v3_v3(d, s->target);\r
+                       BKE_sim_debug_data_add_vector(clmd->debug_data, x_kl, d, 0.8, 0.8, 0.2, "target", hash_vertex(7982, s->kl));\r
+                       \r
+//                     copy_v3_v3(d, s->target_ij);\r
+//                     BKE_sim_debug_data_add_vector(clmd->debug_data, x, d, 1, 0.4, 0.4, "target", hash_vertex(7983, s->kl));\r
+               }\r
+#endif\r
+       }\r
+}\r
+\r
+static void hair_get_boundbox(ClothModifierData *clmd, float gmin[3], float gmax[3])\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       Implicit_Data *data = cloth->implicit;\r
+       unsigned int numverts = cloth->numverts;\r
+       int i;\r
+       \r
+       INIT_MINMAX(gmin, gmax);\r
+       for (i = 0; i < numverts; i++) {\r
+               float x[3];\r
+               BPH_mass_spring_get_motion_state(data, i, x, NULL);\r
+               DO_MINMAX(x, gmin, gmax);\r
+       }\r
+}\r
+\r
+static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time)\r
+{\r
+       /* Collect forces and derivatives:  F, dFdX, dFdV */\r
+       Cloth *cloth = clmd->clothObject;\r
+       Implicit_Data *data = cloth->implicit;\r
+       unsigned int i  = 0;\r
+       float           drag    = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */\r
+       float           gravity[3] = {0.0f, 0.0f, 0.0f};\r
+       MFace           *mfaces         = cloth->mfaces;\r
+       unsigned int numverts = cloth->numverts;\r
+       ClothVertex *vert;\r
+       \r
+#ifdef CLOTH_FORCE_GRAVITY\r
+       /* global acceleration (gravitation) */\r
+       if (clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {\r
+               /* scale gravity force */\r
+               mul_v3_v3fl(gravity, clmd->scene->physics_settings.gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity);\r
+       }\r
+       vert = cloth->verts;\r
+       for (i = 0; i < cloth->numverts; i++, vert++) {\r
+               BPH_mass_spring_force_gravity(data, i, vert->mass, gravity);\r
+       }\r
+#endif\r
+\r
+       /* cloth_calc_volume_force(clmd); */\r
+\r
+#ifdef CLOTH_FORCE_DRAG\r
+       BPH_mass_spring_force_drag(data, drag);\r
+#endif\r
+       \r
+       /* handle external forces like wind */\r
+       if (effectors) {\r
+               /* cache per-vertex forces to avoid redundant calculation */\r
+               float (*winvec)[3] = (float (*)[3])MEM_callocN(sizeof(float) * 3 * numverts, "effector forces");\r
+               for (i = 0; i < cloth->numverts; i++) {\r
+                       float x[3], v[3];\r
+                       EffectedPoint epoint;\r
+                       \r
+                       BPH_mass_spring_get_motion_state(data, i, x, v);\r
+                       pd_point_from_loc(clmd->scene, x, v, i, &epoint);\r
+                       pdDoEffectors(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL);\r
+               }\r
+               \r
+               for (i = 0; i < cloth->numfaces; i++) {\r
+                       MFace *mf = &mfaces[i];\r
+                       BPH_mass_spring_force_face_wind(data, mf->v1, mf->v2, mf->v3, mf->v4, winvec);\r
+               }\r
+\r
+               /* Hair has only edges */\r
+               if (cloth->numfaces == 0) {\r
+                       for (LinkNode *link = cloth->springs; link; link = link->next) {\r
+                               ClothSpring *spring = (ClothSpring *)link->link;\r
+                               if (spring->type == CLOTH_SPRING_TYPE_STRUCTURAL)\r
+                                       BPH_mass_spring_force_edge_wind(data, spring->ij, spring->kl, winvec);\r
+                       }\r
+               }\r
+\r
+               MEM_freeN(winvec);\r
+       }\r
+       \r
+       // calculate spring forces\r
+       for (LinkNode *link = cloth->springs; link; link = link->next) {\r
+               ClothSpring *spring = (ClothSpring *)link->link;\r
+               // only handle active springs\r
+               if (!(spring->flags & CLOTH_SPRING_FLAG_DEACTIVATE))\r
+                       cloth_calc_spring_force(clmd, spring, time);\r
+       }\r
+}\r
+\r
+#if 0\r
+static void cloth_calc_volume_force(ClothModifierData *clmd)\r
+{\r
+       ClothSimSettings *parms = clmd->sim_parms;\r
+       Cloth *cloth = clmd->clothObject;\r
+       Implicit_Data *data = cloth->implicit;\r
+       int numverts = cloth->numverts;\r
+       ClothVertex *vert;\r
+       \r
+       /* 2.0f is an experimental value that seems to give good results */\r
+       float smoothfac = 2.0f * parms->velocity_smooth;\r
+       float collfac = 2.0f * parms->collider_friction;\r
+       float pressfac = parms->pressure;\r
+       float minpress = parms->pressure_threshold;\r
+       float gmin[3], gmax[3];\r
+       int i;\r
+       \r
+       hair_get_boundbox(clmd, gmin, gmax);\r
+       \r
+       /* gather velocities & density */\r
+       if (smoothfac > 0.0f || pressfac > 0.0f) {\r
+               HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);\r
+               \r
+               vert = cloth->verts;\r
+               for (i = 0; i < numverts; i++, vert++) {\r
+                       float x[3], v[3];\r
+                       \r
+                       if (vert->solver_index < 0) {\r
+                               copy_v3_v3(x, vert->x);\r
+                               copy_v3_v3(v, vert->v);\r
+                       }\r
+                       else {\r
+                               BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);\r
+                       }\r
+                       BPH_hair_volume_add_vertex(vertex_grid, x, v);\r
+               }\r
+               BPH_hair_volume_normalize_vertex_grid(vertex_grid);\r
+               \r
+#if 0\r
+               /* apply velocity filter */\r
+               BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);\r
+#endif\r
+               \r
+               vert = cloth->verts;\r
+               for (i = 0; i < numverts; i++, vert++) {\r
+                       float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3];\r
+                       \r
+                       if (vert->solver_index < 0)\r
+                               continue;\r
+                       \r
+                       /* calculate volumetric forces */\r
+                       BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v);\r
+                       BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv);\r
+                       /* apply on hair data */\r
+                       BPH_mass_spring_force_extern(data, vert->solver_index, f, dfdx, dfdv);\r
+               }\r
+               \r
+               BPH_hair_volume_free_vertex_grid(vertex_grid);\r
+       }\r
+}\r
+#endif\r
+\r
+static void cloth_continuum_fill_grid(HairVertexGrid *grid, Cloth *cloth)\r
+{\r
+       Implicit_Data *data = cloth->implicit;\r
+       int numverts = cloth->numverts;\r
+       ClothVertex *vert;\r
+       int i;\r
+       \r
+       for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {\r
+               float x[3], v[3];\r
+               \r
+               BPH_mass_spring_get_position(data, i, x);\r
+               BPH_mass_spring_get_new_velocity(data, i, v);\r
+               BPH_hair_volume_add_vertex(grid, x, v);\r
+       }\r
+       BPH_hair_volume_normalize_vertex_grid(grid);\r
+}\r
+\r
+static void cloth_continuum_step(ClothModifierData *clmd)\r
+{\r
+       ClothSimSettings *parms = clmd->sim_parms;\r
+       Cloth *cloth = clmd->clothObject;\r
+       Implicit_Data *data = cloth->implicit;\r
+       int numverts = cloth->numverts;\r
+       ClothVertex *vert;\r
+       \r
+       const float fluid_factor = 0.95f; /* blend between PIC and FLIP methods */\r
+       float smoothfac = parms->velocity_smooth;\r
+       float pressfac = parms->pressure;\r
+       float minpress = parms->pressure_threshold;\r
+       float gmin[3], gmax[3];\r
+       int i;\r
+       \r
+       /* clear grid info */\r
+       zero_v3_int(clmd->hair_grid_res);\r
+       zero_v3(clmd->hair_grid_min);\r
+       zero_v3(clmd->hair_grid_max);\r
+       \r
+       hair_get_boundbox(clmd, gmin, gmax);\r
+       \r
+       /* gather velocities & density */\r
+       if (smoothfac > 0.0f || pressfac > 0.0f) {\r
+               HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);\r
+               cloth_continuum_fill_grid(vertex_grid, cloth);\r
+               \r
+#if 0\r
+               /* apply velocity filter */\r
+               BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size);\r
+#endif\r
+               \r
+               for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {\r
+                       float x[3], v[3], nv[3];\r
+                       \r
+                       /* calculate volumetric velocity influence */\r
+                       BPH_mass_spring_get_position(data, i, x);\r
+                       BPH_mass_spring_get_new_velocity(data, i, v);\r
+                       \r
+                       BPH_hair_volume_grid_velocity(vertex_grid, x, v, fluid_factor, nv);\r
+                       \r
+                       interp_v3_v3v3(nv, v, nv, smoothfac);\r
+                       \r
+                       /* apply on hair data */\r
+                       BPH_mass_spring_set_new_velocity(data, i, nv);\r
+               }\r
+               \r
+               /* store basic grid info in the modifier data */\r
+               BPH_hair_volume_grid_geometry(vertex_grid, NULL, clmd->hair_grid_res, clmd->hair_grid_min, clmd->hair_grid_max);\r
+               \r
+#if 0 /* DEBUG hair velocity vector field */\r
+               {\r
+                       const int size = 64;\r
+                       int i, j;\r
+                       float offset[3], a[3], b[3];\r
+                       const int axis = 0;\r
+                       const float shift = 0.45f;\r
+                       \r
+                       copy_v3_v3(offset, clmd->hair_grid_min);\r
+                       zero_v3(a);\r
+                       zero_v3(b);\r
+                       \r
+                       offset[axis] = interpf(clmd->hair_grid_max[axis], clmd->hair_grid_min[axis], shift);\r
+                       a[(axis+1) % 3] = clmd->hair_grid_max[(axis+1) % 3] - clmd->hair_grid_min[(axis+1) % 3];\r
+                       b[(axis+2) % 3] = clmd->hair_grid_max[(axis+2) % 3] - clmd->hair_grid_min[(axis+2) % 3];\r
+                       \r
+                       for (j = 0; j < size; ++j) {\r
+                               for (i = 0; i < size; ++i) {\r
+                                       float x[3], v[3], nv[3];\r
+                                       madd_v3_v3v3fl(x, offset, a, (float)i / (float)(size-1));\r
+                                       madd_v3_v3fl(x, b, (float)j / (float)(size-1));\r
+                                       zero_v3(v);\r
+                                       \r
+                                       BPH_hair_volume_grid_velocity(vertex_grid, x, v, 0.0f, nv);\r
+                                       \r
+                                       BKE_sim_debug_data_add_vector(clmd->debug_data, x, nv, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));\r
+                               }\r
+                       }\r
+               }\r
+#endif\r
+               \r
+               BPH_hair_volume_free_vertex_grid(vertex_grid);\r
+       }\r
+}\r
+\r
+static void cloth_clear_result(ClothModifierData *clmd)\r
+{\r
+       ClothSolverResult *sres = clmd->solver_result;\r
+       \r
+       sres->status = 0;\r
+       sres->max_error = sres->min_error = sres->avg_error = 0.0f;\r
+       sres->max_iterations = sres->min_iterations = 0;\r
+       sres->avg_iterations = 0.0f;\r
+}\r
+\r
+static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *result, int steps)\r
+{\r
+       ClothSolverResult *sres = clmd->solver_result;\r
+       \r
+       if (sres->status) { /* already initialized ? */\r
+               /* error only makes sense for successful iterations */\r
+               if (result->status == BPH_SOLVER_SUCCESS) {\r
+                       sres->min_error = min_ff(sres->min_error, result->error);\r
+                       sres->max_error = max_ff(sres->max_error, result->error);\r
+                       sres->avg_error += result->error / (float)steps;\r
+               }\r
+               \r
+               sres->min_iterations = min_ii(sres->min_iterations, result->iterations);\r
+               sres->max_iterations = max_ii(sres->max_iterations, result->iterations);\r
+               sres->avg_iterations += (float)result->iterations / (float)steps;\r
+       }\r
+       else {\r
+               /* error only makes sense for successful iterations */\r
+               if (result->status == BPH_SOLVER_SUCCESS) {\r
+                       sres->min_error = sres->max_error = result->error;\r
+                       sres->avg_error += result->error / (float)steps;\r
+               }\r
+               \r
+               sres->min_iterations = sres->max_iterations  = result->iterations;\r
+               sres->avg_iterations += (float)result->iterations / (float)steps;\r
+       }\r
+       \r
+       sres->status |= result->status;\r
+}\r
+\r
+int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors)\r
+{\r
+       unsigned int i=0;\r
+       float step=0.0f, tf=clmd->sim_parms->timescale;\r
+       Cloth *cloth = clmd->clothObject;\r
+       ClothVertex *verts = cloth->verts/*, *cv*/;\r
+       unsigned int numverts = cloth->numverts;\r
+       float dt = clmd->sim_parms->timescale / clmd->sim_parms->stepsPerFrame;\r
+       Implicit_Data *id = cloth->implicit;\r
+       ColliderContacts *contacts = NULL;\r
+       int totcolliders = 0;\r
+       \r
+       BPH_mass_spring_solver_debug_data(id, clmd->debug_data);\r
+       \r
+       BKE_sim_debug_data_clear_category(clmd->debug_data, "collision");\r
+       \r
+       if (!clmd->solver_result)\r
+               clmd->solver_result = (ClothSolverResult *)MEM_callocN(sizeof(ClothSolverResult), "cloth solver result");\r
+       cloth_clear_result(clmd);\r
+       \r
+       if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) { /* do goal stuff */\r
+               for (i = 0; i < numverts; i++) {\r
+                       // update velocities with constrained velocities from pinned verts\r
+                       if (verts [i].flags & CLOTH_VERT_FLAG_PINNED) {\r
+                               float v[3];\r
+                               sub_v3_v3v3(v, verts[i].xconst, verts[i].xold);\r
+                               // mul_v3_fl(v, clmd->sim_parms->stepsPerFrame);\r
+                               BPH_mass_spring_set_velocity(id, i, v);\r
+                       }\r
+               }\r
+       }\r
+       \r
+       if (clmd->debug_data) {\r
+               for (i = 0; i < numverts; i++) {\r
+//                     BKE_sim_debug_data_add_dot(clmd->debug_data, verts[i].x, 1.0f, 0.1f, 1.0f, "points", hash_vertex(583, i));\r
+               }\r
+       }\r
+       \r
+       while (step < tf) {\r
+               ImplicitSolverResult result;\r
+               \r
+               /* initialize forces to zero */\r
+               BPH_mass_spring_clear_forces(id);\r
+               BPH_mass_spring_clear_constraints(id);\r
+               \r
+               /* copy velocities for collision */\r
+               for (i = 0; i < numverts; i++) {\r
+                       BPH_mass_spring_get_motion_state(id, i, NULL, verts[i].tv);\r
+                       copy_v3_v3(verts[i].v, verts[i].tv);\r
+               }\r
+               \r
+               /* determine contact points */\r
+               if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED) {\r
+                       if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_POINTS) {\r
+                               cloth_find_point_contacts(ob, clmd, 0.0f, tf, &contacts, &totcolliders);\r
+                       }\r
+               }\r
+               \r
+               /* setup vertex constraints for pinned vertices and contacts */\r
+               cloth_setup_constraints(clmd, contacts, totcolliders, dt);\r
+               \r
+               // damping velocity for artistic reasons\r
+               // this is a bad way to do it, should be removed imo - lukas_t\r
+               if (clmd->sim_parms->vel_damping != 1.0f) {\r
+                       for (i = 0; i < numverts; i++) {\r
+                               float v[3];\r
+                               BPH_mass_spring_get_motion_state(id, i, NULL, v);\r
+                               mul_v3_fl(v, clmd->sim_parms->vel_damping);\r
+                               BPH_mass_spring_set_velocity(id, i, v);\r
+                       }\r
+               }\r
+               \r
+               // calculate forces\r
+               cloth_calc_force(clmd, frame, effectors, step);\r
+               \r
+               // calculate new velocity and position\r
+               BPH_mass_spring_solve_velocities(id, dt, &result);\r
+               cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame);\r
+               \r
+               cloth_continuum_step(clmd);\r
+               \r
+               BPH_mass_spring_solve_positions(id, dt);\r
+               \r
+               BPH_mass_spring_apply_result(id);\r
+               \r
+               /* move pinned verts to correct position */\r
+               for (i = 0; i < numverts; i++) {\r
+                       if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) {\r
+                               if (verts[i].flags & CLOTH_VERT_FLAG_PINNED) {\r
+                                       float x[3];\r
+                                       interp_v3_v3v3(x, verts[i].xold, verts[i].xconst, step + dt);\r
+                                       BPH_mass_spring_set_position(id, i, x);\r
+                               }\r
+                       }\r
+                       \r
+                       BPH_mass_spring_get_motion_state(id, i, verts[i].txold, NULL);\r
+                       \r
+//                     if (!(verts[i].flags & CLOTH_VERT_FLAG_PINNED) && i > 0) {\r
+//                             BKE_sim_debug_data_add_line(clmd->debug_data, id->X[i], id->X[i-1], 0.6, 0.3, 0.3, "hair", hash_vertex(4892, i));\r
+//                             BKE_sim_debug_data_add_line(clmd->debug_data, id->Xnew[i], id->Xnew[i-1], 1, 0.5, 0.5, "hair", hash_vertex(4893, i));\r
+//                     }\r
+//                     BKE_sim_debug_data_add_vector(clmd->debug_data, id->X[i], id->V[i], 0, 0, 1, "velocity", hash_vertex(3158, i));\r
+               }\r
+               \r
+               /* free contact points */\r
+               if (contacts) {\r
+                       cloth_free_contacts(contacts, totcolliders);\r
+               }\r
+               \r
+               step += dt;\r
+       }\r
+       \r
+       /* copy results back to cloth data */\r
+       for (i = 0; i < numverts; i++) {\r
+               BPH_mass_spring_get_motion_state(id, i, verts[i].x, verts[i].v);\r
+               copy_v3_v3(verts[i].txold, verts[i].x);\r
+       }\r
+       \r
+       BPH_mass_spring_solver_debug_data(id, NULL);\r
+       \r
+       return 1;\r
+}\r
+\r
+bool BPH_cloth_solver_get_texture_data(Object *UNUSED(ob), ClothModifierData *clmd, VoxelData *vd)\r
+{\r
+       Cloth *cloth = clmd->clothObject;\r
+       HairVertexGrid *grid;\r
+       float gmin[3], gmax[3];\r
+       \r
+       if (!clmd->clothObject || !clmd->clothObject->implicit)\r
+               return false;\r
+       \r
+       hair_get_boundbox(clmd, gmin, gmax);\r
+       \r
+       grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);\r
+       cloth_continuum_fill_grid(grid, cloth);\r
+       \r
+       BPH_hair_volume_get_texture_data(grid, vd);\r
+       \r
+       BPH_hair_volume_free_vertex_grid(grid);\r
+       \r
+       return true;\r
+}\r
index b5c418802681d2633ad3eb96ab1a66355a3b03fd..988bf4b9024cf7cbf44a3346fd9b0e8adc4961d5 100644 (file)
@@ -34,6 +34,8 @@
 #include "BLI_math.h"
 #include "BLI_utildefines.h"
 
+#include "DNA_texture_types.h"
+
 #include "implicit.h"
 
 /* ================ Volumetric Hair Interaction ================
@@ -509,30 +511,16 @@ static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd,
 }
 #endif
 
-#if 0
-bool implicit_hair_volume_get_texture_data(Object *UNUSED(ob), ClothModifierData *clmd, ListBase *UNUSED(effectors), VoxelData *vd)
+bool BPH_hair_volume_get_texture_data(HairVertexGrid *grid, VoxelData *vd)
 {
-       lfVector *lX, *lV;
-       HairGridVert *hairgrid/*, *collgrid*/;
-       int numverts;
        int totres, i;
        int depth;
 
-       if (!clmd->clothObject || !clmd->clothObject->implicit)
-               return false;
-
-       lX = clmd->clothObject->implicit->X;
-       lV = clmd->clothObject->implicit->V;
-       numverts = clmd->clothObject->numverts;
-
-       hairgrid = hair_volume_create_hair_grid(clmd, lX, lV, numverts);
-//     collgrid = hair_volume_create_collision_grid(clmd, lX, numverts);
-
-       vd->resol[0] = hair_grid_res;
-       vd->resol[1] = hair_grid_res;
-       vd->resol[2] = hair_grid_res;
+       vd->resol[0] = grid->res;
+       vd->resol[1] = grid->res;
+       vd->resol[2] = grid->res;
        
-       totres = hair_grid_size(hair_grid_res);
+       totres = hair_grid_size(grid->res);
        
        if (vd->hair_type == TEX_VD_HAIRVELOCITY) {
                depth = 4;
@@ -549,20 +537,21 @@ bool implicit_hair_volume_get_texture_data(Object *UNUSED(ob), ClothModifierData
                for (i = 0; i < totres; ++i) {
                        switch (vd->hair_type) {
                                case TEX_VD_HAIRDENSITY:
-                                       vd->dataset[i] = hairgrid[i].density;
+                                       vd->dataset[i] = grid->verts[i].density;
                                        break;
                                
                                case TEX_VD_HAIRRESTDENSITY:
                                        vd->dataset[i] = 0.0f; // TODO
                                        break;
                                
-                               case TEX_VD_HAIRVELOCITY:
-                                       vd->dataset[i + 0*totres] = hairgrid[i].velocity[0];
-                                       vd->dataset[i + 1*totres] = hairgrid[i].velocity[1];
-                                       vd->dataset[i + 2*totres] = hairgrid[i].velocity[2];
-                                       vd->dataset[i + 3*totres] = len_v3(hairgrid[i].velocity);
+                               case TEX_VD_HAIRVELOCITY: {
+                                       float tmp = grid->verts[i].velocity[1];
+                                       vd->dataset[i + 0*totres] = grid->verts[i].velocity[0];
+                                       vd->dataset[i + 1*totres] = grid->verts[i].velocity[1];
+                                       vd->dataset[i + 2*totres] = grid->verts[i].velocity[2];
+                                       vd->dataset[i + 3*totres] = len_v3(grid->verts[i].velocity);
                                        break;
-                               
+                               }
                                case TEX_VD_HAIRENERGY:
                                        vd->dataset[i] = 0.0f; // TODO
                                        break;
@@ -573,10 +562,5 @@ bool implicit_hair_volume_get_texture_data(Object *UNUSED(ob), ClothModifierData
                vd->dataset = NULL;
        }
        
-       MEM_freeN(hairgrid);
-//     MEM_freeN(collgrid);
-       
        return true;
 }
-
-#endif
index 9ec7e763415e3dd31aea6d0f4e7056b6d4e166ce..b7eeea146bc7ccc2fe1eeb81f02725244f4c524b 100644 (file)
@@ -168,6 +168,9 @@ bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, const
 struct HairVertexGrid;
 struct HairColliderGrid;
 
+struct Object;
+struct VoxelData;
+
 struct HairVertexGrid *BPH_hair_volume_create_vertex_grid(int res, const float gmin[3], const float gmax[3]);
 void BPH_hair_volume_free_vertex_grid(struct HairVertexGrid *grid);
 void BPH_hair_volume_grid_geometry(struct HairVertexGrid *grid, float cellsize[3], int res[3], float gmin[3], float gmax[3]);
@@ -195,6 +198,8 @@ void BPH_hair_volume_vertex_grid_forces(struct HairVertexGrid *grid, const float
                                         float smoothfac, float pressurefac, float minpressure,
                                         float f[3], float dfdx[3][3], float dfdv[3][3]);
 
+bool BPH_hair_volume_get_texture_data(struct HairVertexGrid *grid, struct VoxelData *vd);
+
 #ifdef __cplusplus
 }
 #endif
index e516c95473703cb27492ce2b326a3fbf529efaa4..36b9f8ae362527ef36b2ba45c432702e72964157 100644 (file)
@@ -33,6 +33,7 @@ set(INC
        ../imbuf
        ../makesdna
        ../makesrna
+       ../physics
        ../../../intern/guardedalloc
        ../../../intern/mikktspace
        ../../../intern/smoke/extern
index ebd131082e4591cda251f3d33f7548ff18461514..0cbf188525ba2b93cf50636589142403bb3f8517 100644 (file)
@@ -40,6 +40,7 @@ incs = [
     '../imbuf',
     '../makesdna',
     '../makesrna',
+    '../physics',
     '../../../intern/mikktspace',
     '../../../intern/smoke/extern',
     ]
index 6f3c50bca62a61cdf134361d5d4ca10498dfd358..8bc5c7a1c8f0969fd2bba49f9ac0fdddabab6a87 100644 (file)
@@ -58,6 +58,7 @@
 #include "BKE_modifier.h"
 
 #include "smoke_API.h"
+#include "BPH_mass_spring.h"
 
 #include "DNA_texture_types.h"
 #include "DNA_object_force.h"
@@ -380,8 +381,7 @@ static void init_frame_hair(VoxelData *vd, int UNUSED(cfra))
                ParticleSystemModifierData *pmd = (ParticleSystemModifierData *)md;
                
                if (pmd->psys && pmd->psys->clmd) {
-                       // XXX TODO was moved into own subfolder, figure out how to handle this (perhaps make a wrapper in BKE)
-//                     found |= implicit_hair_volume_get_texture_data(ob, pmd->psys->clmd, NULL, vd);
+                       found |= BPH_cloth_solver_get_texture_data(ob, pmd->psys->clmd, vd);
                }
        }