0b78ead13c8023b9e264bec15ae1413352f07d90
[blender.git] / source / blender / blenkernel / intern / softbody.c
1 /*  softbody.c      
2  * 
3  * $Id: 
4  *
5  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
6  *
7  * This program is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License
9  * as published by the Free Software Foundation; either version 2
10  * of the License, or (at your option) any later version. The Blender
11  * Foundation also sells licenses for use in proprietary software under
12  * the Blender License.  See http://www.blender.org/BL/ for information
13  * about this.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program; if not, write to the Free Software Foundation,
22  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23  *
24  * The Original Code is Copyright (C) Blender Foundation
25  * All rights reserved.
26  *
27  * The Original Code is: all of this file.
28  *
29  * Contributor(s): none yet.
30  *
31  * ***** END GPL/BL DUAL LICENSE BLOCK *****
32  */
33
34 /*
35 ******
36 variables on the UI for now
37
38         float mediafrict;  friction to env 
39         float nodemass;   softbody mass of *vertex* 
40         float grav;        softbody amount of gravitaion to apply 
41         
42         float goalspring;  softbody goal springs 
43         float goalfrict;   softbody goal springs friction 
44         float mingoal;     quick limits for goal 
45         float maxgoal;
46
47         float inspring;   softbody inner springs 
48         float infrict;     softbody inner springs friction 
49
50 *****
51 */
52
53
54 #include <math.h>
55 #include <stdlib.h>
56 #include <string.h>
57
58 #include "MEM_guardedalloc.h"
59
60 /* types */
61 #include "DNA_curve_types.h"
62 #include "DNA_object_types.h"
63 #include "DNA_mesh_types.h"
64 #include "DNA_meshdata_types.h"
65 #include "DNA_lattice_types.h"
66 #include "DNA_scene_types.h"
67
68 #include "BLI_blenlib.h"
69 #include "BLI_arithb.h"
70
71 #include "BKE_displist.h"
72 #include "BKE_effect.h"
73 #include "BKE_global.h"
74 #include "BKE_object.h"
75 #include "BKE_softbody.h"
76 #include "BKE_utildefines.h"
77
78 #include  "BIF_editdeform.h"
79
80 extern bDeformGroup *get_named_vertexgroup(Object *ob, char *name);
81 extern int  get_defgroup_num (Object *ob, bDeformGroup        *dg);
82
83
84 /* ********** soft body engine ******* */
85 #define SOFTGOALSNAP  0.999f 
86 // if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
87 // removes *unnecessary* stiffnes from ODE system
88 #define HEUNWARNLIMIT 1 // 50 would be fine i think for detecting severe *stiff* stuff
89
90
91 float SoftHeunTol = 1.0f; // humm .. this should be calculated from sb parameters and sizes
92
93 /* local prototypes */
94 static void free_softbody_intern(SoftBody *sb);
95
96
97 /*+++ frame based timing +++*/
98
99 //physical unit of force is [kg * m / sec^2]
100
101 static float sb_grav_force_scale(Object *ob)
102 // since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm
103 // put it to a function here, so we can add user options later without touching simulation code
104 {
105         return (0.001f);
106 }
107
108 static float sb_fric_force_scale(Object *ob)
109 // rescaling unit of drag [1 / sec] to somehow reasonable
110 // put it to a function here, so we can add user options later without touching simulation code
111 {
112         return (0.01f);
113 }
114
115 static float sb_time_scale(Object *ob)
116 // defining the frames to *real* time relation
117 {
118         SoftBody *sb= ob->soft; // is supposed to be there
119         if (sb){
120                 return(sb->physics_speed); //hrms .. this could be IPO as well :) 
121                 // estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)]
122                 // 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions]  has period 65 frames
123         // theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM) 
124         }
125         return (1.0f);
126         /* 
127         this would be frames/sec independant timing assuming 25 fps is default
128         but does not work very well with NLA
129                 return (25.0f/G.scene->r.frs_sec)
130         */
131 }
132 /*--- frame based timing ---*/
133
134
135 static int count_mesh_quads(Mesh *me)
136 {
137         int a,result = 0;
138         MFace *mface= me->mface;
139         
140         if(mface) {
141                 for(a=me->totface; a>0; a--, mface++) {
142                         if(mface->v4) result++;
143                 }
144         }       
145         return result;
146 }
147
148 static void add_mesh_quad_diag_springs(Object *ob)
149 {
150         Mesh *me= ob->data;
151         MFace *mface= me->mface;
152         BodyPoint *bp;
153         BodySpring *bs, *bs_new;
154         int a ;
155         
156         if (ob->soft){
157                 int nofquads;
158                 
159                 nofquads = count_mesh_quads(me);
160                 if (nofquads) {
161                         /* resize spring-array to hold additional quad springs */
162                         bs_new= MEM_callocN( (ob->soft->totspring + nofquads *2 )*sizeof(BodySpring), "bodyspring");
163                         memcpy(bs_new,ob->soft->bspring,(ob->soft->totspring )*sizeof(BodySpring));
164                         
165                         if(ob->soft->bspring)
166                                 MEM_freeN(ob->soft->bspring); /* do this before reassigning the pointer  or have a 1st class memory leak */
167                         ob->soft->bspring = bs_new; 
168                         
169                         /* fill the tail */
170                         a = 0;
171                         bs = bs_new+ob->soft->totspring;
172                         bp= ob->soft->bpoint;
173                         if(mface ) {
174                                 for(a=me->totface; a>0; a--, mface++) {
175                                         if(mface->v4) {
176                                                 bs->v1= mface->v1;
177                                                 bs->v2= mface->v3;
178                                                 bs->strength= 1.0;
179                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
180                                                 bs++;
181                                                 bs->v1= mface->v2;
182                                                 bs->v2= mface->v4;
183                                                 bs->strength= 1.0;
184                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
185                                                 bs++;
186                                                 
187                                         }
188                                 }       
189                         }
190                         
191             /* now we can announce new springs */
192                         ob->soft->totspring += nofquads *2;
193                 }
194         }
195 }
196
197
198 static void add_bp_springlist(BodyPoint *bp,int springID)
199 {
200         int *newlist;
201         
202         if (bp->springs == NULL) {
203                 bp->springs = MEM_callocN( sizeof(int), "bpsprings");
204                 bp->springs[0] = springID;
205                 bp->nofsprings = 1;
206         }
207         else {
208                 bp->nofsprings++;
209                 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings");
210                 memcpy(newlist,bp->springs,(bp->nofsprings-1)* sizeof(int));
211                 MEM_freeN(bp->springs);
212                 bp->springs = newlist;
213                 bp->springs[bp->nofsprings-1] = springID;
214         }
215 }
216
217 /* do this once when sb is build
218 it is O(N^2) so scanning for springs every iteration is too expensive
219 */
220 static void build_bps_springlist(Object *ob)
221 {
222         SoftBody *sb= ob->soft; // is supposed to be there
223         BodyPoint *bp;  
224         BodySpring *bs; 
225         int a,b;
226         
227         if (sb==NULL) return; // paranoya check
228         
229         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
230                 /* scan for attached inner springs */   
231                 for(b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) {
232                         if (( (sb->totpoint-a) == bs->v1) ){ 
233                                 add_bp_springlist(bp,sb->totspring -b);
234                         }
235                         if (( (sb->totpoint-a) == bs->v2) ){ 
236                                 add_bp_springlist(bp,sb->totspring -b);
237                         }
238                 }//for springs
239                 // if (bp->nofsprings) printf(" node %d has %d spring links\n",a,bp->nofsprings);
240         }//for bp               
241 }
242
243
244 /* creates new softbody if didn't exist yet, makes new points and springs arrays */
245 /* called in mesh_to_softbody */
246 static void renew_softbody(Object *ob, int totpoint, int totspring)  
247 {
248         SoftBody *sb;
249         
250         if(ob->soft==NULL) ob->soft= sbNew();
251         else free_softbody_intern(ob->soft);
252         sb= ob->soft;
253            
254         if(totpoint) {
255                 sb->totpoint= totpoint;
256                 sb->totspring= totspring;
257                 
258                 sb->bpoint= MEM_mallocN( totpoint*sizeof(BodyPoint), "bodypoint");
259                 if(totspring) 
260                         sb->bspring= MEM_mallocN( totspring*sizeof(BodySpring), "bodyspring");
261         }
262 }
263
264 /* only frees internal data */
265 static void free_softbody_intern(SoftBody *sb)
266 {
267         if(sb) {
268                 int a;
269                 BodyPoint *bp;
270                 
271                 if(sb->bpoint){
272                         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
273                                 /* free spring list */ 
274                                 if (bp->springs != NULL) {
275                                         MEM_freeN(bp->springs);
276                                 }
277                         }
278                         MEM_freeN(sb->bpoint);
279                 }
280                 
281                 if(sb->bspring) MEM_freeN(sb->bspring);
282                 
283                 sb->totpoint= sb->totspring= 0;
284                 sb->bpoint= NULL;
285                 sb->bspring= NULL;
286         }
287 }
288
289
290 /* ************ dynamics ********** */
291
292 /* aye this belongs to arith.c */
293 static void Vec3PlusStVec(float *v, float s, float *v1)
294 {
295         v[0] += s*v1[0];
296         v[1] += s*v1[1];
297         v[2] += s*v1[2];
298 }
299
300 static int sb_deflect_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *force,float *cf ,float *bounce)
301 {
302         int deflected;
303         float s_actpos[3], s_futurepos[3];
304         VECCOPY(s_actpos,actpos);
305         if(futurepos)
306         VECCOPY(s_futurepos,futurepos);
307         if (bounce) *bounce *= 1.5f;
308                                 
309                                 
310         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
311                                         facenormal, cf, force , 1,
312                                         G.scene->r.cfra, ob->lay, ob);
313         return(deflected);
314                                 
315 }
316
317 /* for future use (BM)
318 static int sb_deflect_edge_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *slip ,float *bounce)
319 {
320         int deflected;
321         float dummy[3],s_actpos[3], s_futurepos[3];
322         SoftBody *sb= ob->soft; // is supposed to be there
323         VECCOPY(s_actpos,actpos);
324         VECCOPY(s_futurepos,futurepos);
325         if (slip)   *slip   *= 0.98f;
326         if (bounce) *bounce *= 1.5f;
327                                 
328                                 
329         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
330                                         facenormal, dummy, dummy , 2,
331                                         G.scene->r.cfra, ob->lay, ob);
332         return(deflected);
333                                 
334 }
335 */
336 // some functions removed here .. to help HOS on next merge (BM)
337
338 #define USES_FIELD              1
339 #define USES_DEFLECT    2
340 static int is_there_deflection(unsigned int layer)
341 {
342         Base *base;
343         int retval= 0;
344         
345         for(base = G.scene->base.first; base; base= base->next) {
346                 if( (base->lay & layer) && base->object->pd) {
347                         if(base->object->pd->forcefield) retval |= USES_FIELD;
348                         if(base->object->pd->deflect) retval |= USES_DEFLECT;
349                 }
350         }
351         return retval;
352 }
353
354 static void softbody_calc_forces(Object *ob, float forcetime)
355 {
356 /* rule we never alter free variables :bp->vec bp->pos in here ! 
357  * this will ruin adaptive stepsize AKA heun! (BM) 
358  */
359         SoftBody *sb= ob->soft; // is supposed to be there
360         BodyPoint  *bp;
361         BodyPoint *bproot;
362         BodySpring *bs; 
363         float iks, ks, kd, gravity, actspringlen, forcefactor, sd[3],
364         fieldfactor = 1000.0f, 
365         windfactor  = 250.0f;   
366         int a, b, do_effector;
367         
368         /* clear forces */
369         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
370                 bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
371         }
372         
373         gravity = sb->grav * sb_grav_force_scale(ob);   
374         /* check! */
375         do_effector= is_there_deflection(ob->lay);
376         iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
377         bproot= sb->bpoint; /* need this for proper spring addressing */
378         
379         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
380                 if(bp->goal < SOFTGOALSNAP){ // ommit this bp when i snaps
381                         float auxvect[3]; // aux unit vector  
382                         float velgoal[3];
383                         float absvel =0, projvel= 0;
384                         
385                         /* do goal stuff */
386                         if(ob->softflag & OB_SB_GOAL) {
387                                 /* true elastic goal */
388                                 VecSubf(auxvect,bp->origT,bp->pos);
389                                 ks  = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
390                                 bp->force[0]= ks*(auxvect[0]);
391                                 bp->force[1]= ks*(auxvect[1]);
392                                 bp->force[2]= ks*(auxvect[2]);
393                                 /* calulate damping forces generated by goals*/
394                                 VecSubf(velgoal,bp->origS, bp->origE);
395                                 kd =  sb->goalfrict * sb_fric_force_scale(ob) ;
396                                 
397                                 if (forcetime > 0.0 ) { // make sure friction does not become rocket motor on time reversal
398                                         bp->force[0]-= kd * (velgoal[0] + bp->vec[0]);
399                                         bp->force[1]-= kd * (velgoal[1] + bp->vec[1]);
400                                         bp->force[2]-= kd * (velgoal[2] + bp->vec[2]);
401                                 }
402                                 else {
403                                         bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
404                                         bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
405                                         bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
406                                 }
407                         }
408                         /* done goal stuff */
409                         
410                         
411                         /* gravitation */
412                         bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
413                         
414                         /* particle field & vortex */
415                         if(do_effector & USES_FIELD) {
416                                 float force[3]= {0.0f, 0.0f, 0.0f};
417                                 float speed[3]= {0.0f, 0.0f, 0.0f};
418                                 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); // just for calling functio once
419                                 
420                                 pdDoEffector(bp->pos, force, speed, (float)G.scene->r.cfra, ob->lay,PE_WIND_AS_SPEED);
421                                 // note: now we have wind as motion of media, so we can do anisotropic stuff here, 
422                                 // if we had vertex normals here(BM)
423                                 /* apply forcefield*/
424                                 VecMulf(force,fieldfactor* eval_sb_fric_force_scale); 
425                                 VECADD(bp->force, bp->force, force);
426                                 
427                                 /* friction in moving media */  
428                                 kd= sb->mediafrict* eval_sb_fric_force_scale;  
429                                 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
430                                 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
431                                 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
432                                 /* now we'll have nice centrifugal effect for vortex */
433                                 
434                         }
435                         else {
436                                 /* friction in media (not) moving*/
437                                 kd= sb->mediafrict* sb_fric_force_scale(ob);  
438                                 /* assume it to be proportional to actual velocity */
439                                 bp->force[0]-= bp->vec[0]*kd;
440                                 bp->force[1]-= bp->vec[1]*kd;
441                                 bp->force[2]-= bp->vec[2]*kd;
442                                 /* friction in media done */
443                         }
444                         
445                         /*other forces*/
446                         /* this is the place where other forces can be added
447                         yes, constraints and collision stuff should go here too (read baraff papers on that!)
448                         */
449                         /* try to match moving collision targets */
450                         /* master switch to turn collision off (BM)*/
451                         //if(0) {
452                         if(do_effector & USES_DEFLECT) {
453                                 /*sorry for decl. here i'll move 'em up when WIP is done (BM) */
454                                 float defforce[3];
455                                 float collisionpos[3],facenormal[3];
456                                 float cf = 1.0f;
457                                 float bounce = 0.5f;
458                                 kd = 1.0f;
459                                 defforce[0] = 0.0f;
460                                 defforce[1] = 0.0f;
461                                 defforce[2] = 0.0f;
462                                 
463                                 if (sb_deflect_face(ob,bp->pos, bp->pos, collisionpos, facenormal,defforce,&cf,&bounce)){
464                                         bp->force[0] += defforce[0]*kd;
465                                         bp->force[1] += defforce[1]*kd;
466                                         bp->force[2] += defforce[2]*kd;
467                                         bp->contactfrict = cf;
468                                 }
469                                 else{ 
470                                         bp->contactfrict = 0.0f;
471                                 }
472                                 
473                         }
474                         
475                         /*other forces done*/
476                         /* nice things could be done with anisotropic friction
477                         like wind/air resistance in normal direction
478                         --> having a piece of cloth sailing down 
479                         but this needs to have a *valid* vertex normal
480                         *valid* means to be calulated on time axis
481                         hrms .. may be a rough one could be used as well .. let's see 
482                         */
483                         
484                         if(ob->softflag & OB_SB_EDGES) {
485                                 if (sb->bspring){ // spring list exists at all ? 
486                                         for(b=bp->nofsprings;b>0;b--){
487                                                 bs = sb->bspring + bp->springs[b-1];
488                                                 if (( (sb->totpoint-a) == bs->v1) ){ 
489                                                         actspringlen= VecLenf( (bproot+bs->v2)->pos, bp->pos);
490                                                         VecSubf(sd,(bproot+bs->v2)->pos, bp->pos);
491                                                         Normalise(sd);
492                                                         
493                                                         // friction stuff V1
494                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v2)->vec);
495                                                         kd = sb->infrict * sb_fric_force_scale(ob);
496                                                         absvel  = Normalise(velgoal);
497                                                         projvel = ABS(Inpf(sd,velgoal));
498                                                         kd *= absvel * projvel;
499                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
500                                                         
501                                                         if(bs->len > 0.0) /* check for degenerated springs */
502                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
503                                                         else
504                                                                 forcefactor = actspringlen * iks;
505                                                         
506                                                         Vec3PlusStVec(bp->force,-forcefactor,sd);
507                                                         
508                                                 }
509                                                 
510                                                 if (( (sb->totpoint-a) == bs->v2) ){ 
511                                                         actspringlen= VecLenf( (bproot+bs->v1)->pos, bp->pos);
512                                                         VecSubf(sd,bp->pos,(bproot+bs->v1)->pos);
513                                                         Normalise(sd);
514                                                         
515                                                         // friction stuff V2
516                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v1)->vec);
517                                                         kd = sb->infrict * sb_fric_force_scale(ob);
518                                                         absvel  = Normalise(velgoal);
519                                                         projvel = ABS(Inpf(sd,velgoal));
520                                                         kd *= absvel * projvel;
521                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
522                                                         
523                                                         if(bs->len > 0.0)
524                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
525                                                         else
526                                                                 forcefactor = actspringlen * iks;
527                                                         Vec3PlusStVec(bp->force,+forcefactor,sd);                                                       
528                                                 }
529                                         }/* loop springs */
530                                 }/* existing spring list */ 
531                         }/*any edges*/
532                 }/*omit on snap */
533         }/*loop all bp's*/
534 }
535
536 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err)
537 {
538         /* time evolution */
539         /* actually does an explicit euler step mode == 0 */
540         /* or heun ~ 2nd order runge-kutta steps, mode 1,2 */
541         SoftBody *sb= ob->soft; // is supposed to be there
542         BodyPoint *bp;
543         float dx[3],dv[3];
544         float timeovermass;
545         float maxerr = 0.0;
546         int a, do_effector;
547
548     forcetime *= sb_time_scale(ob);
549         /* check! */
550         do_effector= is_there_deflection(ob->lay);
551     
552         // claim a minimum mass for vertex 
553         if (sb->nodemass > 0.09999f) timeovermass = forcetime/sb->nodemass;
554         else timeovermass = forcetime/0.09999f;
555         
556         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
557                 if(bp->goal < SOFTGOALSNAP){
558                         
559                         /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces  + more forces*/
560                         /* the ( ... )' operator denotes derivate respective time */
561                         /* the euler step for velocity then becomes */
562                         /* v(t + dt) = v(t) + a(t) * dt */ 
563                         bp->force[0]*= timeovermass; /* individual mass of node here */ 
564                         bp->force[1]*= timeovermass;
565                         bp->force[2]*= timeovermass;
566                         /* some nasty if's to have heun in here too */
567                         VECCOPY(dv,bp->force); 
568                         if (mode == 1){
569                                 VECCOPY(bp->prevvec, bp->vec);
570                                 VECCOPY(bp->prevdv, dv);
571                         }
572                         if (mode ==2){
573                                 /* be optimistic and execute step */
574                                 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]);
575                                 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]);
576                                 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]);
577                                 /* compare euler to heun to estimate error for step sizing */
578                                 maxerr = MAX2(maxerr,ABS(dv[0] - bp->prevdv[0]));
579                                 maxerr = MAX2(maxerr,ABS(dv[1] - bp->prevdv[1]));
580                                 maxerr = MAX2(maxerr,ABS(dv[2] - bp->prevdv[2]));
581                         }
582                         else {VECADD(bp->vec, bp->vec, bp->force);}
583
584                         /* so here is (x)'= v(elocity) */
585                         /* the euler step for location then becomes */
586                         /* x(t + dt) = x(t) + v(t) * dt */ 
587                         
588                         VECCOPY(dx,bp->vec);
589                         dx[0]*=forcetime ; 
590                         dx[1]*=forcetime ; 
591                         dx[2]*=forcetime ; 
592                         
593                         /* again some nasty if's to have heun in here too */
594                         if (mode ==1){
595                                 VECCOPY(bp->prevpos,bp->pos);
596                                 VECCOPY(bp->prevdx ,dx);
597                         }
598                         
599                         if (mode ==2){
600                                 bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]);
601                                 bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]);
602                                 bp->pos[2] = bp->prevpos[2] + 0.5f* ( dx[2] + bp->prevdx[2]);
603                                 maxerr = MAX2(maxerr,ABS(dx[0] - bp->prevdx[0]));
604                                 maxerr = MAX2(maxerr,ABS(dx[1] - bp->prevdx[1]));
605                                 maxerr = MAX2(maxerr,ABS(dx[2] - bp->prevdx[2]));
606 /* kind of hack .. while inside collision target .. make movement more *viscous* */
607                                 if (bp->contactfrict > 0.0f){
608                                         bp->vec[0] *= (1.0 - bp->contactfrict);
609                                         bp->vec[1] *= (1.0 - bp->contactfrict);
610                                         bp->vec[2] *= (1.0 - bp->contactfrict);
611                                 }
612                         }
613                         else { VECADD(bp->pos, bp->pos, dx);}
614 // experimental particle collision suff was here .. just to help HOS on next merge (BM)
615                 }//snap
616         } //for
617         if (err){ /* so step size will be controlled by biggest difference in slope */
618                 *err = maxerr;
619         }
620 }
621
622 /* used by heun when it overshoots */
623 static void softbody_restore_prev_step(Object *ob)
624 {
625         SoftBody *sb= ob->soft; // is supposed to be there
626         BodyPoint *bp;
627         int a;
628         
629         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
630                 VECCOPY(bp->vec, bp->prevvec);
631                 VECCOPY(bp->pos, bp->prevpos);
632         }
633 }
634
635
636 static void softbody_apply_goalsnap(Object *ob)
637 {
638         SoftBody *sb= ob->soft; // is supposed to be there
639         BodyPoint *bp;
640         int a;
641         
642         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
643                 if (bp->goal >= SOFTGOALSNAP){
644                         VECCOPY(bp->prevpos,bp->pos);
645                         VECCOPY(bp->pos,bp->origT);
646                 }               
647         }
648 }
649
650 /* unused */
651 #if 0
652 static void softbody_force_goal(Object *ob)
653 {
654         SoftBody *sb= ob->soft; // is supposed to be there
655         BodyPoint *bp;
656         int a;
657         
658         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {           
659                 VECCOPY(bp->pos,bp->origT);
660                 bp->vec[0] = bp->origE[0] - bp->origS[0];
661                 bp->vec[1] = bp->origE[1] - bp->origS[1];
662                 bp->vec[2] = bp->origE[2] - bp->origS[2];               
663         }
664 }
665 #endif
666
667 /* expects full initialized softbody */
668 static void interpolate_exciter(Object *ob, int timescale, int time)
669 {
670         SoftBody *sb= ob->soft;
671         BodyPoint *bp;
672         float f;
673         int a;
674
675         // note: i removed Mesh usage here, softbody should remain generic! (ton)
676         
677         f = (float)time/(float)timescale;
678         
679         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {   
680                 bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]); 
681                 bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]); 
682                 bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]); 
683                 if (bp->goal >= SOFTGOALSNAP){
684                         bp->vec[0] = bp->origE[0] - bp->origS[0];
685                         bp->vec[1] = bp->origE[1] - bp->origS[1];
686                         bp->vec[2] = bp->origE[2] - bp->origS[2];
687                 }
688         }
689         
690         if(ob->softflag & OB_SB_EDGES) {
691                 /* hrms .. do springs alter their lenght ?
692                 bs= ob->soft->bspring;
693                 bp= ob->soft->bpoint;
694                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, bs++) {
695                         bs->len= VecLenf( (bp+bs->v1)->origT, (bp+bs->v2)->origT);
696                 }
697                 */
698         }
699 }
700
701
702 /* ************ convertors ********** */
703
704 /*  for each object type we need;
705     - xxxx_to_softbody(Object *ob)      : a full (new) copy
706         - xxxx_update_softbody(Object *ob)  : update refreshes current positions
707     - softbody_to_xxxx(Object *ob)      : after simulation, copy vertex locations back
708 */
709
710 /* helper  call */
711 static int object_has_edges(Object *ob) 
712 {
713         if(ob->type==OB_MESH) {
714                 Mesh *me= ob->data;
715                 if(me->medge) return 1;
716         }
717         else if(ob->type==OB_LATTICE) {
718                 ;
719         }
720         
721         return 0;
722 }
723
724 /* helper  call */
725 static void set_body_point(Object *ob, BodyPoint *bp, float *vec)
726 {
727         
728         VECCOPY(bp->pos, vec);
729         Mat4MulVecfl(ob->obmat, bp->pos);  // yep, sofbody is global coords
730         VECCOPY(bp->origS, bp->pos);
731         VECCOPY(bp->origE, bp->pos);
732         VECCOPY(bp->origT, bp->pos);
733         
734         bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0;
735         bp->weight= 1.0;
736         bp->goal= ob->soft->defgoal;
737         
738         bp->nofsprings= 0;
739         bp->springs= NULL;
740         bp->contactfrict = 0.0f;
741 }
742
743
744 /* copy original (new) situation in softbody, as result of matrices or deform */
745 /* is assumed to enter function with ob->soft, but can be without points */
746 static void mesh_update_softbody(Object *ob)
747 {
748         Mesh *me= ob->data;
749         MVert *mvert= me->mvert;
750 /*      MEdge *medge= me->medge;  */ /*unused*/
751         BodyPoint *bp;
752         int a;
753         
754         /* possible after a file read... */
755         if(ob->soft->totpoint!=me->totvert) sbObjectToSoftbody(ob);
756         
757         if(me->totvert) {
758         
759                 bp= ob->soft->bpoint;
760                 for(a=0; a<me->totvert; a++, mvert++, bp++) {
761                         VECCOPY(bp->origS, bp->origE);
762                         VECCOPY(bp->origE, mvert->co);
763                         Mat4MulVecfl(ob->obmat, bp->origE);
764                         VECCOPY(bp->origT, bp->origE);
765                 }
766                 
767                 if(ob->softflag & OB_SB_EDGES) {
768                         
769                         /* happens when in UI edges was set */
770                         if(ob->soft->bspring==NULL) 
771                                 if(object_has_edges(ob)) sbObjectToSoftbody(ob);
772                 
773                         /* hrms .. do springs alter their lenght ? (yes, mesh keys would (ton))
774                         if(medge) {
775                                 bs= ob->soft->bspring;
776                                 bp= ob->soft->bpoint;
777                                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, medge++, bs++) { 
778                                         bs->len= VecLenf( (bp+bs->v1)->origE, (bp+bs->v2)->origE);
779                                 }
780                         }
781                         */
782                 }
783         }
784 }
785
786
787 static void get_scalar_from_vertexgroup(Object *ob, int vertID, short groupindex, float *target)
788 /* result 0 on success, else indicates error number
789 -- kind of *inverse* result defintion,
790 -- but this way we can signal error condition to caller  
791 -- and yes this function must not be here but in a *vertex group module*
792 */
793 {
794         MDeformVert *dv;
795         int i;
796         
797         /* spot the vert in deform vert list at mesh */
798         if(ob->type==OB_MESH) {
799                 if (((Mesh *)ob->data)->dvert) {
800                         dv = ((Mesh*)ob->data)->dvert + vertID; 
801                         /* Lets see if this vert is in the weight group */
802                         for (i=0; i<dv->totweight; i++){
803                                 if (dv->dw[i].def_nr == groupindex){
804                                         *target= dv->dw[i].weight; /* got it ! */
805                                         break;
806                                 }
807                         }
808                 }
809         }
810
811
812 /* makes totally fresh start situation */
813 static void mesh_to_softbody(Object *ob)
814 {
815         SoftBody *sb;
816         Mesh *me= ob->data;
817         MVert *mvert= me->mvert;
818         MEdge *medge= me->medge;
819         BodyPoint *bp;
820         BodySpring *bs;
821         float goalfac;
822         int a, totedge;
823         
824         if (ob->softflag & OB_SB_EDGES) totedge= me->totedge;
825         else totedge= 0;
826         
827         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
828         renew_softbody(ob, me->totvert, totedge);
829                 
830         /* we always make body points */
831         sb= ob->soft;   
832         bp= sb->bpoint;
833         goalfac= ABS(sb->maxgoal - sb->mingoal);
834         
835         for(a=me->totvert; a>0; a--, mvert++, bp++) {
836                 
837                 set_body_point(ob, bp, mvert->co);
838                 
839                 /* get scalar values needed  *per vertex* from vertex group functions,
840                 so we can *paint* them nicly .. 
841                 they are normalized [0.0..1.0] so may be we need amplitude for scale
842                 which can be done by caller but still .. i'd like it to go this way 
843                 */ 
844
845                 if(sb->vertgroup) {
846                         get_scalar_from_vertexgroup(ob, me->totvert - a, sb->vertgroup-1, &bp->goal);
847                         // do this always, regardless successfull read from vertex group
848                         bp->goal= sb->mingoal + bp->goal*goalfac;
849                 }
850                 /* a little ad hoc changing the goal control to be less *sharp* */
851                 bp->goal = (float)pow(bp->goal, 4.0f);
852                         
853                 /* to proove the concept
854                 this would enable per vertex *mass painting*
855                 strcpy(name,"SOFTMASS");
856                 error = get_scalar_from_named_vertexgroup(ob,name,me->totvert - a,&temp);
857                 if (!error) bp->mass = temp * ob->rangeofmass;
858                 */
859         }
860
861         /* but we only optionally add body edge springs */
862         if (ob->softflag & OB_SB_EDGES) {
863                 if(medge) {
864                         bs= sb->bspring;
865                         bp= sb->bpoint;
866                         for(a=me->totedge; a>0; a--, medge++, bs++) {
867                                 bs->v1= medge->v1;
868                                 bs->v2= medge->v2;
869                                 bs->strength= 1.0;
870                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
871                         }
872
873                 
874                         /* insert *diagonal* springs in quads if desired */
875                         if (ob->softflag & OB_SB_QUADS) {
876                                 add_mesh_quad_diag_springs(ob);
877                         }
878
879                         build_bps_springlist(ob); /* big mesh optimization */
880                 }
881         }
882         
883 }
884
885 /* copies current sofbody position in mesh, so do this within modifier stacks! */
886 static void softbody_to_mesh(Object *ob)
887 {
888         Mesh *me= ob->data;
889         MVert *mvert;
890         BodyPoint *bp;
891         int a;
892
893         bp= ob->soft->bpoint;
894         mvert= me->mvert;
895         for(a=me->totvert; a>0; a--, mvert++, bp++) {
896                 VECCOPY(mvert->co, bp->pos);
897                 Mat4MulVecfl(ob->imat, mvert->co);      // softbody is in global coords
898         }
899
900 }
901
902 /* makes totally fresh start situation */
903 static void lattice_to_softbody(Object *ob)
904 {
905         SoftBody *sb;
906         Lattice *lt= ob->data;
907         BodyPoint *bop;
908         BPoint *bp;
909         int a, totvert;
910         
911         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
912         
913         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
914         renew_softbody(ob, totvert, 0);
915         
916         /* we always make body points */
917         sb= ob->soft;   
918         
919         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
920                 set_body_point(ob, bop, bp->vec);
921         }
922 }
923
924 /* copies current sofbody position */
925 static void softbody_to_lattice(Object *ob)
926 {
927         SoftBody *sb;
928         Lattice *lt= ob->data;
929         BodyPoint *bop;
930         BPoint *bp;
931         int a, totvert;
932         
933         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
934         sb= ob->soft;   
935         
936         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
937                 VECCOPY(bp->vec, bop->pos);
938                 Mat4MulVecfl(ob->imat, bp->vec);        // softbody is in global coords
939         }
940 }
941
942 /* copy original (new) situation in softbody, as result of matrices or deform */
943 /* is assumed to enter function with ob->soft, but can be without points */
944 static void lattice_update_softbody(Object *ob)
945 {
946         Lattice *lt= ob->data;
947         BodyPoint *bop;
948         BPoint *bp;
949         int a, totvert;
950         
951         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
952         
953         /* possible after a file read... */
954         if(ob->soft->totpoint!=totvert) sbObjectToSoftbody(ob);
955         
956         for(a= totvert, bp= lt->def, bop= ob->soft->bpoint; a>0; a--, bp++, bop++) {
957                 VECCOPY(bop->origS, bop->origE);
958                 VECCOPY(bop->origE, bp->vec);
959                 Mat4MulVecfl(ob->obmat, bop->origE);
960                 VECCOPY(bop->origT, bop->origE);
961         }       
962 }
963
964
965 /* copies softbody result back in object */
966 /* only used in sbObjectStep() */
967 static void softbody_to_object(Object *ob)
968 {
969         
970         if(ob->soft==NULL) return;
971         
972         /* inverse matrix is not uptodate... */
973         Mat4Invert(ob->imat, ob->obmat);
974         
975         switch(ob->type) {
976         case OB_MESH:
977                 softbody_to_mesh(ob);
978                 break;
979         case OB_LATTICE:
980                 softbody_to_lattice(ob);
981                 break;
982         }
983 }
984
985 /* copy original (new) situation in softbody, as result of matrices or deform */
986 /* used in sbObjectStep() and sbObjectReset() */
987 /* assumes to have ob->soft, but can be entered without points */
988 static void object_update_softbody(Object *ob)
989 {
990         
991         switch(ob->type) {
992         case OB_MESH:
993                 mesh_update_softbody(ob);
994                 break;
995         case OB_LATTICE:
996                 lattice_update_softbody(ob);
997                 break;
998         }
999         
1000 }
1001
1002
1003
1004 /* ************ Object level, exported functions *************** */
1005
1006 /* allocates and initializes general main data */
1007 SoftBody *sbNew(void)
1008 {
1009         SoftBody *sb;
1010         
1011         sb= MEM_callocN(sizeof(SoftBody), "softbody");
1012         
1013         sb->mediafrict= 0.5; 
1014         sb->nodemass= 1.0;
1015         sb->grav= 0.0; 
1016         sb->physics_speed= 1.0;
1017         sb->rklimit= 0.1;
1018
1019         sb->goalspring= 0.5; 
1020         sb->goalfrict= 0.0; 
1021         sb->mingoal= 0.0;  
1022         sb->maxgoal= 1.0;
1023         sb->defgoal= 0.7;
1024         
1025         sb->inspring= 0.5;
1026         sb->infrict= 0.5; 
1027         
1028         return sb;
1029 }
1030
1031 /* frees all */
1032 void sbFree(SoftBody *sb)
1033 {
1034         free_softbody_intern(sb);
1035         MEM_freeN(sb);
1036 }
1037
1038
1039 /* makes totally fresh start situation */
1040 void sbObjectToSoftbody(Object *ob)
1041 {
1042
1043         switch(ob->type) {
1044         case OB_MESH:
1045                 mesh_to_softbody(ob);
1046                 break;
1047         case OB_LATTICE:
1048                 lattice_to_softbody(ob);
1049                 break;
1050         }
1051         
1052         if(ob->soft) ob->soft->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1053         ob->softflag &= ~OB_SB_REDO;
1054 }
1055
1056 /* reset all motion */
1057 void sbObjectReset(Object *ob)
1058 {
1059         SoftBody *sb= ob->soft;
1060         BodyPoint *bp;
1061         int a;
1062         
1063         if(sb==NULL) return;
1064         sb->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1065         
1066         object_update_softbody(ob);
1067         
1068         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
1069                 // origS is previous timestep
1070                 VECCOPY(bp->origS, bp->origE);
1071                 VECCOPY(bp->pos, bp->origE);
1072                 VECCOPY(bp->origT, bp->origE);
1073                 bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f;
1074
1075                 // no idea about the Heun stuff! (ton)
1076                 VECCOPY(bp->prevpos, bp->pos);
1077                 VECCOPY(bp->prevvec, bp->vec);
1078                 VECCOPY(bp->prevdx, bp->vec);
1079                 VECCOPY(bp->prevdv, bp->vec);
1080         }
1081 }
1082
1083
1084 /* simulates one step. framenr is in frames */
1085 /* copies result back to object, displist */
1086 void sbObjectStep(Object *ob, float framenr)
1087 {
1088         SoftBody *sb;
1089         Base *base;
1090         float dtime;
1091         int timescale,t;
1092         float ctime, forcetime;
1093         float err;
1094
1095         /* this variable is set while transform(). with lattices also having a softbody, 
1096            it calls lattice_modifier() all the time... has no displist yet. Is temporal
1097            hack which should be resolved with proper depgraph usage + storage of deformed
1098            vertices in lattice (ton) */
1099         if(G.moving) return;
1100         
1101         /* remake softbody if: */
1102         if( (ob->softflag & OB_SB_REDO) ||              // signal after weightpainting
1103                 (ob->soft==NULL) ||                                     // just to be nice we allow full init
1104                 (ob->soft->totpoint==0) )                       // after reading new file, or acceptable as signal to refresh
1105                         sbObjectToSoftbody(ob);
1106         
1107         sb= ob->soft;
1108
1109         /* still no points? go away */
1110         if(sb->totpoint==0) return;
1111         
1112         /* reset deflector cache, sumohandle is free, but its still sorta abuse... (ton) */
1113         for(base= G.scene->base.first; base; base= base->next) {
1114                 base->object->sumohandle= NULL;
1115         }
1116         
1117         /* checking time: */
1118         ctime= bsystem_time(ob, NULL, framenr, 0.0);
1119         dtime= ctime - sb->ctime;
1120                 // bail out for negative or for large steps
1121         if(dtime<0.0 || dtime >= 9.9*G.scene->r.framelen) { // note: what is G.scene->r.framelen for ? (BM)
1122                 sbObjectReset(ob);
1123                 return;
1124         }
1125         
1126         /* the simulator */
1127         
1128         if(dtime > 0.0) {       // note: what does this mean now? (ton)
1129                 //answer (BM) :
1130                 //dtime is still in [frames]
1131                 //we made sure dtime is >= 0.0
1132                 //but still need to handle dtime == 0.0 -> just return sb as is, just to be nice
1133                 object_update_softbody(ob);
1134                 
1135                 if (TRUE) {     // RSOL1 always true now (ton)
1136                         /* special case of 2nd order Runge-Kutta type AKA Heun */
1137                         float timedone =0.0; // how far did we get without violating error condition
1138                         /* loops = counter for emergency brake
1139                          * we don't want to lock up the system if physics fail
1140                          */
1141                         int loops =0 ; 
1142                         SoftHeunTol = sb->rklimit; // humm .. this should be calculated from sb parameters and sizes
1143
1144                         forcetime = dtime; /* hope for integrating in one step */
1145                         while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) )
1146                         {
1147                                 if (ABS(dtime) > 9.0 ){
1148                                         if(G.f & G_DEBUG) printf("SB_STEPSIZE \n");
1149                                         break; // sorry but i must assume goal movement can't be interpolated any more
1150                                 }
1151                                 //set goals in time
1152                                 interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
1153                                 // do predictive euler step
1154                                 softbody_calc_forces(ob, forcetime);
1155                                 softbody_apply_forces(ob, forcetime, 1, NULL);
1156                                 // crop new slope values to do averaged slope step
1157                                 softbody_calc_forces(ob, forcetime);
1158                                 softbody_apply_forces(ob, forcetime, 2, &err);
1159                                 softbody_apply_goalsnap(ob);
1160
1161                                 if (err > SoftHeunTol){ // error needs to be scaled to some quantity
1162                                         softbody_restore_prev_step(ob);
1163                                         forcetime /= 2.0;
1164                                 }
1165                                 else {
1166                                         float newtime = forcetime * 1.1f; // hope for 1.1 times better conditions in next step
1167                                         
1168                                         if (err > SoftHeunTol/2.0){ // stay with this stepsize unless err really small
1169                                                 newtime = forcetime;
1170                                         }
1171                                         timedone += forcetime;
1172                                         if (forcetime > 0.0)
1173                                                 forcetime = MIN2(dtime - timedone,newtime);
1174                                         else 
1175                                                 forcetime = MAX2(dtime - timedone,newtime);
1176                                 }
1177                                 loops++;
1178                         }
1179                         // move snapped to final position
1180                         interpolate_exciter(ob, 2, 2);
1181                         softbody_apply_goalsnap(ob);
1182                         
1183                         if(G.f & G_DEBUG) {
1184                                 if (loops > HEUNWARNLIMIT) /* monitor high loop counts say 1000 after testing */
1185                                         printf("%d heun integration loops/frame \n",loops);
1186                         }
1187                 }
1188                 else{
1189                         /* do brute force explicit euler */
1190                         /* inner intagration loop */
1191                         /* */
1192                         // loop n times so that n*h = duration of one frame := 1
1193                         // x(t+h) = x(t) + h*v(t);
1194                         // v(t+h) = v(t) + h*f(x(t),t);
1195                         timescale = (int)(sb->rklimit * ABS(dtime)); 
1196                         for(t=1 ; t <= timescale; t++) {
1197                                 if (ABS(dtime) > 15 ) break;
1198                                 
1199                                 /* the *goal* mesh must use the n*h timing too !
1200                                 use *cheap* linear intepolation for that  */
1201                                 interpolate_exciter(ob,timescale,t);                    
1202                                 if (timescale > 0 ) {
1203                                         forcetime = dtime/timescale;
1204                                         
1205                                         /* does not fit the concept sloving ODEs :) */
1206                                         /*                      softbody_apply_goal(ob,forcetime );  */
1207                                         
1208                                         /* explicit Euler integration */
1209                                         /* we are not controling a nuclear power plant! 
1210                                         so rought *almost* physical behaviour is acceptable.
1211                                         in cases of *mild* stiffnes cranking up timscale -> decreasing stepsize *h*
1212                                         avoids instability      */
1213                                         softbody_calc_forces(ob,forcetime);
1214                                         softbody_apply_forces(ob,forcetime,0, NULL);
1215                                         softbody_apply_goalsnap(ob);
1216                                         
1217                                         //      if (0){
1218                                         /* ok here comes the ├╝berhammer
1219                                         use a semi implicit euler integration to tackle *all* stiff conditions 
1220                                         but i doubt the cost/benifit holds for most of the cases
1221                                         -- to be coded*/
1222                                         //      }
1223                                         
1224                                 }
1225                         }
1226                 }
1227                 
1228                 /* and apply to vertices */
1229                  softbody_to_object(ob);
1230                 
1231                 sb->ctime= ctime;
1232         } // if(ABS(dtime) > 0.0) 
1233         else {
1234                 // rule : you have asked for the current state of the softobject 
1235                 // since dtime= ctime - ob->soft->ctime== 0.0;
1236                 // and we were not notifified about any other time changes 
1237                 // so here it is !
1238                 softbody_to_object(ob);
1239         }
1240
1241         /* reset deflector cache */
1242         for(base= G.scene->base.first; base; base= base->next) {
1243                 if(base->object->sumohandle) {
1244                         MEM_freeN(base->object->sumohandle);
1245                         base->object->sumohandle= NULL;
1246                 }
1247         }
1248         
1249 }
1250