 |
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
|
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
00001 /************************************************************************ 00002 * MechSys - Open Library for Mechanical Systems * 00003 * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand * 00004 * Copyright (C) 2009 Sergio Galindo * 00005 * * 00006 * This program is free software: you can redistribute it and/or modify * 00007 * it under the terms of the GNU General Public License as published by * 00008 * the Free Software Foundation, either version 3 of the License, or * 00009 * any later version. * 00010 * * 00011 * This program is distributed in the hope that it will be useful, * 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00014 * GNU General Public License for more details. * 00015 * * 00016 * You should have received a copy of the GNU General Public License * 00017 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00018 ************************************************************************/ 00019 00020 #ifndef MECHSYS_FEM_UWPSOLVER_H 00021 #define MECHSYS_FEM_UWPSOLVER_H 00022 00023 // MechSys 00024 #include <mechsys/fem/solver.h> 00025 #include <mechsys/fem/uwpelem.h> 00026 00027 using std::cout; 00028 using std::endl; 00029 00030 namespace FEM 00031 { 00032 00033 class UWPSolver : public Solver 00034 { 00035 public: 00036 // Constructor 00037 UWPSolver (Domain & Dom, SDPair const & Flags, pOutFun OutFun=NULL, void * OutDat=NULL, 00038 pOutFun DbgFun=NULL, void * DbgDat=NULL); 00039 00040 // Methods 00041 String Name () const { return "UWPSolver"; } 00042 void Initialize (); 00043 void Assemble (double Time, Vec_t const & V, Vec_t const & Ww, Vec_t const & Pw); 00044 void RKFunc (double Time, Vec_t const & V, Vec_t const & Ww, Vec_t const & Pw); 00045 void Solve (double tf, double dt, double dtOut, char const * FileKey=NULL); 00046 00047 // Data 00048 Array<long> pEqU, pEqW, pEqP; 00049 Array<bool> pU, pW, pP; 00050 Vec_t U, V, Ww, Pw; 00051 Vec_t V_fe, Ww_fe, Pw_fe; 00052 Vec_t V_me, Ww_me, Pw_me; 00053 Vec_t V_er, Ww_er, Pw_er; 00054 Vec_t DV, DWw, DPw; 00055 Vec_t Rv, Rww; 00056 Vec_t dVdt, dWwdt, dPwdt; 00057 Vec_t f, fw, hw, cw, ew; 00058 Sparse::Triplet<double,int> M, Mb, Mw, Hw; 00059 bool WithInfo; 00060 double STOL, mMin, mMax, dTini, dTlast; 00061 size_t MaxSS; 00062 }; 00063 00064 00066 00067 00068 inline UWPSolver::UWPSolver (Domain & Dom, SDPair const & Flags, pOutFun OutFun, void * OutDat, pOutFun DbgFun, void * DbgDat) 00069 : Solver (Dom, Flags, OutFun, OutDat, DbgFun, DbgDat), WithInfo(true), 00070 STOL(1.0e-2), mMin(0.2), mMax(2.0), dTini(0.01), dTlast(-1.0), MaxSS(2000) 00071 { 00072 } 00073 00074 inline void UWPSolver::Initialize () 00075 { 00076 // assign equation numbers 00077 size_t NEqU = 0; 00078 size_t NEqW = 0; 00079 size_t NEqP = 0; 00080 for (size_t i=0; i<Dom.ActNods.Size(); ++i) 00081 for (size_t j=0; j<Dom.ActNods[i]->NDOF(); ++j) 00082 { 00083 if (Dom.ActNods[i]->UKey(j)=="wwx" || Dom.ActNods[i]->UKey(j)=="wwy" || Dom.ActNods[i]->UKey(j)=="wwz") 00084 { 00085 Dom.ActNods[i]->Eq(j) = NEqW; 00086 NEqW++; 00087 } 00088 else if (Dom.ActNods[i]->UKey(j)=="pw") 00089 { 00090 Dom.ActNods[i]->Eq(j) = NEqP; 00091 NEqP++; 00092 } 00093 else 00094 { 00095 Dom.ActNods[i]->Eq(j) = NEqU; 00096 NEqU++; 00097 } 00098 } 00099 00100 // prescribed equations 00101 pEqU.Resize(0); pU.Resize(NEqU); pU.SetValues(false); 00102 pEqW.Resize(0); pW.Resize(NEqW); pW.SetValues(false); 00103 pEqP.Resize(0); pP.Resize(NEqP); pP.SetValues(false); 00104 for (size_t i=0; i<Dom.NodsWithPU.Size(); ++i) 00105 { 00106 Node const * nod = Dom.NodsWithPU[i]; 00107 for (size_t j=0; j<nod->NPU(); ++j) 00108 { 00109 int eq = nod->EqPU(j); 00110 if (nod->PUKey(j)=="wwx" || nod->PUKey(j)=="wwy" || nod->PUKey(j)=="wwz") { pEqW.Push(eq); pW[eq]=true; } 00111 else if (nod->PUKey(j)=="pw") { pEqP.Push(eq); pP[eq]=true; } 00112 else { pEqU.Push(eq); pU[eq]=true; } 00113 } 00114 } 00115 //cout << "pEqU = " << pEqU << endl; 00116 //cout << "pEqW = " << pEqW << endl; 00117 //cout << "pEqP = " << pEqP << endl; 00118 00119 // find total number of non-zero entries, including duplicates 00120 size_t Mb_size = 0; 00121 size_t Mw_size = 0; 00122 size_t Hw_size = 0; 00123 for (size_t k=0; k<Dom.ActEles.Size(); ++k) 00124 { 00125 if (Dom.ActEles[k]->IsUWP) 00126 { 00127 Dom.ActEles[k]->GetLoc (); 00128 for (size_t i=0; i<UWPElem::LocU .Size(); ++i) 00129 for (size_t j=0; j<UWPElem::LocU .Size(); ++j) Mb_size++; 00130 for (size_t i=0; i<UWPElem::LocWw.Size(); ++i) 00131 for (size_t j=0; j<UWPElem::LocWw.Size(); ++j) Mw_size++; 00132 for (size_t i=0; i<UWPElem::LocPw.Size(); ++i) 00133 for (size_t j=0; j<UWPElem::LocPw.Size(); ++j) Hw_size++; 00134 } 00135 else throw new Fatal("UWPElem::Initialize: this method does not work with elements that are not UWP yet"); 00136 } 00137 00138 // resize vectors and matrices 00139 U . change_dim (NEqU); 00140 V . change_dim (NEqU); 00141 Ww . change_dim (NEqW); 00142 Pw . change_dim (NEqP); 00143 V_fe . change_dim (NEqU); 00144 Ww_fe . change_dim (NEqW); 00145 Pw_fe . change_dim (NEqP); 00146 V_me . change_dim (NEqU); 00147 Ww_me . change_dim (NEqW); 00148 Pw_me . change_dim (NEqP); 00149 V_er . change_dim (NEqU); 00150 Ww_er . change_dim (NEqW); 00151 Pw_er . change_dim (NEqP); 00152 DV . change_dim (NEqU); 00153 DWw . change_dim (NEqW); 00154 DPw . change_dim (NEqP); 00155 Rv . change_dim (NEqU); 00156 Rww . change_dim (NEqW); 00157 dVdt . change_dim (NEqU); 00158 dWwdt . change_dim (NEqW); 00159 dPwdt . change_dim (NEqP); 00160 f . change_dim (NEqU); 00161 fw . change_dim (NEqW); 00162 hw . change_dim (NEqW); 00163 cw . change_dim (NEqU); 00164 ew . change_dim (NEqP); 00165 M . AllocSpace (NEqU,NEqU,Mb_size); 00166 Mb . AllocSpace (NEqU,NEqU,Mb_size); 00167 Mw . AllocSpace (NEqW,NEqW,Mw_size); 00168 Hw . AllocSpace (NEqP,NEqP,Hw_size); 00169 00170 // set initial values (from previous stage) 00171 for (size_t i=0; i<Dom.ActNods.Size(); ++i) 00172 for (size_t j=0; j<Dom.ActNods[i]->NDOF(); ++j) 00173 { 00174 if (Dom.ActNods[i]->UKey(j)=="wwx" || Dom.ActNods[i]->UKey(j)=="wwy" || Dom.ActNods[i]->UKey(j)=="wwz") 00175 { 00176 Ww(Dom.ActNods[i]->Eq(j)) = Dom.ActNods[i]->U(j); 00177 } 00178 else if (Dom.ActNods[i]->UKey(j)=="pw") 00179 { 00180 Pw(Dom.ActNods[i]->Eq(j)) = Dom.ActNods[i]->U(j); 00181 } 00182 else 00183 { 00184 U(Dom.ActNods[i]->Eq(j)) = Dom.ActNods[i]->U(j); 00185 V(Dom.ActNods[i]->Eq(j)) = Dom.ActNods[i]->V(j); 00186 } 00187 } 00188 } 00189 00190 inline void UWPSolver::Assemble (double Time, Vec_t const & TheV, Vec_t const & TheWw, Vec_t const & ThePw) 00191 { 00192 M .ResetTop (); 00193 Mb.ResetTop (); 00194 Mw.ResetTop (); 00195 Hw.ResetTop (); 00196 set_to_zero (f); 00197 set_to_zero (fw); 00198 set_to_zero (hw); 00199 set_to_zero (cw); 00200 set_to_zero (ew); 00201 for (size_t k=0; k<Dom.ActEles.Size(); ++k) 00202 { 00203 if (Dom.ActEles[k]->IsUWP) 00204 { 00205 Dom.ActEles[k]->GetLoc (); 00206 Dom.ActEles[k]->ElemEqs (Time, TheV, TheWw, ThePw); 00207 00208 for (size_t i=0; i<UWPElem::LocU .Size(); ++i) 00209 { 00210 f (UWPElem::LocU[i]) += UWPElem::f (i); 00211 cw(UWPElem::LocU[i]) += UWPElem::cw(i); 00212 for (size_t j=0; j<UWPElem::LocU.Size(); ++j) 00213 { 00214 M.PushEntry (UWPElem::LocU[i], UWPElem::LocU[j], UWPElem::M(i,j)); 00215 //if (!pU[UWPElem::LocU[i]] && !pU[UWPElem::LocU[j]]) Mb.PushEntry (UWPElem::LocU[i], UWPElem::LocU[j], UWPElem::M(i,j) - UWPElem::Mw(i,j)); 00216 Mb.PushEntry (UWPElem::LocU[i], UWPElem::LocU[j], UWPElem::M(i,j) - UWPElem::Mw(i,j)); 00217 } 00218 } 00219 00220 for (size_t i=0; i<UWPElem::LocWw.Size(); ++i) 00221 { 00222 fw(UWPElem::LocWw[i]) += UWPElem::fw(i); 00223 hw(UWPElem::LocWw[i]) += UWPElem::hw(i); 00224 for (size_t j=0; j<UWPElem::LocWw.Size(); ++j) 00225 { 00226 //if (!pW[UWPElem::LocWw[i]] && !pW[UWPElem::LocWw[j]]) Mw.PushEntry (UWPElem::LocWw[i], UWPElem::LocWw[j], UWPElem::Mw(i,j)); 00227 Mw.PushEntry (UWPElem::LocWw[i], UWPElem::LocWw[j], UWPElem::Mw(i,j)); 00228 } 00229 } 00230 00231 for (size_t i=0; i<UWPElem::LocPw.Size(); ++i) 00232 { 00233 ew(UWPElem::LocPw[i]) += UWPElem::ew(i); 00234 for (size_t j=0; j<UWPElem::LocPw.Size(); ++j) 00235 { 00236 //if (!pP[UWPElem::LocPw[i]] && !pP[UWPElem::LocPw[j]]) Hw.PushEntry (UWPElem::LocPw[i], UWPElem::LocPw[j], UWPElem::Hw(i,j)); 00237 Hw.PushEntry (UWPElem::LocPw[i], UWPElem::LocPw[j], UWPElem::Hw(i,j)); 00238 } 00239 } 00240 00241 //cout << "M = \n" << PrintMatrix(UWPElem::M); 00242 //cout << "Mw = \n" << PrintMatrix(UWPElem::Mw); 00243 //cout << "Hw = \n" << PrintMatrix(UWPElem::Hw); 00244 //cout << "f = \n" << PrintVector(UWPElem::f); 00245 //cout << "fw = \n" << PrintVector(UWPElem::fw); 00246 //cout << "hw = \n" << PrintVector(UWPElem::hw); 00247 //cout << "cw = \n" << PrintVector(UWPElem::cw); 00248 //cout << "ew = \n" << PrintVector(UWPElem::ew); 00249 } 00250 else 00251 { 00252 throw new Fatal("UWPElem::Assembly: this method does not work with elements that are not u-w-p yet"); 00253 Array<size_t> loc; 00254 Mat_t M; 00255 Dom.ActEles[k]->GetLoc (loc); 00256 Dom.ActEles[k]->CalcM (M); 00257 } 00258 } 00259 00260 //for (size_t i=0; i<pEqU.Size(); ++i) Mb.PushEntry (pEqU[i],pEqU[i], 1.0); 00261 //for (size_t i=0; i<pEqW.Size(); ++i) Mw.PushEntry (pEqW[i],pEqW[i], 1.0); 00262 //for (size_t i=0; i<pEqP.Size(); ++i) Hw.PushEntry (pEqP[i],pEqP[i], 1.0); 00263 00264 /* 00265 M .WriteSMAT("M"); 00266 Mb.WriteSMAT("Mb"); 00267 Mw.WriteSMAT("Mw"); 00268 Hw.WriteSMAT("Hw"); 00269 Sparse::Matrix<double,int> sM(M), sMb(Mb), sMw(Mw), sHw(Hw); 00270 Mat_t Mmat, Mbmat, Mwmat, Hwmat, Mi, Mbi, Mwi, Hwi; 00271 sM .GetDense (Mmat); 00272 sMb.GetDense (Mbmat); 00273 sMw.GetDense (Mwmat); 00274 sHw.GetDense (Hwmat); 00275 cout << "M = \n" << PrintMatrix(Mmat, "%10.5f"); 00276 cout << "Mb = \n" << PrintMatrix(Mbmat, "%10.5f"); 00277 cout << "Mw = \n" << PrintMatrix(Mwmat, "%10.5f"); 00278 cout << "Hw = \n" << PrintMatrix(Hwmat, "%10.5f"); 00279 Inv (Mmat, Mi); 00280 Inv (Mbmat, Mbi); 00281 Inv (Mwmat, Mwi); 00282 Inv (Hwmat, Hwi); 00283 cout << "Inv(M) = \n" << PrintMatrix(Mi, "%10.5f"); 00284 cout << "Inv(Mbi) = \n" << PrintMatrix(Mbi, "%10.5f"); 00285 cout << "Inv(Mw) = \n" << PrintMatrix(Mwi, "%10.5f"); 00286 cout << "Inv(Hw) = \n" << PrintMatrix(Hwi, "%10.1f"); 00287 cout << "det(M) = " << UMFPACK::Det (sM) << endl; 00288 cout << "det(Mb) = " << UMFPACK::Det (sMb) << endl; 00289 cout << "det(Mw) = " << UMFPACK::Det (sMw) << endl; 00290 cout << "det(Hw) = " << UMFPACK::Det (sHw) << endl; 00291 */ 00292 } 00293 00294 inline void UWPSolver::RKFunc (double Time, Vec_t const & TheV, Vec_t const & TheWw, Vec_t const & ThePw) 00295 { 00296 Assemble (Time, TheV, TheWw, ThePw); 00297 00298 Rv = f + hw - fw; 00299 Rww = f - cw; 00300 00301 /* 00302 cout << "f = " << PrintVector(f, "%10.5f"); 00303 cout << "fw = " << PrintVector(fw, "%10.5f"); 00304 cout << "hw = " << PrintVector(hw, "%10.5f"); 00305 cout << "cw = " << PrintVector(cw, "%10.5f"); 00306 cout << "ew = " << PrintVector(ew, "%10.5f"); 00307 cout << "Rv = " << PrintVector(Rv, "%10.5f"); 00308 cout << "Rww = " << PrintVector(Rww, "%10.5f"); 00309 */ 00310 00311 UMFPACK::Solve (Mb, Rv, dVdt); // dVdt = inv(Mb)*Rv 00312 00313 for (size_t i=0; i<pEqU.Size(); ++i) dVdt(pEqU[i]) = 0.0; 00314 00315 Sparse::SubMult (M, dVdt, Rww); // Rww -= M*dVdt 00316 UMFPACK::Solve (Mw, Rww, dWwdt); // dWwdt = inv(Mw)*Rww 00317 00318 for (size_t i=0; i<pEqW.Size(); ++i) dWwdt(pEqW[i]) = 0.0; 00319 00320 UMFPACK::Solve (Hw, ew, dPwdt); // dPwdt = inv(Hw)*ew 00321 00322 // set prescribed U, V and A 00323 for (size_t i=0; i<Dom.NodsWithPU.Size(); ++i) 00324 { 00325 Node * const nod = Dom.NodsWithPU[i]; 00326 for (size_t j=0; j<nod->NPU(); ++j) 00327 { 00328 if (nod->PUKey(j)=="pw") 00329 { 00330 int eq = nod->EqPU(j); 00331 dPwdt(eq) = nod->PV(j, Time); 00332 } 00333 } 00334 } 00335 00336 /* 00337 cout << "Rv (after) = " << PrintVector(Rv, "%10.5f"); 00338 cout << "Rww (after) = " << PrintVector(Rww, "%10.5f"); 00339 cout << "dVdt = " << PrintVector(dVdt, "%10.5f"); 00340 cout << "dWwdt = " << PrintVector(dWwdt, "%10.5f"); 00341 cout << "dPwdt = " << PrintVector(dPwdt, "%10.5f"); 00342 cout << "V = " << PrintVector(V, "%10.5f"); 00343 cout << "Ww = " << PrintVector(Ww, "%10.5f"); 00344 cout << "Pw = " << PrintVector(Pw, "%10.5f") << "\n"; 00345 */ 00346 } 00347 00348 inline void UWPSolver::Solve (double tf, double, double dtOut, char const * FileKey) 00349 { 00350 // info 00351 Util::Stopwatch stopwatch(WithInfo); 00352 00353 // initialize global matrices and vectors 00354 Initialize (); 00355 00356 // output initial state 00357 printf ("\n%s--- Stage solution --- u-w-p Solver -- Runge-Kutta --------------------------%s\n",TERM_CLR1,TERM_RST); 00358 printf ("%s%10s %4s%s\n",TERM_CLR2,"Time","NSS",TERM_RST); 00359 printf ("%10.6f %4d\n",Dom.Time,0); 00360 if (Dom.IdxOut==0) 00361 { 00362 Dom.OutResults (FileKey); 00363 if (OutFun!=NULL) (*OutFun) (this, OutDat); 00364 if (DbgFun!=NULL) (*DbgFun) (this, DbgDat); 00365 } 00366 00367 // time for output 00368 double tout = Dom.Time + dtOut; 00369 while (Dom.Time<tf) 00370 { 00371 // update to tout 00372 size_t stp = 0; 00373 double T = 0.0; 00374 double dT = (dTlast>0.0 ? dTlast : dTini); 00375 double Dt = tout - Dom.Time; 00376 for (stp=0; stp<MaxSS; ++stp) 00377 { 00378 // exit 00379 if (T>=1.0) break; 00380 00381 // time increment 00382 double dt = dT*Dt; 00383 00384 // FE state 00385 RKFunc (Dom.Time, V, Ww, Pw); 00386 //DU = V * dt; 00387 DV = dVdt * dt; 00388 DWw = dWwdt * dt; 00389 DPw = dPwdt * dt; 00390 //U_fe = U + DU; 00391 V_fe = V + DV; 00392 Ww_fe = Ww + DWw; 00393 Pw_fe = Pw + DPw; 00394 for (size_t k=0; k<Dom.ActEles.Size(); ++k) Dom.ActEles[k]->Update (0, V, dPwdt, dt); // should update U here as well (because of large displacements shape functions) 00395 00396 // ME State 00397 RKFunc (Dom.Time+dt, V_fe, Ww_fe, Pw_fe); 00398 V_me = V + 0.5*DV + (0.5*dt)*dVdt; 00399 Ww_me = Ww + 0.5*DWw + (0.5*dt)*dWwdt; 00400 Pw_me = Pw + 0.5*DPw + (0.5*dt)*dPwdt; 00401 for (size_t k=0; k<Dom.ActEles.Size(); ++k) Dom.ActEles[k]->Update (1, V_fe, dPwdt, dt); 00402 00403 // error estimate 00404 V_er = V_me - V_fe; 00405 Ww_er = Ww_me - Ww_fe; 00406 Pw_er = Pw_me - Pw_fe; 00407 double V_error = Norm(V_er) / (1.0+Norm(V_me)); 00408 double Ww_error = Norm(Ww_er) / (1.0+Norm(Ww_me)); 00409 double Pw_error = Norm(Pw_er) / (1.0+Norm(Pw_me)); 00410 double error = Util::Max (V_error, Ww_error, Pw_error); 00411 00412 // step multiplier 00413 double m = (error>0.0 ? 0.9*sqrt(STOL/error) : mMax); 00414 00415 //cout << "error = " << error << ", elems_error = " << elems_error << endl; 00416 //cout << "error = " << error << endl; 00417 //break; 00418 00419 // update 00420 if (error<STOL) 00421 { 00422 T += dT; 00423 Dom.Time += dt; 00424 U += V_me*dt; 00425 V = V_me; 00426 Ww = Ww_me; 00427 Pw = Pw_me; 00428 if (m>mMax) m = mMax; 00429 } 00430 else 00431 { 00432 for (size_t k=0; k<Dom.ActEles.Size(); ++k) Dom.ActEles[k]->Restore (); 00433 if (m<mMin) m = mMin; 00434 } 00435 dT = m * dT; 00436 if (dT>1.0-T) dT = 1.0-T; 00437 else dTlast = dT; 00438 } 00439 if (stp==MaxSS) throw new Fatal("UWPSolver::Solve: Runge-Kutta (2nd order / ME) did not converge after %zd substeps",stp); 00440 00441 // update nodes to tout 00442 for (size_t i=0; i<Dom.ActNods.Size(); ++i) 00443 for (size_t j=0; j<Dom.ActNods[i]->NDOF(); ++j) 00444 { 00445 if (Dom.ActNods[i]->UKey(j)=="wwx" || Dom.ActNods[i]->UKey(j)=="wwy" || Dom.ActNods[i]->UKey(j)=="wwz") 00446 { 00447 Dom.ActNods[i]->U(j) = Ww(Dom.ActNods[i]->Eq(j)); 00448 } 00449 else if (Dom.ActNods[i]->UKey(j)=="pw") 00450 { 00451 Dom.ActNods[i]->U(j) = Pw(Dom.ActNods[i]->Eq(j)); 00452 } 00453 else 00454 { 00455 Dom.ActNods[i]->U(j) = U(Dom.ActNods[i]->Eq(j)); 00456 Dom.ActNods[i]->V(j) = V(Dom.ActNods[i]->Eq(j)); 00457 } 00458 } 00459 00460 // output 00461 printf ("%10.6f %4zd\n",Dom.Time,stp); 00462 Dom.OutResults (FileKey); 00463 if (OutFun!=NULL) (*OutFun) (this, OutDat); 00464 00465 // next tout 00466 tout = Dom.Time + dtOut; 00467 } 00468 } 00469 00470 00472 00473 00474 Solver * UWPSolverMaker(Domain & Dom, SDPair const & Flags, Solver::pOutFun OutFun, void * OutDat, Solver::pOutFun DbgFun, void * DbgDat) 00475 { 00476 return new UWPSolver (Dom, Flags, OutFun, OutDat, DbgFun, DbgDat); 00477 } 00478 00479 int UWPSolverRegister() 00480 { 00481 SolverFactory["UWPSolver"] = UWPSolverMaker; 00482 return 0; 00483 } 00484 00485 int __UWPSolver_dummy_int = UWPSolverRegister(); 00486 00487 }; 00488 00489 #endif // MECHSYS_FEM_UWPSOLVER_H
1.7.6.1
