00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "slipsurf.h"
00004 #include "global.h"
00005 #include "gtopology.h"
00006 #include "mechmat.h"
00007 #include "mechtop.h"
00008 #include "element.h"
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 long detect_plastic_zones(long *plast_ip, double gamma_min)
00028 {
00029 long i, j, k;
00030 long max_gamma_ip;
00031 long max_gamma_eid;
00032 double max_gamma = 0.0;
00033 long *adjelem;
00034 long nadjelem;
00035 long new_nadjelem;
00036 long *gpl_adjelem;
00037
00038
00039
00040
00041 long *plast_adjelem;
00042
00043 long nplast_ip;
00044 long zone_id = 0;
00045
00046
00047 for (i=0; i<Mm->tnip; i++)
00048 plast_ip[i] = -1;
00049
00050
00051
00052 Gtm->adjacelem(NULL);
00053 gpl_adjelem = new long[Mt->ne];
00054 do
00055 {
00056
00057 max_gamma = detect_max_gamma(plast_ip, max_gamma_ip, max_gamma_eid, gamma_min);
00058
00059 if (max_gamma_ip == -1)
00060 break;
00061
00062 zone_id++;
00063 nadjelem = Gtm->nadjelel[max_gamma_eid];
00064 adjelem = new long [nadjelem];
00065 memcpy(adjelem, Gtm->adjelel[max_gamma_eid], sizeof(*Gtm->adjelel[max_gamma_eid])*nadjelem);
00066 plast_adjelem = new long [nadjelem];
00067 memset(plast_adjelem, 0, sizeof(*plast_adjelem)*nadjelem);
00068 memset(gpl_adjelem, 0, sizeof(*gpl_adjelem)*Mt->ne);
00069 do
00070 {
00071 nplast_ip = detect_plast_ip(nadjelem, adjelem, plast_ip, zone_id, plast_adjelem);
00072 if (nplast_ip == 0)
00073 {
00074 nadjelem = 0;
00075 delete [] adjelem;
00076 delete [] plast_adjelem;
00077 }
00078 else
00079 {
00080 for(i=0; i<nadjelem; i++)
00081 gpl_adjelem[adjelem[i]]=-1;
00082
00083 new_nadjelem = 0;
00084 for(i=0; i<nadjelem; i++)
00085 {
00086 if (plast_adjelem[i])
00087 {
00088 for (j=0; j<Gtm->nadjelel[adjelem[i]]; j++)
00089 {
00090 k = Gtm->adjelel[adjelem[i]][j];
00091 if (gpl_adjelem[k] == 0)
00092 {
00093 new_nadjelem++;
00094 gpl_adjelem[k] = 1;
00095 }
00096 }
00097 }
00098 }
00099 if (new_nadjelem > nadjelem)
00100 {
00101 delete [] plast_adjelem;
00102 plast_adjelem = new long[new_nadjelem];
00103 delete [] adjelem;
00104 adjelem = new long[new_nadjelem];
00105 }
00106 j=0;
00107 for (i=0; i<Mt->ne; i++)
00108 {
00109 if (gpl_adjelem[i] == 1)
00110 {
00111 adjelem[j] = i;
00112 j++;
00113 }
00114 }
00115 nadjelem = new_nadjelem;
00116 memset(plast_adjelem, 0, sizeof(*plast_adjelem)*nadjelem);
00117 }
00118
00119 }while (nadjelem);
00120
00121 }while (1);
00122
00123 delete [] gpl_adjelem;
00124 return zone_id;
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 long detect_plast_ip(long nadjelem, long *adjelem, long *plast_ip, long zone_id, long *plast_adjelem)
00149 {
00150 long i, j, nplast_ip, tnip, ipp;
00151
00152 nplast_ip = 0;
00153 for(i=0; i<nadjelem; i++)
00154 {
00155 if (Gtm->leso[adjelem[i]]==0)
00156 continue;
00157
00158 tnip = Mt->give_tnip(adjelem[i]);
00159 ipp = Mt->elements[adjelem[i]].ipp[0][0];
00160 for(j=0; j<tnip; j++)
00161 {
00162 if (plast_ip[ipp] == 0)
00163 {
00164 plast_ip[ipp] = zone_id;
00165 plast_adjelem[i] = 1;
00166 nplast_ip++;
00167 }
00168 ipp++;
00169 }
00170 }
00171 return nplast_ip;
00172 }
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 double detect_max_gamma(long *plast_ip, long &max_gamma_ip, long &max_gamma_eid, double gamma_min)
00197 {
00198 long ipp = 0;
00199 long nplast_ip = 0;
00200 long i, j;
00201 double gamma;
00202 double max_gamma = gamma_min;
00203 max_gamma_ip = -1;
00204 max_gamma_eid = -1;
00205
00206 for(i=0; i<Mt->ne; i++)
00207 {
00208 if (Gtm->leso[i]==1)
00209 {
00210
00211 ipp = Mt->elements[i].ipp[0][0];
00212 for (j=0; j<Mt->give_tnip(i); j++)
00213 {
00214 if (plast_ip[ipp] < 1)
00215 {
00216 gamma = Mm->give_consparam(ipp);
00217 if (gamma > max_gamma)
00218 {
00219 max_gamma = gamma;
00220 max_gamma_ip = ipp;
00221 max_gamma_eid = i;
00222 }
00223 if (gamma > gamma_min )
00224 {
00225 plast_ip[ipp] = 0;
00226 nplast_ip++;
00227 }
00228 }
00229 ipp++;
00230 }
00231 }
00232 }
00233
00234 return max_gamma;
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 long approx_slip_surf(long id, long *plast_ip, double &a, double &b)
00254 {
00255 long i,j, ipp, npip, tnip;
00256 double sw2_x_lny = 0.0;
00257 double sw2_x = 0.0;
00258 double sw2_lny = 0.0;
00259 double sw2_x2 = 0.0;
00260 double sw2 = 0.0;
00261 double w;
00262 matrix coord;
00263
00264
00265
00266
00267
00268 npip=0;
00269 for(i=0; i<Mt->ne; i++)
00270 {
00271 if (Gtm->leso[i]==1)
00272 {
00273
00274 ipp = Mt->elements[i].ipp[0][0];
00275 tnip = Mt->give_tnip(i);
00276 reallocm(tnip, 3, coord);
00277 Mt->give_ipcoord_elem(i, coord);
00278 for (j=0; j<tnip; j++)
00279 {
00280 if (plast_ip[ipp] == id)
00281 {
00282 npip++;
00283 w = Mm->give_consparam(ipp);
00284 sw2_x_lny += w*w*coord[j][0]*log(coord[j][1]);
00285 sw2_x += w*w*coord[j][0];
00286 sw2_lny += w*w*log(coord[j][1]);
00287 sw2_x2 += w*w*coord[j][0]*coord[j][0];
00288 sw2 += w*w;
00289 }
00290 ipp++;
00291 }
00292 }
00293 }
00294 if (npip < 2)
00295 {
00296 print_err("Cannot determine the slip surface for the given plastic zone %ld\n."
00297 "At least two different points in the plastic state must be given.", __FILE__, __LINE__, __func__, id);
00298 return 1;
00299 }
00300
00301 a = (sw2_x_lny*sw2_x - sw2_lny*sw2_x2)/(sw2_x*sw2_x-sw2*sw2_x2);
00302 a = exp(a);
00303 b = (sw2_lny*sw2_x - sw2_x_lny*sw2)/(sw2_x*sw2_x-sw2*sw2_x2);
00304
00305
00306 return 0;
00307 }