1 //#************************************************************** 2 //# 3 //# filename: ANCFBeam3D.h 4 //# 5 //# author: Gerstmayr Johannes 6 //# 7 //# generated: 17.October 2004 8 //# description: 3D Element Library 9 //# remarks: 10 //# 11 //# Copyright (c) 2003-2013 Johannes Gerstmayr, Linz Center of Mechatronics GmbH, Austrian 12 //# Center of Competence in Mechatronics GmbH, Institute of Technical Mechanics at the 13 //# Johannes Kepler Universitaet Linz, Austria. All rights reserved. 14 //# 15 //# This file is part of HotInt. 16 //# HotInt is free software: you can redistribute it and/or modify it under the terms of 17 //# the HOTINT license. See folder 'licenses' for more details. 18 //# 19 //# bug reports are welcome!!! 20 //# WWW: www.hotint.org 21 //# email: bug_reports@hotint.org or support@hotint.org 22 //#*************************************************************************************** 23 24 #ifndef ANCFBEAM3D__H 25 #define ANCFBEAM3D__H 26 27 #include "body3d.h" 28 29 30 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 31 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 32 // ANCFBEAM3D ANCFBEAM3D ANCFBEAM3D ANCFBEAM3D ANCFBEAM3D ANCFBEAM3D ANCFBEAM3D 33 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 34 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 35 36 37 const int ANCFBeamMaxIP=45;//45; 38 //rigid cube 39 class ANCFBeam3D: public Body3D 40 { 41 public: 42 //Body3D():Element() {mbs = NULL;}; ANCFBeam3D(MBS * mbsi)43 ANCFBeam3D(MBS* mbsi):Body3D(mbsi), massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 44 T1(), T2(), T(), xg(), xgd(), e0() {}; ANCFBeam3D(const ANCFBeam3D & e)45 ANCFBeam3D(const ANCFBeam3D& e):Body3D(e.mbs),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 46 T1(), T2(), T(), xg(), xgd(), e0() {CopyFrom(e);}; 47 //nodal coordinates of first and second point (x1,x2, x1.Length()==12) 48 ANCFBeam3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, double rhoi, double Emi, double nui, 49 const Vector3D& si, const Vector3D& coli, int beammodel = 0); 50 51 //Element shares nodes with other elements, n1 and n2 are nodenumbers; element sets initial conditions for nodes 52 ANCFBeam3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, int n1i, int n2i, double rhoi, double Emi, double nui, 53 const Vector3D& si, const Vector3D& coli, int beammodel = 0); 54 55 //Element shares nodes with other elements, n1 and n2 are nodenumbers; element sets initial conditions for nodes 56 ANCFBeam3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, const Vector& vc1, const Vector& vc2, int n1i, int n2i, double rhoi, double Emi, double nui, 57 const Vector3D& si, const Vector3D& coli, int beammodel = 0); 58 GetCopy()59 virtual Element* GetCopy() 60 { 61 Element* ec = new ANCFBeam3D(*this); 62 return ec; 63 } CopyFrom(const Element & e)64 virtual void CopyFrom(const Element& e) 65 { 66 Body3D::CopyFrom(e); 67 const ANCFBeam3D& ce = (const ANCFBeam3D&)e; 68 lx = ce.lx; 69 ly = ce.ly; 70 lz = ce.lz; 71 Em = ce.Em; 72 nu = ce.nu; 73 //size = ce.size; 74 xg = ce.xg; 75 xgd = ce.xgd; 76 massmatrix = ce.massmatrix; 77 Hmatrix = ce.Hmatrix; 78 SV = ce.SV; 79 DS = ce.DS; 80 elasticforce_beam = ce.elasticforce_beam; 81 82 //integration points 83 x1 = ce.x1; x2 = ce.x2; x3 = ce.x3; 84 w1 = ce.w1; w2 = ce.w2; w3 = ce.w3; 85 orderx = ce.orderx; 86 orderyz = ce.orderyz; 87 88 for (int i=0; i < ANCFBeamMaxIP; i++) 89 { 90 grad[i] = ce.grad[i]; 91 jacdet[i] = ce.jacdet[i]; 92 } 93 T1 = ce.T1; 94 T2 = ce.T2; 95 T = ce.T; 96 e0 = ce.e0; 97 jac1 = ce.jac1; 98 jac2 = ce.jac2; 99 100 sos2 = ce.sos2; 101 n1 = ce.n1; n2 = ce.n2; 102 concentratedmass1 = ce.concentratedmass1; 103 concentratedmass2 = ce.concentratedmass2; 104 105 beamEA = ce.beamEA; 106 beamEIy = ce.beamEIy; 107 beamEIz = ce.beamEIz; 108 beamGJkx = ce.beamGJkx; 109 beamGAky = ce.beamGAky; 110 beamGAkz = ce.beamGAkz; 111 112 orderWl = ce.orderWl; 113 orderWs = ce.orderWs; 114 orderWsHR = ce.orderWsHR; 115 orderWt = ce.orderWt; 116 orderWb = ce.orderWb; 117 118 Wlepsx0 = ce.Wlepsx0 ; 119 Wlepsy0 = ce.Wlepsy0 ; 120 Wlepsz0 = ce.Wlepsz0 ; 121 Wlgamyz0= ce.Wlgamyz0; 122 WsHRk10 = ce.WsHRk10 ; 123 WsHRk20 = ce.WsHRk20 ; 124 Wtkapx0 = ce.Wtkapx0 ; 125 Wbkapy0 = ce.Wbkapy0 ; 126 Wbkapz0 = ce.Wbkapz0 ; 127 128 factstiffWl = ce.factstiffWl; 129 130 rho= ce.rho; //DR 2013-02-04 deleted rho from class element 131 } 132 Initialize()133 virtual void Initialize() 134 { 135 Body3D::Initialize(); 136 } 137 virtual void LinkToElements(); 138 virtual void BuildDSMatrices(); InitConstructor()139 virtual void InitConstructor() //initialize while constructor is called 140 { 141 concentratedmass1 = 0; 142 concentratedmass2 = 0; 143 } 144 GetElementSpec()145 virtual const char* GetElementSpec() const {return "ANCFBeam";} 146 virtual int CheckConsistency(mystr& errorstr); //rv==0 --> OK, rv==1 --> can not compute, rv==2 --> can not draw and not compute 147 virtual void GetElementData(ElementDataContainer& edc); //fill in all element data 148 virtual int SetElementData(ElementDataContainer& edc); //set element data acc. to ElementDataContainer, return 0 if failed or values invalid! 149 SOS()150 virtual int SOS() const {return 24;}; //size of K and M SOSowned()151 virtual int SOSowned() const {return sos2;}; //len(u) ES()152 virtual int ES() const {return 0;}; //size of first order explicit equations IS()153 virtual int IS() const {return 0;}; //implicit (algebraic) size 154 NNodes()155 virtual int NNodes() const {if (SOSowned() == 0) return 2; else return 0;}; NodeNum(int i)156 virtual const int& NodeNum(int i) const 157 { 158 if (i == 1) return n1; else return n2; 159 } NodeNum(int i)160 virtual int& NodeNum(int i) 161 { 162 if (i == 1) return n1; else return n2; 163 } GetNodeLocPos(int i)164 virtual Vector3D GetNodeLocPos(int i) const 165 { 166 if (i==1) return Vector3D(-0.5*lx,0.,0.); 167 else if (i==2) return Vector3D(0.5*lx,0.,0.); 168 169 return Vector3D(0.); 170 } IsRigid()171 virtual int IsRigid() const {return 0;} SetConcentratedMass(double cm1,double cm2)172 virtual void SetConcentratedMass(double cm1, double cm2) {concentratedmass1=cm1; concentratedmass2=cm2;} 173 174 virtual void SetInitialCurvatures(); //set initial vectors for curvature, stretch, shear deformation for Schwab-Meijaard model 175 virtual void SetRectangularBeamParameters(); //set beam parameters for rectangular cross-section, Schwab-Meijaard model 176 177 virtual void SetBeamParameters(double EA, double EIw, double EIv, double GJ, double GA); 178 179 virtual void SetBeamParameters2(double beamEAi, double beamEIyi, double beamEIzi, double beamGJkxi, 180 double beamGAkyi, double beamGAkzi); 181 182 // (AD) changed () to .Get() XGP(int iloc)183 virtual const double& XGP(int iloc) const {return GetXact(ltg.Get(iloc+24));} XGPD(int iloc)184 virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg.Get(iloc+24));} 185 // virtual const double& XGP(int iloc) const {return GetXact(ltg(iloc+24));} 186 // virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg(iloc+24));} 187 NS()188 virtual int NS() const {return 8;} 189 190 //compute v = T*v 191 virtual void ApplyT(Vector& v) const; 192 virtual void ApplyTD(Vector& v) const; 193 194 virtual void ApplyTtp(Vector& v) const; 195 196 virtual void ApplyTtp(Matrix& m) const; 197 198 virtual void GetS0(Vector& sf, const Vector3D& ploc) const; 199 200 virtual void GetDSMatrix0(const Vector3D& ploc, Matrix& sf) const; 201 202 virtual void GetS0x(Vector& sfx, const Vector3D& ploc) const; 203 virtual void GetS0y(Vector& sfx, const Vector3D& ploc) const; 204 virtual void GetS0z(Vector& sfx, const Vector3D& ploc) const; 205 virtual void GetS0yx(Vector& sfx, const Vector3D& ploc) const; 206 virtual void GetS0zx(Vector& sfx, const Vector3D& ploc) const; 207 virtual void GetS0xx(Vector& sfx, const Vector3D& ploc) const; 208 virtual void GetRot(Matrix3D& rot, const Vector& xg) const; 209 210 virtual double GetTau(const Vector& xg) const; 211 virtual void GetDeltaTau(const Vector& xg, Vector& deltatau) const; 212 virtual void GetDeltaKappa(const double& x, const Vector& xg, Vector& dkappa, double& kappa) const; //test only 213 GetSMatrix(const Vector & sf,Matrix & sm)214 virtual void GetSMatrix(const Vector& sf, Matrix& sm) const 215 { 216 for (int j = 1; j <= 3; j++) 217 for (int i = 1; i<=sf.Length(); i++) 218 sm(j,i) = sf(i); 219 } 220 GetCoordinates(Vector & dc)221 virtual void GetCoordinates(Vector& dc) const 222 { 223 for (int i = 1; i <= SOS(); i++) 224 dc(i) = XG(i); 225 } 226 GetCoordinatesP(Vector & dc)227 virtual void GetCoordinatesP(Vector& dc) const 228 { 229 for (int i = 1; i <= SOS(); i++) 230 dc(i) = XGP(i); 231 } 232 GetDrawCoordinates(Vector & dc)233 virtual void GetDrawCoordinates(Vector& dc) const 234 { 235 for (int i = 1; i <= SOS(); i++) 236 dc(i) = XGD(i); 237 } 238 GetDrawCoordinatesP(Vector & dc)239 virtual void GetDrawCoordinatesP(Vector& dc) const 240 { 241 for (int i = 1; i <= SOS(); i++) 242 dc(i) = XGPD(i); 243 } 244 245 virtual Vector3D GetPos(const Vector3D& p_loc) const; 246 247 virtual Vector3D GetVel(const Vector3D& p_loc) const; 248 249 virtual Vector3D GetPosD(const Vector3D& p_loc) const; 250 251 //in reference element coordinates (-1..1) 252 virtual Vector3D GetPos0D(const Vector3D& p_loc, double def_scale = 1) const; 253 254 virtual Vector3D GetVelD(const Vector3D& p_loc) const; 255 256 virtual Vector3D GetPosx(const Vector3D& p_loc) const; 257 258 virtual Vector3D GetDisplacement(const Vector3D& p_loc) const; 259 virtual Vector3D GetDisplacementD(const Vector3D& p_loc) const; 260 SetComputeCoordinates()261 virtual void SetComputeCoordinates() 262 { 263 for (int i = 1; i <= SOS(); i++) 264 xg(i) = XG(i); 265 } 266 267 virtual void GetH(Matrix& H); 268 269 virtual void EvalM(Matrix& m, double t); 270 GetJacobi(Matrix3D & jac,const Vector3D & p,const Matrix & DS,const Vector & x0)271 virtual void GetJacobi(Matrix3D& jac, const Vector3D& p, const Matrix& DS, const Vector& x0) const 272 { 273 int ns = NS(); 274 for (int j = 1; j <= 3; j++) 275 { 276 for (int i = 1; i <= 3; i++) 277 { 278 jac(i,j) = 0; 279 for (int k=1; k <= ns; k++) 280 { 281 jac(i,j) += DS(j,k)*x0((k-1)*3+i); 282 } 283 } 284 } 285 //global_uo << "jac=" << jac << "\n"; 286 } 287 GetDMatrix(Matrix & D,double nu,double em)288 void GetDMatrix(Matrix& D, double nu, double em) const 289 { 290 //double nu = 1./2.*la/(la+mu); 291 //double em = mu*(3.*la+2.*mu)/(la+mu); 292 D(1,1)=em*(-1.+nu)/(1.+nu)/(-1.+2.*nu);D(1,2)=em*nu/(1.+nu)/(1-2*nu);D(1,3)=em*nu/(1.+nu)/(1-2*nu);D(1,4)=0;D(1,5)=0;D(1,6)=0; 293 D(2,1)=em*nu/(1.+nu)/(1-2*nu);D(2,2)=em*(-1.+nu)/(1.+nu)/(-1.+2.*nu);D(2,3)=em*nu/(1.+nu)/(1-2*nu);D(2,4)=0;D(2,5)=0;D(2,6)=0; 294 D(3,1)=em*nu/(1.+nu)/(1-2*nu);D(3,2)=em*nu/(1.+nu)/(1-2*nu);D(3,3)=em*(-1.+nu)/(1.+nu)/(-1.+2.*nu);D(3,4)=0;D(3,5)=0;D(3,6)=0; 295 D(4,1)=0;D(4,2)=0;D(4,3)=0;D(4,4)=.5000000000*em/(1.+nu);D(4,5)=0;D(4,6)=0; 296 D(5,1)=0;D(5,2)=0;D(5,3)=0;D(5,4)=0;D(5,5)=.5000000000*em/(1.+nu);D(5,6)=0; 297 D(6,1)=0;D(6,2)=0;D(6,3)=0;D(6,4)=0;D(6,5)=0;D(6,6)=.5000000000*em/(1.+nu); 298 } 299 300 virtual void EvalF2(Vector& f, double t); 301 302 virtual int FastStiffnessMatrix() const; //1==Position derivatives done analytically, damping added with numerical differentiation 303 virtual void StiffnessMatrix(Matrix& m); //fill in sos x sos components, m might be larger 304 305 306 virtual Matrix3D GetRotMatrix(const Vector3D& ploc) const; 307 virtual Matrix3D GetRotMatrixP(const Vector3D& ploc) const; 308 virtual Matrix3D GetRotMatrixD(const Vector3D& ploc) const; 309 310 virtual Matrix3D ComputeTangentFrame(double x, const Vector& xg) const; 311 virtual Matrix3D ComputeTangentFrameP(double x, const Vector& xg, const Vector& xgp) const; 312 virtual void GetTangentFramedRotvdqT(const Vector3D& vloc, double x, const Vector& xg, Matrix& d); 313 314 virtual Matrix3D ComputeCrosssectionFrame(double x, const Vector& q) const; 315 virtual Matrix3D ComputeCrosssectionFrameP(double x, const Vector& q, const Vector& qp) const; 316 virtual void GetCrosssectionFramedRotvdqT(const Vector3D& vloc, double x, const Vector& q, Matrix& dAvdq); 317 318 GetAngularVel(const Vector3D & p_loc)319 virtual Vector3D GetAngularVel(const Vector3D& p_loc) const 320 { 321 Matrix3D omega_skew = GetRotMatrixP(p_loc)*GetRotMatrix(p_loc).GetTp(); 322 323 return Vector3D(omega_skew(2,3),omega_skew(1,3),-omega_skew(1,2)); 324 } 325 326 327 //ploc -1 ... +1 Gradu(const Vector3D & ploc,const Vector & u,Matrix3D & gradu)328 virtual void Gradu(const Vector3D& ploc, const Vector& u, Matrix3D& gradu) const 329 { 330 static Matrix DS; 331 GetDSMatrix0(ploc,DS); 332 333 Matrix3D jac, jacinv; 334 GetJacobi(jac,ploc,DS,e0); 335 336 jac.GetInverse(jacinv); 337 jacinv = jacinv.GetTp(); 338 339 static Matrix grad; 340 grad.SetSize(3,NS()); 341 Mult(jacinv, DS, grad); 342 343 gradu.SetAll(0); 344 int dim = Dim(); 345 int l; 346 for (int j = 1; j <= dim; j++) 347 { 348 for (int i = 1; i <= NS(); i++) 349 { 350 l = (i-1)*dim+j; 351 for (int k = 1; k <= dim; k++) 352 { 353 gradu(j,k) += grad(k,i)*u(l); 354 } 355 } 356 } 357 } 358 359 //for body loads: 360 //Computes f = d p_ref/d q * x ApplyDprefdq(Vector & f,const Vector3D & x)361 virtual void ApplyDprefdq(Vector& f, const Vector3D& x) 362 { 363 //fill in, f.Length is already set 364 UO() << "Not yet implemented\n"; 365 366 } 367 //Computes f = d rot_ref/d q * x, rot bedeutet rotation um x, y, und z-Achse ApplyDrotrefdq(Vector & f,const Vector3D & x)368 virtual void ApplyDrotrefdq(Vector& f, const Vector3D& x) 369 { 370 //fill in, f.Length is already set 371 UO() << "Not yet implemented\n"; 372 } ApplyDrotdq(Vector & f,const Vector3D & x)373 virtual void ApplyDrotdq(Vector& f, const Vector3D& x) 374 { 375 //fill in, f.Length is already set 376 UO() << "Not yet implemented\n"; 377 } 378 //only displacements, rotations makes no sense, even in rigid body 379 //->only for volumeloads (gravity ...) GetIntDuDq(Matrix & dudq)380 virtual void GetIntDuDq(Matrix& dudq) 381 { 382 GetH(dudq); 383 //UO() << "Not yet implemented\n"; 384 } 385 virtual void GetdRotvdqT(const Vector3D& vloc, const Vector3D& ploc, Matrix& d); 386 387 virtual void GetdPosdqT(const Vector3D& ploc, Matrix& d); 388 389 virtual void GetdRotdqT(const Vector3D& ploc, Matrix& d); 390 GetdPosdx(const Vector3D & ploc,Vector3D & dpdx)391 virtual void GetdPosdx(const Vector3D& ploc, Vector3D& dpdx) 392 { 393 //optimierungspotenzial 500% !!!!!!!!!!!!!!!!!!! 394 double p0 = ploc.X()/(0.5*lx); 395 396 SV.SetLen(NS()); 397 GetS0x(SV,p0); 398 399 static Vector xg; 400 xg.SetLen(SOS()); 401 GetCoordinates(xg); 402 ApplyT(xg); 403 dpdx = 0; 404 for (int i=1; i <= 3; i++) 405 { 406 for (int j=1; j <= NS(); j++) 407 { 408 dpdx(i) += SV(j) * xg((j-1)*3+i); 409 } 410 } 411 412 }; 413 414 virtual void ComputeCorotationalFrame(Vector3D& pref, Matrix3D& Aref); 415 virtual void ComputeCorotationalFrameD(Vector3D& pref, Matrix3D& Aref); 416 virtual void ComputeCorotationalFrameDAvdq(Vector3D& v, Matrix& dAvdq); 417 virtual void ComputeCorotationalFrameDATvdq(Vector3D& v, Matrix& dATvdq); 418 419 //ploc -1 ... +1 420 virtual void GraduD(const Vector3D& ploc, const Vector& u, Matrix3D& gradu) const; 421 422 // variables, available for post-processing and sensing 423 virtual void GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables); 424 // computation of the variables 425 virtual double GetFieldVariableValue(const FieldVariableDescriptor & fvd, const Vector3D & local_position, bool flagD); 426 427 virtual void DrawElement(); 428 429 protected: 430 //mechanical: 431 double lx, ly, lz, Em, nu; 432 433 int n1, n2, sos2; 434 double concentratedmass1, concentratedmass2; 435 436 //temporary storage for acceleration: 437 Vector xg, xgd; 438 Matrix massmatrix, Hmatrix; 439 Matrix DS; //3*8 440 Vector SV; //8 441 Vector temp; 442 int elasticforce_beam; 443 444 //integration points 445 Vector x1,x2,x3,w1,w2,w3; 446 int orderx, orderyz; 447 Matrix grad[ANCFBeamMaxIP]; //DS transformed, grad u = grad[i]*e 448 double jacdet[ANCFBeamMaxIP]; 449 450 Matrix T1, T2, T; //slope discontinuities 451 Matrix3D jac1, jac2; //slope discontinuities 452 Vector e0; //initial vector in e, not in p 453 454 Matrix K0; //stiffness matrix at initial configuration 455 456 double beamEA; //for beam formulation 457 double beamEIy; 458 double beamEIz; 459 double beamGJkx; 460 double beamGAky; 461 double beamGAkz; 462 463 double factstiffWl; //reduce stiffness of thickness modes, standard==1 464 465 int orderWl; //5; 5 and 6 coincide with 8 up to 10 digits, 7=8 466 int orderWs; //2; more than 2 gives locking, 1 is rather inaccurate 467 int orderWsHR; //4; for Hellinger Reissner, 4 is exact 468 int orderWt; //2; 1=2=4, 1 does not work in L-shape!!! 469 int orderWb; //3; 3=4=6, 2 gives shear locking 470 471 472 Vector Wlepsx0; 473 Vector Wlepsy0; 474 Vector Wlepsz0; 475 Vector Wlgamyz0; 476 Vector WsHRk10; 477 Vector WsHRk20; 478 Vector Wtkapx0; 479 Vector Wbkapy0; 480 Vector Wbkapz0; 481 482 double rho; //DR 2013-02-04 deleted rho from class element 483 }; 484 485 486 487 488 #endif 489 490