 |
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 * 00004 * * 00005 * This program is free software: you can redistribute it and/or modify * 00006 * it under the terms of the GNU General Public License as published by * 00007 * the Free Software Foundation, either version 3 of the License, or * 00008 * any later version. * 00009 * * 00010 * This program is distributed in the hope that it will be useful, * 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00013 * GNU General Public License for more details. * 00014 * * 00015 * You should have received a copy of the GNU General Public License * 00016 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00017 ************************************************************************/ 00018 00019 /* MPoints2D - Copyright (C) 2007 Dorival de Moraes Pedroso */ 00020 00021 #ifndef MPM_MPOINTS2D_H 00022 #define MPM_MPOINTS2D_H 00023 00024 // STL 00025 #include <cmath> // for sqrt, pow, etc. 00026 #include <cfloat> // for DBL_EPSILON 00027 00028 // MechSys 00029 #include <mechsys/util/array.h> 00030 #include <mechsys/models/model.h> 00031 #include <mechsys/models/equilibstate.h> 00032 00033 // Local 00034 #include <mechsys/mpm/defs.h> 00035 #include <mechsys/mpm/grid2d.h> 00036 #include <mechsys/mpm/tensors.h> 00037 00039 namespace MPM { 00040 00041 class MPoints2D 00042 { 00043 public: 00044 static const int NDIM = 2; 00045 typedef std::map<int,Model*> Models_t; 00046 00048 MPoints2D (size_t nPcell, Grid2D * G, 00049 ptIsPointInGeom pIsPointInGeom, 00050 ptDensity pDensity, 00051 ptModelData pModelData, 00052 ptIniVelocity pIniVelocity, 00053 ptHasTraction pHasTraction, 00054 ptHasAppDisp pHasAppDisp, 00055 ptAppDisp pAppDisp=NULL); 00056 00057 // Destructor 00058 ~MPoints2D (); 00059 00060 // Methods 00061 size_t nPoints () const { return _p.Size(); } 00062 size_t nPointsOnBry () const { return _ipbry.Size(); } 00063 void SetGrid (Grid2D * G) { _g = G; } 00064 void ReInit () { _initialize_points(); } 00065 00066 // Update method 00067 bool TimeUpdate (bool FwEuler, double Dt, ptB pB, ptLdM pLdM, double & t, double & DsE); 00068 00069 // Set Methods 00070 void SetCPDI (bool UseCPDI=true) { _cpdi = UseCPDI; } 00071 void SetMPM (bool UseMPM=true) { _mpm = UseMPM; } 00072 void SetUSF (bool DoUSF =true) { _usf = DoUSF; } 00073 00074 // Access methods 00075 Array<ATensor2> const & Fs () const { return _F; } 00076 Array<Vector3D> const & Ps () const { return _p; } 00077 Array<Vector3D> const & v () const { return _v; } 00078 Array<Vector3D> const & u () const { return _u; } 00079 Vector3D const & P (size_t i) const { return _p[i]; } 00080 Vector3D const & B (size_t i) const { return _p[_ipbry[i]]; } 00081 double const & m (size_t i) const { return _m[i]; } 00082 Vector3D const & f (size_t i) const { return _f[i]; } 00083 Vector3D const & t (size_t i) const { return _t[i]; } 00084 Vector3D const & fu (size_t i) const { return _fu[i]; } 00085 Vector3D const & u (size_t i) const { return _u[i]; } 00086 Vector3D const & v (size_t i) const { return _v[i]; } 00087 Vec_t const & s (size_t i) const { return _sta[i]->Sig; } 00088 Vec_t const & e (size_t i) const { return _sta[i]->Eps; } 00089 bool HasT (size_t i) const { return _has_t[i]; } 00090 bool HasU (size_t i) const { return _has_fu[i]; } 00091 int Clr (size_t i) const { return _clr[i]; } 00092 Array<EquilibState*> const & Sta () const { return _sta; } 00093 00094 // Access corners 00095 mutable Array<Vector3D> C; 00096 void CalcC (size_t p) const 00097 { 00098 if (NRVEC==2) 00099 { 00100 C[0] = _p[p] - _sg[p].R[0] - _sg[p].R[1]; 00101 C[1] = _p[p] + _sg[p].R[0] - _sg[p].R[1]; 00102 C[2] = _p[p] + _sg[p].R[0] + _sg[p].R[1]; 00103 C[3] = _p[p] - _sg[p].R[0] + _sg[p].R[1]; 00104 } 00105 else for (size_t i=0; i<4; ++i) C[i] = _p[p] + _sg[p].R[i]; 00106 } 00107 00108 // Limit methods 00109 void MinMaxS (int iComp, double & MinS, double & MaxS) const; 00110 void MinMaxE (int iComp, double & MinE, double & MaxE) const; 00111 void MinMaxCamp (double & Minp, double & Maxp) const; 00112 void MinMaxCamq (double & Minq, double & Maxq) const; 00113 void MinMaxVelN (double & MinNv, double & MaxNv) const; 00114 00115 // Constants 00116 double MINMASS; 00117 00118 // Data 00119 Models_t Mdls; 00120 00121 private: 00122 // Data 00123 bool _cpdi; 00124 bool _mpm; 00125 bool _usf; 00126 size_t _npcell; 00127 Grid2D * _g; 00128 ptIsPointInGeom _is_point_in_geom; 00129 ptDensity _density; 00130 ptModelData _model_data; 00131 ptIniVelocity _ini_velocity; 00132 ptHasTraction _has_traction; 00133 ptHasAppDisp _has_app_disp; 00134 ptAppDisp _app_disp; 00135 00136 // Geometry 00137 Array<Vector3D> _p; 00138 Array<int> _ipbry; 00139 00140 // Data (size == _np) 00141 Array<Vector3D> _p0; 00142 Array<double> _V0; 00143 Array<double> _V; 00144 Array<double> _m; 00145 Array<Vector3D> _f; 00146 Array<Vector3D> _u; 00147 Array<Vector3D> _v; 00148 Array<ATensor2> _F; 00149 Array<EquilibState*> _sta; 00150 Array<ShapeAndGrads2D> _sg; 00151 Array<Model*> _mdl; 00152 Array<int> _clr; 00153 00154 // Boundary conditions (size == np) 00155 Array<bool> _has_t; 00156 Array<bool> _has_fu; 00157 Array<Vector3D> _t; 00158 Array<Vector3D> _fu; 00159 00160 // Scratch pad 00161 mutable STensor2 _sp0; 00162 mutable Vector3D _vn; 00163 mutable ATensor2 _Lp; 00164 mutable ATensor2 _DF; 00165 00166 // Update functions 00167 void _initialize_points (); 00168 void _calc_veloc_grad (size_t p); 00169 void _stress_update (double Dt, double & DsE); 00170 bool _explicit_update (double Dt, Vector3D const & B, double LdM, Array<Vector3D> * DqDt, double & DsE); 00171 00172 // Interpolation functions 00173 double _mpm_snp (int n, int p, int icomp) const; 00174 double _mpm_gnp (int n, int p, int icomp) const; 00175 double _gmpm_snp (int n, int p, int icomp) const; 00176 double _gmpm_gnp (int n, int p, int icomp) const; 00177 void _Nnp (int n, int p, double & N) const; 00178 double _S (int n, Vector3D const & X) const; 00179 void _Snp (int n, int p, double & S) const; 00180 void _Gnp (int n, int p, Vector3D & G) const; 00181 }; 00182 00183 00185 00186 00187 /* public */ 00188 00189 inline MPoints2D::MPoints2D(size_t nPcell, Grid2D * G, ptIsPointInGeom pIsPointInGeom, ptDensity pDensity, ptModelData pModelData, ptIniVelocity pIniVelocity, ptHasTraction pHasTraction, ptHasAppDisp pHasAppDisp, ptAppDisp pAppDisp) 00190 : MINMASS (DBL_EPSILON*10.0), 00191 _cpdi (true), 00192 _mpm (false), 00193 _usf (true), 00194 _npcell (nPcell), 00195 _g (G), 00196 _is_point_in_geom (pIsPointInGeom), 00197 _density (pDensity), 00198 _model_data (pModelData), 00199 _ini_velocity (pIniVelocity), 00200 _has_traction (pHasTraction), 00201 _has_app_disp (pHasAppDisp), 00202 _app_disp (pAppDisp) 00203 { 00204 C.Resize(4); 00205 _initialize_points(); 00206 } 00207 00208 inline MPoints2D::~MPoints2D () 00209 { 00210 for (size_t i=0; i<_sta.Size(); ++i) if (_sta[i]!=NULL) delete _sta[i]; 00211 for (Models_t::iterator it=Mdls.begin(); it!=Mdls.end(); ++it) delete it->second; 00212 } 00213 00214 inline bool MPoints2D::TimeUpdate (bool FwEuler, double Dt, ptB pB, ptLdM pM, double & t, double & DsE) 00215 { 00216 // Check if access to grid was set 00217 if (_g==NULL) return false; // failed 00218 00219 // Variables 00220 double dsek; 00221 Vector3D bk; 00222 double Mk; 00223 00224 // Update 00225 if (FwEuler) 00226 { 00227 double nincs = 1; // Number of sub-increments 00228 double dtk = Dt/nincs; // Sub-divide time-step 00229 for (int k=0; k<nincs; ++k) 00230 { 00231 // Increments 00232 (*pB) (t+dtk, bk); // Body forces 00233 (*pM) (t+dtk, Mk); // Multiplier for external forces 00234 00235 // Update points 00236 bool ok = _explicit_update (dtk, bk, Mk, &_g->dqdt(), dsek); 00237 00238 // Update time and strain energy 00239 t += dtk; 00240 DsE += dsek; 00241 00242 // Check 00243 if (!ok) return false; 00244 } 00245 return true; 00246 } 00247 else return false; // not implemented yet 00248 } 00249 00250 inline void MPoints2D::MinMaxS (int iComp, double & MinS, double & MaxS) const 00251 { 00252 if (_sta.Size()<1) { MinS=0; MaxS=0; return; } 00253 MinS = _sta[0]->Sig(iComp); 00254 MaxS = _sta[0]->Sig(iComp); 00255 for (size_t i=1; i<_p.Size(); ++i) 00256 { 00257 if (_sta[i]->Sig(iComp)<MinS) MinS = _sta[i]->Sig(iComp); 00258 if (_sta[i]->Sig(iComp)>MaxS) MaxS = _sta[i]->Sig(iComp); 00259 } 00260 } 00261 00262 inline void MPoints2D::MinMaxE (int iComp, double & MinE, double & MaxE) const 00263 { 00264 if (_sta.Size()<1) { MinE=0; MaxE=0; return; } 00265 MinE = _sta[0]->Eps(iComp); 00266 MaxE = _sta[0]->Eps(iComp); 00267 for (size_t i=1; i<_p.Size(); ++i) 00268 { 00269 if (_sta[i]->Eps(iComp)<MinE) MinE = _sta[i]->Eps(iComp); 00270 if (_sta[i]->Eps(iComp)>MaxE) MaxE = _sta[i]->Eps(iComp); 00271 } 00272 } 00273 00274 inline void MPoints2D::MinMaxCamp (double & Minp, double & Maxp) const 00275 { 00276 if (_sta.Size()<1) { Minp=0; Maxp=0; return; } 00277 Minp = Calc_pcam(_sta[0]->Sig); 00278 Maxp = Minp; 00279 for (size_t i=1; i<_p.Size(); ++i) 00280 { 00281 double p = Calc_pcam(_sta[i]->Sig); 00282 if (p<Minp) Minp = p; 00283 if (p>Maxp) Maxp = p; 00284 } 00285 } 00286 00287 inline void MPoints2D::MinMaxCamq (double & Minq, double & Maxq) const 00288 { 00289 if (_sta.Size()<1) { Minq=0; Maxq=0; return; } 00290 Minq = Calc_qcam(_sta[0]->Sig); 00291 Maxq = Minq; 00292 for (size_t i=1; i<_p.Size(); ++i) 00293 { 00294 double q = Calc_qcam(_sta[i]->Sig); 00295 if (q<Minq) Minq = q; 00296 if (q>Maxq) Maxq = q; 00297 } 00298 } 00299 00300 inline void MPoints2D::MinMaxVelN (double & MinNv, double & MaxNv) const 00301 { 00302 if (_v.Size()<1) { MinNv=0; MaxNv=0; return; } 00303 MinNv = MPM::NormV(_v[0]); 00304 MaxNv = MinNv; 00305 for (size_t i=1; i<_v.Size(); ++i) 00306 { 00307 double nv = MPM::NormV(_v[i]); 00308 if (nv<MinNv) MinNv = nv; 00309 if (nv>MaxNv) MaxNv = nv; 00310 } 00311 } 00312 00313 00314 /* update functions */ 00315 00316 inline void MPoints2D::_initialize_points () 00317 { 00318 // Check input 00319 if (!(_npcell==1 || _npcell==4 || _npcell==9 || _npcell==16)) throw new Fatal("MPoints2D::_initialize_points: Number of points per cell (NPcell=%d) must be 1, 4, 9 or 16.",_npcell); 00320 00321 /* x sqrt(x) 2*sqrt(x) log2(2*sqrt(x)) 1/sqrt(x) 1-1/sqrt(x) 00322 * 1 1 2 1 1 0.0 00323 * 4 2 4 2 0.5 0.5 00324 * 16 4 8 3 0.25 0.75 00325 * 00326 * 1 4 9 16 00327 * 1_ _______ 1_ _______ 1_ ________ 1_ _______ 00328 * | | 0.5_| | _ |_3L/4 | | | | 0.75>|_|_|_|_|<7L/8 00329 * 0_| _ |_L/2 0_|___|___| 0_|--|--|--|_L/2 0_|_|_|_|_| 00330 * | | | | _ |_L/4 |--|--|--| |_|_|_|_| 00331 * -1_|_______| -1_|___|___| 1_|__|__|__| -1_|_|_|_|_|<L/8 */ 00332 00333 // Gauss points 00334 double sqnpc = sqrt(static_cast<double>(_npcell)); 00335 double gv = 1.0-1.0/sqnpc; // Gauss point value 00336 double ksi[ 4] = {-1,1,1,-1}; // Natural coordinates 00337 double eta[ 4] = {-1,-1,1,1}; // Natural coordinates 00338 double kgp[16] = {-gv, gv,gv,-gv, -gv,-gv, -gv/3.,-gv/3.,-gv/3.,-gv/3., gv/3.,gv/3.,gv/3.,gv/3., gv,gv }; // ksi Gauss point 00339 double egp[16] = {-gv,-gv,gv, gv, -gv/3.,gv/3., -gv,-gv/3.,gv/3.,gv, -gv,-gv/3.,gv/3.,gv, -gv/3.,gv/3. }; // eta Gauss point 00340 if (_npcell==9) 00341 { 00342 gv = (1.0-(-1.0))/3.0; // one third of cell length in nat cooords 00343 kgp[0]=-gv; kgp[1]=0.0; kgp[2]=gv; kgp[3]=-gv; kgp[4]=0.0; kgp[5]=gv; kgp[6]=-gv; kgp[7]=0.0; kgp[8]=gv; 00344 egp[0]=-gv; egp[1]=-gv; egp[2]=-gv; egp[3]=0.0; egp[4]=0.0; egp[5]=0.0; egp[6]=gv; egp[7]=gv; egp[8]=gv; 00345 } 00346 00347 // Loop over cells 00348 _p .Resize(0); 00349 _clr .Resize(0); 00350 _t .Resize(0); 00351 _has_t .Resize(0); 00352 _fu .Resize(0); 00353 _has_fu.Resize(0); 00354 for (size_t c=0; c<_g->nCells(); ++c) 00355 { 00356 // For each cell c, generate _npcell mat points 00357 for (size_t j=0; j<_npcell; ++j) 00358 { 00359 // For each node j, sum the contribution of node c 00360 Vector3D p; p = 0.0; 00361 Array<Vector3D> nodes(4); 00362 for (size_t l=0; l<4; ++l) 00363 { 00364 nodes[l] = _g->X(_g->C(c)(l)); 00365 double shape = (1.0+ksi[l]*kgp[j])*(1.0+eta[l]*egp[j])/4.0; 00366 p(0) += shape*nodes[l](0); 00367 p(1) += shape*nodes[l](1); 00368 } 00369 00370 // Check if the mat point is inside the geometry 00371 int tag = -1; 00372 if ((*_is_point_in_geom)(p, tag)) 00373 { 00374 // Coords 00375 _p .Push (p); // coordinates 00376 _clr.Push (tag); // clr index 00377 00378 // Boundary conds - tractions 00379 _t .Push(0.0); 00380 _has_t.Push((*_has_traction) (_p[_p.Size()-1], nodes, _g->L(), _npcell, _t[_p.Size()-1])); 00381 00382 // Boundary conds - displacements 00383 _fu .Push(0.0); 00384 _has_fu.Push((*_has_app_disp) (_p[_p.Size()-1], nodes, _g->L(), _npcell, _fu[_p.Size()-1])); 00385 } 00386 } 00387 } 00388 00389 // Delete previous EquilibState 00390 for (size_t i=0; i<_sta.Size(); ++i) if (_sta[i]!=NULL) delete _sta[i]; 00391 00392 // Delete previous Models 00393 for (Models_t::iterator it=Mdls.begin(); it!=Mdls.end(); ++it) delete it->second; 00394 Mdls.clear(); 00395 00396 // Data for models 00397 String name, scheme; 00398 SDPair prms, inis; 00399 00400 // Initialize arrays 00401 _p0 .Resize(_p.Size()); 00402 _V0 .Resize(_p.Size()); 00403 _V .Resize(_p.Size()); 00404 _m .Resize(_p.Size()); 00405 _f .Resize(_p.Size()); 00406 _u .Resize(_p.Size()); 00407 _v .Resize(_p.Size()); 00408 _F .Resize(_p.Size()); 00409 _sta.Resize(_p.Size()); 00410 _sg .Resize(_p.Size()); 00411 _mdl.Resize(_p.Size()); 00412 double mt = 0.0; // total mass 00413 for (size_t p=0; p<_p.Size(); ++p) 00414 { 00415 double lx = 0.5*_g->L(0)/sqnpc; // half-length 00416 double ly = 0.5*_g->L(1)/sqnpc; // half-length 00417 _p0 [p] = _p[p]; 00418 _V0 [p] = 4.0*lx*ly; 00419 _V [p] = _V0[p]; 00420 _m [p] = _V[p] * (*_density) (_p[p]); 00421 _f [p] = 0.0, 0.0, 0.0; 00422 _u [p] = 0.0, 0.0, 0.0; 00423 _v [p] = 0.0, 0.0, 0.0; 00424 (*_ini_velocity) (_p[p], _v[p]); 00425 _F [p] = 1.0,1.0,1.0, 0.0,0.0,0.0, 0.0,0.0,0.0; 00426 if (NRVEC==2) 00427 { 00428 _sg[p].R0[0] = lx, 0., 0.; 00429 _sg[p].R0[1] = 0., ly, 0.; 00430 } 00431 else 00432 { 00433 _sg[p].R0[0] = -lx, -ly, 0.0; 00434 _sg[p].R0[1] = lx, -ly, 0.0; 00435 _sg[p].R0[2] = lx, ly, 0.0; 00436 _sg[p].R0[3] = -lx, ly, 0.0; 00437 } 00438 for (size_t k=0; k<NRVEC; ++k) _sg[p].R[k] = _sg[p].R0[k]; 00439 mt += _m[p]; 00440 00441 // allocate models and set initial state 00442 int tag = -1; 00443 (*_model_data) (_p[p], name, tag, NDIM, prms, inis, scheme); 00444 Models_t::iterator it = Mdls.find(tag); 00445 if (it==Mdls.end()) 00446 { 00447 Mdls[tag] = AllocModel (name, NDIM, prms, NULL); 00448 it = Mdls.find(tag); 00449 if (scheme!="") it->second->SUp.SetScheme (scheme); 00450 } 00451 _mdl[p] = it->second; 00452 _sta[p] = new EquilibState (NDIM); 00453 _mdl[p]->InitIvs (inis, _sta[p]); 00454 } 00455 std::cout << "MPoints2D::_initialize_points: total mass = " << mt << std::endl; 00456 } 00457 00458 inline void MPoints2D::_calc_veloc_grad (size_t p) 00459 { 00460 _Lp = 0.,0.,0., 0.,0.,0., 0.,0.,0.; 00461 for (n2i_it it=_sg[p].n2i.begin(); it!=_sg[p].n2i.end(); ++it) 00462 { 00463 int n = it->first; // grid node number 00464 int i = it->second; // index in N,S,G 00465 _vn = 0., 0., 0.; // node velocity 00466 if (_g->m(n)>MINMASS) _vn = _g->q(n)/_g->m(n); // node velocity 00467 if (_g->IsFixed(n,0)) _vn(0) = 0.0; // Fix nodes: x-direction 00468 if (_g->IsFixed(n,1)) _vn(1) = 0.0; // Fix nodes: y-direction 00469 DyadUp (_vn, _sg[p].G[i], _Lp); // Lp += vp dyad G 00470 } 00471 } 00472 00473 inline void MPoints2D::_stress_update (double Dt, double & DsE) 00474 { 00475 DsE = 0.0; 00476 for (size_t p=0; p<_p.Size(); ++p) 00477 { 00478 // Compute velocity gradient: _Lp 00479 _calc_veloc_grad (p); 00480 00481 // Update F and stress state 00482 DsE += _mdl[p]->SUp.Update (Dt, _Lp, _F[p], _sta[p]) * _V[p]; 00483 00484 // Update volume 00485 _V[p] = Det(_F[p]) * _V0[p]; 00486 00487 // Update centre-to-corner vectors 00488 for (size_t i=0; i<NRVEC; ++i) Dot (_F[p], _sg[p].R0[i], _sg[p].R[i]); 00489 } 00490 } 00491 00492 inline bool MPoints2D::_explicit_update (double Dt, Vector3D const & B, double LdM, Array<Vector3D> * DqDt, double & DsE) 00493 { 00494 // 1) Discard previous grid 00495 _g->ClearState(); 00496 00497 // 2) Compute interpolation values 00498 if (_cpdi) 00499 { 00500 for (size_t p=0; p<_p.Size(); ++p) 00501 { 00502 CalcC(p); // calculate corners 00503 _sg[p].n2i.clear(); 00504 int k = 0.0; 00505 for (size_t i=0; i<4; ++i) // for each corner i 00506 { 00507 int l = static_cast<int>((C[i](0)-_g->xMin())/_g->L(0)); // left-grid-index 00508 int b = static_cast<int>((C[i](1)-_g->yMin())/_g->L(1)); // bottom-grid-index 00509 int r = l+1; // right-grid-index 00510 int t = b+1; // top-grid-index 00511 int nlb = l+b*_g->nCol(); 00512 int nrb = r+b*_g->nCol(); 00513 int nlt = l+t*_g->nCol(); 00514 int nrt = r+t*_g->nCol(); 00515 _sg[p].n2i[nlb]=k; _Nnp(nlb,p,_sg[p].N[k]); _Snp(nlb,p,_sg[p].S[k]); _Gnp(nlb,p,_sg[p].G[k]); k++; 00516 _sg[p].n2i[nrb]=k; _Nnp(nrb,p,_sg[p].N[k]); _Snp(nrb,p,_sg[p].S[k]); _Gnp(nrb,p,_sg[p].G[k]); k++; 00517 _sg[p].n2i[nlt]=k; _Nnp(nlt,p,_sg[p].N[k]); _Snp(nlt,p,_sg[p].S[k]); _Gnp(nlt,p,_sg[p].G[k]); k++; 00518 _sg[p].n2i[nrt]=k; _Nnp(nrt,p,_sg[p].N[k]); _Snp(nrt,p,_sg[p].S[k]); _Gnp(nrt,p,_sg[p].G[k]); k++; 00519 } 00520 } 00521 } 00522 else 00523 { 00524 int nctb = (_mpm ? 2 : 4); // number of contributions 00525 for (size_t p=0; p<_p.Size(); ++p) 00526 { 00527 _sg[p].n2i.clear(); 00528 int k = 0; 00529 int lbn = _g->LBN(_p[p], _mpm); // Compute reference node 00530 for (int i=0; i<nctb; ++i) 00531 for (int j=0; j<nctb; ++j) 00532 { 00533 int n = lbn+i+j*_g->nCol(); // grid node number 00534 _sg[p].n2i[n] = k; // set map: node => index in S,N,G 00535 _Nnp (n,p, _sg[p].N[k]); // calc interpolation function (C0) 00536 _Snp (n,p, _sg[p].S[k]); // calc interpolation function (C1) 00537 _Gnp (n,p, _sg[p].G[k]); // calc deriv. interp. function 00538 k++; 00539 } 00540 } 00541 } 00542 00543 // 3) Initialize grid state (mass and momentum) 00544 for (size_t p=0; p<_p.Size(); ++p) 00545 { 00546 for (n2i_it it=_sg[p].n2i.begin(); it!=_sg[p].n2i.end(); ++it) 00547 { 00548 int n = it->first; // grid node number 00549 int i = it->second; // index in N,S,G 00550 _g->m(n) += _sg[p].S[i]*_m[p]; // node mass 00551 _g->q(n) += _sg[p].S[i]*_m[p]*_v[p]; // node momentum 00552 } 00553 } 00554 00555 // 4) (USF) Update strain and stress 00556 if (_usf) _stress_update (Dt, DsE); 00557 00558 // 5) Compute internal and external forces 00559 for (size_t p=0; p<_p.Size(); ++p) 00560 { 00561 for (n2i_it it=_sg[p].n2i.begin(); it!=_sg[p].n2i.end(); ++it) 00562 { 00563 int n = it->first; // grid node number 00564 int i = it->second; // index in N,S,G 00565 ScDotUp (_V[p],_sta[p]->Sig,_sg[p].G[i], _g->fi(n)); // fi += V_p*Sig.G 00566 _g->fe(n) += (_sg[p].S[i]*_m[p])*B; // fe += mp*b*S 00567 if (_has_t[p]) _g->fe(n) += LdM*_sg[p].N[i]*_t[p]; // fe += int(S*t)dA 00568 } 00569 } 00570 00571 // 6) Compute rate of momentum and momentum 00572 for (size_t n=0; n<_g->nNodes(); ++n) 00573 { 00574 (*DqDt)[n] = _g->fe(n) - _g->fi(n); // Nodes rate of momentum 00575 _g->q(n) += (*DqDt)[n]*Dt; // Update nodes momentum 00576 if (_g->IsFixed(n,0)) { (*DqDt)[n](0)=0.0; _g->q(n)(0)=0.0; } // Fix nodes: x-direction 00577 if (_g->IsFixed(n,1)) { (*DqDt)[n](1)=0.0; _g->q(n)(1)=0.0; } // Fix nodes: y-direction 00578 } 00579 00580 // 7) Update material points 00581 for (size_t p=0; p<_p.Size(); ++p) 00582 { 00583 // Update position and velocity 00584 for (n2i_it it=_sg[p].n2i.begin(); it!=_sg[p].n2i.end(); ++it) 00585 { 00586 int n = it->first; // grid node number 00587 int i = it->second; // index in N,S,G 00588 if (_g->m(n)>MINMASS) 00589 { 00590 _p[p] += (Dt*_sg[p].S[i]/_g->m(n))*_g->q(n); // point position 00591 _v[p] += (Dt*_sg[p].S[i]/_g->m(n))*(*DqDt)[n]; // point velocity 00592 _f[p] += _sg[p].S[i]*_g->fe(n); // point force 00593 } 00594 } 00595 00596 // Update displacements 00597 _u[p] = _p[p] - _p0[p]; 00598 00599 // Check 00600 if (!_g->IsInside(_p[p], /*WithBorder*/true)) 00601 { 00602 std::cout << " [1;31mError:[0m Material point moved outside the valid region inside the grid (with borders)" << std::endl; 00603 return false; 00604 } 00605 } 00606 00607 // 8) (USL) Update strain and stress 00608 if (!_usf) _stress_update (Dt, DsE); 00609 00610 // end 00611 return true; 00612 } 00613 00614 00615 /* interpolation functions */ 00616 00617 inline double MPoints2D::_mpm_snp (int n, int p, int icomp) const 00618 { 00619 double d = _p[p](icomp)-_g->x(n)(icomp); 00620 double L = _g->L(icomp); 00621 return (d<=-L ? 0.0 : 00622 (d> L ? 0.0 : 00623 (d<=0.0 ? 1.0+d/L : 00624 1.0-d/L ))); 00625 } 00626 00627 inline double MPoints2D::_mpm_gnp (int n, int p, int icomp) const 00628 { 00629 double d = _p[p](icomp)-_g->x(n)(icomp); 00630 double L = _g->L(icomp); 00631 return (d<=-L ? 0.0 : 00632 (d> L ? 0.0 : 00633 (d<= 0.0 ? 1.0/L : 00634 -1.0/L ))); 00635 } 00636 00637 inline double MPoints2D::_gmpm_snp (int n, int p, int icomp) const 00638 { 00639 double d = _p[p](icomp)-_g->x(n)(icomp); 00640 double L = _g->L(icomp); 00641 double lp = fabs(_sg[p].R0[0](icomp)); 00642 return (d<=-L-lp ? 0.0 : 00643 (d> L+lp ? 0.0 : 00644 (d<=-L+lp ? pow(L+lp+d,2.0)/(4.0*L*lp) : 00645 (d<= -lp ? 1.0+d/L : 00646 (d<= lp ? 1.0-(d*d+lp*lp)/(2.0*L*lp) : 00647 (d<= L-lp ? 1.0-d/L : 00648 pow(L+lp-d,2.0)/(4.0*L*lp) )))))); 00649 } 00650 00651 inline double MPoints2D::_gmpm_gnp (int n, int p, int icomp) const 00652 { 00653 double d = _p[p](icomp)-_g->x(n)(icomp); 00654 double L = _g->L(icomp); 00655 double lp = fabs(_sg[p].R0[0](icomp)); 00656 return (d<=-L-lp ? 0.0 : 00657 (d> L+lp ? 0.0 : 00658 (d<=-L+lp ? (L+d+lp)/(2.0*L*lp) : 00659 (d<= -lp ? 1.0/L : 00660 (d<= lp ? -d/(L*lp) : 00661 (d<= L-lp ? -1.0/L : 00662 (-L+d-lp)/(2.0*L*lp) )))))); 00663 } 00664 00665 inline void MPoints2D::_Nnp (int n, int p, double & N) const 00666 { 00667 N = _mpm_snp(n,p,0)*_mpm_snp(n,p,1); 00668 } 00669 00670 inline double MPoints2D::_S (int n, Vector3D const & X) const 00671 { 00672 double dx = X(0)-_g->x(n)(0); 00673 double dy = X(1)-_g->x(n)(1); 00674 double sx = ( dx<-_g->L(0) || dx>_g->L(0) ? 0.0 : (dx<0.0 ? 1.0+dx/_g->L(0) : 1.0-dx/_g->L(0)) ); 00675 double sy = ( dy<-_g->L(1) || dy>_g->L(1) ? 0.0 : (dy<0.0 ? 1.0+dy/_g->L(1) : 1.0-dy/_g->L(1)) ); 00676 return sx*sy; 00677 } 00678 00679 inline void MPoints2D::_Snp (int n, int p, double & S) const 00680 { 00681 if (_cpdi) 00682 { 00683 // remember to calculate corners first !!! 00684 S = (_S(n,C[0])+_S(n,C[1])+_S(n,C[2])+_S(n,C[3]))/4.0; 00685 } 00686 else 00687 { 00688 if (_mpm) S = _mpm_snp (n,p,0)*_mpm_snp (n,p,1); 00689 else S = _gmpm_snp(n,p,0)*_gmpm_snp(n,p,1); 00690 } 00691 } 00692 00693 inline void MPoints2D::_Gnp (int n, int p, Vector3D & G) const 00694 { 00695 if (_cpdi) 00696 { 00697 // remember to calculate corners first !!! 00698 double a = _S(n,C[0]) - _S(n,C[2]); 00699 double b = _S(n,C[1]) - _S(n,C[3]); 00700 G = (a*(_sg[p].R[0](1)-_sg[p].R[1](1))+b*(_sg[p].R[0](1)+_sg[p].R[1](1)))/_V[p], 00701 (a*(_sg[p].R[1](0)-_sg[p].R[0](0))-b*(_sg[p].R[1](0)+_sg[p].R[0](0)))/_V[p], 0.0; 00702 } 00703 else 00704 { 00705 if (_mpm) G = _mpm_snp (n,p,1)*_mpm_gnp (n,p,0), _mpm_snp (n,p,0)*_mpm_gnp (n,p,1), 0.0; 00706 else G = _gmpm_snp(n,p,1)*_gmpm_gnp(n,p,0), _gmpm_snp(n,p,0)*_gmpm_gnp(n,p,1), 0.0; 00707 } 00708 } 00709 00710 }; // namespace MPM 00711 00712 #endif // MPM_MPOINTS2D_H
1.7.6.1
