1 //#************************************************************** 2 //# 3 //# filename: ANCFCable3D.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 ANCFCABLE3D__H 25 #define ANCFCABLE3D__H 26 27 #include "body3d.h" 28 29 30 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 31 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 32 // ANCFCable3D ANCFCable3D ANCFCable3D ANCFCable3D ANCFCable3D ANCFCable3D ANCFCable3D 33 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 34 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 35 36 37 const int ANCFCableMaxIP=45; 38 //rigid cube 39 class ANCFCable3D: public Body3D //$EDC$[beginclass,classname=ANCFCable3D,parentclassname=Body3D] 40 { 41 public: 42 //Body3D():Element() {mbs = NULL;}; ANCFCable3D(MBS * mbsi)43 ANCFCable3D(MBS* mbsi):Body3D(mbsi), massmatrix(), Hmatrix(), SV(), x1(), w1(), 44 xg(), xgd(), e0() { SetElementName("ANCFCable3D"); }; ANCFCable3D(const ANCFCable3D & e)45 ANCFCable3D(const ANCFCable3D& e):Body3D(e.mbs),massmatrix(), Hmatrix(), SV(), x1(), w1(), 46 xg(), xgd(), e0() {CopyFrom(e);}; 47 //nodal coordinates of first and second point (x1,x2, x1.Length()==12) 48 ANCFCable3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, double rhoi, double Emi, 49 const Vector3D& si, const Vector3D& coli, int beammodel = 0, const Vector3D& directori=Vector3D(0,0,0)): Body3D(mbsi)50 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), x1(), w1(), 51 xg(), xgd(), e0() 52 { 53 SetElementName("ANCFCable3D"); 54 n1=0; n2=0; sos2=SOS(); 55 size = si; 56 lx = si.X(); ly = si.Y(); lz = si.Z(); 57 58 mass = lx*ly*lz*rhoi; 59 //lx=1;ly=1;lz=1; 60 61 x_init = xc1.Append(xc2); 62 xg = xc1.Append(xc2); 63 64 Em = Emi; 65 //rho = rhoi; //$ DR 2013-02-04 deleted rho from class element, do not use it here! 66 col = coli; 67 concentratedmass1 = 0; 68 concentratedmass2 = 0; 69 70 BuildDSMatrices(); 71 e0 = x_init; 72 x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!! 73 74 if (directori.X() != 0 || directori.Y() != 0 || directori.Z() != 0) 75 { 76 director.Set(directori.X(),directori.Y(),directori.Z()); 77 use_director = true; 78 } 79 else 80 { 81 use_director = false; 82 } 83 }; 84 85 //Element shares nodes with other elements, n1 and n2 are nodenumbers; element sets initial conditions for nodes 86 ANCFCable3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, int n1i, int n2i, double rhoi, double Emi, 87 const Vector3D& si, const Vector3D& coli, int beammodel = 0, const Vector3D& directori=Vector3D(0,0,0)): Body3D(mbsi)88 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), x1(), w1(), 89 xg(), xgd(), e0() 90 { 91 SetElementName("ANCFCable3D"); 92 n1=n1i; n2=n2i; sos2=0; 93 size = si; 94 lx = si.X(); ly = si.Y(); lz = si.Z(); 95 96 mass = lx*ly*lz*rhoi; 97 //lx=1;ly=1;lz=1; 98 99 x_init = xc1.Append(xc2); 100 xg = xc1.Append(xc2); 101 102 Em = Emi; 103 //rho = rhoi; //$ DR 2013-02-04 deleted rho from class element, do not use it here! 104 col = coli; 105 concentratedmass1 = 0; 106 concentratedmass2 = 0; 107 108 //UO() << "Cable: rho=" << rho << ", E=" << Em << ", size=" << size << ", n1=" << n1 << ", n2=" << n2 << "\n"; 109 110 BuildDSMatrices(); 111 e0 = x_init; 112 x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!! 113 114 if (directori.X() != 0 || directori.Y() != 0 || directori.Z() != 0) 115 { 116 director.Set(directori.X(),directori.Y(),directori.Z()); 117 use_director = true; 118 } 119 else 120 { 121 use_director = false; 122 } 123 }; 124 125 //Element shares nodes with other elements, n1 and n2 are nodenumbers; element sets initial conditions for nodes 126 ANCFCable3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, const Vector& vc1, const Vector& vc2, 127 int n1i, int n2i, double rhoi, double Emi, 128 const Vector3D& si, const Vector3D& coli, int beammodel = 0, const Vector3D& directori=Vector3D(0,0,0)): Body3D(mbsi)129 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), x1(), w1(), 130 xg(), xgd(), e0() 131 { 132 SetElementName("ANCFCable3D"); 133 n1=n1i; n2=n2i; sos2=0; 134 size = si; 135 lx = si.X(); ly = si.Y(); lz = si.Z(); 136 137 mass = lx*ly*lz*rhoi; 138 //lx=1;ly=1;lz=1; 139 140 x_init = xc1.Append(xc2); 141 xg = xc1.Append(xc2); 142 Vector v_init = vc1.Append(vc2); 143 144 Em = Emi; 145 rho = rhoi; 146 col = coli; 147 concentratedmass1 = 0; 148 concentratedmass2 = 0; 149 150 BuildDSMatrices(); 151 e0 = x_init; 152 x_init = x_init.Append(v_init); //Velocity initial conditions can also be transformed by Tinv!!! 153 154 if (directori.X() != 0 || directori.Y() != 0 || directori.Z() != 0) 155 { 156 director.Set(directori.X(),directori.Y(),directori.Z()); 157 use_director = true; 158 } 159 else 160 { 161 use_director = false; 162 } 163 }; 164 165 //To be overwritten in derived class: GetCopy()166 virtual Element* GetCopy() 167 { 168 Element* ec = new ANCFCable3D(*this); 169 return ec; 170 } 171 //To be overwritten in derived class: CopyFrom(const Element & e)172 virtual void CopyFrom(const Element& e) 173 { 174 Body3D::CopyFrom(e); 175 const ANCFCable3D& ce = (const ANCFCable3D&)e; 176 lx = ce.lx; 177 ly = ce.ly; 178 lz = ce.lz; 179 Em = ce.Em; 180 xg = ce.xg; 181 xgd = ce.xgd; 182 massmatrix = ce.massmatrix; 183 Hmatrix = ce.Hmatrix; 184 SV = ce.SV; 185 186 rho= ce.rho; //DR 2013-02-04 deleted rho from class element 187 188 //integration points 189 x1 = ce.x1; 190 w1 = ce.w1; 191 orderx = ce.orderx; 192 193 e0 = ce.e0; 194 195 sos2 = ce.sos2; 196 n1 = ce.n1; n2 = ce.n2; 197 concentratedmass1 = ce.concentratedmass1; 198 concentratedmass2 = ce.concentratedmass2; 199 200 use_director = ce.use_director; 201 if (use_director) 202 { 203 director = ce.director; 204 } 205 } 206 Initialize()207 virtual void Initialize() 208 { 209 Body3D::Initialize(); 210 } 211 virtual void LinkToElements(); 212 virtual void BuildDSMatrices(); 213 214 virtual void GetElementDataAuto(ElementDataContainer& edc); //fill in all element data 215 virtual int SetElementDataAuto(ElementDataContainer& edc); //set element data acc. to ElementDataContainer, return 0 if failed or values invalid! 216 virtual int ReadSingleElementDataAuto(ReadWriteElementDataVariableType& RWdata); //automatically generated function from EDC_converter 217 virtual int WriteSingleElementDataAuto(const ReadWriteElementDataVariableType& RWdata); //automatically generated function from EDC_converter 218 virtual int GetAvailableSpecialValuesAuto(TArrayDynamic<ReadWriteElementDataVariableType>& available_variables); //automatically generated function from EDC_converter 219 GetElementSpec()220 virtual const char* GetElementSpec() const { return "ANCFCable3D"; } 221 SOS()222 virtual int SOS() const {return 12;}; //size of K and M SOSowned()223 virtual int SOSowned() const {return sos2;}; //len(u) ES()224 virtual int ES() const {return 0;}; //size of first order explicit equations IS()225 virtual int IS() const {return 0;}; //implicit (algebraic) size 226 IsRigid()227 virtual int IsRigid() const {return 0;} 228 SetConcentratedMass(double cm1,double cm2)229 virtual void SetConcentratedMass(double cm1, double cm2) {concentratedmass1=cm1; concentratedmass2=cm2;} 230 231 // (AD) changed () to .Get() XGP(int iloc)232 virtual const double& XGP(int iloc) const {return GetXact(ltg.Get(iloc+12));} XGPD(int iloc)233 virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg.Get(iloc+12));} 234 // virtual const double& XGP(int iloc) const {return GetXact(ltg(iloc+12));} 235 // virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg(iloc+12));} 236 NS()237 virtual int NS() const {return 4;} NNodes()238 virtual int NNodes() const {if (SOSowned() == 0) return 2; else return 0;}; NodeNum(int i)239 virtual const int& NodeNum(int i) const 240 { 241 if (i == 1) return n1; else return n2; 242 } NodeNum(int i)243 virtual int& NodeNum(int i) 244 { 245 if (i == 1) return n1; else return n2; 246 } 247 Rho()248 double& Rho() {return rho;} Rho()249 const double& Rho() const {return rho;} 250 251 252 virtual void GetS0(Vector& sf, const double& ploc) const; 253 virtual void GetS0x(Vector& sfx, const double& ploc) const; 254 virtual void GetS0xx(Vector& sfx, const double& ploc) const; 255 GetCoordinates(Vector & dc)256 virtual void GetCoordinates(Vector& dc) const 257 { 258 for (int i = 1; i <= SOS(); i++) 259 dc(i) = XG(i); 260 } 261 GetCoordinatesP(Vector & dc)262 virtual void GetCoordinatesP(Vector& dc) const 263 { 264 for (int i = 1; i <= SOS(); i++) 265 dc(i) = XGP(i); 266 } 267 GetDrawCoordinates(Vector & dc)268 virtual void GetDrawCoordinates(Vector& dc) const 269 { 270 for (int i = 1; i <= SOS(); i++) 271 dc(i) = XGD(i); 272 } 273 GetDrawCoordinatesP(Vector & dc)274 virtual void GetDrawCoordinatesP(Vector& dc) const 275 { 276 for (int i = 1; i <= SOS(); i++) 277 dc(i) = XGPD(i); 278 } 279 280 virtual Vector3D GetPos(const double& p_loc) const; 281 GetPosx0(const double & p_loc,const Vector & xg)282 virtual Vector3D GetPosx0(const double& p_loc, const Vector& xg) const 283 { 284 double p0=p_loc; 285 static Vector SV; 286 GetS0x(SV, p0); 287 Vector3D p(0.,0.,0.); 288 for (int i = 1; i <= Dim(); i++) 289 { 290 for (int j = 1; j <= NS(); j++) 291 { 292 p(i) += SV(j)*xg((j-1)*3+i); 293 } 294 } 295 return p; 296 } 297 GetPosx0D(const double & p_loc)298 virtual Vector3D GetPosx0D(const double& p_loc) const 299 { 300 double p0=p_loc; 301 static Vector SV; 302 GetS0x(SV, p0); 303 static Vector xgd; 304 xgd.SetLen(SOS()); 305 GetDrawCoordinates(xgd); 306 307 Vector3D p(0.,0.,0.); 308 for (int i = 1; i <= Dim(); i++) 309 { 310 for (int j = 1; j <= NS(); j++) 311 { 312 p(i) += SV(j)*xgd((j-1)*3+i); 313 } 314 } 315 return p; 316 } 317 GetPosxx0(const double & p_loc,const Vector & xg)318 virtual Vector3D GetPosxx0(const double& p_loc, const Vector& xg) const 319 { 320 double p0=p_loc; 321 static Vector SV; 322 GetS0xx(SV, p0); 323 Vector3D p(0.,0.,0.); 324 for (int i = 1; i <= Dim(); i++) 325 { 326 for (int j = 1; j <= NS(); j++) 327 { 328 p(i) += SV(j)*xg((j-1)*3+i); 329 } 330 } 331 return p; 332 } 333 334 virtual Vector3D GetVel(const double& p_loc) const; 335 virtual Vector3D GetPosD(const double& p_loc) const; GetRefPosD()336 virtual Vector3D GetRefPosD() const { return GetPosD(0); } //PG & KN GetDOFPosD(int idof)337 virtual Vector3D GetDOFPosD(int idof) const 338 { 339 //returns postion of i-th DOF 340 if (idof < 7) 341 { 342 return GetPosD(-lx/2.); 343 } 344 return GetPosD(lx/2.); 345 } 346 GetDOFDirD(int idof)347 virtual Vector3D GetDOFDirD(int idof) const //returns direction of action of i-th DOF 348 { 349 Matrix3D rot; 350 if (idof < 7) 351 rot = Matrix3D(1.);//GetInitRotMatrix3D(-lx*.5); 352 else 353 rot = Matrix3D(1.);//GetInitRotMatrix3D(lx*.5); 354 355 switch(idof) 356 { 357 case 1: case 7: return Vector3D(rot(1,1),rot(2,1),rot(3,1)); break; 358 case 2: case 8: return Vector3D(rot(1,2),rot(2,2),rot(3,2)); break; 359 case 3: case 9: return Vector3D(rot(1,3),rot(2,3),rot(3,3)); break; 360 } 361 return Vector3D(0.,0.,0.); 362 } 363 364 //in reference element coordinates (-1..1) 365 virtual Vector3D GetPos0D(const double& p_loc) const; 366 367 virtual Vector3D GetVelD(const double& p_loc) const; 368 GetPos(const Vector3D & p_loc)369 virtual Vector3D GetPos(const Vector3D& p_loc) const 370 { 371 if (use_director) 372 { 373 return GetPos(p_loc.X()) + GetRotMatrix(p_loc.X())*Vector3D(0.,p_loc.Y(),p_loc.Z()); 374 } 375 return GetPos(p_loc.X()); 376 } GetPosD(const Vector3D & p_loc)377 virtual Vector3D GetPosD(const Vector3D& p_loc) const 378 { 379 if (use_director) 380 { 381 return GetPosD(p_loc.X()) + GetRotMatrixD(p_loc.X())*Vector3D(0.,p_loc.Y(),p_loc.Z()); 382 } 383 return GetPosD(p_loc.X()); 384 } GetVel(const Vector3D & p_loc)385 virtual Vector3D GetVel(const Vector3D& p_loc) const {return GetVel(p_loc.X());} GetVelD(const Vector3D & p_loc)386 virtual Vector3D GetVelD(const Vector3D& p_loc) const {return GetVelD(p_loc.X());} 387 SetComputeCoordinates()388 virtual void SetComputeCoordinates() 389 { 390 for (int i = 1; i <= SOS(); i++) 391 xg(i) = XG(i); 392 } 393 394 virtual void GetH(Matrix& H); 395 396 virtual void EvalM(Matrix& m, double t); 397 398 virtual void EvalF2(Vector& f, double t); 399 400 virtual double GetKappa(const double& x, const Vector& xg) const; 401 virtual double GetEpsAxial(const double& x, const Vector& xg) const; 402 virtual void GetDeltaKappa(const double& x, const Vector& xg, Vector& dkappa, double& kappa) const; 403 virtual void GetDeltaEpsAxial(const double& x, const Vector& xg, Vector& depsaxial) const; 404 405 //for body loads: 406 //Computes f = d p_ref/d q * x ApplyDprefdq(Vector & f,const Vector3D & x)407 virtual void ApplyDprefdq(Vector& f, const Vector3D& x) 408 { 409 //fill in, f.Length is already set 410 UO() << "Not yet implemented\n"; 411 412 } 413 //Computes f = d rot_ref/d q * x, rot bedeutet rotation um x, y, und z-Achse ApplyDrotrefdq(Vector & f,const Vector3D & x)414 virtual void ApplyDrotrefdq(Vector& f, const Vector3D& x) 415 { 416 //fill in, f.Length is already set 417 UO() << "Not yet implemented\n"; 418 } 419 //only displacements, rotations makes no sense, even in rigid body 420 //->only for volumeloads (gravity ...) GetIntDuDq(Matrix & dudq)421 virtual void GetIntDuDq(Matrix& dudq) 422 { 423 GetH(dudq); 424 } GetRotMatv(const Vector3D & ploc,const Vector3D & v)425 virtual Vector3D GetRotMatv(const Vector3D& ploc, const Vector3D& v) 426 { 427 double diffpar = 1e-8; 428 Vector3D vrot = (1./diffpar)*(GetPos(ploc + diffpar*v) - GetPos(ploc)); 429 return vrot; 430 431 /* 432 Vector3D vrot = GetPos(ploc + diffpar*v) - GetPos(ploc); 433 vrot.Normalize(); 434 return vrot*v.Norm();*/ 435 } GetRotMatPv(const Vector3D & ploc,const Vector3D & v)436 virtual Vector3D GetRotMatPv(const Vector3D& ploc, const Vector3D& v) 437 { 438 double diffpar = 1e-8; 439 Vector3D vrot = (1./diffpar)*(GetVel(ploc + diffpar*v) - GetVel(ploc)); 440 return vrot; 441 //vrot.Normalize(); 442 //return vrot*v.Norm(); 443 } GetdRotvdqT(const Vector3D & vloc,const Vector3D & ploc,Matrix & d)444 virtual void GetdRotvdqT(const Vector3D& vloc, const Vector3D& ploc, Matrix& d) 445 { 446 //Only vloc = (1,0,0) is possible!!!!!!!!!!!!!!!!! 447 448 d.SetSize(NS()*Dim(),3); 449 static Matrix d2; 450 d2.SetSize(NS()*Dim(),3); 451 double diffpar = 1e-8; 452 GetdPosdqT(ploc+diffpar*vloc, d2); 453 d = d2; 454 //UO() << "d=" << d << "\n"; 455 GetdPosdqT(ploc,d2); 456 d -= d2; 457 d *= 1./diffpar; 458 } 459 GetdPosdqT(const Vector3D & ploc,Matrix & d)460 virtual void GetdPosdqT(const Vector3D& ploc, Matrix& d) 461 { 462 //p = S(p.x,p.y,p.z)*q; d(p)/dq 463 static Vector SV; 464 GetS0(SV, ploc.X()/(0.5*lx)); 465 d.SetSize(NS()*Dim(),3); 466 d.FillWithZeros(); 467 for (int i = 1; i <= NS(); i++) 468 { 469 d((i-1)*3+1,1)=SV(i); 470 d((i-1)*3+2,2)=SV(i); 471 d((i-1)*3+3,3)=SV(i); 472 } 473 } 474 GetdPosdx(const Vector3D & ploc,Vector3D & dpdx)475 virtual void GetdPosdx(const Vector3D& ploc, Vector3D& dpdx) 476 { 477 //optimierungspotenzial 500% !!!!!!!!!!!!!!!!!!! 478 double p0 = ploc.X()/(0.5*lx); 479 480 SV.SetLen(NS()); 481 GetS0x(SV,p0); 482 483 static Vector xg; 484 xg.SetLen(SOS()); 485 GetCoordinates(xg); 486 dpdx = 0; 487 for (int i=1; i <= 3; i++) 488 { 489 for (int j=1; j <= NS(); j++) 490 { 491 dpdx(i) += SV(j) * xg((j-1)*3+i); 492 } 493 } 494 }; 495 GetDirector()496 virtual const Vector3D& GetDirector() const {return director;} GetDirector()497 virtual Vector3D& GetDirector() {return director;} GetRotMatrix(const double p_loc)498 virtual Matrix3D GetRotMatrix(const double p_loc) const 499 { 500 // use only, if director is defined (see constructor) 501 assert (use_director); 502 503 // p_loc in [-lx/2,lx/2] 504 // p0 in [-1,1] 505 double p0=p_loc/(0.5*lx); 506 507 static Vector xg(SOS()); 508 for (int j=1; j<=SOS(); j++) 509 { 510 xg(j) = XG(j); 511 } 512 513 Vector3D rx = GetPosx0(p0, xg); 514 rx.Normalize(); 515 516 Vector3D e20 = GetDirector(); 517 rx.GramSchmidt(e20); //d is projected into plane normal to rx (GramSchmidt() already returns a normalized Vector3D) 518 519 Vector3D e30 = rx.Cross(e20); //length=1 520 521 Matrix3D A(rx.X(),e20.X(),e30.X(), 522 rx.Y(),e20.Y(),e30.Y(), 523 rx.Z(),e20.Z(),e30.Z()); 524 525 return A; 526 } GetRotMatrixD(const double p_loc)527 virtual Matrix3D GetRotMatrixD(const double p_loc) const 528 { 529 // use only, if director is defined (see constructor) 530 assert (use_director); 531 532 // p_loc in [-lx/2,lx/2] 533 // p0 in [-1,1] 534 double p0=p_loc/(0.5*lx); 535 536 Vector3D rx = GetPosx0D(p0); 537 rx.Normalize(); 538 539 Vector3D e20 = GetDirector(); 540 rx.GramSchmidt(e20); //d is projected into plane normal to rx (GramSchmidt() already returns a normalized Vector3D) 541 542 Vector3D e30 = rx.Cross(e20); //length=1 543 544 Matrix3D A(rx.X(),e20.X(),e30.X(), 545 rx.Y(),e20.Y(),e30.Y(), 546 rx.Z(),e20.Z(),e30.Z()); 547 548 return A; 549 } 550 551 virtual void DrawElement(); 552 553 protected: 554 //mechanical: 555 double lx, ly, lz, Em; 556 557 int n1, n2, sos2; 558 559 double concentratedmass1, concentratedmass2; 560 //temporary storage for acceleration: 561 Vector xg, xgd; 562 Matrix massmatrix, Hmatrix; 563 Vector SV; //4 564 Vector temp; 565 566 //integration points 567 Vector x1,w1; 568 int orderx; 569 570 Vector e0; //initial vector in e, not in p 571 572 Vector3D director; 573 bool use_director; 574 575 double rho; //DR 2013-02-04 deleted rho from class element 576 577 }; //$EDC$[endclass,ANCFCable3D] 578 579 580 581 582 #endif 583