00001 #include "lssolver.h" 00002 #include "global.h" 00003 #include "globmat.h" 00004 #include "loadcase.h" 00005 #include "gmatrix.h" 00006 #include "mechprint.h" 00007 #include "node.h" 00008 #include "elemswitch.h" 00009 #include <string.h> 00010 00011 /** 00012 function solves linear statics problems 00013 00014 function contains the following steps: 00015 1. stiffness matrix assembling 00016 2. right hand side assembling (load vectors assembling) 00017 3. solution of equation systems by choosen method 00018 4. strain computation (eligible) 00019 5. stress computation (eligible) 00020 6. reaction computation (eligible) 00021 7. results printing 00022 00023 JK 00024 */ 00025 void solve_linear_statics () 00026 { 00027 long i,n; 00028 double *lhs,*rhs,*fl,*fi; 00029 00030 // stiffness matrix assembling 00031 stiffness_matrix (0); 00032 00033 /* 00034 n=Smat->n; 00035 //fprintf (Out,"\n\n jouda \n\n"); 00036 //fprintf (Out,"A=["); 00037 //for (long ijk=0;ijk<n;ijk++){ 00038 // for (long kji=0;kji<n;kji++){ 00039 // fprintf (Out," %le",Smat->dm->a[ijk*n+kji]); 00040 // } 00041 // fprintf (Out,";"); 00042 //} 00043 //fprintf (Out,"]"); 00044 */ 00045 00046 00047 // array is allocated in Lsrs 00048 lhs = Lsrs->give_lhs (0); 00049 rhs = Lsrs->give_rhs (0); 00050 00051 // number of mechanical degrees of freedom 00052 n=Ndofm; 00053 // vector of prescribed force loads, it does not contain forces caused by temperature, etc. 00054 fl = new double [n]; 00055 nullv (fl,n); 00056 // vector of resulting internal forces 00057 fi = new double [n]; 00058 nullv (fi,n); 00059 00060 // loop over load cases 00061 for (i=0;i<Lsrs->nlc;i++){ 00062 // right hand side of system 00063 mefel_right_hand_side (i,rhs); 00064 } 00065 00066 //Smat->printmat (Out); 00067 00068 /* 00069 FILE *out; 00070 out = fopen ("vektor.txt","w"); 00071 for (i=0;i<Ndofm;i++){ 00072 fprintf (out,"%20.15le\n",rhs[i]); 00073 } 00074 fclose (out); 00075 */ 00076 00077 /* 00078 Smat->printmat (Out); 00079 00080 out = fopen ("dof.txt","w"); 00081 for (i=0;i<Gtm->nn;i++){ 00082 fprintf (out,"%6ld %6ld %6ld %6ld\n",i,Gtm->gnodes[i].cn[0],Gtm->gnodes[i].cn[1],Gtm->gnodes[i].cn[2]); 00083 } 00084 fclose (out); 00085 00086 out = fopen ("sit.txt","w"); 00087 fprintf (out,"%6ld\n",Gtm->nn); 00088 for (i=0;i<Gtm->nn;i++){ 00089 fprintf (out,"%6ld %le %le %le\n",i,Gtm->gnodes[i].x,Gtm->gnodes[i].y,Gtm->gnodes[i].z); 00090 } 00091 fprintf (out,"%6ld\n",Gtm->ne); 00092 for (i=0;i<Gtm->ne;i++){ 00093 fprintf (out,"%6ld %6ld %6ld %6ld %6ld\n",i,Gtm->gelements[i].nodes[0],Gtm->gelements[i].nodes[1],Gtm->gelements[i].nodes[2],Gtm->gelements[i].nodes[3]); 00094 } 00095 fclose (out); 00096 */ 00097 00098 /* 00099 FILE *out; 00100 out = fopen ("matice_pred.m","w"); 00101 Smat->printmat (out); 00102 fclose (out); 00103 */ 00104 00105 00106 // solution of equation system 00107 for (i=0;i<Lsrs->nlc;i++){ 00108 Mp->ssle->solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out); 00109 } 00110 00111 /* 00112 out = fopen ("reseni.txt","w"); 00113 for (i=0;i<Ndofm;i++){ 00114 fprintf (out,"%8ld % 20.15le\n",i,lhs[i]); 00115 } 00116 fclose (out); 00117 */ 00118 00119 /* 00120 FILE *out; 00121 out = fopen ("posuny.txt","w"); 00122 for (i=0;i<Gtm->nn;i++){ 00123 fprintf (out,"%6ld",i+1); 00124 00125 if (Gtm->gnodes[i].cn[0]==-1) 00126 fprintf (out," % 12.9le",-0.008); 00127 if (Gtm->gnodes[i].cn[0]==0) 00128 fprintf (out," % 12.9le",0.0); 00129 if (Gtm->gnodes[i].cn[0]>0) 00130 fprintf (out," % 12.9le",Lsrs->lhs[Gtm->gnodes[i].cn[0]-1]); 00131 00132 if (Gtm->gnodes[i].cn[1]==-1) 00133 fprintf (out," % 12.9le",-0.008); 00134 if (Gtm->gnodes[i].cn[1]==0) 00135 fprintf (out," % 12.9le",0.0); 00136 if (Gtm->gnodes[i].cn[1]>0) 00137 fprintf (out," % 12.9le",Lsrs->lhs[Gtm->gnodes[i].cn[1]-1]); 00138 00139 if (Gtm->gnodes[i].cn[2]==-1) 00140 fprintf (out," % 12.9le",-0.008); 00141 if (Gtm->gnodes[i].cn[2]==0) 00142 fprintf (out," % 12.9le",0.0); 00143 if (Gtm->gnodes[i].cn[2]>0) 00144 fprintf (out," % 12.9le",Lsrs->lhs[Gtm->gnodes[i].cn[2]-1]); 00145 00146 fprintf (out,"\n"); 00147 } 00148 fclose (out); 00149 */ 00150 00151 /* 00152 out = fopen ("matice_po.m","w"); 00153 Smat->printmat (out); 00154 fclose (out); 00155 */ 00156 00157 00158 print_init(-1, "wt"); 00159 for (i=0;i<Lsrs->nlc;i++){ 00160 // indicator of strain computation 00161 // it means, strains at integration points have not been computed 00162 Mp->strainstate=0; 00163 // indicator of stress computation 00164 // it means, stresses at integration points have not been computed 00165 Mp->stressstate=0; 00166 // indicator of computation of other array 00167 // it means, stresses at integration points have not been computed 00168 Mp->otherstate=0; 00169 // prints required quantities 00170 // computes required quantities 00171 compute_req_val (i); 00172 00173 if (Outdm->nog.selnforce.st != sel_no) // nodal forces are required 00174 { 00175 // rhs was rewritten by the solver call -> rebuild load vector 00176 mefel_right_hand_side (i,Lsrs->give_rhs(i),fl); 00177 } 00178 00179 // in artificial time 0.0 the nodal forces contains only the load vector fl 00180 print_step(i, 0, 0.0, fl); 00181 00182 // if the nodal forces are required for graphical output 00183 // reaction + load vector are printed in the artificial time 1.0 00184 if (Outdm->nog.selnforce.st != sel_no) 00185 { 00186 // indicator of strain computation 00187 // it means, strains at integration points have not been computed 00188 Mp->strainstate=0; 00189 // indicator of stress computation 00190 // it means, stresses at integration points have not been computed 00191 Mp->stressstate=0; 00192 // indicator of computation of other array 00193 // it means, stresses at integration points have not been computed 00194 Mp->otherstate=0; 00195 // calculate nodal forces and reactions 00196 internal_forces(i, fi); 00197 // in artificial time 1.0 the nodal forces contains the load vector + reactions 00198 print_step(i, 1, 1.0, fi); 00199 } 00200 } 00201 print_close(); 00202 00203 00204 if (Mp->adaptivityflag) 00205 Ada->run (2,2,0); 00206 00207 delete [] fl; 00208 } 00209