1 /*
2 Copyright (C) 1998 Dennis Roddeman
3 email: dennis.roddeman@feat.nl
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software Foundation
17 59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
18 */
19
20 #include "tochnog.h"
21
22 #define EPS_SIG 1.e-5
23 #define EPS_DIR 1.e-10
24 #define EPS_PRISCO_1 1.e-8
25 #define EPS_PRISCO_2 1.e-5
26
plasti_rule(long int element,long int gr,long int plasti_on_boundary,double user_data[],double new_unknowns[],double new_grad_new_unknowns[],double old_hisv[],double new_hisv[],double old_epp[],double inc_epp[],double inc_ept[],long int task,long int & plasti_type,double sig[],double & f,double & new_f,double dir[])27 void plasti_rule( long int element, long int gr,
28 long int plasti_on_boundary, double user_data[],
29 double new_unknowns[], double new_grad_new_unknowns[],
30 double old_hisv[], double new_hisv[], double old_epp[],
31 double inc_epp[], double inc_ept[], long int task, long int &plasti_type,
32 double sig[], double &f, double &new_f, double dir[] )
33
34 {
35
36 long int i=0, swit=0, idim=0, jdim=0, ind=0, tmp_plasti_type=-1, in_apex=0,
37 length=0, group_materi_plasti_user=-NO,
38 test1=0, test2=0, test3=0, test4=0, test5=0,
39 test6=0, test7=0, test8=0, test9=0, ldum=0, icontrol=0,
40 tension_cutoff=-NO, tension_limit=-NO,
41 options_skip_plasticity=-NO, idum[1];
42 double sig_yield=0., sigm=0., signorm=0., K=0., K_flow=0., alpha=0.,
43 dev=0., devp=0., q=0., p=0., C=0., kappa_0=0., n=0., options_nonlocal=0.,
44 f_flow_right=0., f_flow_left=0., phi=0., c=0., sig0=0.,
45 phi_flow=0., alpha_flow=0., f_yield=NO_YIELD_F, f_flow=NO_YIELD_F,
46 sig_eq=0., tmp_sig_eq=0., void_fraction=0., tmp=0., tmp1=0., tmp2=0.,
47 q1=0., q2=0., q3=0., m1=0., m20=0., m21=0., m22=0., m23=0.,
48 prisco_gamma=0., prisco_j2eta=0., prisco_j3eta=0.,
49 prisco_r=0., prisco_betaf=0., prisco_rc=0., prisco_cp=0.,
50 prisco_kinc=0., prisco_bp=0., prisco_ksi=0., prisco_ksi_c=0.,
51 prisco_ksi_e=0., prisco_thetahat=1.,
52 prisco_theta=0., prisco_thetahat_c=0., prisco_thetahat_e=0.,
53 prisco_tp=0., prisco_betafhat=0.,
54 prisco_r0=0., prisco_st0=0., prisco_b4=0., prisco_b5=0., prisco_b7=0.,
55 prisco_beta=0., prisco_betaf0=0., prisco_r1tmp=0.,
56 prisco_r2tmp=0., prisco_rt=0.,
57 pa=0., pb=0., R=0., t=0., epp_vol=0.,
58 phi0=0., c0=0., phi0_flow=0., phi1=0., c1=0.,
59 phi1_flow=0., kappa=0., kappa_crit=0., I1=0., I2=0., I3=0.,
60 plasti_on_boundary_factor=BOUNDARY_REDUCTION_FACTOR,
61 m=0., lambda=0., e=0., N=0., de=0., p0=0., dp0=0., rdum=0.,
62 prisco_rv[MDIM][MDIM], prisco_st[MDIM][MDIM], prisco_rr[MDIM][MDIM],
63 prisco_st1[MDIM][MDIM], prisco_st2[MDIM][MDIM], prisco_sv0[MDIM][MDIM],
64 prisco_chi[MDIM*MDIM], prisco_chihat[MDIM*MDIM],
65 prisco_ss[MDIM*MDIM], prisco_etas[MDIM*MDIM], prisco_epinc[MDIM*MDIM],
66 prisco_epp_dev[MDIM*MDIM], prisco_incepp_dev[MDIM*MDIM],
67 sig_princ_apex[MDIM], sig_dev[MDIM*MDIM], sig_tmp[MDIM*MDIM],
68 plasti_data[DATA_ITEM_SIZE], sig_princ[MDIM], ddum[MDIM*MDIM],
69 e_vec[MDIM], f_vec[MDIM], work[MDIM*MDIM],
70 geo_sig[MDIM*MDIM], geo_inc_epp[MDIM*MDIM],
71 geo_old_epp[MDIM*MDIM], geo_inc_ept[MDIM*MDIM];
72
73 db( OPTIONS_SKIP_PLASTICITY, 0, &options_skip_plasticity, ddum, ldum,
74 VERSION_NORMAL, GET_IF_EXISTS );
75 db( ICONTROL, 0, &icontrol, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
76 db( CONTROL_OPTIONS_SKIP_PLASTICITY, icontrol, &options_skip_plasticity, ddum, ldum,
77 VERSION_NORMAL, GET_IF_EXISTS );
78 if ( options_skip_plasticity==-YES ) return;
79
80 swit = set_swit(element,-1,"plasti_rule");
81 if ( task==GET_FLOW_RULE ) swit = 0;
82 if ( swit ) pri( "In PLASTI_RULE" );
83
84 if ( swit ) {
85 if ( task==GET_YIELD_RULE )
86 pri( "GET_YIELD_RULE" );
87 else if ( task==GET_FLOW_RULE ) {
88 pri( "GET_FLOW_RULE" );
89 pri( "plasti_type", plasti_type );
90 }
91 else if ( task==GET_FLOW_RULE_GRAD ) {
92 pri( "GET_FLOW_RULE_GRAD" );
93 pri( "plasti_type", plasti_type );
94 }
95 pri( "sig", sig, MDIM, MDIM );
96 }
97
98 db( OPTIONS_NONLOCAL, 0, idum, &options_nonlocal, ldum,
99 VERSION_NORMAL, GET_IF_EXISTS );
100 db( GROUP_MATERI_PLASTI_USER, 0, &group_materi_plasti_user, ddum, ldum,
101 VERSION_NORMAL, GET_IF_EXISTS );
102 db( GROUP_MATERI_PLASTI_BOUNDARY_FACTOR, gr, idum, &plasti_on_boundary_factor,
103 ldum, VERSION_NORMAL, GET_IF_EXISTS );
104
105 // mean stress
106 sigm = ( sig[0] + sig[4] + sig[8] ) / 3.;
107 if ( swit ) pri( "sigm", sigm );
108
109 // deviatoric strains and stresses
110 array_move( sig, sig_dev, MDIM*MDIM );
111 for ( idim=0; idim<MDIM; idim++ ) sig_dev[idim*MDIM+idim] -= sigm;
112
113 // equivalent stress
114 sig_eq = sqrt(scalar_dabs(array_inproduct(sig_dev,sig_dev,MDIM*MDIM))/2.);
115 if ( swit ) pri( "sig_eq", sig_eq );
116
117 // geoetechnical convention
118 array_multiply( sig, geo_sig, -1., MDIM*MDIM );
119 array_multiply( inc_ept, geo_inc_ept, -1., MDIM*MDIM );
120 array_multiply( inc_epp, geo_inc_epp, -1., MDIM*MDIM );
121 if ( materi_strain_plasti ) array_multiply( old_epp, geo_old_epp, -1., MDIM*MDIM );
122
123 f = NO_YIELD_F;
124 if ( get_group_data( GROUP_MATERI_PLASTI_CAMCLAY, gr, element, new_unknowns,
125 plasti_data, ldum, GET_IF_EXISTS ) ) {
126 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
127 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_CAMCLAY;
128 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_CAMCLAY;
129 if ( test1 || test2 || test3 ) {
130 if ( swit ) pri( "plasti_camclay" );
131 if ( materi_history_variables<2 ) {
132 pri( "Error: GROUP_MATERI_PLASTI_CAMCLAY needs 'materi_history_variables 2'" );
133 exit(1);
134 }
135 m = plasti_data[0];
136 kappa = plasti_data[1];
137 lambda = plasti_data[2];
138 N = plasti_data[3];
139 if ( scalar_dabs(lambda-kappa)==0. ) db_error( GROUP_MATERI_PLASTI_CAMCLAY, gr );
140 e = old_hisv[0];
141 //p0 = old_hisv[1];
142 p = -sigm;
143 p0 = exp((N-kappa*log(p)-(1+e))/(lambda-kappa));
144 if ( swit ) {
145 pri( "e", e );
146 pri( "p0", p0 );
147 }
148 dev = (geo_inc_ept[0]+geo_inc_ept[4]+geo_inc_ept[8]);
149 de = -dev*(1.+e);
150 devp = (geo_inc_epp[0]+geo_inc_epp[4]+geo_inc_epp[8]);
151 dp0 = devp * p0 * (1.+e)/(lambda-kappa);
152 if ( swit ) {
153 pri( "dev", dev );
154 pri( "devp", devp );
155 pri( "de", de );
156 pri( "dp0", dp0 );
157 }
158 e += de;
159 p0 += dp0;
160 q = sqrt(
161 0.5*( scalar_square(geo_sig[0]-geo_sig[4])+
162 scalar_square(geo_sig[4]-geo_sig[8])+
163 scalar_square(geo_sig[0]-geo_sig[8]) ) +
164 3.*( scalar_square(geo_sig[1]) +
165 scalar_square(geo_sig[2]) +
166 scalar_square(geo_sig[5]) ) );
167 f_yield = q*q - m*m*(p*(p0-p));
168 f_flow = f_yield;
169 if ( task==GET_YIELD_RULE ) {
170 if ( f_yield>f ) {
171 f = f_yield;
172 tmp_plasti_type = GROUP_MATERI_PLASTI_CAMCLAY;
173 }
174 if ( swit ) pri( "f_yield", f_yield );
175
176 new_hisv[0] = e;
177 new_hisv[1] = p0;
178 }
179 else {
180 assert( task==GET_FLOW_RULE );
181 f = f_flow;
182 if ( swit ) pri( "f_flow", f_flow );
183 }
184 }
185 }
186 if ( get_group_data( GROUP_MATERI_PLASTI_CAP, gr, element, new_unknowns,
187 plasti_data, length, GET_IF_EXISTS ) ) {
188 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
189 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_CAP;
190 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_CAP;
191 if ( test1 || test2 || test3 ) {
192 if ( swit ) pri( "check plasti_cap" );
193 /* get data */
194 c = plasti_data[0];
195 phi = plasti_data[1];
196 alpha = plasti_data[2];
197 R = plasti_data[3];
198 epp_vol = scalar_dabs( new_unknowns[epp_indx+0] + new_unknowns[epp_indx+3] +
199 new_unknowns[epp_indx+5] );
200 /* epp_vol-pb table */
201 length -= 4;
202 if ( !table_xy( &plasti_data[4], "GROUP_MATERI_PLASTI_CAP",
203 length, epp_vol, pb ) ) db_error( GROUP_MATERI_PLASTI_CAP, gr );
204 /* invariants */
205 p = -sigm;
206 t = sqrt(3.) * sig_eq;
207 /* pa */
208 tmp = (1.+R*tan(phi));
209 if ( tmp==0. ) db_error( GROUP_MATERI_PLASTI_CAP, gr );
210 pa = (pb-R*c)/tmp;
211 /* yield rule and flow rule */
212 if ( cos(phi)==0. ) db_error( GROUP_MATERI_PLASTI_CAP, gr );
213 tmp = 1. + alpha - (alpha/cos(phi));
214 if ( tmp==0. ) db_error( GROUP_MATERI_PLASTI_CAP, gr );
215 f_yield = sqrt( scalar_square(p-pa) + scalar_square(R*t/tmp) ) -
216 R*(c+pa*tan(phi));
217 f_flow = f_yield;
218 if ( task==GET_YIELD_RULE ) {
219 if ( f_yield>f ) {
220 f = f_yield;
221 tmp_plasti_type = GROUP_MATERI_PLASTI_CAP;
222 }
223 if ( swit ) pri( "f_yield", f_yield );
224 }
225 else {
226 assert( task==GET_FLOW_RULE );
227 f = f_flow;
228 if ( swit ) pri( "f_flow", f_flow );
229 }
230 }
231 }
232 if ( get_group_data( GROUP_MATERI_PLASTI_COMPRESSION, gr, element, new_unknowns,
233 plasti_data, ldum, GET_IF_EXISTS ) ) {
234 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
235 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_COMPRESSION;
236 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_COMPRESSION;
237 if ( test1 || test2 || test3 ) {
238 if ( swit ) pri( "check plasti_compression" );
239 sig_yield = plasti_data[0];
240 if ( swit ) pri( "sig_yield", sig_yield );
241 matrix_eigenvalues( sig, sig_princ );
242 if ( swit ) pri( "sig_princ", sig_princ, MDIM );
243 tmp_sig_eq = 0.;
244 for ( idim=0; idim<MDIM; idim++ ) {
245 if ( sig_princ[idim]<0. ) {
246 tmp_sig_eq += sig_princ[idim] * sig_princ[idim];
247 }
248 }
249 tmp_sig_eq = sqrt( scalar_dabs(tmp_sig_eq) );
250 f_yield = tmp_sig_eq - sig_yield;
251 f_flow = tmp_sig_eq - sig_yield;
252 if ( task==GET_YIELD_RULE ) {
253 if ( f_yield>f ) {
254 f = f_yield;
255 tmp_plasti_type = GROUP_MATERI_PLASTI_COMPRESSION;
256 }
257 if ( swit ) pri( "f_yield", f_yield );
258 }
259 else {
260 assert( task==GET_FLOW_RULE );
261 f = f_flow;
262 if ( swit ) pri( "f_flow", f_flow );
263 }
264 }
265 }
266 if ( get_group_data( GROUP_MATERI_PLASTI_DIPRISCO, gr, element, new_unknowns,
267 plasti_data, ldum, GET_IF_EXISTS ) ) {
268 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
269 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_DIPRISCO;
270 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_DIPRISCO;
271 if ( test1 || test2 || test3 ) {
272 get_group_data( GROUP_MATERI_PLASTI_DIPRISCO_RT, gr, element, new_unknowns,
273 &prisco_rt, ldum, GET_IF_EXISTS );
274 // geotechnical notation
275 array_multiply( sig, geo_sig, -1., MDIM*MDIM );
276 array_multiply( inc_epp, geo_inc_epp, -1., MDIM*MDIM );
277 array_multiply( old_epp, geo_old_epp, -1., MDIM*MDIM );
278 if ( swit ) pri( "plasti_diprisco" );
279 assert( materi_history_variables>=11 );
280 sigm = ( geo_sig[0] + geo_sig[4] + geo_sig[8] ) / 3.;
281 array_move( geo_sig, sig_dev, MDIM*MDIM );
282 for ( idim=0; idim<MDIM; idim++ ) sig_dev[idim*MDIM+idim] -= sigm;
283 if ( swit ) pri( "sig_dev", sig_dev, MDIM, MDIM );
284 // history variables
285 array_move( old_hisv, prisco_chi, MDIM*MDIM );
286 tmp = array_size( prisco_chi, MDIM*MDIM );
287 assert( tmp!=0. );
288 array_multiply( prisco_chi, prisco_chi, 1./tmp, MDIM*MDIM );
289 prisco_beta = old_hisv[MDIM*MDIM];
290 prisco_rc = - old_hisv[MDIM*MDIM+1];
291 // user data
292 prisco_gamma = plasti_data[0];
293 prisco_betafhat = plasti_data[1];
294 prisco_bp = plasti_data[2];
295 prisco_cp = plasti_data[3];
296 prisco_tp = plasti_data[4];
297 prisco_thetahat_c = plasti_data[5];
298 prisco_thetahat_e = plasti_data[6];
299 prisco_ksi_c = plasti_data[7];
300 prisco_ksi_e = plasti_data[8];
301 prisco_betaf0 = plasti_data[9];
302 // fixed data
303 prisco_r = array_inproduct( geo_sig, prisco_chi, MDIM*MDIM );
304 if ( prisco_bp<=0. ) db_error( GROUP_MATERI_PLASTI_DIPRISCO, gr );
305 if ( swit ) {
306 pri( "prisco_chi", prisco_chi, MDIM*MDIM );
307 pri( "prisco_betaf", prisco_betaf );
308 pri( "prisco_rc", prisco_rc );
309 pri( "prisco_gamma", prisco_gamma );
310 pri( "prisco_betafhat", prisco_betafhat );
311 pri( "prisco_bp", prisco_bp );
312 pri( "prisco_cp", prisco_cp );
313 pri( "prisco_tp", prisco_tp );
314 pri( "prisco_thetahat_c", prisco_thetahat_c );
315 pri( "prisco_thetahat_e", prisco_thetahat_e );
316 pri( "prisco_ksi_c", prisco_ksi_c );
317 pri( "prisco_ksi_e", prisco_ksi_e );
318 pri( "prisco_r", prisco_r );
319 }
320 for ( idim=0; idim<MDIM; idim++ ) {
321 for ( jdim=0; jdim<MDIM; jdim++ ) {
322 ind = idim*MDIM+jdim;
323 prisco_ss[ind] = geo_sig[ind] - prisco_r * prisco_chi[ind];
324 if ( scalar_dabs(prisco_r)<EPS_PRISCO_1 ) {
325 prisco_etas[ind] = 0.;
326 }
327 else
328 prisco_etas[ind] = prisco_ss[ind] * sqrt(3.) / (prisco_r+prisco_rt);
329 }
330 }
331 if ( swit ) {
332 pri( "prisco_ss", prisco_ss, MDIM, MDIM );
333 pri( "prisco_etas", prisco_etas, MDIM, MDIM );
334 }
335 prisco_j2eta = array_inproduct( prisco_etas, prisco_etas, MDIM*MDIM );
336 matrix_ab( prisco_etas, prisco_etas, work, MDIM, MDIM, MDIM );
337 prisco_j3eta = array_inproduct( work, prisco_etas, MDIM*MDIM );
338 if ( swit ) {
339 pri( "prisco_j2eta", prisco_j2eta );
340 pri( "prisco_j3eta", prisco_j3eta );
341 }
342 // betaf calculation
343 prisco_betaf = prisco_betafhat + ( prisco_betaf0 - prisco_betafhat ) *
344 exp( - prisco_beta * prisco_tp );
345 // second yield function
346 prisco_r1tmp = (prisco_rc/sqrt(3.))*EPS_PRISCO_2;
347 prisco_r2tmp = 0.;
348 for ( idim=0; idim<MDIM; idim++ ) {
349 for ( jdim=0; jdim<MDIM; jdim++ ) {
350 ind = idim*MDIM+jdim;
351 if ( prisco_etas[ind]>prisco_r2tmp ) prisco_r2tmp = prisco_etas[ind];
352 }
353 }
354 if ( prisco_r<prisco_r1tmp && task==GET_YIELD_RULE && prisco_rt==0. ) {
355 f_yield = NO_YIELD_F;
356 f_flow = NO_YIELD_F;
357 }
358 else if ( prisco_r2tmp >= 2.7 && task==GET_YIELD_RULE ) {
359 f_yield = NO_YIELD_F;
360 f_flow = NO_YIELD_F;
361 }
362 else {
363 // normal yield rule and flow rule
364 tmp = (prisco_r+prisco_rt)/(prisco_rc+prisco_rt);
365 f_yield = 3. * prisco_betaf * ( prisco_gamma - 3. ) *
366 log ( scalar_dabs(tmp) ) -
367 prisco_gamma * prisco_j3eta + (9./4.) * ( prisco_gamma - 1. ) * prisco_j2eta;
368 f_flow = 9. * ( prisco_gamma - 3. ) *
369 log ( scalar_dabs(tmp) ) -
370 prisco_gamma * prisco_j3eta + (9./4.) * ( prisco_gamma - 1. ) * prisco_j2eta;
371 }
372 tmp = ( prisco_etas[0] + prisco_etas[4] + prisco_etas[8] ) ;
373 // load angle
374 array_set( &prisco_st1[0][0], 0., MDIM*MDIM );
375 array_set( &prisco_st2[0][0], 0., MDIM*MDIM );
376 prisco_st1[0][0] = 0.8164965809;
377 prisco_st1[1][1] = - prisco_st1[0][0] / 2.;
378 prisco_st1[2][2] = - prisco_st1[0][0] / 2.;
379 tmp = array_size( geo_sig, MDIM*MDIM );
380 if ( tmp!=0. )
381 array_multiply( geo_sig, &prisco_sv0[0][0], 1./tmp, MDIM*MDIM );
382 else
383 array_set( &prisco_sv0[0][0], 0., MDIM*MDIM );
384
385 prisco_r0 = prisco_chi[0]+prisco_chi[4]+prisco_chi[8];
386 prisco_st0 = prisco_sv0[0][0]+prisco_sv0[1][1]+
387 prisco_sv0[2][2];
388
389 for ( idim=0; idim<MDIM; idim++ ) {
390 for ( jdim=0; jdim<MDIM; jdim++ ) {
391 ind = idim*MDIM+jdim;
392 prisco_rv[idim][jdim] = prisco_chi[ind];
393 prisco_st[idim][jdim] = prisco_sv0[idim][jdim];
394 if ( idim==jdim ) {
395 prisco_rv[idim][jdim] -= prisco_r0/3.;
396 prisco_st[idim][jdim] -= prisco_st0/3.;
397 }
398 prisco_rr[idim][jdim] = prisco_st[idim][jdim] - prisco_rv[idim][jdim];
399 }
400 }
401 if ( swit ) {
402 pri( "prisco_betaf", prisco_betaf );
403 pri( "prisco_r0", prisco_r0 );
404 pri( "prisco_st0", prisco_st0 );
405 pri( "prisco_sv0", &prisco_sv0[0][0], MDIM*MDIM );
406 pri( "prisco_rv", &prisco_rv[0][0], MDIM*MDIM );
407 pri( "prisco_st", &prisco_st[0][0], MDIM*MDIM );
408 pri( "prisco_rr", &prisco_rr[0][0], MDIM*MDIM );
409 }
410 tmp = array_size( &prisco_rr[0][0], MDIM*MDIM );
411 if ( tmp<EPS_PRISCO_1 ) tmp = EPS_PRISCO_1;
412 array_multiply( &prisco_rr[0][0], &prisco_rr[0][0], 1./tmp, MDIM*MDIM );
413 // alternative method
414 // matrix_eigenvalues( &prisco_rr[0][0], work );
415 // prisco_st2[0][0] = work[0];
416 // prisco_st2[1][1] = work[1];
417 // prisco_st2[2][2] = work[2];
418 tmp1 = (prisco_rr[0][0]+prisco_rr[1][1])/2.;
419 tmp2 = (prisco_rr[0][0]-prisco_rr[1][1])/2.;
420 tmp2 = tmp2*tmp2+prisco_rr[0][1]*prisco_rr[0][1];
421 tmp2 = sqrt(tmp2);
422 prisco_st2[0][0] = tmp1+tmp2;
423 prisco_st2[1][1] = tmp1-tmp2;
424 prisco_st2[2][2] = prisco_rr[2][2];
425
426 prisco_b4 = array_inproduct( &prisco_st2[0][0], &prisco_st1[0][0], MDIM*MDIM );
427 prisco_b7 = scalar_dabs(prisco_b4);
428 if ( prisco_b7>1. ) prisco_b7 = 1.;
429 if ( prisco_b4>1. ) prisco_b4 = 1.;
430 if ( prisco_b4<0. ) prisco_b4 = -prisco_b7;
431 prisco_b5 = acos(prisco_b4);
432 if ( prisco_b5<(1./3)*PIRAD ) {
433 prisco_b5 = scalar_dabs(prisco_b5);
434 }
435 else {
436 if ( prisco_b5<(3./3.)*PIRAD ) {
437 prisco_b5 = prisco_b5 - (2./3.)*PIRAD;
438 prisco_b5 = scalar_dabs(prisco_b5);
439 }
440 else {
441 if ( prisco_b5<(5./3.)*PIRAD ) {
442 prisco_b5 = prisco_b5 - (4./3.)*PIRAD;
443 prisco_b5 = scalar_dabs(prisco_b5);
444 }
445 }
446 }
447 // ksi, theta
448 prisco_ksi = 0.5*(prisco_ksi_c+prisco_ksi_e) +
449 0.5*(prisco_ksi_c-prisco_ksi_e)*(2.*prisco_b5*prisco_b5-4.*prisco_b5+1.);
450 prisco_thetahat = 0.5*(prisco_thetahat_c+prisco_thetahat_e) +
451 0.5*(prisco_thetahat_c-prisco_thetahat_e)*(2.*prisco_b5*prisco_b5-4.*prisco_b5+1.);
452 prisco_theta = prisco_thetahat;
453 if ( swit ) {
454 pri( "prisco_b4", prisco_b4 );
455 pri( "prisco_b5", prisco_b5 );
456 pri( "prisco_b7", prisco_b7 );
457 pri( "prisco_kinc", prisco_kinc );
458 pri( "prisco_ksi", prisco_ksi );
459 pri( "prisco_thetahat", prisco_thetahat );
460 pri( "prisco_theta", prisco_theta );
461 }
462 tmp1 = array_inproduct( geo_inc_epp, prisco_chi, MDIM*MDIM );
463 for ( idim=0; idim<MDIM; idim++ ) {
464 for ( jdim=0; jdim<MDIM; jdim++ ) {
465 ind = idim*MDIM+jdim;
466 prisco_epinc[ind] = geo_inc_epp[ind] - tmp1 * prisco_chi[ind];
467 }
468 }
469 tmp2 = prisco_ksi * array_size( prisco_epinc, MDIM*MDIM );
470 prisco_rc += (prisco_rc/prisco_bp) * ( tmp1 + tmp2 );
471 if (prisco_rc <= EPS_PRISCO_2) prisco_rc = EPS_PRISCO_2;
472 if ( swit ) {
473 pri( "prisco_epinc", prisco_epinc, MDIM, MDIM );
474 pri( "tmp1", tmp1 );
475 pri( "tmp2", tmp2 );
476 pri( "prisco_rc", prisco_rc );
477 }
478 for ( idim=0; idim<MDIM; idim++ ) {
479 for ( jdim=0; jdim<MDIM; jdim++ ) {
480 ind = idim*MDIM+jdim;
481 prisco_chihat[ind] = prisco_rr[idim][jdim] * sin( prisco_theta );
482 if ( idim==jdim ) prisco_chihat[ind] += cos( prisco_theta ) / sqrt( 3. );
483 }
484 }
485 prisco_kinc = prisco_cp * array_size( geo_inc_epp, MDIM*MDIM );
486 for ( idim=0; idim<MDIM; idim++ ) {
487 for ( jdim=0; jdim<MDIM; jdim++ ) {
488 ind = idim*MDIM+jdim;
489 prisco_chi[ind] += prisco_kinc * ( prisco_chihat[ind] - prisco_chi[ind] );
490 }
491 }
492 array_normalize( prisco_chi, MDIM*MDIM );
493 if ( swit ) {
494 pri( "prisco_chihat", prisco_chihat, MDIM, MDIM );
495 pri( "prisco_chi", prisco_chi, MDIM, MDIM );
496 }
497 array_move( geo_old_epp, work, MDIM*MDIM );
498 array_add( work, geo_inc_epp, work, MDIM*MDIM );
499 tmp1 = (work[0] + work[4] + work[8])/3. ;
500 tmp2 = (geo_inc_epp[0] + geo_inc_epp[4] + geo_inc_epp[8])/3. ;
501 array_move( work, prisco_epp_dev, MDIM*MDIM );
502 array_move( geo_inc_epp, prisco_incepp_dev, MDIM*MDIM );
503 for ( idim=0; idim<MDIM; idim++ ) {
504 ind = idim*MDIM + idim;
505 prisco_epp_dev[ind] -= tmp1;
506 prisco_incepp_dev[ind] -= tmp2;
507 }
508
509 if ( swit ) {
510 pri( "prisco_epp_dev", prisco_epp_dev, MDIM*MDIM );
511 pri( "prisco_incepp_dev", prisco_incepp_dev, MDIM*MDIM );
512 }
513 tmp1 = array_inproduct( prisco_epp_dev, prisco_incepp_dev, MDIM*MDIM );
514 tmp2 = array_inproduct( prisco_epp_dev, prisco_epp_dev, MDIM*MDIM );
515 tmp2 = sqrt( scalar_dabs( tmp2 ) );
516 if ( tmp2<EPS_PRISCO_1 ) tmp2 = EPS_PRISCO_1;
517 tmp = tmp1 / tmp2;
518 prisco_beta += tmp;
519 if ( swit ) pri( "prisco_beta new", prisco_beta );
520 array_move( prisco_chi, new_hisv, MDIM*MDIM );
521 new_hisv[MDIM*MDIM] = prisco_beta;
522 new_hisv[MDIM*MDIM+1] = - prisco_rc;
523 if ( task==GET_YIELD_RULE ) {
524 if ( f_yield>f ) {
525 f = f_yield;
526 tmp_plasti_type = GROUP_MATERI_PLASTI_DIPRISCO;
527 }
528 if ( swit ) {
529 pri( "f_yield", f_yield );
530 pri( "new_hisv", new_hisv, materi_history_variables );
531 }
532 }
533 else {
534 assert( task==GET_FLOW_RULE );
535 f = f_flow;
536 if ( swit ) pri( "f_flow", f_flow );
537 }
538 }
539 }
540 if ( get_group_data( GROUP_MATERI_PLASTI_DRUCKPRAG, gr, element, new_unknowns,
541 plasti_data, ldum, GET_IF_EXISTS ) ) {
542 tension_cutoff = -NO;
543 db( GROUP_MATERI_PLASTI_DRUCKPRAG_TENSIONCUTOFF, 0, &tension_cutoff, ddum, ldum,
544 VERSION_NORMAL, GET_IF_EXISTS );
545 tension_limit = -NO;
546 db( GROUP_MATERI_PLASTI_DRUCKPRAG_TENSIONLIMIT, 0, &tension_limit, ddum, ldum,
547 VERSION_NORMAL, GET_IF_EXISTS );
548 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
549 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG;
550 test3 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG_APEX;
551 test4 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG;
552 test5 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG_APEX;
553 if ( test1 || test2 || test3 || test4 || test5 ) {
554 if ( swit ) pri( "check plasti_druckprag" );
555 phi = plasti_data[0];
556 phi_flow = plasti_data[2];
557 if ( plasti_on_boundary ) {
558 phi *= plasti_on_boundary_factor;
559 phi_flow *= plasti_on_boundary_factor;
560 }
561 c = plasti_data[1];
562 alpha = 2.*sin(phi)/((3.-sin(phi))*sqrt(3.));
563 alpha_flow = 2.*sin(phi_flow)/((3.-sin(phi_flow))*sqrt(3.));
564 K = 6.*c*cos(phi)/((3.-sin(phi))*sqrt(3.));
565 K_flow = 6.*c*cos(phi_flow)/((3.-sin(phi_flow))*sqrt(3.));
566 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE && phi>0.;
567 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG_APEX;
568 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG_APEX;
569 if ( test1 || test2 || test3 ) {
570 // separate yield and flow rule for apex
571 matrix_eigenvalues( sig, sig_princ );
572 array_set( sig_princ_apex, c*cos(phi)/sin(phi), MDIM );
573 array_set( e_vec, 1./sqrt(3.), MDIM );
574 array_subtract( sig_princ, sig_princ_apex, f_vec, MDIM );
575 in_apex = array_inproduct(f_vec,e_vec,MDIM)>0. && tension_cutoff==-YES;
576 if ( in_apex || plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG_APEX ) {
577 f_yield = array_size( f_vec, MDIM );
578 f_flow = array_size( f_vec, MDIM );
579 }
580 else {
581 f_yield = NO_YIELD_F;
582 f_flow = NO_YIELD_F;
583 }
584 if ( task==GET_YIELD_RULE ) {
585 if ( f_yield>f ) {
586 f = f_yield;
587 tmp_plasti_type = GROUP_MATERI_PLASTI_DRUCKPRAG_APEX;
588 }
589 if ( swit ) pri( "f_yield", f_yield );
590 }
591 else {
592 assert( task==GET_FLOW_RULE );
593 f = f_flow;
594 if ( swit ) pri( "f_flow", f_flow );
595 }
596 }
597 if ( !in_apex ) {
598 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
599 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG;
600 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_DRUCKPRAG;
601 if ( test1 || test2 || test3 ) {
602 if ( tension_limit==-YES && sigm>0. ) {
603 f_yield = sig_eq - K;
604 f_flow = sig_eq - K_flow;
605 }
606 else {
607 f_yield = 3.*alpha*sigm + sig_eq - K;
608 f_flow = 3.*alpha_flow*sigm + sig_eq - K_flow;
609 }
610 if ( task==GET_YIELD_RULE ) {
611 if ( f_yield>f ) {
612 f = f_yield;
613 tmp_plasti_type = GROUP_MATERI_PLASTI_DRUCKPRAG;
614 }
615 if ( swit ) pri( "f_yield", f_yield );
616 }
617 else {
618 assert( task==GET_FLOW_RULE );
619 f = f_flow;
620 if ( swit ) pri( "f_flow", f_flow );
621 }
622 }
623 }
624 }
625 }
626 if ( get_group_data( GROUP_MATERI_PLASTI_GURSON, gr, element, new_unknowns,
627 plasti_data, ldum, GET_IF_EXISTS ) ) {
628 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
629 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_GURSON;
630 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_GURSON;
631 if ( test1 || test2 || test3 ) {
632 if ( swit ) pri( "plasti_gurson" );
633 void_fraction = new_unknowns[void_indx];
634 sig_yield = plasti_data[0];
635 q1 = plasti_data[1];
636 q2 = plasti_data[2];
637 q3 = plasti_data[3];
638 f_yield = (3.*sig_eq*sig_eq)/(sig_yield*sig_yield) +
639 2.*q1*void_fraction*cosh(q2*3.*sigm/(2.*sig_yield)) -
640 (1.+(q3*void_fraction)*(q3*void_fraction) );
641 f_flow = f_yield;
642 if ( task==GET_YIELD_RULE ) {
643 if ( f_yield>f ) {
644 f = f_yield;
645 tmp_plasti_type = GROUP_MATERI_PLASTI_GURSON;
646 }
647 if ( swit ) pri( "f_yield", f_yield );
648 }
649 else {
650 assert( task==GET_FLOW_RULE );
651 f = f_flow;
652 if ( swit ) pri( "f_flow", f_flow );
653 }
654 }
655 }
656 if ( get_group_data( GROUP_MATERI_PLASTI_HLC, gr, element, new_unknowns,
657 plasti_data, ldum, GET_IF_EXISTS ) ) {
658 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
659 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_HLC;
660 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_HLC;
661 if ( test1 || test2 || test3 ) {
662 if ( swit ) pri( "plasti_hlc" );
663 void_fraction = new_unknowns[void_indx];
664 sig_yield = plasti_data[0];
665 m21 = plasti_data[1];
666 m22 = plasti_data[2];
667 m23 = plasti_data[3];
668 m1 = plasti_data[4];
669 m20 = (m21 - m22 * void_fraction) * exp(m23*sigm/sig_yield);
670 f_yield = (3.*sig_eq*sig_eq)/(sig_yield*sig_yield) +
671 (1+1/m20)*void_fraction*m1*exp(3.*sigm/(2.*sig_yield)) - 1.;
672 f_flow = f_yield;
673 if ( task==GET_YIELD_RULE ) {
674 if ( f_yield>f ) {
675 f = f_yield;
676 tmp_plasti_type = GROUP_MATERI_PLASTI_HLC;
677 }
678 if ( swit ) pri( "f_yield", f_yield );
679 }
680 else {
681 assert( task==GET_FLOW_RULE );
682 f = f_flow;
683 if ( swit ) pri( "f_flow", f_flow );
684 }
685 }
686 }
687 if ( get_group_data( GROUP_MATERI_PLASTI_MATSUOKANAKAI, gr, element, new_unknowns,
688 plasti_data, ldum, GET_IF_EXISTS ) ) {
689 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
690 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI;
691 test3 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX;
692 test4 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI;
693 test5 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX;
694 if ( test1 || test2 || test3 || test4 || test5 ) {
695 if ( swit ) pri( "check plasti_matsuokanakai" );
696 phi = plasti_data[0];
697 c = plasti_data[1];
698 phi_flow = plasti_data[2];
699 if ( plasti_on_boundary ) {
700 phi *= plasti_on_boundary_factor;
701 phi_flow *= plasti_on_boundary_factor;
702 }
703 alpha_flow = 2.*sin(phi_flow)/((3.-sin(phi_flow))*sqrt(3.));
704 if ( phi>0. )
705 sig0 = c * cos(phi) / sin(phi);
706 else
707 sig0 = 0.;
708 array_move( sig, sig_tmp, MDIM*MDIM );
709 for ( idim=0; idim<MDIM; idim++ ) sig_tmp[idim*MDIM+idim] -= sig0;
710 matrix_eigenvalues( sig_tmp, sig_princ );
711 I1 = sig_princ[0] + sig_princ[1] + sig_princ[2];
712 I2 = -sig_princ[0]*sig_princ[1] - sig_princ[1]*sig_princ[2] -
713 sig_princ[2]*sig_princ[0];
714 I3 = sig_princ[0] * sig_princ[1] * sig_princ[2];
715 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE && phi>0.;
716 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX;
717 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX;
718 if ( test1 || test2 || test3 ) {
719 // separate yield and flow rule for apex
720 array_set( sig_princ_apex, 0., MDIM );
721 array_set( e_vec, 1./sqrt(3.), MDIM );
722 array_subtract( sig_princ, sig_princ_apex, f_vec, MDIM );
723 in_apex = array_inproduct(f_vec,e_vec,MDIM)>0. && tension_cutoff==-YES;
724 if ( in_apex || plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX ) {
725 f_yield = array_size( f_vec, MDIM );
726 f_flow = array_size( f_vec, MDIM );
727 }
728 else {
729 f_yield = NO_YIELD_F;
730 f_flow = NO_YIELD_F;
731 }
732 if ( task==GET_YIELD_RULE ) {
733 if ( f_yield>f ) {
734 f = f_yield;
735 tmp_plasti_type = GROUP_MATERI_PLASTI_MATSUOKANAKAI_APEX;
736 }
737 if ( swit ) pri( "f_yield", f_yield );
738 }
739 else {
740 assert( task==GET_FLOW_RULE );
741 f = f_flow;
742 if ( swit ) pri( "f_flow", f_flow );
743 }
744 }
745 if ( !in_apex ) {
746 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
747 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI;
748 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MATSUOKANAKAI;
749 if ( test1 || test2 || test3 ) {
750 f_yield = I3 + scalar_square(cos(phi))/(9.-scalar_square(sin(phi))) * I1*I2;
751 f_flow = 3.*alpha_flow*sigm + sig_eq;
752 if ( task==GET_YIELD_RULE ) {
753 if ( f_yield>f ) {
754 f = f_yield;
755 tmp_plasti_type = GROUP_MATERI_PLASTI_MATSUOKANAKAI;
756 }
757 if ( swit ) pri( "f_yield", f_yield );
758 }
759 else {
760 assert( task==GET_FLOW_RULE );
761 f = f_flow;
762 if ( swit ) pri( "f_flow", f_flow );
763 }
764 }
765 }
766 }
767 }
768 if ( get_group_data( GROUP_MATERI_PLASTI_MOHRCOUL, gr, element, new_unknowns,
769 plasti_data, ldum, GET_IF_EXISTS ) ||
770 get_group_data( GROUP_MATERI_PLASTI_MOHRCOUL_SOFTENING, gr, element, new_unknowns,
771 plasti_data, ldum, GET_IF_EXISTS ) ) {
772 db( GROUP_MATERI_PLASTI_MOHRCOUL_TENSIONCUTOFF, 0, &tension_cutoff, ddum, ldum,
773 VERSION_NORMAL, GET_IF_EXISTS );
774 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
775 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_01;
776 test3 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_12;
777 test4 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_20;
778 test5 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_APEX;
779 test6 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_01;
780 test7 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_12;
781 test8 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_20;
782 test9 = GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_APEX;
783 if ( test1 || test2 || test3 || test4 || test5 || test6 || test7 || test8 || test9 ) {
784 if ( swit ) pri( "check plasti_mohrcoul" );
785 if ( get_group_data( GROUP_MATERI_PLASTI_MOHRCOUL_SOFTENING, gr, element,
786 new_unknowns, plasti_data, ldum, GET_IF_EXISTS ) ) {
787 phi0 = plasti_data[0];
788 c0 = plasti_data[1];
789 phi0_flow = plasti_data[2];
790 phi1 = plasti_data[3];
791 c1 = plasti_data[4];
792 phi1_flow = plasti_data[5];
793 kappa_crit = plasti_data[6];
794 kappa = new_unknowns[kap_indx];
795 if ( kappa>=kappa_crit ) {
796 phi = phi1;
797 c = c1;
798 phi_flow = phi1_flow;
799 }
800 else {
801 tmp = kappa/kappa_crit;
802 phi = phi0 + tmp*(phi1-phi0);
803 c = c0 + tmp*(c1-c0);
804 phi_flow = phi0_flow + tmp*(phi1_flow-phi0_flow);
805 }
806 }
807 else {
808 phi = plasti_data[0];
809 c = plasti_data[1];
810 phi_flow = plasti_data[2];
811 }
812 if ( plasti_on_boundary ) {
813 phi *= plasti_on_boundary_factor;
814 phi_flow *= plasti_on_boundary_factor;
815 }
816 matrix_eigenvalues( sig, sig_princ );
817 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE && phi>0.;
818 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_APEX;
819 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_APEX;
820 if ( test1 || test2 || test3 ) {
821 // separate yield and flow rule for apex
822 array_set( sig_princ_apex, c*cos(phi)/sin(phi), MDIM );
823 array_set( e_vec, 1./sqrt(3.), MDIM );
824 array_subtract( sig_princ, sig_princ_apex, f_vec, MDIM );
825 in_apex = array_inproduct(f_vec,e_vec,MDIM)>0. && tension_cutoff==-YES;
826 if ( in_apex || plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_APEX ) {
827 f_yield = array_size( f_vec, MDIM );
828 f_flow = array_size( f_vec, MDIM );
829 }
830 else {
831 f_yield = NO_YIELD_F;
832 f_flow = NO_YIELD_F;
833 }
834 if ( task==GET_YIELD_RULE ) {
835 if ( f_yield>f ) {
836 f = f_yield;
837 tmp_plasti_type = GROUP_MATERI_PLASTI_MOHRCOUL_APEX;
838 }
839 if ( swit ) pri( "f_yield", f_yield );
840 }
841 else {
842 assert( task==GET_FLOW_RULE );
843 f = f_flow;
844 if ( swit ) pri( "f_flow", f_flow );
845 }
846 }
847 if ( !in_apex ) {
848 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
849 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_01;
850 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_01;
851 if ( test1 || test2 || test3 ) {
852 f_yield = 0.5*scalar_dabs(sig_princ[0]-sig_princ[1]) +
853 0.5*(sig_princ[0]+sig_princ[1])*sin(phi) - c*cos(phi);
854 f_flow = 0.5*scalar_dabs(sig_princ[0]-sig_princ[1]) +
855 0.5*(sig_princ[0]+sig_princ[1])*sin(phi_flow) - c*cos(phi_flow);
856 if ( task==GET_YIELD_RULE ) {
857 if ( f_yield>f ) {
858 f = f_yield;
859 tmp_plasti_type = GROUP_MATERI_PLASTI_MOHRCOUL_01;
860 if ( swit ) pri( "f_yield", f_yield );
861 }
862 }
863 else {
864 assert( task==GET_FLOW_RULE );
865 f = f_flow;
866 if ( swit ) pri( "f_flow", f_flow );
867 }
868 }
869 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
870 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_12;
871 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_12;
872 if ( test1 || test2 || test3 ) {
873 f_yield = 0.5*scalar_dabs(sig_princ[1]-sig_princ[2]) +
874 0.5*(sig_princ[1]+sig_princ[2])*sin(phi) - c*cos(phi);
875 f_flow = 0.5*scalar_dabs(sig_princ[1]-sig_princ[2]) +
876 0.5*(sig_princ[1]+sig_princ[2])*sin(phi_flow) - c*cos(phi_flow);
877 if ( task==GET_YIELD_RULE ) {
878 if ( f_yield>f ) {
879 f = f_yield;
880 tmp_plasti_type = GROUP_MATERI_PLASTI_MOHRCOUL_12;
881 if ( swit ) pri( "f_yield", f_yield );
882 }
883 }
884 else {
885 assert( task==GET_FLOW_RULE );
886 f = f_flow;
887 if ( swit ) pri( "f_flow", f_flow );
888 }
889 }
890 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
891 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_20;
892 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_MOHRCOUL_20;
893 if ( test1 || test2 || test3 ) {
894 f_yield = 0.5*scalar_dabs(sig_princ[2]-sig_princ[0]) +
895 0.5*(sig_princ[2]+sig_princ[0])*sin(phi) - c*cos(phi);
896 f_flow = 0.5*scalar_dabs(sig_princ[2]-sig_princ[0]) +
897 0.5*(sig_princ[2]+sig_princ[0])*sin(phi_flow) - c*cos(phi_flow);
898 if ( task==GET_YIELD_RULE ) {
899 if ( f_yield>f ) {
900 f = f_yield;
901 tmp_plasti_type = GROUP_MATERI_PLASTI_MOHRCOUL_20;
902 if ( swit ) pri( "f_yield", f_yield );
903 }
904 }
905 else {
906 assert( task==GET_FLOW_RULE );
907 f = f_flow;
908 if ( swit ) pri( "f_flow", f_flow );
909 }
910 }
911 }
912 }
913 }
914 if ( get_group_data( GROUP_MATERI_PLASTI_STRESS, gr, element, new_unknowns,
915 plasti_data, ldum, GET_IF_EXISTS ) ) {
916 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
917 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_STRESS;
918 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_STRESS;
919 if ( test1 || test2 || test3 ) {
920 if ( swit ) pri( "check plasti_stress" );
921 sig_yield = plasti_data[0];
922 if ( swit ) pri( "sig_yield", sig_yield );
923 matrix_eigenvalues( sig, sig_princ );
924 if ( swit ) pri( "sig_princ", sig_princ, MDIM );
925 tmp_sig_eq = 0.;
926 for ( idim=0; idim<MDIM; idim++ ) {
927 tmp_sig_eq += sig_princ[idim] * sig_princ[idim];
928 }
929 tmp_sig_eq = sqrt( scalar_dabs(tmp_sig_eq) );
930 f_yield = tmp_sig_eq - sig_yield;
931 f_flow = tmp_sig_eq - sig_yield;
932 if ( task==GET_YIELD_RULE ) {
933 if ( f_yield>f ) {
934 f = f_yield;
935 tmp_plasti_type = GROUP_MATERI_PLASTI_STRESS;
936 }
937 if ( swit ) pri( "f_yield", f_yield );
938 }
939 else {
940 assert( task==GET_FLOW_RULE );
941 f = f_flow;
942 if ( swit ) pri( "f_flow", f_flow );
943 }
944 }
945 }
946 if ( get_group_data( GROUP_MATERI_PLASTI_TENSION, gr, element, new_unknowns,
947 plasti_data, ldum, GET_IF_EXISTS ) ) {
948 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
949 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_TENSION;
950 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_TENSION;
951 if ( test1 || test2 || test3 ) {
952 if ( swit ) pri( "check plasti_tension" );
953 sig_yield = plasti_data[0];
954 if ( swit ) pri( "sig_yield", sig_yield );
955 matrix_eigenvalues( sig, sig_princ );
956 if ( swit ) pri( "sig_princ", sig_princ, MDIM );
957 tmp_sig_eq = 0.;
958 for ( idim=0; idim<MDIM; idim++ ) {
959 if ( sig_princ[idim]>0. ) {
960 tmp_sig_eq += sig_princ[idim] * sig_princ[idim];
961 }
962 }
963 tmp_sig_eq = sqrt( scalar_dabs(tmp_sig_eq) );
964 f_yield = tmp_sig_eq - sig_yield;
965 f_flow = tmp_sig_eq - sig_yield;
966 if ( task==GET_YIELD_RULE ) {
967 if ( f_yield>f ) {
968 f = f_yield;
969 tmp_plasti_type = GROUP_MATERI_PLASTI_TENSION;
970 }
971 if ( swit ) pri( "f_yield", f_yield );
972 }
973 else {
974 assert( task==GET_FLOW_RULE );
975 f = f_flow;
976 if ( swit ) pri( "f_flow", f_flow );
977 }
978 }
979 }
980 if ( get_group_data( GROUP_MATERI_PLASTI_VONMISES, gr, element, new_unknowns,
981 plasti_data, ldum, GET_IF_EXISTS ) ) {
982 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
983 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_VONMISES;
984 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_VONMISES;
985 if ( test1 || test2 || test3 ) {
986 if ( swit ) pri( "plasti_vonmises" );
987 if ( get_group_data( GROUP_MATERI_PLASTI_VONMISES_NADAI, gr, element,
988 new_unknowns, work, ldum, GET_IF_EXISTS ) ) {
989 C = work[0];
990 kappa_0 = work[1];
991 n = work[2];
992 sig_yield = plasti_data[0] + C * scalar_power(kappa_0+new_unknowns[kap_indx],n);
993 }
994 else
995 sig_yield = plasti_data[0];
996 f_yield = sig_eq*sqrt(3.) - sig_yield;
997 f_flow = f_yield;
998 if ( task==GET_YIELD_RULE ) {
999 if ( f_yield>f ) {
1000 f = f_yield;
1001 tmp_plasti_type = GROUP_MATERI_PLASTI_VONMISES;
1002 }
1003 if ( swit ) pri( "f_yield", f_yield );
1004 }
1005 else {
1006 assert( task==GET_FLOW_RULE );
1007 f = f_flow;
1008 if ( swit ) pri( "f_flow", f_flow );
1009 }
1010 }
1011 }
1012 test1 = task==GET_YIELD_RULE&&plasti_type==-NONE;
1013 test2 = task==GET_YIELD_RULE&&plasti_type==GROUP_MATERI_PLASTI_USER;
1014 test3 = task==GET_FLOW_RULE&&plasti_type==GROUP_MATERI_PLASTI_USER;
1015 if ( test1 || test2 || test3 ) {
1016 if ( swit ) pri( "check plasti_user" );
1017 if ( task==GET_YIELD_RULE ) {
1018 f_yield = NO_YIELD_F;
1019 if ( group_materi_plasti_user==-YES )
1020 user_plasti( 1, user_data, new_unknowns, old_hisv, new_hisv, sig, f_yield );
1021 if ( f_yield>f ) {
1022 f = f_yield;
1023 tmp_plasti_type = GROUP_MATERI_PLASTI_USER;
1024 }
1025 if ( swit ) pri( "f_yield", f_yield );
1026 }
1027 else {
1028 assert( task==GET_FLOW_RULE );
1029 if ( group_materi_plasti_user==-YES )
1030 user_plasti( 2, user_data, new_unknowns, old_hisv, new_hisv, sig, f_flow );
1031 f = f_flow;
1032 if ( swit ) pri( "f_flow", f_flow );
1033 }
1034 }
1035 new_f = f;
1036 if ( materi_plasti_f && derivatives ) {
1037 for ( idim=0; idim<ndim; idim++ ) {
1038 tmp = new_grad_new_unknowns[idim*nuknwn+f_indx+idim*nder];
1039 f += options_nonlocal * tmp;
1040 }
1041 }
1042 if ( task==GET_YIELD_RULE && plasti_type==-NONE ) plasti_type = tmp_plasti_type;
1043 if ( task==GET_YIELD_RULE || task==GET_FLOW_RULE ) {
1044 if ( swit ) {
1045 pri( "new_f", new_f );
1046 pri( "f", f );
1047 pri( "plasti_type", plasti_type );
1048 }
1049 }
1050
1051 // use central differences to determine flow direction
1052 signorm = array_size( sig, MDIM*MDIM );
1053 if ( task==GET_FLOW_RULE_GRAD ) {
1054 array_move( sig, sig_tmp, MDIM*MDIM );
1055 for ( i=0; i<MDIM*MDIM; i++ ) {
1056 sig_tmp[i] += EPS_SIG * signorm;
1057 plasti_rule( element, gr, plasti_on_boundary,
1058 user_data, new_unknowns, new_grad_new_unknowns,
1059 old_hisv, new_hisv, old_epp, inc_epp, inc_ept,
1060 GET_FLOW_RULE, plasti_type, sig_tmp, f_flow_right, rdum, ddum );
1061 sig_tmp[i] -= 2. * EPS_SIG * signorm;
1062 plasti_rule( element, gr, plasti_on_boundary,
1063 user_data, new_unknowns, new_grad_new_unknowns,
1064 old_hisv, new_hisv, old_epp, inc_epp, inc_ept,
1065 GET_FLOW_RULE, plasti_type, sig_tmp, f_flow_left, rdum, ddum );
1066 dir[i] = (f_flow_right-f_flow_left)/(2.*EPS_SIG*signorm);
1067 sig_tmp[i] += EPS_SIG * signorm;
1068 }
1069 }
1070
1071 if ( swit ) pri( "Out PLASTI_RULE" );
1072 }
1073