00001 #include "elemmat.h"
00002 #include <math.h>
00003 #include <string.h>
00004 #include <stdlib.h>
00005
00006 elemmat::elemmat (void)
00007 {
00008
00009 n=0;
00010
00011 ne=0;
00012
00013 decompid=0;
00014
00015 mem=0;
00016
00017 maxarr=0;
00018
00019
00020 a=NULL;
00021
00022 acn=NULL;
00023
00024 neme = NULL;
00025
00026 andofe = NULL;
00027 }
00028
00029 elemmat::~elemmat (void)
00030 {
00031 long i;
00032
00033 for (i=0;i<ne;i++){
00034 delete [] a[i];
00035 delete [] acn[i];
00036 }
00037 delete [] a;
00038 delete [] acn;
00039
00040 delete [] neme;
00041 delete [] andofe;
00042 }
00043
00044 void elemmat::alloc (gtopology *top)
00045 {
00046 long i,ndofe;
00047
00048
00049 ne=top->ne;
00050
00051 neme = new long [ne];
00052
00053 andofe = new long [ne];
00054
00055 mem=0; maxarr=0;
00056 for (i=0;i<ne;i++){
00057 ndofe = top->give_ndofe (i);
00058 andofe[i] = ndofe;
00059 neme[i] = ndofe*ndofe;
00060 mem += neme[i];
00061 if (ndofe>maxarr) maxarr=ndofe;
00062 }
00063
00064
00065 a = new double* [ne];
00066
00067 acn = new long* [ne];
00068 for (i=0;i<ne;i++){
00069 a[i] = new double [neme[i]];
00070 acn[i] = new long [andofe[i]];
00071 }
00072 }
00073
00074 double** elemmat::status ()
00075 {
00076 return a;
00077 }
00078
00079 long elemmat::decomp ()
00080 {
00081 return decompid;
00082 }
00083
00084 void elemmat::changedecomp ()
00085 {
00086 if (decompid==0) decompid=1;
00087 else decompid=0;
00088 }
00089
00090 void elemmat::nullmat ()
00091 {
00092 long i,j;
00093 for (i=0;i<ne;i++){
00094 for (j=0;j<neme[i];j++){
00095 a[i][j]=0.0;
00096 }
00097 }
00098 }
00099
00100 void elemmat::localize (matrix &b,long *cn,long eid)
00101 {
00102 long i,j,k;
00103
00104 k=0;
00105 for (i=0;i<b.m;i++){
00106 acn[eid][i]=cn[i];
00107 for (j=0;j<b.m;j++){
00108 a[eid][k]=b[i][j];
00109 k++;
00110 }
00111 }
00112 }
00113
00114 void elemmat::localized (double *b,long *cn,long eid,long m)
00115 {
00116 long i,j,k;
00117
00118 k=0;
00119 for (i=0;i<m;i++){
00120 acn[eid][i]=cn[i];
00121 for (j=0;j<m;j++){
00122 a[eid][k]=b[i*m+j];
00123 k++;
00124 }
00125 }
00126 }
00127
00128 void elemmat::initiate (gtopology *top,long ndof,long mespr)
00129 {
00130 if (status ()==NULL){
00131 alloc (top);
00132 }
00133 else{
00134 nullmat ();
00135 }
00136
00137 n=ndof;
00138
00139
00140 }
00141
00142
00143 void elemmat::mxv_em (double *b,double *c)
00144 {
00145 long i,j,ii,uj;
00146 double *u,*v;
00147
00148 u = new double [maxarr];
00149 v = new double [maxarr];
00150
00151 nullv (c,n);
00152
00153 for (i=0;i<ne;i++){
00154 uj=andofe[i];
00155
00156 nullv (u,maxarr);
00157 for (j=0;j<uj;j++){
00158 ii=acn[i][j];
00159 if (ii>0)
00160 u[j]=b[ii-1];
00161 }
00162
00163 mxv (a[i],u,v,uj,uj);
00164
00165 for (j=0;j<uj;j++){
00166 ii=acn[i][j];
00167 if (ii>0)
00168 c[ii-1]+=v[j];
00169 }
00170 }
00171
00172 delete [] u;
00173 delete [] v;
00174 }
00175
00176 void elemmat::addmat_em (double c,elemmat &dm)
00177 {
00178 long i,j;
00179 for (i=0;i<ne;i++){
00180 for (j=0;j<neme[i];j++){
00181 a[i][j]+=c*dm.a[i][j];
00182 }
00183 }
00184 }
00185
00186 void elemmat::scalmat_em (double c)
00187 {
00188 long i,j;
00189 for (i=0;i<ne;i++){
00190 for (j=0;j<neme[i];j++){
00191 a[i][j]*=c;
00192 }
00193 }
00194 }
00195
00196 void elemmat::printmat (FILE *out)
00197 {
00198 long i,j;
00199
00200 fprintf (out,"\n\n");
00201 for (i=0;i<ne;i++){
00202 fprintf (out,"\n%ld",i);
00203 for (j=0;j<neme[i];j++){
00204 fprintf (out," %5.2f",a[i][j]);
00205 }
00206 }
00207 }
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227 void elemmat::cg (double *x,double *y,
00228 long ni,double err,long &ani,double &ares,double zero,long iv)
00229 {
00230 long i,j;
00231 double nom,denom,nory,alpha,beta;
00232 double *d,*r,*p;
00233
00234 d = new double [n];
00235 r = new double [n];
00236 p = new double [n];
00237
00238
00239 if (iv==0){
00240 for (i=0;i<n;i++)
00241 x[i]=0.0;
00242 }
00243 mxv_em (x,p);
00244
00245 nory=0.0; nom=0.0;
00246 for (i=0;i<n;i++){
00247 nory+=y[i]*y[i];
00248 r[i]=p[i]-y[i];
00249 nom+=r[i]*r[i];
00250 d[i]=-1.0*r[i];
00251 }
00252
00253 if (nory<zero){
00254 fprintf (stderr,"\n\n norm of right hand side in conjugate gradient method is smaller than %e",zero);
00255 fprintf (stderr,"\n see file %s, line %d.\n",__FILE__,__LINE__);
00256 ares=nory; ani=0;
00257 return;
00258 }
00259
00260
00261
00262 for (i=0;i<ni;i++){
00263
00264
00265
00266
00267 mxv_em (d,p);
00268
00269 denom = ss (d,p,n);
00270 if (fabs(denom)<zero){
00271 fprintf (stdout,"\n there is zero denominator in alpha computation in conjugate gradient method (cr.cpp)\n");
00272 break;
00273 }
00274
00275 alpha = nom/denom;
00276
00277
00278 for (j=0;j<n;j++){
00279 x[j]+=alpha*d[j];
00280 r[j]+=alpha*p[j];
00281 }
00282
00283 denom=nom;
00284
00285 nom = ss (r,r,n);
00286
00287 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/nory);
00288
00289 if (nom/nory<err) break;
00290
00291
00292
00293 beta = nom/denom;
00294
00295
00296 for (j=0;j<n;j++){
00297 d[j]=beta*d[j]-r[j];
00298 }
00299 }
00300
00301 ani=i; ares=nom;
00302
00303 delete [] p; delete [] r; delete [] d;
00304 }
00305