00001 #include "saddpoint.h"
00002
00003 saddpoint::saddpoint ()
00004 {
00005
00006 nrdof=0;
00007
00008 nm=0;
00009
00010
00011 dm = new densemat ();
00012
00013 c=NULL;
00014 }
00015
00016 saddpoint::~saddpoint ()
00017 {
00018 delete [] c;
00019 delete dm;
00020 }
00021
00022
00023
00024
00025 void saddpoint::initiate (seqselnodes *selnodfeti,gtopology *top,FILE *out)
00026 {
00027 long i,j,k,l,ndofn,dof,nid;
00028 long *red;
00029 red = new long [ns+1];
00030
00031
00032
00033
00034 ncdofd = new long [ns];
00035 for (i=0;i<ns;i++){
00036
00037 }
00038
00039
00040
00041
00042 edofs = new long* [ns];
00043 for (i=0;i<ns;i++){
00044 edofs[i] = new long [ncdofd[i]];
00045 }
00046
00047 for (i=0;i<ns;i++){
00048 for (j=0;j<ncdofd[i];j++){
00049
00050 }
00051 }
00052
00053 red[0]=0;
00054 nid=0;
00055 for (i=0;i<ns;i++){
00056 red[i+1]=0;
00057 for (j=0;j<top->stop->nnsd[i];j++){
00058 ndofn = top->give_ndofn (nid);
00059 for (k=0;k<ndofn;k++){
00060 dof = top->give_dof (nid,k);
00061 if (dof>red[i+1]){
00062 red[i+1]=dof;
00063 }
00064 }
00065 nid++;
00066 }
00067 }
00068
00069 for (i=0;i<ns;i++){
00070 for (j=0;j<ncdofd[i];j++){
00071 if (edofs[i][j]>0)
00072 edofs[i][j] -= red[i];
00073 else
00074 edofs[i][j] += red[i];
00075 }
00076 }
00077
00078 fprintf (out,"\n\n\n osel \n");
00079 for (i=0;i<ns;i++){
00080 fprintf (out,"\n");
00081 for (j=0;j<ncdofd[i];j++){
00082 fprintf (out,"\n edofs %ld",edofs[i][j]);
00083 }
00084 }
00085 fprintf (out,"\n\n\n osel \n");
00086
00087
00088
00089
00090
00091
00092
00093
00094 ccn = new long* [ns];
00095 for (i=0;i<ns;i++){
00096 ccn[i] = new long [ncdofd[i]];
00097 }
00098
00099 for (i=0;i<ns;i++){
00100 for (j=0;j<ncdofd[i];j++){
00101
00102 }
00103 }
00104
00105 delete [] red;
00106
00107
00108
00109
00110
00111
00112 nsid = new long [top->nn];
00113 l=0;
00114 for (i=0;i<ns;i++){
00115 for (j=0;j<top->stop->nnsd[i];j++){
00116 nsid[l]=i;
00117 l++;
00118 }
00119 }
00120
00121
00122 ndofdom = new long [ns];
00123 for (i=0;i<ns;i++){
00124 ndofdom[i]=0;
00125 }
00126
00127 for (i=0;i<top->nn;i++){
00128 ndofn=top->give_ndofn (i);
00129 l=nsid[i];
00130 for (j=0;j<ndofn;j++){
00131 k=top->give_dof (i,j);
00132 if (k>0)
00133 ndofdom[l]++;
00134 }
00135 }
00136
00137 cndom = new long* [ns];
00138 for (i=0;i<ns;i++){
00139 cndom[i] = new long [ndofdom[i]];
00140 ndofdom[i]=0;
00141 }
00142
00143 for (i=0;i<top->nn;i++){
00144 ndofn=top->give_ndofn (i);
00145 l=nsid[i];
00146 for (j=0;j<ndofn;j++){
00147 k=top->give_dof (i,j);
00148 if (k>0){
00149 cndom[l][ndofdom[l]]=k;
00150 ndofdom[l]++;
00151 }
00152 }
00153 }
00154
00155
00156 }
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 void saddpoint::matrix_e ()
00185 {
00186
00187 }
00188
00189
00190 void saddpoint::get_jumps (long nc,double *jumps)
00191 {
00192 long i;
00193
00194 if (c==NULL)
00195 c = new double [nc];
00196
00197 for (i=0;i<nc;i++){
00198 c[i]=jumps[i];
00199
00200
00201
00202 }
00203 }
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 void saddpoint::solve_system (gtopology *top,gmatrix *gm,double *lhs,double *rhs,FILE *out)
00217 {
00218 long i,j;
00219 double limit=1.0e-7;
00220 double *condmat,*condvect,*x,*y;
00221
00222
00223 nrdof=top->nbdof;
00224
00225 nidof=top->nidof;
00226
00227
00228 condmat = new double [nrdof*nrdof];
00229
00230 condvect = new double [nrdof];
00231
00232
00233 gm->condense (top,condmat,condvect,lhs,rhs,nrdof,1,out);
00234
00235
00236
00237 n=nrdof+nrdof/2;
00238
00239
00240 nm=nrdof/2;
00241
00242
00243 dm->alloc (n);
00244 x = new double [n];
00245 y = new double [n];
00246
00247 for (i=0;i<nrdof;i++){
00248 for (j=0;j<nrdof;j++){
00249 dm->a[i*n+j]=condmat[i*nrdof+j];
00250 }
00251 }
00252
00253
00254 for (i=0;i<nm;i++){
00255 dm->a[i*n+nrdof+i]=1.0;
00256 }
00257 j=nm;
00258 for (i=0;i<nm;i++){
00259 dm->a[j*n+nrdof+i]=-1.0;
00260 j++;
00261 }
00262
00263 j=nrdof;
00264 for (i=0;i<nm;i++){
00265 dm->a[j*n+i]=1.0;
00266 j++;
00267 }
00268
00269 j=nrdof;
00270 for (i=0;i<nm;i++){
00271 dm->a[j*n+nm+i]=-1.0;
00272 j++;
00273 }
00274
00275 j=nidof;
00276 for (i=0;i<nrdof;i++){
00277 y[i]=rhs[j];
00278 j++;
00279 }
00280
00281 j=nrdof;
00282 for (i=0;i<nm;i++){
00283 y[j]=c[i];
00284 j++;
00285 }
00286
00287
00288 dm->gemp (x,y,1,limit,1);
00289
00290 gm->condense (top,condmat,x,lhs,rhs,nrdof,2,out);
00291
00292
00293 dm->dealloc();
00294 delete [] x;
00295 delete [] y;
00296 delete [] condmat;
00297 delete [] condvect;
00298 }