00001 #include "mdsolver.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 <string.h> 00009 00010 /** 00011 function solves problems of molecular dynamics 00012 00013 JK, 19.6.2005 00014 */ 00015 void solve_molecular_dynamics () 00016 { 00017 00018 00019 00020 00021 00022 00023 00024 00025 00026 00027 00028 00029 00030 00031 00032 00033 00034 00035 00036 00037 00038 00039 long i; 00040 double *lhs,*rhs; 00041 00042 // stiffness matrix assembling 00043 stiffness_matrix (0); 00044 00045 // array is allocated in Lsrs 00046 lhs = Lsrs->give_lhs (0); 00047 rhs = Lsrs->give_rhs (0); 00048 00049 // loop over load cases 00050 for (i=0;i<Lsrs->nlc;i++){ 00051 // right hand side of system 00052 mefel_right_hand_side (i,rhs); 00053 } 00054 00055 /* 00056 Smat->printmat(Out); 00057 FILE *p; 00058 p = fopen ("vektor.txt","w"); 00059 fprintf (p,"%ld\n",Ndofm); 00060 for (i=0;i<Ndofm;i++){ 00061 fprintf (p,"%le\n",rhs[i]); 00062 } 00063 fclose (p); 00064 */ 00065 00066 // solution of equation system 00067 00068 for (i=0;i<Lsrs->nlc;i++){ 00069 Smat->solve_system (Lsrs->give_lhs(i),Lsrs->give_rhs(i)); 00070 } 00071 00072 // zacatek Risoviny 00073 /* 00074 long j; 00075 long n=6; 00076 double *condmat,*condvect,*condsolv; 00077 condsolv = new double [n]; 00078 condvect = new double [n]; 00079 condmat = new double [n*n]; 00080 nullvr (condmat,n*n); 00081 nullvr (condvect,n); 00082 nullvr (condsolv,n); 00083 00084 //Smat->condense (Gtm,condmat,lhs,rhs,3,1); 00085 Smat->condense (Gtm,condmat,lhs,rhs,n,1); 00086 00087 fprintf (Out,"\n\n kondenzovany vektor"); 00088 j=0; 00089 for (i=Ndofm-n;i<Ndofm;i++){ 00090 condvect[j]=rhs[i]; 00091 fprintf (Out,"\n %e",condvect[j]); 00092 j++; 00093 } 00094 00095 fprintf (Out,"\n\n kondenzovana matice"); 00096 for (i=0;i<n;i++){ 00097 fprintf (Out,"\n"); 00098 for (j=0;j<n;j++){ 00099 fprintf (Out," %e",condmat[i*n+j]); 00100 } 00101 } 00102 00103 gemp (condmat,condsolv,condvect,n,1,1.0e-15,1); 00104 00105 fprintf (Out,"\n\n vyresena cast vysledku"); 00106 j=0; 00107 for (i=Ndofm-n;i<Ndofm;i++){ 00108 lhs[i]=condsolv[j]; 00109 fprintf (Out,"\n %e",condsolv[j]); 00110 j++; 00111 } 00112 00113 // back substitution on subdomains 00114 //Smat->condense (Gtm,condmat,lhs,rhs,3,2); 00115 Smat->condense (Gtm,condmat,lhs,rhs,n,2); 00116 00117 delete [] condmat; 00118 delete [] condvect; 00119 delete [] condsolv; 00120 */ 00121 // konec Risoviny 00122 00123 00124 /* 00125 // zacatek JK testu 00126 long nse=6,dim; 00127 long *se; 00128 double *rbm; 00129 00130 se = new long [nse]; 00131 rbm = new double [nse*Ndofm]; 00132 00133 Smat->kernel (rbm,dim,se,nse,1.0e-3,3); 00134 00135 fprintf (Out,"\n\n\n RIGID BODY MOTIONS \n"); 00136 fprintf (Out,"\n pocet RBM %ld\n",dim); 00137 for (i=0;i<dim;i++){ 00138 fprintf (Out," %ld",se[i]); 00139 } 00140 for (i=0;i<dim;i++){ 00141 fprintf (Out,"\n\n RBM number %ld",i); 00142 for (j=0;j<Ndofm;j++){ 00143 fprintf (Out,"\n %4ld %4ld %e",i,j,rbm[i*Ndofm+j]); 00144 } 00145 } 00146 // konec JK testu 00147 */ 00148 /* 00149 // loop over load cases 00150 for (i=0;i<Lsrs->nlc;i++){ 00151 00152 // strains 00153 if (Mp->straincomp==1){ 00154 Mm->computestrains (i); 00155 Mm->stra.transformvalues (1); 00156 } 00157 00158 // stresses 00159 if (Mp->stresscomp==1){ 00160 Mm->computestresses (i); 00161 Mm->stre.transformvalues(0); 00162 } 00163 00164 // reactions 00165 if (Mp->reactcomp==1){ 00166 Mb->lc[i].compute_reactions (i); 00167 } 00168 00169 00170 } 00171 */ 00172 00173 print_init(-1, "wt"); 00174 for (i=0;i<Lsrs->nlc;i++){ 00175 // computes and prints required quantities 00176 Mm->allipstrains (i); 00177 Mm->allipstresses (i); 00178 compute_req_val (i); 00179 print_step(i, 0, 0.0, NULL); 00180 } 00181 print_close(); 00182 00183 00184 /* 00185 for (i=0;i<Ndofm;i++){ 00186 printf ("\n posun %20.10lf",Lsrs->lhs[i]); 00187 } 00188 */ 00189 00190 } 00191