00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <math.h>
00004 #include "splas1d.h"
00005 #include "global.h"
00006 #include "vecttens.h"
00007 #include "intpoints.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "stochdriver.h"
00011
00012
00013 splas1d::splas1d (void)
00014 {
00015 fs=0.0; k=0.0;
00016 }
00017
00018 splas1d::~splas1d (void)
00019 {
00020
00021 }
00022
00023 void splas1d::read (FILE *in)
00024 {
00025 fscanf (in,"%lf %lf",&fs,&k);
00026 sra.read (in);
00027 }
00028
00029 double splas1d::yieldfunction (matrix &sig,vector &q)
00030
00031
00032 {
00033 double f;
00034 f = fabs(sig[0][0])-(fs+k*q[0]);
00035 return f;
00036 }
00037
00038
00039
00040
00041
00042
00043
00044 void splas1d::deryieldfsigma (matrix &sig,matrix &dfds)
00045 {
00046 if (sig[0][0]<0.0) dfds[0][0]=-1.0;
00047 else dfds[0][0]=1.0;
00048 }
00049
00050
00051
00052
00053
00054
00055
00056 void splas1d::deryieldfq (vector &dq)
00057 {
00058 dq[0]=k*(-1.0);
00059 }
00060
00061 void splas1d::plasmod (matrix &h)
00062
00063
00064 {
00065 h[0][0]=k;
00066 }
00067
00068 double splas1d::plasmodscalar(vector &qtr)
00069 {
00070 double ret;
00071 vector dfq(qtr.n), hp(qtr.n);
00072 matrix h(qtr.n, qtr.n);
00073
00074 deryieldfq(dfq);
00075 plasmod (h);
00076 mxv (h,dfq,hp);
00077 scprd (hp,dfq,ret);
00078 return ret;
00079 }
00080
00081 void splas1d::updateq(double dgamma, vector &q)
00082 {
00083 long j;
00084 vector dfq(q.n), hp(q.n);
00085 matrix h(q.n, q.n);
00086
00087 deryieldfq(dfq);
00088 plasmod (h);
00089 mxv (h,dfq,hp);
00090 for (j=0;j<q.n;j++)
00091 q[j]-=dgamma*hp[j];
00092 }
00093
00094 void splas1d::nlstresses (long ipp, long im, long ido)
00095 {
00096 long ni;
00097 double err;
00098 vector epsn(1),epsp(1),q(1);
00099 double gamma;
00100
00101
00102 epsn[0] = Mm->ip[ipp].strain[0];
00103 gamma = Mm->ip[ipp].eqother[ido+0];
00104 epsp[0] = Mm->ip[ipp].eqother[ido+1];
00105 q[0] = Mm->ip[ipp].eqother[ido+2];
00106
00107
00108 if (sra.give_tsra () == cp){
00109 ni=sra.give_ni ();
00110 err=sra.give_err ();
00111
00112 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00113 }
00114 else{
00115 fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
00116 abort ();
00117 }
00118
00119
00120 Mm->ip[ipp].other[ido+0]=gamma;
00121 Mm->ip[ipp].other[ido+1]=epsp[0];
00122 Mm->ip[ipp].other[ido+2]=q[0];
00123 }
00124
00125 void splas1d::nonloc_nlstresses (long ipp, long im, long ido)
00126 {
00127 long ni;
00128 double err;
00129 vector epsn(1),epsp(1),q(1);
00130 double gamma;
00131
00132
00133 epsn[0] = Mm->ip[ipp].strain[0];
00134 gamma = Mm->ip[ipp].eqother[ido+0];
00135 epsp[0] = Mm->ip[ipp].nonloc[ido+1];
00136 q[0] = Mm->ip[ipp].eqother[ido+2];
00137
00138
00139 if (sra.give_tsra () == cp){
00140 ni=sra.give_ni ();
00141 err=sra.give_err ();
00142
00143 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00144 }
00145 else{
00146 fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
00147 abort ();
00148 }
00149
00150
00151 Mm->ip[ipp].other[ido+0]=gamma;
00152 Mm->ip[ipp].other[ido+1]=epsp[0];
00153 Mm->ip[ipp].other[ido+2]=q[0];
00154 }
00155
00156 void splas1d::matstiff (matrix &d,long ipp,long ido)
00157 {
00158 double denom1,denom2;
00159 vector sig(1),dq(1),u(1);
00160 matrix de(1,1),h(1,1),g(1,1),sigt(3,3);
00161
00162 if (Mp->nlman->stmat==0){
00163 Mm->elmatstiff (d,ipp);
00164 }
00165 else{
00166
00167
00168 Mm->elmatstiff (de,ipp);
00169
00170
00171 sig[0]=Mm->ip[ipp].stress[0];
00172 if (Mm->ip[ipp].eqother[ido+1]<Mp->zero){
00173 Mm->elmatstiff (d,ipp);
00174 }
00175 else{
00176
00177
00178 deryieldfsigma (sigt,sigt);
00179 tensor_vector (sig,sigt,Mm->ip[ipp].ssst,stress);
00180
00181 vxv (sig,sig,h);
00182 mxm (de,h,g);
00183 mxm (g,de,h);
00184
00185 mxv (de,sig,u);
00186 scprd (u,sig,denom1);
00187
00188 deryieldfq (dq);
00189 scprd (dq,dq,denom2);
00190
00191 cmulm (1.0/(denom1+denom2),h);
00192
00193 subm (de,h,d);
00194 }
00195 }
00196 }
00197
00198 void splas1d::updateval (long ipp, long ido)
00199 {
00200 Mm->ip[ipp].eqother[ido+0]=Mm->ip[ipp].other[ido+0];
00201 Mm->ip[ipp].eqother[ido+1]=Mm->ip[ipp].other[ido+1];
00202 Mm->ip[ipp].eqother[ido+2]=Mm->ip[ipp].other[ido+2];
00203 }
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 void splas1d::giveirrstrains (long ipp, long ido, vector &epsp)
00217 {
00218 epsp[0] = Mm->ip[ipp].eqother[ido+1];
00219 }
00220
00221
00222
00223 void splas1d::changeparam (atsel &atm,vector &val)
00224 {
00225 long i;
00226
00227 for (i=0;i<atm.num;i++){
00228 switch (atm.atrib[i]){
00229 case 0:{
00230 fs=val[i];
00231 break;
00232 }
00233 case 1:{
00234 k=val[i];
00235 break;
00236 }
00237 default:{
00238 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (%s, line %d).\n",__FILE__,__LINE__);
00239 }
00240 }
00241 }
00242 }
00243
00244 double splas1d::give_consparam (long ipp, long ido)
00245 {
00246 double gamma;
00247 gamma = Mm->ip[ipp].other[ido+0];
00248 return gamma;
00249 }
00250
00251 long splas1d::give_num_interparam ()
00252 {
00253 return 1;
00254 }
00255
00256 void splas1d::give_interparam (long ipp,long ido,vector &q)
00257 {
00258 long ncompstr=Mm->ip[ipp].ncompstr;
00259
00260 q[0]=Mm->ip[ipp].eqother[ido+ncompstr+1];
00261 }