tcontactmodels.cuh - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tcontactmodels.cuh (21251B)
---
1 #ifndef CONTACTMODELS_CUH_
2 #define CONTACTMODELS_CUH_
3
4 // contactmodels.cuh
5 // Functions governing repulsive forces between contacts
6
7
8 // Linear viscoelastic contact model for particle-wall interactions
9 // with tangential friction and rolling resistance
10 __device__ Float contactLinear_wall(
11 Float3* F,
12 Float3* T,
13 Float* es_dot,
14 Float* ev_dot,
15 Float* p,
16 const unsigned int idx_a,
17 const Float radius_a,
18 const Float4* __restrict__ dev_vel_sorted,
19 const Float4* __restrict__ dev_angvel_sorted,
20 const Float3 n,
21 const Float delta,
22 const Float wvel)
23 {
24 // Fetch particle velocities from global memory
25 Float4 vel_tmp = dev_vel_sorted[idx_a];
26 Float4 angvel_tmp = dev_angvel_sorted[idx_a];
27
28 // Convert velocities to three-component vectors
29 Float3 vel_linear = MAKE_FLOAT3(
30 vel_tmp.x,
31 vel_tmp.y,
32 vel_tmp.z);
33 Float3 angvel = MAKE_FLOAT3(
34 angvel_tmp.x,
35 angvel_tmp.y,
36 angvel_tmp.z);
37
38 // Store the length of the angular velocity for later use
39 Float angvel_length = length(angvel);
40
41 // Contact velocity is the sum of the linear and
42 // rotational components
43 Float3 vel = vel_linear + (radius_a + delta/2.0) * cross(n, angvel) + wvel;
44
45 // Normal component of the contact velocity
46 Float vel_n = dot(vel_linear, n);
47
48 // The tangential velocity is the contact velocity
49 // with the normal component subtracted
50 const Float3 vel_t = vel - n * (dot(n, vel));
51 const Float vel_t_length = length(vel_t);
52
53 // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
54 // 2013) based on Young's modulus if params.E > 0.0.
55 // Use the calculated stiffness for normal and tangential components.
56 Float k_n = devC_params.k_n;
57 if (devC_params.E > .001) {
58 k_n = 3.141592654/2.0 * devC_params.E * radius_a;
59 }
60
61 // Normal force component: Elastic - viscous damping
62 Float3 f_n = fmax(0.0, -k_n*delta
63 - devC_params.gamma_wn*vel_n) * n;
64 const Float f_n_length = length(f_n); // Save length for later use
65
66 // Store the energy lost by viscous damping. See derivation in
67 // contactLinear()
68 *ev_dot += devC_params.gamma_wn * vel_n * vel_n;
69
70 // Initialize vectors
71 Float3 f_t = MAKE_FLOAT3(0.0, 0.0, 0.0);
72
73 // Check that the tangential velocity is high enough to avoid
74 // divide by zero (producing a NaN)
75 if (vel_t_length > 0.0 && devC_params.gamma_wt > 0.0) {
76
77 // Tangential force by viscous model
78 const Float3 f_t_visc = devC_params.gamma_wt * vel_t;
79 const Float f_t_visc_length = length(f_t_visc);
80
81 // Determine max. friction
82 Float f_t_limit;
83 if (vel_t_length > 0.0) { // Dynamic
84 f_t_limit = devC_params.mu_wd * f_n_length;
85 } else { // Static
86 f_t_limit = devC_params.mu_ws * f_n_length;
87 }
88
89 // If the shear force component exceeds the friction,
90 // the particle slips and energy is dissipated
91 if (f_t_visc_length < f_t_limit) {
92 f_t = -1.0*f_t_visc;
93
94 } else { // Dynamic friction, friction failure
95 f_t = -f_t_limit * vel_t/vel_t_length;
96
97 // Shear energy production rate [W]
98 //*es_dot += -dot(vel_t, f_t);
99 *es_dot += length(length(f_t) * vel_t * devC_dt) / devC_dt;
100 }
101 }
102
103 // Total force from wall
104 *F += f_n + f_t;
105
106 // Total torque from wall
107 *T += cross(-(radius_a + delta*0.5)*n, f_t);
108
109 // Pressure excerted onto particle from this contact
110 *p += f_n_length / (4.0f * PI * radius_a*radius_a);
111
112 // Return force excerted onto the wall
113 return dot(f_n, n);
114 }
115
116
117 // Linear vicoelastic contact model for particle-particle interactions
118 // with tangential friction and rolling resistance
119 __device__ void contactLinearViscous(
120 Float3* F,
121 Float3* T,
122 Float* es_dot,
123 Float* ev_dot,
124 Float* p,
125 const unsigned int idx_a,
126 const unsigned int idx_b,
127 const Float4* __restrict__ dev_vel_sorted,
128 const Float4* __restrict__ dev_angvel_sorted,
129 const Float radius_a,
130 const Float radius_b,
131 const Float3 x_ab,
132 const Float x_ab_length,
133 const Float delta_ab,
134 const Float kappa)
135 {
136
137 // Allocate variables and fetch missing time=t values for particle A and B
138 Float4 vel_a = dev_vel_sorted[idx_a];
139 Float4 vel_b = dev_vel_sorted[idx_b];
140 Float4 angvel4_a = dev_angvel_sorted[idx_a];
141 Float4 angvel4_b = dev_angvel_sorted[idx_b];
142
143 // Convert to Float3's
144 Float3 angvel_a = MAKE_FLOAT3(angvel4_a.x, angvel4_a.y, angvel4_a.z);
145 Float3 angvel_b = MAKE_FLOAT3(angvel4_b.x, angvel4_b.y, angvel4_b.z);
146
147 // Force between grain pair decomposed into normal- and tangential part
148 Float3 f_n, f_t, f_c;
149
150 // Normal vector of contact
151 Float3 n_ab = x_ab/x_ab_length;
152
153 // Relative contact interface velocity, w/o rolling
154 Float3 vel_ab_linear = MAKE_FLOAT3(vel_a.x - vel_b.x,
155 vel_a.y - vel_b.y,
156 vel_a.z - vel_b.z);
157
158 // Relative contact interface velocity of particle surfaces at
159 // the contact, with rolling (Hinrichsen and Wolf 2004, eq. 13.10)
160 Float3 vel_ab = vel_ab_linear
161 + (radius_a + delta_ab/2.0f) * cross(n_ab, angvel_a)
162 + (radius_b + delta_ab/2.0f) * cross(n_ab, angvel_b);
163
164 // Relative contact interface rolling velocity
165 Float3 angvel_ab = angvel_a - angvel_b;
166 Float angvel_ab_length = length(angvel_ab);
167
168 // Normal component of the relative contact interface velocity
169 Float vel_n_ab = dot(vel_ab_linear, n_ab);
170
171 // Tangential component of the relative contact interface velocity
172 // Hinrichsen and Wolf 2004, eq. 13.9
173 Float3 vel_t_ab = vel_ab - (n_ab * dot(vel_ab, n_ab));
174 Float vel_t_ab_length = length(vel_t_ab);
175
176 // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
177 // 2013) based on Young's modulus if params.E > 0.0.
178 // Use the calculated stiffness for normal and tangential components.
179 Float k_n = devC_params.k_n;
180 if (devC_params.E > .001) {
181 k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
182 }
183
184 // Normal force component: Elastic - viscous damping
185 f_n = fmax(0.0, -k_n * delta_ab
186 - devC_params.gamma_n * vel_n_ab) * n_ab;
187
188 // Make sure the viscous damping doesn't exceed the elastic component,
189 // i.e. the damping factor doesn't exceed the critical damping, 2*sqrt(m*k_n)
190 if (dot(f_n, n_ab) < 0.0f)
191 f_n = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
192
193 Float f_n_length = length(f_n);
194
195 // Add max. capillary force
196 f_c = -kappa * sqrtf(radius_a * radius_b) * n_ab;
197
198 // Initialize force vectors to zero
199 f_t = MAKE_FLOAT3(0.0, 0.0, 0.0);
200
201 // Shear force component: Nonlinear relation
202 // Coulomb's law of friction limits the tangential force to less or equal
203 // to the normal force
204 if (vel_t_ab_length > 0.0) {
205
206 // Tangential force by viscous model
207 Float f_t_visc = devC_params.gamma_t * vel_t_ab_length;
208
209 // Determine max. friction
210 Float f_t_limit;
211 if (vel_t_ab_length > 0.001f) { // Dynamic
212 f_t_limit = devC_params.mu_d * length(f_n-f_c);
213 } else { // Static
214 f_t_limit = devC_params.mu_s * length(f_n-f_c);
215 }
216
217 // If the shear force component exceeds the friction,
218 // the particle slips and energy is dissipated
219 if (f_t_visc < f_t_limit) { // Static
220 f_t = -f_t_visc * vel_t_ab/vel_t_ab_length;
221
222 } else { // Dynamic, friction failure
223 f_t = -f_t_limit * vel_t_ab/vel_t_ab_length;
224
225 // Shear friction production rate [W]
226 //*es_dot += -dot(vel_t_ab, f_t);
227 }
228 }
229
230 // Add force components from this collision to total force for particle
231 *F += f_n + f_t + f_c;
232 *T += -(radius_a + delta_ab/2.0f) * cross(n_ab, f_t);
233
234 // Pressure excerted onto the particle from this contact
235 *p += f_n_length / (4.0f * PI * radius_a*radius_a);
236
237 } // End of contactLinearViscous()
238
239
240 // Linear elastic contact model for particle-particle interactions
241 __device__ void contactLinear(
242 Float3* F,
243 Float3* T,
244 Float* es_dot,
245 Float* ev_dot,
246 Float* p,
247 const unsigned int idx_a_orig,
248 const unsigned int idx_b_orig,
249 const Float4 vel_a,
250 const Float4* __restrict__ dev_vel,
251 const Float3 angvel_a,
252 const Float4* __restrict__ dev_angvel,
253 const Float radius_a,
254 const Float radius_b,
255 const Float3 x,
256 const Float x_length,
257 const Float delta,
258 Float4* __restrict__ dev_delta_t,
259 const unsigned int mempos)
260 {
261
262 // Allocate variables and fetch missing time=t values for particle A and B
263 Float4 vel_b = dev_vel[idx_b_orig];
264 Float4 angvel4_b = dev_angvel[idx_b_orig];
265
266 // Fetch previous sum of shear displacement for the contact pair
267 Float4 delta_t0_4 = dev_delta_t[mempos];
268
269 Float3 delta_t0_uncor = MAKE_FLOAT3(
270 delta_t0_4.x,
271 delta_t0_4.y,
272 delta_t0_4.z);
273
274 // Convert to Float3
275 Float3 angvel_b = MAKE_FLOAT3(
276 angvel4_b.x,
277 angvel4_b.y,
278 angvel4_b.z);
279
280 // Force between grain pair decomposed into normal- and tangential part
281 Float3 f_n, f_t, f_c;
282
283 // Normal vector of contact
284 Float3 n = x / x_length;
285
286 // Relative contact interface velocity, w/o rolling
287 Float3 vel_linear = MAKE_FLOAT3(
288 vel_a.x - vel_b.x,
289 vel_a.y - vel_b.y,
290 vel_a.z - vel_b.z);
291
292 // Relative contact interface velocity of particle surfaces at
293 // the contact, with rolling (Hinrichsen and Wolf 2004, eq. 13.10,
294 // or Luding 2008, eq. 10)
295 Float3 vel = vel_linear
296 + (radius_a + delta/2.0) * cross(n, angvel_a)
297 + (radius_b + delta/2.0) * cross(n, angvel_b);
298
299 // Relative contact interface rolling velocity
300 Float3 angvel = angvel_a - angvel_b;
301 Float angvel_length = length(angvel);
302
303 // Normal component of the relative contact interface velocity
304 Float vel_n = -dot(vel_linear, n);
305
306 // Tangential component of the relative contact interface velocity
307 // Hinrichsen and Wolf 2004, eq. 13.9
308 Float3 vel_t = vel - n * dot(n, vel);
309 Float vel_t_length = length(vel_t);
310
311 // Correct tangential displacement vector, which is
312 // necessary if the tangential plane rotated
313 Float3 delta_t0 = delta_t0_uncor - (n * dot(n, delta_t0_uncor));
314 Float delta_t0_length = length(delta_t0);
315
316 // New tangential displacement vector
317 Float3 delta_t;
318
319 // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
320 // 2013) based on Young's modulus if params.E > 0.0.
321 // Use the calculated stiffness for normal and tangential components.
322 Float k_n = devC_params.k_n;
323 Float k_t = devC_params.k_t;
324 if (devC_params.E > .001) {
325 k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
326 k_t = k_n;
327 }
328
329 // Normal force component: Elastic - viscous damping
330 f_n = fmax(0.0, -k_n*delta + devC_params.gamma_n * vel_n) * n;
331 Float f_n_length = length(f_n);
332
333 // Store energy dissipated in normal viscous component
334 // watt = gamma_n * vel_n * dx_n / dt
335 // watt = gamma_n * vel_n * vel_n * dt / dt
336 // watt = gamma_n * vel_n * vel_n
337 // watt = N*m/s = N*s/m * m/s * m/s * s / s
338 *ev_dot += 0.5 * devC_params.gamma_n * vel_n * vel_n;
339
340 // Add max. capillary force
341 f_c = -devC_params.kappa * sqrtf(radius_a * radius_b) * n;
342
343 // Initialize force vectors to zero
344 f_t = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
345
346 // Apply a tangential force if the previous tangential displacement
347 // is non-zero, or the current sliding velocity is non-zero.
348 if (delta_t0_length > 0.0 || vel_t_length > 0.0) {
349
350 // Add tangential displacement to total tangential displacement
351 delta_t = delta_t0 + vel_t * devC_dt;
352
353 // Tangential force: Visco-Elastic, before limitation criterion
354 Float3 f_t_elast = -k_t * delta_t;
355 Float3 f_t_visc = -devC_params.gamma_t * vel_t;
356 f_t = f_t_elast + f_t_visc;
357 Float f_t_length = length(f_t);
358
359 // Static frictional limit
360 Float f_t_limit = devC_params.mu_s * length(f_n-f_c);
361
362 // Store energy dissipated in tangential viscous component
363 *ev_dot += 0.5 * devC_params.gamma_t * vel_t_length * vel_t_length;
364
365 // If failure criterion is not met, contact is viscous-linear elastic.
366 // If failure criterion is met, contact force is limited,
367 // resulting in a slip and energy dissipation
368 if (f_t_length > f_t_limit) { // Static friciton exceeded: Dynamic case
369
370 // tangential vector
371 Float3 t = f_t/length(f_t);
372
373 // Frictional force is reduced to equal the dynamic limit
374 f_t = f_t_limit * t;
375
376 // In the sliding friction case, the tangential spring is adjusted
377 // to a length consistent with Coulombs (dynamic) condition (Luding
378 // 2008)
379 delta_t = -1.0/k_t
380 * (devC_params.mu_d * length(f_n-f_c) * t
381 + devC_params.gamma_t * vel_t);
382
383 // Shear friction heat production rate:
384 // The energy lost from the tangential spring is dissipated as heat
385 *es_dot += 0.5*length(length(f_t) * vel_t * devC_dt) / devC_dt; // Seen in ESyS-Particle
386
387 }
388 }
389
390
391 // Add force components from this collision to total force for particle
392 *F += f_n + f_t + f_c;
393
394 // Add torque components from this collision to total torque for particle
395 // Comment out the line below to disable rotation
396 *T += cross(-(radius_a + delta*0.5) * n, f_t);
397
398 // Pressure excerted onto the particle from this contact
399 *p += f_n_length / (4.0f * PI * radius_a*radius_a);
400
401 // Store sum of tangential displacements
402 dev_delta_t[mempos] = MAKE_FLOAT4(
403 delta_t.x,
404 delta_t.y,
405 delta_t.z,
406 0.0f);
407
408 } // End of contactLinear()
409
410
411 // Non-linear contact model for particle-particle interactions
412 // Based on Hertzian and Mindlin contact theories (e.g. Hertz, 1882, Mindlin and
413 // Deresiewicz, 1953, Johnson, 1985). See Yohannes et al 2012 for example.
414 __device__ void contactHertz(
415 Float3* F,
416 Float3* T,
417 Float* es_dot,
418 Float* ev_dot,
419 Float* p,
420 const unsigned int idx_a_orig,
421 const unsigned int idx_b_orig,
422 const Float4 vel_a,
423 const Float4* __restrict__ dev_vel,
424 const Float3 angvel_a,
425 const Float4* __restrict__ dev_angvel,
426 const Float radius_a,
427 const Float radius_b,
428 const Float3 x_ab,
429 const Float x_ab_length,
430 const Float delta_ab,
431 Float4* __restrict__ dev_delta_t,
432 const unsigned int mempos)
433 {
434
435 // Allocate variables and fetch missing time=t values for particle A and B
436 Float4 vel_b = dev_vel[idx_b_orig];
437 Float4 angvel4_b = dev_angvel[idx_b_orig];
438
439 // Fetch previous sum of shear displacement for the contact pair
440 Float4 delta_t0_4 = dev_delta_t[mempos];
441
442 Float3 delta_t0_uncor = MAKE_FLOAT3(delta_t0_4.x,
443 delta_t0_4.y,
444 delta_t0_4.z);
445
446 // Convert to Float3
447 Float3 angvel_b = MAKE_FLOAT3(angvel4_b.x, angvel4_b.y, angvel4_b.z);
448
449 // Force between grain pair decomposed into normal- and tangential part
450 Float3 f_n, f_t, f_c, T_res;
451
452 // Normal vector of contact
453 Float3 n_ab = x_ab/x_ab_length;
454
455 // Relative contact interface velocity, w/o rolling
456 Float3 vel_ab_linear = MAKE_FLOAT3(vel_a.x - vel_b.x,
457 vel_a.y - vel_b.y,
458 vel_a.z - vel_b.z);
459
460 // Relative contact interface velocity of particle surfaces at
461 // the contact, with rolling (Hinrichsen and Wolf 2004, eq. 13.10)
462 Float3 vel_ab = vel_ab_linear
463 + (radius_a + delta_ab/2.0f) * cross(n_ab, angvel_a)
464 + (radius_b + delta_ab/2.0f) * cross(n_ab, angvel_b);
465
466 // Relative contact interface rolling velocity
467 Float3 angvel_ab = angvel_a - angvel_b;
468 Float angvel_ab_length = length(angvel_ab);
469
470 // Normal component of the relative contact interface velocity
471 Float vel_n_ab = dot(vel_ab_linear, n_ab);
472
473 // Tangential component of the relative contact interface velocity
474 // Hinrichsen and Wolf 2004, eq. 13.9
475 Float3 vel_t_ab = vel_ab - (n_ab * dot(vel_ab, n_ab));
476 Float vel_t_ab_length = length(vel_t_ab);
477
478 // Correct tangential displacement vector, which is
479 // necessary if the tangential plane rotated
480 Float3 delta_t0 = delta_t0_uncor - (n_ab * dot(delta_t0_uncor, n_ab));
481 Float delta_t0_length = length(delta_t0);
482
483 // New tangential displacement vector
484 Float3 delta_t;
485
486 // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
487 // 2013) based on Young's modulus if params.E > 0.0.
488 // Use the calculated stiffness for normal and tangential components.
489 Float k_n = devC_params.k_n;
490 Float k_t = devC_params.k_t;
491 if (devC_params.E > .001) {
492 k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
493 k_t = k_n;
494 }
495
496 // Normal force component
497 f_n = (-k_n * powf(delta_ab, 3.0f/2.0f) -devC_params.gamma_n *
498 powf(delta_ab, 1.0f/4.0f) * vel_n_ab)
499 * n_ab;
500
501 // Store energy dissipated in normal viscous component
502 // watt = gamma_n * vel_n * dx_n / dt
503 // watt = gamma_n * vel_n * vel_n * dt / dt
504 // watt = gamma_n * vel_n * vel_n
505 // watt = N*m/s = N*s/m * m/s * m/s * s / s
506 *ev_dot += devC_params.gamma_n * vel_n_ab * vel_n_ab;
507
508 // Make sure the viscous damping doesn't exceed the elastic component,
509 // i.e. the damping factor doesn't exceed the critical damping, 2*sqrt(m*k_n)
510 if (dot(f_n, n_ab) < 0.0f)
511 f_n = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
512
513 Float f_n_length = length(f_n);
514
515 // Add max. capillary force
516 f_c = -devC_params.kappa * sqrtf(radius_a * radius_b) * n_ab;
517
518 // Initialize force vectors to zero
519 f_t = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
520 T_res = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
521
522 // Apply a tangential force if the previous tangential displacement
523 // is non-zero, or the current sliding velocity is non-zero.
524 if (delta_t0_length > 0.f || vel_t_ab_length > 0.f) {
525
526 // Shear force: Visco-Elastic, limited by Coulomb friction
527 Float3 f_t_elast = -k_t * powf(delta_ab, 1.0f/2.0f) * delta_t0;
528 Float3 f_t_visc = -devC_params.gamma_t * powf(delta_ab, 1.0f/4.0f) * vel_t_ab;
529 Float f_t_limit;
530
531 if (vel_t_ab_length > 0.001f) { // Dynamic friciton
532 f_t_limit = devC_params.mu_d * length(f_n-f_c);
533 } else { // Static friction
534 f_t_limit = devC_params.mu_s * length(f_n-f_c);
535 }
536
537 // Tangential force before friction limit correction
538 f_t = f_t_elast + f_t_visc;
539 Float f_t_length = length(f_t);
540
541 // If failure criterion is not met, contact is viscous-linear elastic.
542 // If failure criterion is met, contact force is limited,
543 // resulting in a slip and energy dissipation
544 if (f_t_length > f_t_limit) { // Dynamic case
545
546 // Frictional force is reduced to equal the limit
547 f_t *= f_t_limit/f_t_length;
548
549 // A slip event zeros the displacement vector
550 //delta_t = make_Float3(0.0f, 0.0f, 0.0f);
551
552 // In a slip event, the tangential spring is adjusted to a
553 // length which is consistent with Coulomb's equation
554 // (Hinrichsen and Wolf, 2004)
555 delta_t = (f_t + devC_params.gamma_t * powf(delta_ab, 1.0f/4.0f) * vel_t_ab)
556 / (-k_t * powf(delta_ab, 1.0f/2.0f));
557
558 // Shear friction heat production rate:
559 // The energy lost from the tangential spring is dissipated as heat
560 *es_dot += length(delta_t0 - delta_t) * k_t / devC_dt; // Seen in EsyS-Particle
561
562 } else { // Static case
563
564 // No correction of f_t is required
565
566 // Add tangential displacement to total tangential displacement
567 delta_t = delta_t0 + vel_t_ab * devC_dt;
568 }
569 }
570
571 if (angvel_ab_length > 0.f) {
572 // Apply rolling resistance (Zhou et al. 1999)
573 //T_res = -angvel_ab/angvel_ab_length * devC_params.mu_r * R_bar * length(f_n);
574
575 // New rolling resistance model
576 T_res = -1.0f * fmin(devC_params.gamma_r * radius_a * angvel_ab_length,
577 devC_params.mu_r * radius_a * f_n_length)
578 * angvel_ab/angvel_ab_length;
579 }
580
581 // Add force components from this collision to total force for particle
582 *F += f_n + f_t + f_c;
583 *T += -(radius_a + delta_ab/2.0f) * cross(n_ab, f_t) + T_res;
584
585 // Pressure excerted onto the particle from this contact
586 *p += f_n_length / (4.0f * PI * radius_a*radius_a);
587
588 // Store sum of tangential displacements
589 dev_delta_t[mempos] = MAKE_FLOAT4(delta_t.x, delta_t.y, delta_t.z, 0.0f);
590
591 } // End of contactHertz()
592
593
594 #endif
595 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4