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_H 1.e-8
23
materi(long int element,long int gr,long int nnol,long int npoint,long int nodes[],long int plasti_on_boundary,double coord_ip[],double old_coord[],double h[],double new_d[],double new_b[],double volume,double old_unknowns[],double new_unknowns[],double old_grad_old_unknowns[],double old_grad_new_unknowns[],double new_grad_new_unknowns[],double element_lhside[],double element_matrix[],double element_rhside[],double element_residue[],double tendon_element_rhside[])24 void materi( long int element, long int gr, long int nnol,
25 long int npoint, long int nodes[], long int plasti_on_boundary,
26 double coord_ip[], double old_coord[],
27 double h[], double new_d[], double new_b[],
28 double volume, double old_unknowns[],
29 double new_unknowns[], double old_grad_old_unknowns[],
30 double old_grad_new_unknowns[], double new_grad_new_unknowns[],
31 double element_lhside[], double element_matrix[],
32 double element_rhside[], double element_residue[],
33 double tendon_element_rhside[])
34
35 {
36 long int i=0, j=0, idim=0, jdim=0, kdim=0,
37 inol=0, jnol=0, m=0, n=0, indx=0, ipuknwn=0, iuknwn=0, jpuknwn=0,
38 swit=0, indxi=0, indxj=0, indx1=0, indx2=0, memory=-UPDATED,
39 ind_ddsdde=0, ldum=0, idum[1];
40 double rdum=0., dens=0., dtime=0., materi_expansion_linear=0.,
41 materi_expansion_volume=0., temp=0., tmp=0., damping=0., fac=0,
42 plasti_heatgeneration=0., viscosity_heatgeneration=0.,
43 viscosity=0., old_damage=0., new_damage=0., old_kappa=0., new_kappa=0.,
44 old_f=0., new_f=0., void_fraction=0., new_pres=0., old_substeps=0., new_substeps=0.,
45 softvar_nonl=0, softvar_l=0,
46 static_pressure=0., total_pressure=0., location=0.,
47 J=0., ddum[1], *force_gravity=NULL,
48 *old_deften=NULL, *new_deften=NULL, *inv_deften=NULL,
49 *old_epe=NULL, *inc_epe=NULL, *new_epe=NULL,
50 *old_epp=NULL, *inc_epp=NULL,
51 *old_ept=NULL, *inc_ept=NULL, *new_ept=NULL,
52 *old_rot=NULL, *inc_rot=NULL, *new_rot=NULL,
53 *inv_rot=NULL, *old_sig=NULL, *rot=NULL,
54 *new_sig=NULL, *total_new_sig=NULL,
55 *rotated_old_sig=NULL, *rotated_new_sig=NULL,
56 *rotated_old_epi=NULL, *rotated_new_epi=NULL,
57 *rotated_old_rho=NULL, *rotated_new_rho=NULL,
58 *rotated_old_msig=NULL, *rotated_new_msig=NULL,
59 *old_msig=NULL, *new_msig=NULL,
60 *old_rho=NULL, *new_rho=NULL,
61 *old_epi=NULL, *new_epi=NULL,
62 *old_hisv=NULL, *new_hisv=NULL,
63 *ddsdde=NULL, *ddsdde_tendon=NULL,
64 *ddsdde_total=NULL, *sigvec=NULL,
65 *work=NULL, *stiffness=NULL, *force=NULL,
66 *dbl_array=NULL,
67 *new_sig_nonrot=NULL;
68
69 swit = set_swit(element,-1,"materi");
70 if ( swit ) pri( "In routine MATERI." );
71
72 n = nnol*MDIM*nnol*MDIM +
73 nnol*ndim*nnol*ndim + nnol*ndim +
74 MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
75 MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
76 MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
77 MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
78 MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
79 MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
80 MDIM*MDIM + MDIM*MDIM +
81 materi_maxwell_stress*MDIM*MDIM + materi_maxwell_stress*MDIM*MDIM +
82 materi_maxwell_stress*MDIM*MDIM + materi_maxwell_stress*MDIM*MDIM + MDIM*MDIM + MDIM*MDIM +
83 nuknwn + nuknwn + MSTRAIN*MSTRAIN + MSTRAIN*MSTRAIN +
84 MSTRAIN*MSTRAIN + MSTRAIN*MSTRAIN +
85 nnol*MDIM*nnol*MDIM + nnol*ndim*nnol*ndim + nnol*ndim + MDIM*MDIM;
86 dbl_array = get_new_dbl(n);
87 assert( indx<=n );
88
89 indx = 0;
90 work = &dbl_array[indx]; indx += nnol*MDIM*nnol*MDIM;
91 stiffness = &dbl_array[indx]; indx += nnol*ndim*nnol*ndim;
92 force = &dbl_array[indx]; indx += nnol*ndim;
93 force_gravity = &dbl_array[indx]; indx += MDIM;
94 old_deften = &dbl_array[indx]; indx += MDIM*MDIM;
95 new_deften = &dbl_array[indx]; indx += MDIM*MDIM;
96 inv_deften = &dbl_array[indx]; indx += MDIM*MDIM;
97 old_epe = &dbl_array[indx]; indx += MDIM*MDIM;
98 inc_epe = &dbl_array[indx]; indx += MDIM*MDIM;
99 new_epe = &dbl_array[indx]; indx += MDIM*MDIM;
100 old_epp = &dbl_array[indx]; indx += MDIM*MDIM;
101 inc_epp = &dbl_array[indx]; indx += MDIM*MDIM;
102 old_ept = &dbl_array[indx]; indx += MDIM*MDIM;
103 inc_ept = &dbl_array[indx]; indx += MDIM*MDIM;
104 new_ept = &dbl_array[indx]; indx += MDIM*MDIM;
105 old_rot = &dbl_array[indx]; indx += MDIM*MDIM;
106 inc_rot = &dbl_array[indx]; indx += MDIM*MDIM;
107 new_rot = &dbl_array[indx]; indx += MDIM*MDIM;
108 inv_rot = &dbl_array[indx]; indx += MDIM*MDIM;
109 old_sig = &dbl_array[indx]; indx += MDIM*MDIM;
110 rot = &dbl_array[indx]; indx += MDIM*MDIM;
111 new_sig = &dbl_array[indx]; indx += MDIM*MDIM;
112 total_new_sig = &dbl_array[indx]; indx += MDIM*MDIM;
113 rotated_old_sig = &dbl_array[indx]; indx += MDIM*MDIM;
114 rotated_new_sig = &dbl_array[indx]; indx += MDIM*MDIM;
115 rotated_old_epi = &dbl_array[indx]; indx += MDIM*MDIM;
116 rotated_new_epi = &dbl_array[indx]; indx += MDIM*MDIM;
117 rotated_old_rho = &dbl_array[indx]; indx += MDIM*MDIM;
118 rotated_new_rho = &dbl_array[indx]; indx += MDIM*MDIM;
119 rotated_old_msig = &dbl_array[indx]; indx += materi_maxwell_stress*MDIM*MDIM;
120 rotated_new_msig = &dbl_array[indx]; indx += materi_maxwell_stress*MDIM*MDIM;
121 old_msig = &dbl_array[indx]; indx += materi_maxwell_stress*MDIM*MDIM;
122 new_msig = &dbl_array[indx]; indx += materi_maxwell_stress*MDIM*MDIM;
123 old_rho = &dbl_array[indx]; indx += MDIM*MDIM;
124 new_rho = &dbl_array[indx]; indx += MDIM*MDIM;
125 old_epi = &dbl_array[indx]; indx += MDIM*MDIM;
126 new_epi = &dbl_array[indx]; indx += MDIM*MDIM;
127 old_hisv = &dbl_array[indx]; indx += nuknwn;
128 new_hisv = &dbl_array[indx]; indx += nuknwn;
129 ddsdde = &dbl_array[indx]; indx += MSTRAIN*MSTRAIN;
130 ddsdde_tendon = &dbl_array[indx]; indx += MSTRAIN*MSTRAIN;
131 ddsdde_total = &dbl_array[indx]; indx += MSTRAIN*MSTRAIN;
132 sigvec = &dbl_array[indx]; indx += MSTRAIN*MSTRAIN;
133 work = &dbl_array[indx]; indx += nnol*MDIM*nnol*MDIM;
134 stiffness = &dbl_array[indx]; indx += nnol*ndim*nnol*ndim;
135 force = &dbl_array[indx]; indx += nnol*ndim;
136 new_sig_nonrot = &dbl_array[indx]; indx += MDIM*MDIM;
137 assert( indx<=n );
138
139 array_set( inc_epe, 0., MDIM*MDIM );
140 array_set( old_ept, 0., MDIM*MDIM );
141 array_set( old_epp, 0., MDIM*MDIM );
142 array_set( old_epe, 0., MDIM*MDIM );
143 array_set( old_epi, 0., MDIM*MDIM );
144 array_set( old_rot, 0., MDIM*MDIM );
145 array_set( old_sig, 0., MDIM*MDIM );
146 array_set( new_sig, 0., MDIM*MDIM );
147 array_set( new_sig_nonrot, 0., MDIM*MDIM );
148 array_set( old_rho, 0., MDIM*MDIM );
149 array_set( old_deften, 0., MDIM*MDIM );
150 array_set( old_hisv, 0., nuknwn );
151 array_set( new_hisv, 0., nuknwn );
152 array_set( old_msig, 0., materi_maxwell_stress*MDIM*MDIM );
153
154 db( DTIME, 0, idum, &dtime, ldum, VERSION_NEW, GET );
155 get_group_data( GROUP_MATERI_DAMPING, gr, element, new_unknowns,
156 &damping, ldum, GET_IF_EXISTS );
157 get_group_data( GROUP_MATERI_PLASTI_HEATGENERATION, gr, element,
158 new_unknowns, &plasti_heatgeneration, ldum, GET_IF_EXISTS );
159 db( GROUP_MATERI_MEMORY, gr, &memory, ddum, ldum,
160 VERSION_NORMAL, GET_IF_EXISTS );
161 dens = get_materi_density( element, gr, nnol, nodes, new_unknowns );
162 force_gravity_calculate( force_gravity );
163
164 if ( condif_temperature ) {
165 get_group_data( GROUP_MATERI_EXPANSION_VOLUME, gr, element, new_unknowns,
166 &materi_expansion_volume, ldum, GET_IF_EXISTS );
167 get_group_data( GROUP_MATERI_EXPANSION_LINEAR, gr, element, new_unknowns,
168 &materi_expansion_linear, ldum, GET_IF_EXISTS );
169 temp = new_unknowns[temp_indx];
170 dens = (1.-materi_expansion_volume*temp) * dens;
171 }
172
173 // get old stresses, etc.
174 for ( idim=0; idim<MDIM; idim++ ) {
175 for ( jdim=0; jdim<MDIM; jdim++ ) {
176 indx = idim*MDIM + jdim;
177 if ( materi_stress ) {
178 tmp = old_unknowns[stres_indx+stress_indx(idim,jdim)*nder];
179 old_sig[indx] = tmp;
180 new_sig[indx] = tmp;
181 }
182 if ( materi_strain_elasti ) old_epe[indx] =
183 old_unknowns[epe_indx+stress_indx(idim,jdim)*nder];
184 if ( materi_strain_total ) old_ept[indx] =
185 old_unknowns[ept_indx+stress_indx(idim,jdim)*nder];
186 for ( m=0; m<materi_maxwell_stress; m++ ) {
187 indx = m*MDIM*MDIM+idim*MDIM+jdim;
188 old_msig[indx] = old_unknowns[mstres_indx+
189 (m*6+stress_indx(idim,jdim))*nder];
190 }
191 if ( materi_plasti_rho ) old_rho[indx] =
192 old_unknowns[rho_indx+stress_indx(idim,jdim)*nder];
193 if ( materi_strain_intergranular ) old_epi[indx] =
194 old_unknowns[epi_indx+stress_indx(idim,jdim)*nder];
195 if ( materi_strain_plasti ) old_epp[indx] =
196 old_unknowns[epp_indx+stress_indx(idim,jdim)*nder];
197 }
198 }
199
200 if ( materi_history_variables ) {
201 for ( i=0; i<materi_history_variables; i++ ) {
202 iuknwn = hisv_indx + i;
203 old_hisv[i] = old_unknowns[iuknwn];
204 new_hisv[i] = new_unknowns[iuknwn];
205 }
206 }
207 if ( materi_plasti_kappa ) {
208 iuknwn = kap_indx;
209 old_kappa = old_unknowns[iuknwn];
210 new_kappa = new_unknowns[iuknwn];
211 }
212 if ( materi_damage ) {
213 iuknwn = dam_indx;
214 old_damage = old_unknowns[iuknwn];
215 new_damage = new_unknowns[iuknwn];
216 }
217 if ( materi_plasti_f ) {
218 iuknwn = f_indx;
219 old_f = old_unknowns[iuknwn];
220 }
221 if ( materi_plasti_incremental_substeps ) {
222 iuknwn = substeps_indx;
223 old_substeps = old_unknowns[iuknwn];
224 }
225 if ( materi_plasti_softvar_nonlocal ) {
226 iuknwn = svnonloc_indx;
227 softvar_nonl = new_unknowns[iuknwn];
228 }
229 if ( materi_plasti_softvar_local ) {
230 iuknwn = svloc_indx;
231 softvar_l = new_unknowns[iuknwn];
232 }
233
234 set_deften_etc( element, gr, nnol, h, old_coord, old_unknowns,
235 new_unknowns, old_grad_old_unknowns, old_grad_new_unknowns,
236 old_deften, new_deften, old_ept, inc_ept, new_ept,
237 old_rot, inc_rot, new_rot );
238
239 // back rotate to old configuration
240 if ( memory==-TOTAL || memory==-TOTAL_PIOLA ) {
241 if ( !matrix_inverse( old_rot, inv_rot, rdum, MDIM ) ) {
242 pri ("Error detected for element ", element );
243 pri ("Probably too large distortions." );
244 exit_tn_on_error();
245 }
246 if ( materi_stress ) {
247 if ( memory==-TOTAL_PIOLA ) {
248 if ( !matrix_inverse( old_deften, inv_deften, rdum, MDIM ) ) {
249 pri ("Error detected for element ", element );
250 pri ("Probably too large distortions." );
251 exit_tn_on_error();
252 }
253 J = matrix_determinant( old_deften, MDIM );
254 matrix_abat( inv_deften, old_sig, rotated_old_sig, work, MDIM);
255 array_multiply( rotated_old_sig, rotated_old_sig, J, MDIM*MDIM );
256 }
257 else
258 matrix_abat( inv_rot, old_sig, rotated_old_sig, work, MDIM );
259 }
260 if ( materi_plasti_rho ) {
261 matrix_abat( inv_rot, old_rho, rotated_old_rho, work, MDIM );
262 if ( swit ) pri( "rotated_old_rho", rotated_old_rho, MDIM, MDIM );
263 }
264 if ( materi_strain_intergranular ) {
265 matrix_abat( inv_rot, old_epi, rotated_old_epi, work, MDIM );
266 if ( swit ) pri( "rotated_old_epi", rotated_old_epi, MDIM, MDIM );
267 }
268 if ( materi_maxwell_stress ) {
269 for ( m=0; m<materi_maxwell_stress; m++ ) {
270 indx = m*MDIM*MDIM;
271 matrix_abat( inv_rot, &old_msig[indx], &rotated_old_msig[indx], work, MDIM );
272 }
273 }
274 }
275 else if ( memory==-UPDATED || memory==-TOTAL_LINEAR ||
276 memory==-UPDATED_WITHOUT_ROTATION ) {
277 if ( materi_stress )
278 array_move( old_sig, rotated_old_sig, MDIM*MDIM );
279 if ( materi_plasti_rho )
280 array_move( old_rho, rotated_old_rho, MDIM*MDIM );
281 if ( materi_strain_intergranular )
282 array_move( old_epi, rotated_old_epi, MDIM*MDIM );
283 if ( materi_maxwell_stress )
284 array_move( old_msig, rotated_old_msig, materi_maxwell_stress*MDIM*MDIM );
285 }
286 else {
287 array_set( rotated_old_sig, 0., MDIM*MDIM );
288 }
289 if ( swit ) {
290 if ( materi_stress ) pri( "rotated_old_sig", rotated_old_sig, MDIM, MDIM );
291 }
292
293 // set stress, strain
294 if ( materi_stress ) {
295 set_stress( element, gr, plasti_on_boundary, coord_ip,
296 old_unknowns, new_unknowns,
297 old_grad_old_unknowns, new_grad_new_unknowns,
298 rotated_old_sig, new_sig,
299 rotated_old_msig, new_msig, inc_ept, new_ept,
300 old_epe, inc_epe, old_epp, inc_epp, old_rho, new_rho,
301 old_epi, new_epi, old_hisv, new_hisv,
302 old_damage, new_damage, old_kappa, new_kappa, new_f, new_substeps,
303 old_deften, new_deften, inc_rot,
304 ddsdde, viscosity, viscosity_heatgeneration, softvar_nonl, softvar_l);
305 tendons( element, gr, nnol, npoint, volume, new_d, old_unknowns, new_unknowns,
306 new_rot, inc_ept, tendon_element_rhside, ddsdde_tendon );
307 array_add( ddsdde, ddsdde_tendon, ddsdde_total, MSTRAIN*MSTRAIN );
308 if ( swit ) pri( "ddsdde_total", ddsdde_total, MSTRAIN, MSTRAIN );
309 }
310
311 if(find_local_softvar) {
312 for ( inol=0; inol<nnol; inol++ ) {
313 if ( scalar_dabs(h[inol])>EPS_H ) {
314 if ( materi_plasti_softvar_nonlocal && materi_plasti_softvar_local ) {
315 //added for options_element_dof
316 if(options_element_dof==-YES) new_unknowns[svloc_indx/nder] = softvar_l;
317 }
318 }
319 }
320 }
321
322 //Not used when searching for local values of softening variable
323 if(!find_local_softvar) {
324
325 array_move( new_sig, new_sig_nonrot, MDIM*MDIM );
326 // rotate to new configuration
327 if ( memory==-UPDATED || memory==-TOTAL || memory==-TOTAL_PIOLA ) {
328 if ( memory==-UPDATED )
329 array_move( inc_rot, rot, MDIM*MDIM );
330 else {
331 assert( memory==-TOTAL || memory==-TOTAL_PIOLA );
332 array_move( new_rot, rot, MDIM*MDIM );
333 }
334 if ( materi_stress ) {
335 if ( memory==-TOTAL_PIOLA ) {
336 matrix_abat( new_deften, new_sig, rotated_new_sig, work, MDIM );
337 J = matrix_determinant( new_deften, MDIM );
338 if ( J<=0. ) {
339 pri ("Error detected for element ", element );
340 pri ("Non-positive jacobian." );
341 pri ("Probably too large distortions." );
342 exit_tn_on_error();
343 }
344 array_multiply( rotated_new_sig, rotated_new_sig, 1./J, MDIM*MDIM );
345 }
346 else
347 matrix_abat( rot, new_sig, rotated_new_sig, work, MDIM );
348 array_move( rotated_new_sig, new_sig, MDIM*MDIM );
349 if ( swit ) pri( "new_sig", new_sig, MDIM, MDIM );
350 }
351 if ( materi_plasti_rho ) {
352 matrix_abat( rot, new_rho, rotated_new_rho, work, MDIM );
353 array_move( rotated_new_rho, new_rho, MDIM*MDIM );
354 }
355 if ( materi_strain_intergranular ) {
356 matrix_abat( rot, new_epi, rotated_new_epi, work, MDIM );
357 array_move( rotated_new_epi, new_epi, MDIM*MDIM );
358 }
359 if ( materi_maxwell_stress ) {
360 for ( m=0; m<materi_maxwell_stress; m++ ) {
361 indx = m*MDIM*MDIM;
362 matrix_abat( rot, &new_msig[indx], &rotated_new_msig[indx], work, MDIM );
363 }
364 array_move( rotated_new_msig, new_msig, materi_maxwell_stress*MDIM*MDIM );
365 }
366 }
367
368 // forces on nodes
369 if ( materi_stress ) {
370 array_move( new_sig, total_new_sig, MDIM*MDIM );
371 if ( groundflow_pressure ) {
372 new_pres = new_unknowns[pres_indx];
373 if ( groundflow_phreatic_coord( -1, coord_ip, new_unknowns,
374 total_pressure, static_pressure, location ) ) new_pres = total_pressure;
375 for ( idim=0; idim<MDIM; idim++ ) total_new_sig[idim*MDIM+idim] += new_pres;
376 }
377 for ( idim=0; idim<MDIM; idim++ ) {
378 for ( jdim=idim; jdim<MDIM; jdim++ ) {
379 indx = stress_indx(idim,jdim);
380 sigvec[indx] = total_new_sig[idim*MDIM+jdim];
381 }
382 }
383 matrix_atb( new_b, sigvec, force, MSTRAIN, nnol*ndim, 1 );
384 matrix_atba( new_b, ddsdde_total, stiffness, work, MSTRAIN, nnol*ndim );
385 if ( swit ) {
386 pri( "force", force, nnol*ndim );
387 pri( "stiffness", stiffness, nnol*ndim, nnol*ndim );
388 pri( "sigvec", sigvec, MSTRAIN );
389 }
390 }
391
392 // new elastic strains
393 array_add( old_epe, inc_epe, new_epe, MDIM*MDIM );
394
395 // add to right hand side and left hand side
396 fac = ((double)nnol)/2.;
397 for ( inol=0; inol<nnol; inol++ ) {
398
399 // velocity
400 for ( idim=0; idim<ndim; idim++ ) {
401 ipuknwn = vel_indx/nder+idim;
402 indx = inol*npuknwn + ipuknwn;
403 indxi = inol*npuknwn + ipuknwn;
404 iuknwn = vel_indx + idim*nder;
405 // damping
406 tmp = - h[inol] * damping * new_grad_new_unknowns[idim*nuknwn+iuknwn];
407 element_rhside[indx] += volume * tmp;
408 if ( residue ) element_residue[indx] -= tmp;
409 for ( jnol=0; jnol<nnol; jnol++ ) {
410 indxj = jnol*npuknwn + ipuknwn;
411 tmp = volume * h[inol] * damping * new_d[idim*nnol+jnol];
412 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
413 if ( jnol==inol ) element_lhside[indx] += tmp;
414 }
415 // force_gravity
416 tmp = h[inol] * dens * force_gravity[idim];
417 element_rhside[indx] += volume * tmp;
418 if ( residue ) element_residue[indx] -= tmp;
419 if ( materi_stress ) {
420 // stress gradient (rhside with green partial integration)
421 tmp = force[inol*ndim+idim];
422 element_rhside[indx] -= volume * tmp;
423 for ( jdim=0; jdim<ndim; jdim++ ) {
424 iuknwn = stres_indx+stress_indx(idim,jdim)*nder;
425 if ( residue ) element_residue[indx] += h[inol] *
426 new_grad_new_unknowns[jdim*nuknwn+iuknwn];
427 }
428 for ( jnol=0; jnol<nnol; jnol++ ) {
429 // groundflow
430 if ( groundflow_pressure ) {
431 jpuknwn = pres_indx/nder;
432 indxj = jnol*npuknwn + jpuknwn;
433 tmp = volume * new_d[idim*nnol+inol] * h[jnol];
434 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
435 }
436 // temperature stiffness
437 if ( condif_temperature ) {
438 for ( kdim=0; kdim<ndim; kdim++ ) {
439 jpuknwn = temp_indx/nder;
440 indxj = jnol*npuknwn + jpuknwn;
441 i = stress_indx(idim,idim);
442 j = stress_indx(kdim,kdim);
443 ind_ddsdde = i*MSTRAIN+j;
444 tmp = volume * materi_expansion_linear * new_d[idim*nnol+inol] *
445 -h[jnol] * ddsdde[ind_ddsdde];
446 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
447 }
448 }
449 for ( jdim=0; jdim<ndim; jdim++ ) {
450 // stiffness
451 indx1 = inol*ndim + idim;
452 indx2 = jnol*ndim + jdim;
453 jpuknwn = vel_indx/nder + jdim;
454 indxj = jnol*npuknwn + jpuknwn;
455 tmp = volume * dtime * stiffness[indx1*nnol*ndim+indx2];
456 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
457 if ( indxi==indxj ) element_lhside[indx] += fac * tmp;
458 // viscosity
459 jpuknwn = vel_indx/nder + jdim;
460 indxj = jnol*npuknwn + jpuknwn;
461 tmp = volume * new_d[jdim*nnol+inol] * viscosity * new_d[idim*nnol+jnol];
462 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
463 if ( idim==jdim && inol==jnol ) element_lhside[indx] += fac * tmp;
464 jpuknwn = vel_indx/nder + idim;
465 indxj = jnol*npuknwn + jpuknwn;
466 tmp = volume * new_d[jdim*nnol+inol] * viscosity * new_d[jdim*nnol+jnol];
467 element_matrix[indxi*nnol*npuknwn+indxj] += tmp;
468 if ( idim==jdim && inol==jnol ) element_lhside[indx] += fac * tmp;
469 }
470 }
471 }
472 }
473 if ( scalar_dabs(h[inol])>EPS_H ) {
474
475 if ( materi_history_variables ) {
476 for ( i=0; i<materi_history_variables; i++ ) {
477 ipuknwn = hisv_indx/nder + i;
478 indx = inol*npuknwn + ipuknwn;
479 tmp = volume * h[inol] * ( new_hisv[i] - old_hisv[i] ) / dtime;
480 element_rhside[indx] += tmp;
481 //added for options_element_dof
482 if(options_element_dof==-YES) new_unknowns[ipuknwn] = new_hisv[i];
483 ipuknwn++;
484 }
485 }
486
487 if ( materi_damage ) {
488 ipuknwn = dam_indx/nder;
489 indx = inol*npuknwn + ipuknwn;
490 tmp = volume * h[inol] * ( new_damage - old_damage ) / dtime;
491 element_rhside[indx] += tmp;
492 }
493
494 if ( materi_stress ) {
495 ipuknwn = stres_indx/nder;
496 for ( idim=0; idim<MDIM; idim++ ) {
497 for ( jdim=idim; jdim<MDIM; jdim++ ) {
498 indx = inol*npuknwn + ipuknwn;
499 tmp = volume * h[inol] * ( new_sig[idim*MDIM+jdim] -
500 old_sig[idim*MDIM+jdim] ) / dtime;
501 element_rhside[indx] += tmp;
502 //added for options_element_dof
503 if(options_element_dof==-YES) new_unknowns[ipuknwn] = new_sig[idim*MDIM+jdim];// new_sig_nonrot[idim*MDIM+jdim];
504 ipuknwn++;
505 }
506 }
507 }
508 if ( materi_plasti_f ) {
509 ipuknwn = f_indx/nder;
510 indx = inol*npuknwn + ipuknwn;
511 tmp = volume * h[inol] * ( new_f - old_f ) / dtime;
512 element_rhside[indx] += tmp;
513 }
514
515 if ( materi_plasti_incremental_substeps ) {
516 ipuknwn = substeps_indx/nder;
517 indx = inol*npuknwn + ipuknwn;
518 tmp = volume * h[inol] * ( new_substeps - old_substeps ) / dtime;
519 element_rhside[indx] += tmp;
520 }
521
522 if ( materi_plasti_kappa ) {
523 ipuknwn = kap_indx/nder;
524 indx = inol*npuknwn + ipuknwn;
525 tmp = volume * h[inol] * ( new_kappa - old_kappa ) / dtime;
526 element_rhside[indx] += tmp;
527 }
528
529 if ( materi_plasti_rho ) {
530 ipuknwn = rho_indx/nder;
531 for ( idim=0; idim<MDIM; idim++ ) {
532 for ( jdim=idim; jdim<MDIM; jdim++ ) {
533 indx = inol*npuknwn + ipuknwn;
534 tmp = volume * h[inol] * ( new_rho[idim*MDIM+jdim] -
535 old_rho[idim*MDIM+jdim] ) / dtime;
536 element_rhside[indx] += tmp;
537 ipuknwn++;
538 }
539 }
540 }
541 if ( materi_plasti_softvar_local ) {
542 ipuknwn = svloc_indx/nder;
543 indx = inol*npuknwn + ipuknwn;
544 tmp = volume * h[inol] * ( new_unknowns[ipuknwn] - old_unknowns[ipuknwn] ) / dtime;
545 element_rhside[indx] += tmp;
546 }
547 if ( materi_plasti_softvar_nonlocal ) {
548 ipuknwn = svnonloc_indx/nder;
549 indx = inol*npuknwn + ipuknwn;
550 tmp = volume * h[inol] * ( new_unknowns[ipuknwn] - old_unknowns[ipuknwn] ) / dtime;
551 element_rhside[indx] += tmp;
552 }
553 if ( materi_strainenergy ) {
554 ipuknwn = ener_indx/nder;
555 indx = inol*npuknwn + ipuknwn;
556 tmp = 0.5 * array_inproduct( new_sig, new_epe, MDIM*MDIM );
557 tmp = volume * h[inol] * ( tmp - old_unknowns[ener_indx] ) / dtime;
558 element_rhside[indx] += tmp;
559 }
560
561 if ( materi_strain_intergranular ) {
562 ipuknwn = epi_indx/nder;
563 for ( idim=0; idim<MDIM; idim++ ) {
564 for ( jdim=idim; jdim<MDIM; jdim++ ) {
565 indx = inol*npuknwn + ipuknwn;
566 tmp = volume * h[inol] * ( new_epi[idim*MDIM+jdim] -
567 old_epi[idim*MDIM+jdim] ) / dtime;
568 element_rhside[indx] += tmp;
569 //added for options_element_dof
570 if(options_element_dof==-YES)
571 new_unknowns[ipuknwn] = new_epi[idim*MDIM+jdim];
572 ipuknwn++;
573 }
574 }
575 }
576
577 if ( materi_strain_elasti ) {
578 ipuknwn = epe_indx/nder;
579 for ( idim=0; idim<MDIM; idim++ ) {
580 for ( jdim=idim; jdim<MDIM; jdim++ ) {
581 indx = inol*npuknwn + ipuknwn;
582 tmp = volume * h[inol] * inc_epe[idim*MDIM+jdim] / dtime;
583 element_rhside[indx] += tmp;
584 //added for options_element_dof
585 if(options_element_dof==-YES)
586 new_unknowns[ipuknwn] = old_epe[idim*MDIM+jdim] + inc_epe[idim*MDIM+jdim];
587 ipuknwn++;
588 }
589 }
590 }
591
592 if ( materi_strain_plasti ) {
593 ipuknwn = epp_indx/nder;
594 for ( idim=0; idim<MDIM; idim++ ) {
595 for ( jdim=idim; jdim<MDIM; jdim++ ) {
596 indx = inol*npuknwn + ipuknwn;
597 tmp = volume * h[inol] * inc_epp[idim*MDIM+jdim] / dtime;
598 element_rhside[indx] += tmp;
599 // added for options_element_dof
600 if(options_element_dof==-YES)
601 new_unknowns[ipuknwn] = old_epp[idim*MDIM+jdim] + inc_epp[idim*MDIM+jdim];
602 ipuknwn++;
603 }
604 }
605 if ( condif_temperature ) {
606 ipuknwn = temp_indx/nder;
607 indx = inol*npuknwn + ipuknwn;
608 tmp = plasti_heatgeneration * h[inol] *
609 array_inproduct( new_sig, inc_epp, MDIM*MDIM ) / dtime;
610 element_rhside[indx] += tmp;
611 element_residue[indx] -= tmp;
612 }
613 }
614
615 if ( condif_temperature ) {
616 ipuknwn = temp_indx/nder;
617 indx = inol*npuknwn + ipuknwn;
618 tmp = viscosity_heatgeneration * h[inol];
619 element_rhside[indx] += tmp;
620 element_residue[indx] -= tmp;
621 }
622
623 if ( materi_strain_total ) {
624 ipuknwn = ept_indx/nder;
625 for ( idim=0; idim<MDIM; idim++ ) {
626 for ( jdim=idim; jdim<MDIM; jdim++ ) {
627 indx = inol*npuknwn + ipuknwn;
628 tmp = volume * h[inol] * inc_ept[idim*MDIM+jdim] / dtime;
629 element_rhside[indx] += tmp;
630 //added for options_element_dof
631 if(options_element_dof==-YES)
632 new_unknowns[ipuknwn] = old_ept[idim*MDIM+jdim] + inc_ept[idim*MDIM+jdim];
633 ipuknwn++;
634 }
635 }
636 }
637
638 if ( materi_maxwell_stress ) {
639 ipuknwn = mstres_indx/nder;
640 for ( m=0; m<materi_maxwell_stress; m++ ) {
641 for ( idim=0; idim<MDIM; idim++ ) {
642 for ( jdim=idim; jdim<MDIM; jdim++ ) {
643 indx = inol*npuknwn + ipuknwn;
644 tmp = volume * h[inol] * ( new_msig[m*MDIM*MDIM+idim*MDIM+jdim] -
645 old_msig[m*MDIM*MDIM+idim*MDIM+jdim] ) / dtime;
646 element_rhside[indx] += tmp;
647 ipuknwn++;
648 }
649 }
650 }
651 }
652
653 if ( materi_void_fraction && materi_strain_plasti ) {
654 ipuknwn = void_indx/nder;
655 indx = inol*npuknwn + ipuknwn;
656 void_fraction = new_unknowns[void_indx];
657 tmp = volume * h[inol] * ( 1 - void_fraction ) * void_fraction *
658 ( inc_epp[0*MDIM+0] + inc_epp[1*MDIM+1] + inc_epp[2*MDIM+2] );
659 element_rhside[indx] += tmp;
660 }
661
662 if ( materi_work ) {
663 ipuknwn = work_indx/nder;
664 indx = inol*npuknwn + ipuknwn;
665 array_subtract( new_sig, old_sig, work, MDIM*MDIM );
666 tmp = array_inproduct( work, inc_ept, MDIM*MDIM );
667 tmp = volume * h[inol] * ( tmp - old_unknowns[work_indx] ) / dtime;
668 element_rhside[indx] += tmp;
669 }
670
671 }
672 }//Not used when searching for local values of softening variable -- end
673 }
674
675 if ( swit ) {
676 pri( "tendon_element_rhside", tendon_element_rhside, nnol, npuknwn );
677 pri( "element_lhside", element_lhside, nnol, npuknwn );
678 pri( "element_rhside", element_rhside, nnol, npuknwn );
679 if ( residue ) pri( "element_residue", element_residue, nnol, npuknwn );
680 }
681
682 delete[] dbl_array;
683
684 if ( swit ) pri( "Out function MATERI" );
685
686 }
687
set_deften_etc(long int element,long int gr,long int nnol,double h[],double old_coord[],double old_unknowns[],double new_unknowns[],double old_grad_old_unknowns[],double old_grad_new_unknowns[],double old_deften[],double new_deften[],double old_ept[],double inc_ept[],double new_ept[],double old_rot[],double inc_rot[],double new_rot[])688 void set_deften_etc( long int element, long int gr, long int nnol, double h[],
689 double old_coord[], double old_unknowns[], double new_unknowns[],
690 double old_grad_old_unknowns[], double old_grad_new_unknowns[],
691 double old_deften[], double new_deften[],
692 double old_ept[], double inc_ept[], double new_ept[],
693 double old_rot[], double inc_rot[], double new_rot[] )
694
695 {
696 long int idim=0, jdim=0, indx=0, swit=0, memory=-UPDATED,
697 axisymmetric=-NO, ind=0, ldum=0, idum[1];
698 double dtime=0., radius=0., ddum[1],
699 uTu[MDIM*MDIM], old_u[MDIM*MDIM], new_u[MDIM*MDIM],
700 inc_u[MDIM*MDIM], inc_deften[MDIM*MDIM], coord_ip[MDIM];
701
702 swit = set_swit(element,-1,"set_deften_etc");
703 if ( swit ) pri( "In routine SET_DEFTEN_ETC." );
704
705 // initialize
706 db( DTIME, 0, idum, &dtime, ldum, VERSION_NEW, GET );
707 db( GROUP_MATERI_MEMORY, gr, &memory, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
708 db( GROUP_AXISYMMETRIC, gr, &axisymmetric, ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS );
709 array_set( old_deften, 0., MDIM*MDIM );
710 array_set( new_deften, 0., MDIM*MDIM );
711 array_set( inc_deften, 0., MDIM*MDIM );
712 array_set( old_rot, 0., MDIM*MDIM );
713 array_set( new_rot, 0., MDIM*MDIM );
714 array_set( inc_rot, 0., MDIM*MDIM );
715 array_set( inc_ept, 0., MDIM*MDIM );
716 array_set( new_ept, 0., MDIM*MDIM );
717
718 // deformation tensors
719 if ( materi_displacement || materi_velocity_integrated ) {
720 old_deften[0] = old_deften[4] = old_deften[8] = 1.;
721 new_deften[0] = new_deften[4] = new_deften[8] = 1.;
722 for ( idim=0; idim<ndim; idim++ ) {
723 for ( jdim=0; jdim<ndim; jdim++ ) {
724 indx = idim*MDIM+jdim;
725 old_deften[indx] +=
726 old_grad_old_unknowns[jdim*nuknwn+dis_indx+idim*nder];
727 new_deften[indx] +=
728 old_grad_new_unknowns[jdim*nuknwn+dis_indx+idim*nder];
729 }
730 }
731 }
732 inc_deften[0] = inc_deften[4] = inc_deften[8] = 1.;
733 for ( idim=0; idim<ndim; idim++ ) {
734 for ( jdim=0; jdim<ndim; jdim++ ) {
735 indx = idim*MDIM+jdim;
736 ind = jdim*nuknwn+vel_indx+idim*nder;
737 inc_deften[indx] += old_grad_new_unknowns[ind]*dtime;
738 }
739 }
740 if ( axisymmetric==-YES ) {
741 matrix_ab( h, old_coord, coord_ip, 1, nnol, ndim );
742 radius = scalar_dabs(coord_ip[0]);
743 if ( radius!=0. ) {
744 if ( materi_displacement ) {
745 ind = dis_indx;
746 old_deften[8] += old_unknowns[ind]/radius;
747 new_deften[8] += new_unknowns[ind]/radius;
748 }
749 else if ( materi_velocity_integrated ) {
750 ind = veli_indx;
751 old_deften[8] += old_unknowns[ind]/radius;
752 new_deften[8] += new_unknowns[ind]/radius;
753 }
754 ind = vel_indx;
755 inc_deften[8] += new_unknowns[ind]*dtime/radius;
756 }
757 }
758
759 // rotation matrices
760 if ( memory==-UPDATED_WITHOUT_ROTATION || memory==-TOTAL_LINEAR ) {
761 for ( idim=0; idim<MDIM; idim++ ) {
762 old_rot[idim*MDIM+idim] = 1.;
763 new_rot[idim*MDIM+idim] = 1.;
764 inc_rot[idim*MDIM+idim] = 1.;
765 }
766 }
767 else if ( materi_displacement || materi_velocity_integrated ) {
768 set_deften_u_rot( old_deften, old_u, old_rot );
769 set_deften_u_rot( new_deften, new_u, new_rot );
770 }
771 set_deften_u_rot( inc_deften, inc_u, inc_rot );
772
773 // strain matrices
774 if ( memory==-UPDATED_WITHOUT_ROTATION ) {
775 // linear engineering strains
776 for ( idim=0; idim<MDIM; idim++ ) {
777 for ( jdim=0; jdim<MDIM; jdim++ ) inc_ept[idim*MDIM+jdim] =
778 0.5*(inc_deften[idim*MDIM+jdim]+inc_deften[jdim*MDIM+idim]);
779 }
780 for ( idim=0; idim<MDIM; idim++ ) inc_ept[idim*MDIM+idim] -= 1.;
781 array_add( old_ept, inc_ept, new_ept, MDIM*MDIM );
782 }
783 else if ( memory==-UPDATED ) {
784 // U from incremental polar decomposition
785 for ( idim=0; idim<MDIM; idim++ ) {
786 for ( jdim=0; jdim<MDIM; jdim++ ) inc_ept[idim*MDIM+jdim] =
787 0.5*(inc_u[idim*MDIM+jdim]+inc_u[jdim*MDIM+idim]);
788 }
789 for ( idim=0; idim<MDIM; idim++ ) inc_ept[idim*MDIM+idim] -= 1.;
790 array_add( old_ept, inc_ept, new_ept, MDIM*MDIM );
791 }
792 else if ( memory==-TOTAL ) {
793 // U from total polar decomposition
794 for ( idim=0; idim<MDIM; idim++ ) {
795 for ( jdim=0; jdim<MDIM; jdim++ ) new_ept[idim*MDIM+jdim] =
796 0.5*(new_u[idim*MDIM+jdim]+new_u[jdim*MDIM+idim]);
797 }
798 for ( idim=0; idim<MDIM; idim++ ) new_ept[idim*MDIM+idim] -= 1.;
799 array_subtract( new_ept, old_ept, inc_ept, MDIM*MDIM );
800 }
801 else if ( memory==-TOTAL_LINEAR ) {
802 // linear engineering strains
803 for ( idim=0; idim<MDIM; idim++ ) {
804 for ( jdim=0; jdim<MDIM; jdim++ ) new_ept[idim*MDIM+jdim] =
805 0.5*(new_deften[idim*MDIM+jdim]+new_deften[jdim*MDIM+idim]);
806 }
807 for ( idim=0; idim<MDIM; idim++ ) new_ept[idim*MDIM+idim] -= 1.;
808 array_subtract( new_ept, old_ept, inc_ept, MDIM*MDIM );
809 }
810 else if ( memory==-TOTAL_PIOLA ) {
811 // Green-Lagrange strains
812 matrix_atb( new_deften, new_deften, uTu, MDIM, MDIM, MDIM );
813 for ( idim=0; idim<MDIM; idim++ ) {
814 for ( jdim=0; jdim<MDIM; jdim++ ) {
815 new_ept[idim*MDIM+jdim] += 0.5*uTu[idim*MDIM+jdim];
816 if ( idim==jdim ) new_ept[idim*MDIM+jdim] -= 0.5;
817 }
818 }
819 array_subtract( new_ept, old_ept, inc_ept, MDIM*MDIM );
820 }
821 else
822 db_error( GROUP_MATERI_MEMORY, gr );
823
824 if ( swit ) {
825 pri( "old_deften", old_deften, MDIM, MDIM );
826 pri( "new_deften", new_deften, MDIM, MDIM );
827 pri( "inc_deften", inc_deften, MDIM, MDIM );
828 pri( "old_rot", old_rot, MDIM, MDIM );
829 pri( "inc_rot", inc_rot, MDIM, MDIM );
830 pri( "new_rot", new_rot, MDIM, MDIM );
831 pri( "inc_ept", inc_ept, MDIM, MDIM );
832 pri( "new_ept", new_ept, MDIM, MDIM );
833 }
834
835 if ( swit ) pri( "Out function SET_DEFTEN_ETC" );
836 }
837
set_deften_u_rot(double deften[],double u[],double rot[])838 void set_deften_u_rot( double deften[], double u[], double rot[] )
839
840 {
841 long int idim=0, jdim=0, kdim=0, indx=0, idum[1];
842 double rdum=0., uTu[MDIM*MDIM], val[MDIM], dir[MDIM*MDIM], work[MDIM*MDIM];
843
844 matrix_atb( deften, deften, uTu, MDIM, MDIM, MDIM );
845 matrix_jacobi( uTu, MDIM, val, dir, idum );
846 for ( idim=0; idim<MDIM; idim++ ) {
847 for ( jdim=0; jdim<MDIM; jdim++ ) {
848 indx = idim*MDIM+jdim;
849 u[indx] = 0.;
850 for ( kdim=0; kdim<MDIM; kdim++ ) u[indx] +=
851 sqrt(scalar_dabs(val[kdim]))*dir[idim*MDIM+kdim]*dir[jdim*MDIM+kdim];
852 }
853 }
854 if ( matrix_inverse( u, work, rdum, MDIM ) ) {
855 matrix_ab( deften, work, rot, MDIM, MDIM, MDIM );
856 }
857 else {
858 for ( idim=0; idim<MDIM; idim++ ) {
859 for ( jdim=0; jdim<MDIM; jdim++ ) {
860 u[idim*MDIM+jdim] = 0.5*(deften[idim*MDIM+jdim]+deften[jdim*MDIM+idim]);
861 if ( jdim==idim )
862 rot[idim*MDIM+idim] = 1.;
863 else
864 rot[idim*MDIM+idim] = 0.;
865 }
866 }
867 }
868
869 }
870
get_materi_density(long int element,long int element_group,long int nnol,long int nodes[],double new_unknowns[])871 double get_materi_density( long int element, long int element_group, long int nnol, long int nodes[],
872 double new_unknowns[] )
873
874 {
875
876 long int ldum=0, all_below=0, node_phreaticlevel=0, inol=0, inod=0,
877 idum[1];
878 double materi_dens = 0., ddum[1], group_materi_density_groundflow[2];
879
880 if ( materi_density ) {
881 materi_dens = new_unknowns[dens_indx];
882 if ( materi_dens<0. ) materi_dens = 0.;
883 }
884 else if ( db_active_index( GROUP_MATERI_DENSITY_GROUNDFLOW, element_group, VERSION_NORMAL ) ) {
885 db( GROUP_MATERI_DENSITY_GROUNDFLOW, element_group, idum,
886 group_materi_density_groundflow, ldum, VERSION_NORMAL, GET_IF_EXISTS );
887 all_below = 1;
888 for ( inol=0; inol<nnol; inol++ ) {
889 inod = nodes[inol];
890 node_phreaticlevel = -BELOW;
891 db( NODE_PHREATICLEVEL, inod, &node_phreaticlevel, ddum,
892 ldum, VERSION_NORMAL, GET_IF_EXISTS );
893 if ( node_phreaticlevel==-ABOVE ) all_below = 0;
894 }
895 if ( all_below )
896 materi_dens = group_materi_density_groundflow[0]; // wet
897 else
898 materi_dens = group_materi_density_groundflow[1]; // dry
899 }
900 else {
901 get_group_data( GROUP_MATERI_DENSITY, element_group, element, new_unknowns,
902 &materi_dens, ldum, GET_IF_EXISTS );
903 }
904 return materi_dens;
905
906 }
907