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