 |
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_NODE_H 00021 #define MECHSYS_FEM_NODE_H 00022 00023 // Std Lib 00024 #include <sstream> // for istringstream, ostringstream 00025 #include <map> 00026 00027 // MechSys 00028 #include <mechsys/mesh/mesh.h> 00029 #include <mechsys/util/maps.h> 00030 #include <mechsys/util/fatal.h> 00031 #include <mechsys/util/numstreams.h> 00032 #include <mechsys/linalg/matvec.h> 00033 #include <mechsys/linalg/sparse_triplet.h> 00034 00035 namespace FEM 00036 { 00037 00039 class BCFuncs 00040 { 00041 public: 00042 virtual ~BCFuncs () {} 00043 virtual double u (double Time) { return 0; } 00044 virtual double v (double Time) { return 0; } 00045 virtual double a (double Time) { return 0; } 00046 virtual double fm (double Time) { return 0; } 00047 }; 00048 00049 class Node 00050 { 00051 public: 00052 // Constructor 00053 Node (Mesh::Vertex const & TheVert) : Vert(TheVert), NShares(0), _has_incsup(false) {} 00054 00055 // Methods 00056 void AddDOF (char const * StrU, char const * StrF); 00057 size_t NDOF () const { return _eq.Size(); } 00058 double U (char const * UKey ) const { return _U(UKey); } 00059 double F (char const * FKey ) const { return _F(FKey); } 00060 double & U (String const & UKey ) { return _U(UKey); } 00061 double & F (String const & FKey ) { return _F(FKey); } 00062 double & U (char const * UKey ) { return _U(UKey); } 00063 double & F (char const * FKey ) { return _F(FKey); } 00064 double & V (size_t IdxDOF) { return _V(UKey(IdxDOF)); } 00065 double & U (size_t IdxDOF) { return _U(UKey(IdxDOF)); } 00066 double & F (size_t IdxDOF) { return _F(FKey(IdxDOF)); } 00067 double UOrZero (String const & UKey ) const { return _U.ValOrZero(UKey); } 00068 double VOrZero (String const & UKey ) const { return _V.ValOrZero(UKey); } 00069 double FOrZero (String const & FKey ) const { return _F.ValOrZero(FKey); } 00070 void TrySetU (String const & UKey, double Val) { if (_U.HasKey(UKey)) _U(UKey)=Val; } 00071 void TrySetF (String const & FKey, double Val) { if (_F.HasKey(FKey)) _F(FKey)=Val; } 00072 int Eq (char const * UKey ) const { return _eq[_U2IDOF(UKey)]; } 00073 int Eq (String const & UKey ) const { return _eq[_U2IDOF(UKey)]; } 00074 int Eq (size_t IdxDOF) const { return _eq[IdxDOF]; } 00075 int & Eq (size_t IdxDOF) { return _eq[IdxDOF]; } 00076 String const & UKey (size_t IdxDOF) const { return _U.Keys[IdxDOF]; } 00077 String const & FKey (size_t IdxDOF) const { return _F.Keys[IdxDOF]; } 00078 bool pU (size_t IdxDOF) const { return _pU[IdxDOF]; } 00079 void SetUF (Vec_t const & UVec, Vec_t const & FVec); 00080 void SetUVF (Vec_t const & UVec, Vec_t const & VVec, Vec_t const & FVec); 00081 void GetUF (Vec_t & UVec, Vec_t & FVec); 00082 void Clear (); 00083 void ClearU () { _U.SetValues (0.0); } 00084 00085 // Methods to set/access prescribed U values 00086 void SetPU (String const & UKey, double Val, BCFuncs * BCF=NULL); 00087 void SetMPU (size_t IdxPU, BCFuncs * BCF=NULL){ _MPU[IdxPU] = BCF; } 00088 size_t NPU () const { return _PU.Keys.Size(); } 00089 int EqPU (size_t IdxPU) const { return _eq[_U2IDOF(_PU.Keys[IdxPU])]; } 00090 void DelPUs () { _PU.clear(); _MPU.Resize(0); } 00091 String const & PUKey (size_t IdxPU) const { return _PU.Keys[IdxPU]; } 00092 double PU (size_t IdxPU, double Time) const; 00093 double PV (size_t IdxPU, double Time) const; 00094 double PA (size_t IdxPU, double Time) const; 00095 double PUIdxDOF (size_t IdxDOF, double Time) const { return PU(_PU.IdxKey(UKey(IdxDOF)),Time); } 00096 double PVIdxDOF (size_t IdxDOF, double Time) const { return PV(_PU.IdxKey(UKey(IdxDOF)),Time); } 00097 double PAIdxDOF (size_t IdxDOF, double Time) const { return PA(_PU.IdxKey(UKey(IdxDOF)),Time); } 00098 00099 // Methods to set/access prescribed F values 00100 void AddToPF (String const & FKey, double Val, BCFuncs * BCF=NULL); 00101 void SetMPF (size_t IdxPF, BCFuncs * BCF=NULL){ _MPF[IdxPF] = BCF; } 00102 size_t NPF () const { return _PF.Keys.Size(); } 00103 int EqPF (size_t IdxPF) const { return _eq[_F2IDOF(_PF.Keys[IdxPF])]; } 00104 void DelPFs () { _PF.clear(); _MPF.Resize(0); } 00105 String const & PFKey (size_t IdxPF) const { return _PF.Keys[IdxPF]; } 00106 double PF (size_t IdxPF, double Time) const; 00107 double PFOrZero (size_t IdxDOF, double Time) const { int idxpf=_PF.IdxKey(FKey(IdxDOF)); return (idxpf<0 ? 0.0 : PF(idxpf,Time)); } 00108 void AccumPF (); 00109 void Reactions (std::map<String,double> & R) const; 00110 00111 // Pins 00112 void SetPin (Array<Node*> const & ConnectedNodes) { _pin_nodes=ConnectedNodes; } 00113 void SetLagPin (int & EqLag, Sparse::Triplet<double,int> & A) const; 00114 void ClrRPin (int & EqLag, Vec_t & R) const; 00115 00116 // Inclined supports: 2D 00117 void SetIncSup (double Alpha); 00118 void DelIncSup () { _has_incsup=false; } 00119 bool HasIncSup () const { return _has_incsup; } 00120 void SetLagIncSup (int & EqLag, Sparse::Triplet<double,int> & A) const; 00121 void ClrRIncSup (int & EqLag, Vec_t & R) const; 00122 00123 // Data 00124 Mesh::Vertex const & Vert; 00125 long NShares; 00126 00127 private: 00128 // Data at every node 00129 SDPair _V; 00130 SDPair _U; 00131 SDPair _F; 00132 SIPair _U2IDOF; 00133 SIPair _F2IDOF; 00134 Array<int> _eq; 00135 Array<bool> _pU; 00136 00137 // Data at nodes with prescribed values 00138 SDPair _PU; 00139 SDPair _PF; 00140 SDPair _aPF; 00141 Array<BCFuncs*> _MPU; 00142 Array<BCFuncs*> _MPF; 00143 00144 // Data for pins 00145 Array<Node*> _pin_nodes; 00146 00147 // Data at nodes with inclined supports 00148 double _incsup_alpha; 00149 bool _has_incsup; 00150 00151 friend std::ostream & operator<< (std::ostream & os, Node const & N); 00152 }; 00153 00154 00156 00157 00158 inline void Node::AddDOF (char const * StrU, char const * StrF) 00159 { 00160 std::istringstream u_iss(StrU); 00161 std::istringstream f_iss(StrF); 00162 String u_key, f_key; 00163 while ((u_iss>>u_key) && (f_iss>>f_key)) 00164 { 00165 if (!_U.HasKey(u_key)) // add only if not found 00166 { 00167 _V.Set (u_key.CStr(), 0.0); 00168 _U.Set (u_key.CStr(), 0.0); 00169 _F.Set (f_key.CStr(), 0.0); 00170 _U2IDOF.Set (u_key.CStr(), _eq.Size()); 00171 _F2IDOF.Set (f_key.CStr(), _eq.Size()); 00172 _eq.Push (-1); 00173 _pU.Push (false); 00174 } 00175 else 00176 { 00177 if (!_F.HasKey(f_key)) throw new Fatal("FEM::Node: __internal_error__ Node has UKey=%s but not FKey=%s",u_key.CStr(),f_key.CStr()); 00178 } 00179 } 00180 } 00181 00182 inline void Node::Clear () 00183 { 00184 _U .SetValues (0.0); 00185 _V .SetValues (0.0); 00186 _F .SetValues (0.0); 00187 _eq.SetValues (-1); 00188 _pU.SetValues (false); 00189 DelPUs (); 00190 DelPFs (); 00191 _has_incsup = false; 00192 } 00193 00194 inline void Node::SetUF (Vec_t const & UVec, Vec_t const & FVec) 00195 { 00196 for (size_t idof=0; idof<NDOF(); ++idof) 00197 { 00198 int eq = Eq(idof); 00199 _U(UKey(idof)) = UVec(eq); 00200 _F(FKey(idof)) = FVec(eq); 00201 } 00202 } 00203 00204 inline void Node::SetUVF (Vec_t const & UVec, Vec_t const & VVec, Vec_t const & FVec) 00205 { 00206 for (size_t idof=0; idof<NDOF(); ++idof) 00207 { 00208 int eq = Eq(idof); 00209 _U(UKey(idof)) = UVec(eq); 00210 _V(UKey(idof)) = VVec(eq); 00211 _F(FKey(idof)) = FVec(eq); 00212 } 00213 } 00214 00215 inline void Node::GetUF (Vec_t & UVec, Vec_t & FVec) 00216 { 00217 for (size_t idof=0; idof<NDOF(); ++idof) 00218 { 00219 int eq = Eq(idof); 00220 UVec(eq) = _U(UKey(idof)); 00221 FVec(eq) = _F(FKey(idof)); 00222 } 00223 } 00224 00225 inline void Node::SetPU (String const & UKey, double Val, BCFuncs * BCF) 00226 { 00227 size_t idx_pu = _PU.ReSet (UKey.CStr(), Val); 00228 if (idx_pu==_MPU.Size()) _MPU.Push (BCF); 00229 else _MPU[idx_pu] = BCF; 00230 _has_incsup = false; 00231 _pU[_U2IDOF(UKey)] = true; 00232 } 00233 00234 inline void Node::AddToPF (String const & FKey, double Val, BCFuncs * BCF) 00235 { 00236 size_t idx_pf = _PF.AddVal (FKey.CStr(), Val); 00237 if (idx_pf==_MPF.Size()) _MPF.Push (BCF); 00238 else _MPF[idx_pf] = BCF; 00239 } 00240 00241 inline double Node::PU (size_t IdxPU, double Time) const 00242 { 00243 if (_MPU[IdxPU]==NULL) return _PU(_PU.Keys[IdxPU]); 00244 else return _MPU[IdxPU]->u(Time); 00245 } 00246 00247 inline double Node::PV (size_t IdxPU, double Time) const 00248 { 00249 if (_MPU[IdxPU]==NULL) return 0.0; 00250 else return _MPU[IdxPU]->v(Time); 00251 } 00252 00253 inline double Node::PA (size_t IdxPU, double Time) const 00254 { 00255 if (_MPU[IdxPU]==NULL) return 0.0; 00256 else return _MPU[IdxPU]->a(Time); 00257 } 00258 00259 inline double Node::PF (size_t IdxPF, double Time) const 00260 { 00261 if (_MPF[IdxPF]==NULL) return _PF(_PF.Keys[IdxPF]); 00262 else return _PF(_PF.Keys[IdxPF]) * _MPF[IdxPF]->fm(Time); 00263 } 00264 00265 inline void Node::AccumPF () 00266 { 00267 for (StrDbl_t::iterator it=_PF.begin(); it!=_PF.end(); ++it) 00268 _aPF.AddVal (it->first.CStr(), it->second); 00269 } 00270 00271 inline void Node::Reactions (std::map<String,double> & R) const 00272 { 00273 for (size_t i=0; i<_PU.Keys.Size(); ++i) 00274 { 00275 String const & ukey = _PU.Keys[i]; 00276 int idof = _U2IDOF(ukey); 00277 String const & fkey = _F.Keys[idof]; 00278 double reac = _F(fkey); 00279 if (_aPF.HasKey(fkey)) reac -= _aPF(fkey); 00280 R[ukey] = reac; 00281 } 00282 } 00283 00284 inline void Node::SetLagPin (int & EqLag, Sparse::Triplet<double,int> & A) const 00285 { 00286 for (size_t i=0; i<_pin_nodes.Size(); ++i) 00287 { 00288 for (size_t j=0; j<_U.Keys.Size(); ++j) 00289 { 00290 if (_U.Keys[j]=="ux" || _U.Keys[j]=="uy" || _U.Keys[j]=="uz") 00291 { 00292 int eq0 = Eq(_U.Keys[j]); 00293 int eq1 = _pin_nodes[i]->Eq(_U.Keys[j]); 00294 A.PushEntry (eq0,EqLag,1.0); A.PushEntry (eq1,EqLag,-1.0); 00295 A.PushEntry (EqLag,eq0,1.0); A.PushEntry (EqLag,eq1,-1.0); 00296 EqLag++; // each pin adds one equation per DOF 00297 } 00298 } 00299 } 00300 } 00301 00302 inline void Node::ClrRPin (int & EqLag, Vec_t & R) const 00303 { 00304 for (size_t i=0; i<_pin_nodes.Size(); ++i) 00305 { 00306 for (size_t j=0; j<_U.Keys.Size(); ++j) 00307 { 00308 if (_U.Keys[j]=="ux" || _U.Keys[j]=="uy" || _U.Keys[j]=="uz") 00309 { 00310 int eq0 = Eq(_U.Keys[j]); 00311 int eq1 = _pin_nodes[i]->Eq(_U.Keys[j]); 00312 R(eq0) = 0.0; 00313 R(eq1) = 0.0; 00314 EqLag++; // each pin adds one equation per DOF 00315 } 00316 } 00317 } 00318 } 00319 00320 inline void Node::SetIncSup (double Alpha) 00321 { 00322 if (NPU()>0) return; // skip if this node has prescribed U already 00323 _incsup_alpha = Alpha*Util::PI/180.0; 00324 _has_incsup = true; 00325 } 00326 00327 inline void Node::SetLagIncSup (int & EqLag, Sparse::Triplet<double,int> & A) const 00328 { 00329 double s = sin(_incsup_alpha); 00330 double c = cos(_incsup_alpha); 00331 int eq0 = Eq("ux"); 00332 int eq1 = Eq("uy"); 00333 A.PushEntry (EqLag, eq0, s); 00334 A.PushEntry (EqLag, eq1, -c); 00335 A.PushEntry (eq0, EqLag, s); 00336 A.PushEntry (eq1, EqLag, -c); 00337 EqLag++; // each inclined support adds one equation 00338 } 00339 00340 inline void Node::ClrRIncSup (int & EqLag, Vec_t & R) const 00341 { 00342 int eq0 = Eq("ux"); 00343 int eq1 = Eq("uy"); 00344 R(eq0) = 0.0; 00345 R(eq1) = 0.0; 00346 EqLag++; // each inclined support adds one equation 00347 } 00348 00349 std::ostream & operator<< (std::ostream & os, Node const & N) 00350 { 00351 os << Util::_4 << N.Vert.ID << " " << Util::_4 << N.Vert.Tag << " "; 00352 for (size_t i=0; i<N.NDOF(); ++i) 00353 { 00354 os << N.UKey(i) << " "; 00355 os << N.FKey(i) << " "; 00356 if (i!=N.NDOF()-1) os << " "; 00357 } 00358 os << " NShares=" << (N.NShares>0?TERM_GREEN:TERM_RED) << N.NShares << TERM_RST; 00359 os << " EQ=["; 00360 for (size_t i=0; i<N.NDOF(); ++i) 00361 { 00362 os << N.Eq(i); 00363 if (i!=N.NDOF()-1) os << ","; 00364 else os << "]"; 00365 } 00366 os << " ("; 00367 for (int j=0; j<3; ++j) 00368 { 00369 os << Util::_6_3 << N.Vert.C[j]; 00370 if (j==2) os << ")"; 00371 else os << ", "; 00372 } 00373 if (N._pin_nodes.Size()>0) 00374 { 00375 os << TERM_GREEN <<" PIN:["; 00376 for (size_t i=0; i<N._pin_nodes.Size(); ++i) os << N._pin_nodes[i]->Vert.ID << (i==N._pin_nodes.Size()-1?"]":","); 00377 os << TERM_RST; 00378 } 00379 if (N._has_incsup) os << " INCSUP"; 00380 os << TERM_CLR4 << Util::_reset << " PU=["; 00381 for (size_t i=0; i<N.NPU(); ++i) 00382 { 00383 os << N._PU.Keys[i] << "=" << N._PU(N._PU.Keys[i]); 00384 if (N._MPU[i]!=NULL) os << "*"; 00385 if (i!=N.NPU()-1) os << " "; 00386 } 00387 os << "]" << TERM_CLR5 << " PF=["; 00388 for (size_t i=0; i<N.NPF(); ++i) 00389 { 00390 os << N._PF.Keys[i] << "=" << N._PF(N._PF.Keys[i]); 00391 if (N._MPF[i]!=NULL) os << "*"; 00392 if (i!=N.NPF()-1) os << " "; 00393 } 00394 os << "]" << TERM_RST; 00395 return os; 00396 } 00397 00398 }; //namespace FEM 00399 00400 #endif // MECHSYS_FEM_NODE_H
1.7.6.1
