Completed the implementation of bent rest shapes for hair.
authorLukas Tönne <lukas.toenne@gmail.com>
Thu, 25 Sep 2014 21:19:20 +0000 (23:19 +0200)
committerLukas Tönne <lukas.toenne@gmail.com>
Tue, 20 Jan 2015 08:30:03 +0000 (09:30 +0100)
Basically follows the Pixar approach from "Artistic Simulation of Curly
Hair".

source/blender/blenkernel/BKE_cloth.h
source/blender/blenkernel/intern/cloth.c
source/blender/blenkernel/intern/particle_system.c
source/blender/physics/intern/BPH_mass_spring.cpp
source/blender/physics/intern/implicit.h
source/blender/physics/intern/implicit_blender.c

index 8a2477aaa15d1c282dcaafdb74ce05c8b5497518..2a1741d1941e84e375a1b38cedd840d1c2f8a25f 100644 (file)
@@ -149,9 +149,6 @@ typedef struct ClothSpring {
        
        /* angular bending spring target and derivatives */
        float target[3];
-       float dtarget_dxij[3][3];
-       float dtarget_dxkl[3][3];
-       float dtarget_dxmn[3][3];
 }
 ClothSpring;
 
@@ -251,7 +248,7 @@ void cloth_clear_cache (struct Object *ob, struct ClothModifierData *clmd, float
 // needed for cloth.c
 int cloth_add_spring (struct ClothModifierData *clmd, unsigned int indexA, unsigned int indexB, float restlength, int spring_type);
 
-void cloth_parallel_transport_hair_frame(float mat[3][3], float dir_old[3], const float x_cur[3], const float x_new[3]);
+void cloth_parallel_transport_hair_frame(float mat[3][3], const float dir_old[3], const float dir_new[3]);
 
 ////////////////////////////////////////////////
 
index a6df9942158df2c272386d942ffde070cd1e6902..7768149bc1230c89ddcf17d99b161e637944a3b7 100644 (file)
@@ -1106,7 +1106,7 @@ static void cloth_update_bending_targets(ClothModifierData *clmd)
        Cloth *cloth = clmd->clothObject;
        ClothSpring *spring;
        LinkNode *search = NULL;
-       float hair_frame[3][3], dir[3];
+       float hair_frame[3][3], dir_old[3], dir_new[3];
        bool is_root;
        
        /* XXX Note: we need to propagate frames from the root up,
@@ -1119,7 +1119,7 @@ static void cloth_update_bending_targets(ClothModifierData *clmd)
        
        is_root = true;
        for (search = cloth->springs; search; search = search->next) {
-               ClothHairRoot *hair_info;
+               ClothHairRoot *hair_ij, *hair_kl;
                
                spring = search->link;
                if (spring->type != CLOTH_SPRING_TYPE_BENDING_ANG) {
@@ -1127,20 +1127,23 @@ static void cloth_update_bending_targets(ClothModifierData *clmd)
                        continue;
                }
                
-               hair_info = &clmd->roots[spring->kl];
+               hair_ij = &clmd->roots[spring->ij];
+               hair_kl = &clmd->roots[spring->kl];
                if (is_root) {
                        /* initial hair frame from root orientation */
-                       copy_m3_m3(hair_frame, hair_info->rot);
+                       copy_m3_m3(hair_frame, hair_ij->rot);
                        /* surface normal is the initial direction,
                         * parallel transport then keeps it aligned to the hair direction
                         */
-                       copy_v3_v3(dir, hair_frame[2]);
+                       copy_v3_v3(dir_new, hair_frame[2]);
                }
                
-               /* move frame to next hair segment */
-               cloth_parallel_transport_hair_frame(hair_frame, dir, cloth->verts[spring->kl].x, cloth->verts[spring->mn].x);
+               copy_v3_v3(dir_old, dir_new);
+               sub_v3_v3v3(dir_new, cloth->verts[spring->mn].x, cloth->verts[spring->kl].x);
+               normalize_v3(dir_new);
                
-               if (clmd->debug_data) {
+#if 1
+               if (clmd->debug_data && (spring->ij == 0 || spring->ij == 1)) {
                        float a[3], b[3];
                        
                        copy_v3_v3(a, cloth->verts[spring->kl].x);
@@ -1155,13 +1158,69 @@ static void cloth_update_bending_targets(ClothModifierData *clmd)
                        mul_v3_v3fl(b, hair_frame[2], clmd->sim_parms->avg_spring_len);
                        BKE_sim_debug_data_add_vector(clmd->debug_data, a, b, 0, 0, 1, "frames", hash_vertex(8249, hash_int_2d(spring->kl, spring->mn)));
                }
+#endif
+               
+               /* get local targets for kl/mn vertices by putting rest targets into the current frame,
+                * then multiply with the rest length to get the actual goals
+                */
+               
+               mul_v3_m3v3(spring->target, hair_frame, hair_kl->rest_target);
+               mul_v3_fl(spring->target, spring->restlen);
+               
+               /* move frame to next hair segment */
+               cloth_parallel_transport_hair_frame(hair_frame, dir_old, dir_new);
+               
+               is_root = false; /* next bending spring not connected to root */
+       }
+}
+
+static void cloth_update_bending_rest_targets(ClothModifierData *clmd)
+{
+       Cloth *cloth = clmd->clothObject;
+       ClothSpring *spring;
+       LinkNode *search = NULL;
+       float hair_frame[3][3], dir_old[3], dir_new[3];
+       bool is_root;
+       
+       /* XXX Note: we need to propagate frames from the root up,
+        * but structural hair springs are stored in reverse order.
+        * The bending springs however are then inserted in the same
+        * order as vertices again ...
+        * This messy situation can be resolved when solver data is
+        * generated directly from a dedicated hair system.
+        */
+       
+       is_root = true;
+       for (search = cloth->springs; search; search = search->next) {
+               ClothHairRoot *hair_ij, *hair_kl;
+               
+               spring = search->link;
+               if (spring->type != CLOTH_SPRING_TYPE_BENDING_ANG) {
+                       is_root = true; /* next bending spring connects to root */
+                       continue;
+               }
                
-               /* get target direction by putting rest target into the current frame */
-               mul_v3_m3v3(spring->target, hair_frame, hair_info->rest_target);
-               /* XXX TODO */
-               zero_m3(spring->dtarget_dxij);
-               zero_m3(spring->dtarget_dxkl);
-               zero_m3(spring->dtarget_dxmn);
+               hair_ij = &clmd->roots[spring->ij];
+               hair_kl = &clmd->roots[spring->kl];
+               if (is_root) {
+                       /* initial hair frame from root orientation */
+                       copy_m3_m3(hair_frame, hair_ij->rot);
+                       /* surface normal is the initial direction,
+                        * parallel transport then keeps it aligned to the hair direction
+                        */
+                       copy_v3_v3(dir_new, hair_frame[2]);
+               }
+               
+               copy_v3_v3(dir_old, dir_new);
+               sub_v3_v3v3(dir_new, cloth->verts[spring->mn].xrest, cloth->verts[spring->kl].xrest);
+               normalize_v3(dir_new);
+               
+               /* dir expressed in the hair frame defines the rest target direction */
+               copy_v3_v3(hair_kl->rest_target, dir_new);
+               mul_transposed_m3_v3(hair_frame, hair_kl->rest_target);
+               
+               /* move frame to next hair segment */
+               cloth_parallel_transport_hair_frame(hair_frame, dir_old, dir_new);
                
                is_root = false; /* next bending spring not connected to root */
        }
@@ -1232,23 +1291,15 @@ BLI_INLINE void madd_m3_m3fl(float r[3][3], float m[3][3], float f)
        r[2][2] += m[2][2] * f;
 }
 
-void cloth_parallel_transport_hair_frame(float mat[3][3], float dir_old[3], const float x_cur[3], const float x_new[3])
+void cloth_parallel_transport_hair_frame(float mat[3][3], const float dir_old[3], const float dir_new[3])
 {
-       float dir_new[3];
        float rot[3][3];
        
-       /* next segment direction */
-       sub_v3_v3v3(dir_new, x_new, x_cur);
-       normalize_v3(dir_new);
-       
        /* rotation between segments */
        rotation_between_vecs_to_mat3(rot, dir_old, dir_new);
        
        /* rotate the frame */
        mul_m3_m3m3(mat, rot, mat);
-       
-       /* advance old variables */
-       copy_v3_v3(dir_old, dir_new);
 }
 
 static int cloth_build_springs ( ClothModifierData *clmd, DerivedMesh *dm )
@@ -1491,6 +1542,8 @@ static int cloth_build_springs ( ClothModifierData *clmd, DerivedMesh *dm )
                                search2 = search2->next;
                        }
                }
+               
+               cloth_update_bending_rest_targets(clmd);
        }
        
        /* note: the edges may already exist so run reinsert */
index 661e0c5c1be876fe99f9fb3f3c84b24d2907668f..a840b7400680b622f5de2b8df6235a3530c27c4d 100644 (file)
@@ -4083,7 +4083,7 @@ static void do_hair_dynamics(ParticleSimulationData *sim)
        /* make vgroup for pin roots etc.. */
        psys->particles->hair_index = 1;
        LOOP_PARTICLES {
-               float root_mat[4][4], hair_frame[3][3], dir[3];
+               float root_mat[4][4];
                bool use_hair = psys_hair_use_simulation(pa, max_length);
 
                if (p)
@@ -4093,13 +4093,6 @@ static void do_hair_dynamics(ParticleSimulationData *sim)
                mul_m4_m4m4(root_mat, sim->ob->obmat, hairmat);
                normalize_m4(root_mat);
 
-               /* initial hair frame from root orientation */
-               copy_m3_m4(hair_frame, root_mat);
-               /* surface normal is the initial direction,
-                * parallel transport then keeps it aligned to the hair direction
-                */
-               copy_v3_v3(dir, hair_frame[2]);
-
                for (k=0, key=pa->hair; k<pa->totkey; k++,key++) {
                        ClothHairRoot *root;
                        
@@ -4111,10 +4104,6 @@ static void do_hair_dynamics(ParticleSimulationData *sim)
                                copy_v3_v3(root->loc, root_mat[3]);
                                copy_m3_m4(root->rot, root_mat);
                                
-                               /* dir expressed in the hair frame defines the rest target direction */
-                               copy_v3_v3(root->rest_target, dir);
-                               mul_transposed_m3_v3(hair_frame, root->rest_target);
-                               
                                sub_v3_v3v3(temp, key->co, (key+1)->co);
                                copy_v3_v3(mvert->co, key->co);
                                add_v3_v3v3(mvert->co, mvert->co, temp);
@@ -4133,14 +4122,6 @@ static void do_hair_dynamics(ParticleSimulationData *sim)
                        copy_v3_v3(root->loc, root_mat[3]);
                        copy_m3_m4(root->rot, root_mat);
 
-                       if (k < pa->totkey-1)
-                               /* move frame to next hair segment */
-                               cloth_parallel_transport_hair_frame(hair_frame, dir, key->co, (key+1)->co);
-                       
-                       /* dir expressed in the hair frame defines the rest target direction */
-                       copy_v3_v3(root->rest_target, dir);
-                       mul_transposed_m3_v3(hair_frame, root->rest_target);
-
                        copy_v3_v3(mvert->co, key->co);
                        mul_m4_v3(hairmat, mvert->co);
                        mvert++;
index 458a8f27ee7531c95ade425844194e594dfc24be..89e9447f345b572c792dbece4084ad7a3e2e0920 100644 (file)
@@ -425,13 +425,18 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
                cb = kb = 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->matrix_ij_kl, s->matrix_kl_mn, s->matrix_ij_mn, s->restlen, s->restlen, kb, cb);
+               BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->matrix_ij_kl, s->matrix_kl_mn, s->matrix_ij_mn, s->target, kb, cb);
                
                {
                        float x[3], v[3], d[3];
+                       
                        BPH_mass_spring_get_motion_state(data, s->kl, x, v);
-                       mul_v3_v3fl(d, s->target, clmd->sim_parms->avg_spring_len);
-                       BKE_sim_debug_data_add_vector(clmd->debug_data, x, d, 0.4, 0.4, 1, "target", hash_vertex(7982, s->kl));
+                       
+                       copy_v3_v3(d, s->target);
+                       BKE_sim_debug_data_add_vector(clmd->debug_data, x, 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
        }
index dd91d09cbb978e807bf1c54a655c921ec486a151..40b6174e250426eda4ba47fed5b228ef9015aa8b 100644 (file)
@@ -152,7 +152,7 @@ bool BPH_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int
                                           float r_f[3], float r_dfdx[3][3], float r_dfdv[3][3]);
 /* Angular bending force based on local target vectors */
 bool BPH_mass_spring_force_spring_bending_angular(struct Implicit_Data *data, int i, int j, int k, int block_ij, int block_jk, int block_ik,
-                                                  float restlen_ij, float restlen_jk, float stiffness, float damping);
+                                                  const float target[3], float stiffness, float damping);
 /* Global goal spring */
 bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, int spring_index, const float goal_x[3], const float goal_v[3],
                                        float stiffness, float damping,
index 0fd46b79aa0e3603616323c1b002ac675a7b2592..7a4ff5581b073f8ff93916cabe0d92641b79252f 100644 (file)
@@ -1670,7 +1670,8 @@ BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3]
 }
 
 BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
-                                      float restlen_ij, float restlen_jk, float stiffness, float damping,
+                                      const float goal[3],
+                                      float stiffness, float damping,
                                       int p, int q, const float dx[3], const float dv[3],
                                       float r_f[3])
 {
@@ -1679,7 +1680,7 @@ BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
        float vel_ij[3], vel_jk[3], vel_ortho[3];
        float f_bend[3], f_damp[3];
        float fi[3], fj[3], fk[3];
-       float target[3], dist[3];
+       float dist[3];
        
        zero_v3(fi);
        zero_v3(fj);
@@ -1703,33 +1704,14 @@ BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
        if (q == j) sub_v3_v3(vel_jk, dv);
        if (q == k) add_v3_v3(vel_jk, dv);
        
-       /* XXX fi(x) == fk(x) only holds true as long as we assume straight rest shape!
-        * eventually will become a bit more involved since the opposite segment
-        * gets its own target, under condition of having equal torque on both sides.
-        */
-       
        /* bending force */
-       mul_v3_v3fl(target, dir_jk, restlen_ij);
-       sub_v3_v3v3(dist, target, edge_ij);
-       mul_v3_v3fl(f_bend, dist, stiffness);
-       
-       sub_v3_v3(fi, f_bend);
-       add_v3_v3(fj, f_bend);
-       
-       mul_v3_v3fl(target, dir_ij, restlen_jk);
-       sub_v3_v3v3(dist, target, edge_jk);
+       sub_v3_v3v3(dist, goal, edge_jk);
        mul_v3_v3fl(f_bend, dist, stiffness);
        
        sub_v3_v3(fj, f_bend);
        add_v3_v3(fk, f_bend);
        
        /* damping force */
-       madd_v3_v3v3fl(vel_ortho, vel_ij, dir_ij, -dot_v3v3(vel_ij, dir_ij));
-       mul_v3_v3fl(f_damp, vel_ortho, damping);
-       
-       add_v3_v3(fi, f_damp);
-       sub_v3_v3(fj, f_damp);
-       
        madd_v3_v3v3fl(vel_ortho, vel_jk, dir_jk, -dot_v3v3(vel_jk, dir_jk));
        mul_v3_v3fl(f_damp, vel_ortho, damping);
        
@@ -1746,7 +1728,8 @@ BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k,
 
 /* Finite Differences method for estimating the jacobian of the force */
 BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k,
-                                             float restlen_ij, float restlen_jk, float stiffness, float damping,
+                                             const float goal[3],
+                                             float stiffness, float damping,
                                              int q, int p, float dfdx[3][3])
 {
        const float delta = 0.00001f; // TODO find a good heuristic for this
@@ -1760,12 +1743,14 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j,
        copy_m3_m3(dvec_neg, dvec_pos);
        negate_m3(dvec_neg);
        
+       /* XXX TODO offset targets to account for position dependency */
+       
        for (a = 0; a < 3; ++a) {
-               spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping,
+               spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
                                      q, p, dvec_pos[a], dvec_null[a], f);
                copy_v3_v3(dfdx[a], f);
                
-               spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping,
+               spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
                                      q, p, dvec_neg[a], dvec_null[a], f);
                sub_v3_v3(dfdx[a], f);
                
@@ -1777,7 +1762,8 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j,
 
 /* Finite Differences method for estimating the jacobian of the force */
 BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k,
-                                             float restlen_ij, float restlen_jk, float stiffness, float damping,
+                                             const float goal[3],
+                                             float stiffness, float damping,
                                              int q, int p, float dfdv[3][3])
 {
        const float delta = 0.00001f; // TODO find a good heuristic for this
@@ -1791,12 +1777,14 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j,
        copy_m3_m3(dvec_neg, dvec_pos);
        negate_m3(dvec_neg);
        
+       /* XXX TODO offset targets to account for position dependency */
+       
        for (a = 0; a < 3; ++a) {
-               spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping,
+               spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
                                      q, p, dvec_null[a], dvec_pos[a], f);
                copy_v3_v3(dfdv[a], f);
                
-               spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping,
+               spring_angbend_forces(data, i, j, k, goal, stiffness, damping,
                                      q, p, dvec_null[a], dvec_neg[a], f);
                sub_v3_v3(dfdv[a], f);
                
@@ -1810,31 +1798,34 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j,
  * See "Artistic Simulation of Curly Hair" (Pixar technical memo #12-03a)
  */
 bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, int j, int k, int block_ij, int block_jk, int block_ik,
-                                                  float restlen_ij, float restlen_jk, float stiffness, float damping)
+                                                  const float target[3], float stiffness, float damping)
 {
+       float goal[3];
        float fi[3], fj[3], fk[3];
        float dfi_dxi[3][3], dfj_dxi[3][3], dfj_dxj[3][3], dfk_dxi[3][3], dfk_dxj[3][3], dfk_dxk[3][3];
        float dfi_dvi[3][3], dfj_dvi[3][3], dfj_dvj[3][3], dfk_dvi[3][3], dfk_dvj[3][3], dfk_dvk[3][3];
        
        const float vecnull[3] = {0.0f, 0.0f, 0.0f};
        
-       spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, i, -1, vecnull, vecnull, fi);
-       spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, j, -1, vecnull, vecnull, fj);
-       spring_angbend_forces(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, -1, vecnull, vecnull, fk);
-       
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, i, i, dfi_dxi);
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, j, i, dfj_dxi);
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, j, j, dfj_dxj);
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, i, dfk_dxi);
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, j, dfk_dxj);
-       spring_angbend_estimate_dfdx(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, k, dfk_dxk);
-       
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, i, i, dfi_dvi);
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, j, i, dfj_dvi);
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, j, j, dfj_dvj);
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, i, dfk_dvi);
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, j, dfk_dvj);
-       spring_angbend_estimate_dfdv(data, i, j, k, restlen_ij, restlen_jk, stiffness, damping, k, k, dfk_dvk);
+       world_to_root_v3(data, j, goal, target);
+       
+       spring_angbend_forces(data, i, j, k, goal, stiffness, damping, i, -1, vecnull, vecnull, fi);
+       spring_angbend_forces(data, i, j, k, goal, stiffness, damping, j, -1, vecnull, vecnull, fj);
+       spring_angbend_forces(data, i, j, k, goal, stiffness, damping, k, -1, vecnull, vecnull, fk);
+       
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, i, dfi_dxi);
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, i, dfj_dxi);
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, j, dfj_dxj);
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, i, dfk_dxi);
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, j, dfk_dxj);
+       spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, k, dfk_dxk);
+       
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, i, dfi_dvi);
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, i, dfj_dvi);
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, j, dfj_dvj);
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, i, dfk_dvi);
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, j, dfk_dvj);
+       spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, k, dfk_dvk);
        
        /* add forces and jacobians to the solver data */
        add_v3_v3(data->F[i], fi);