1 //#************************************************************** 2 //# 3 //# filename: ANCFPlate3D.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 ANCFPlate3D__H 25 #define ANCFPlate3D__H 26 27 28 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 29 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 30 // ANCFPlate3D ANCFPlate3D ANCFPlate3D ANCFPlate3D ANCFPlate3D ANCFPlate3D 31 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 32 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 33 34 35 const int ANCFPlateMaxIP=120; //7*4*4 or 4*5*5 for plate ?? 36 //rigid cube 37 class ANCFPlate3D: public Body3D 38 { 39 public: 40 //Body3D():Element() {mbs = NULL;}; ANCFPlate3D(MBS * mbsi)41 ANCFPlate3D(MBS* mbsi):Body3D(mbsi),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 42 T(), xg(), xgd(), e0() {}; ANCFPlate3D(const ANCFPlate3D & e)43 ANCFPlate3D(const ANCFPlate3D& e):Body3D(e.mbs),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 44 T(), xg(), xgd(), e0() {CopyFrom(e);}; 45 46 //nodal coordinates of first and second point (x1,x2, x1.Length()==12) 47 //initial velocities!!!!! ANCFPlate3D(MBS * mbsi,const Vector & xc1,const Vector & xc2,const Vector & xc3,const Vector & xc4,const Vector & vc1,const Vector & vc2,const Vector & vc3,const Vector & vc4,double rhoi,double Emi,double nui,const Vector3D & si,const Vector3D & coli)48 ANCFPlate3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, const Vector& xc3, const Vector& xc4, 49 const Vector& vc1, const Vector& vc2, const Vector& vc3, const Vector& vc4, 50 double rhoi, double Emi, double nui, const Vector3D& si, const Vector3D& coli): 51 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 52 T(), xg(), xgd(), e0() 53 { 54 node[0]=0; node[1]=0; node[2]=0; node[3]=0; 55 sos2=SOS(); 56 size = si; 57 lx = si.X(); ly = si.Y(); lz = si.Z(); 58 59 mass = lx*ly*lz*rhoi; 60 61 int sadd = SOS()/NNodes()-4*Dim(); 62 x_init = (((xc1.Append(xc2)).Append(xc3)).Append(xc4)).Append(Vector(sadd*NNodes())); 63 xg = x_init; 64 Vector v_init = (((vc1.Append(vc2)).Append(vc3)).Append(vc4)).Append(Vector(sadd*NNodes())); 65 66 nu = nui; 67 Em = Emi; 68 //rho = rhoi; //$ DR 2013-02-04 deleted rho from class element, do not use it here! 69 col = coli; 70 concentratedmass[0] = concentratedmass[1] = concentratedmass[2] = concentratedmass[3] = 0; 71 72 BuildDSMatrices(); 73 Matrix Tinv = T; 74 if (!Tinv.Invert2()) {UO() << "ERROR: T matrix not invertible!!!\n";} 75 e0 = x_init; 76 x_init = Tinv * x_init; 77 v_init = Tinv * v_init; 78 //x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!! 79 x_init = x_init.Append(v_init); //Velocity initial conditions can also be transformed by Tinv!!! 80 }; 81 82 //Element shares nodes with other elements, n1i, n2i, n3i, n4i are nodenumbers; element sets initial conditions for nodes ANCFPlate3D(MBS * mbsi,const Vector & xc1,const Vector & xc2,const Vector & xc3,const Vector & xc4,int n1i,int n2i,int n3i,int n4i,double rhoi,double Emi,double nui,const Vector3D & si,const Vector3D & coli)83 ANCFPlate3D(MBS* mbsi, const Vector& xc1, const Vector& xc2, const Vector& xc3, const Vector& xc4, 84 int n1i, int n2i, int n3i, int n4i, double rhoi, double Emi, double nui, 85 const Vector3D& si, const Vector3D& coli): 86 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 87 T(), xg(), xgd(), e0() 88 { 89 node[0]=n1i; node[1]=n2i; node[2]=n3i; node[3]=n4i; 90 sos2=(NS()-4*NNodes())*Dim(); 91 size = si; 92 lx = si.X(); ly = si.Y(); lz = si.Z(); 93 94 mass = lx*ly*lz*rhoi; 95 96 int sadd = SOS()/NNodes()-4*Dim(); 97 x_init = (((xc1.Append(xc2)).Append(xc3)).Append(xc4)).Append(Vector(sadd*NNodes())); 98 xg = x_init; 99 100 nu = nui; 101 Em = Emi; 102 //rho = rhoi; //$ DR 2013-02-04 deleted rho from class element, do not use it here! 103 col = coli; 104 concentratedmass[0] = concentratedmass[1] = concentratedmass[2] = concentratedmass[3] = 0; 105 106 BuildDSMatrices(); 107 Matrix Tinv = T; 108 if (!Tinv.Invert2()) {UO() << "ERROR: T matrix not invertible!!!\n";} 109 e0 = x_init; 110 x_init = Tinv * x_init; 111 x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!! 112 }; 113 114 //Element shares nodes with other elements, n1i, n2i, n3i, n4i are nodenumbers; element sets initial conditions for nodes 115 //for a flat plate element, nodal gradients are computed for a flat plate ANCFPlate3D(MBS * mbsi,const Vector3D & xc1,const Vector3D & xc2,const Vector3D & xc3,const Vector3D & xc4,const Vector3D & vc1,const Vector3D & vc2,const Vector3D & vc3,const Vector3D & vc4,int n1i,int n2i,int n3i,int n4i,double rhoi,double Emi,double nui,const Vector3D & si,const Vector3D & coli)116 ANCFPlate3D(MBS* mbsi, const Vector3D& xc1, const Vector3D& xc2, const Vector3D& xc3, const Vector3D& xc4, 117 const Vector3D& vc1, const Vector3D& vc2, const Vector3D& vc3, const Vector3D& vc4, 118 int n1i, int n2i, int n3i, int n4i, double rhoi, double Emi, double nui, 119 const Vector3D& si, const Vector3D& coli): 120 Body3D(mbsi),massmatrix(), Hmatrix(), SV(), DS(), x1(), x2(), x3(), w1(), w2(), w3(), 121 T(), xg(), xgd(), e0() 122 { 123 node[0]=n1i; node[1]=n2i; node[2]=n3i; node[3]=n4i; 124 sos2=(NS()-4*NNodes())*Dim(); 125 size = si; 126 lx = si.X(); ly = si.Y(); lz = si.Z(); 127 128 mass = lx*ly*lz*rhoi; 129 130 int sadd = SOS()/NNodes()-4*Dim(); 131 Vector3D xx[4]; 132 xx[0] = xc1; 133 xx[1] = xc2; 134 xx[2] = xc3; 135 xx[3] = xc4; 136 Vector nc[4]; //nodal coordinate, 12 components 137 for (int i=0; i < 4; i++) 138 { 139 nc[i].Init(); nc[i].SetLen(12); 140 } 141 Vector3D nvx = xx[1]-xx[0]; //compute gradients for straight plate 142 Vector3D nvy = xx[3]-xx[0]; 143 nvx.Normalize(); 144 nvy.Normalize(); 145 Vector3D nvz = nvx.Cross(nvy); 146 for (int i=0; i < 4; i++) 147 { 148 nc[i]( 1) = xx[i](1);nc[i]( 2) = xx[i](2);nc[i]( 3) = xx[i](3); 149 nc[i]( 4) = nvx(1); nc[i]( 5) = nvx(2); nc[i]( 6) = nvx(3); 150 nc[i]( 7) = nvy(1); nc[i]( 8) = nvy(2); nc[i]( 9) = nvy(3); 151 nc[i](10) = nvz(1); nc[i](11) = nvz(2); nc[i](12) = nvz(3); 152 } 153 154 155 Vector nv1(12); //nodal velocities 156 Vector nv2(12); 157 Vector nv3(12); 158 Vector nv4(12); 159 nv1(1) = vc1(1); nv1(2) = vc1(2); nv1(3) = vc1(3); //set initial velocities, no initial gradient derivatives! 160 nv2(1) = vc2(1); nv2(2) = vc2(2); nv2(3) = vc2(3); 161 nv3(1) = vc3(1); nv3(2) = vc3(2); nv3(3) = vc3(3); 162 nv4(1) = vc4(1); nv4(2) = vc4(2); nv4(3) = vc4(3); 163 164 x_init = (((nc[0].Append(nc[1])).Append(nc[2])).Append(nc[3])).Append(Vector(sadd*NNodes())); 165 xg = x_init; 166 167 nu = nui; 168 Em = Emi; 169 rho = rhoi; 170 col = coli; 171 concentratedmass[0] = concentratedmass[1] = concentratedmass[2] = concentratedmass[3] = 0; 172 173 BuildDSMatrices(); 174 Matrix Tinv = T; 175 if (!Tinv.Invert2()) {UO() << "ERROR: T matrix not invertible!!!\n";} 176 e0 = x_init; 177 178 Vector v_init = (((nv1.Append(nv2)).Append(nv3)).Append(nv4)).Append(Vector(sadd*NNodes())); 179 180 x_init = Tinv * x_init; 181 v_init = Tinv * v_init; 182 183 x_init = x_init.Append(v_init); //Velocity initial conditions can also be transformed by Tinv!!! 184 }; 185 186 //To be overwritten in derived class: GetCopy()187 virtual Element* GetCopy() 188 { 189 Element* ec = new ANCFPlate3D(*this); 190 return ec; 191 } 192 //To be overwritten in derived class: CopyFrom(const Element & e)193 virtual void CopyFrom(const Element& e) 194 { 195 Body3D::CopyFrom(e); 196 const ANCFPlate3D& ce = (const ANCFPlate3D&)e; 197 lx = ce.lx; 198 ly = ce.ly; 199 lz = ce.lz; 200 Em = ce.Em; 201 nu = ce.nu; 202 //size = ce.size; 203 xg = ce.xg; 204 xgd = ce.xgd; 205 massmatrix = ce.massmatrix; 206 Hmatrix = ce.Hmatrix; 207 temp = ce.temp; 208 SV = ce.SV; 209 DS = ce.DS; 210 211 //integration points 212 x1 = ce.x1; x2 = ce.x2; x3 = ce.x3; 213 w1 = ce.w1; w2 = ce.w2; w3 = ce.w3; 214 orderxy = ce.orderxy; 215 orderz = ce.orderz; 216 217 for (int i=0; i < ANCFPlateMaxIP; i++) 218 { 219 grad[i] = ce.grad[i]; 220 jacdet[i] = ce.jacdet[i]; 221 } 222 T = ce.T; 223 e0 = ce.e0; 224 225 sos2 = ce.sos2; 226 for (int i = 0; i < NNodes(); i++) 227 { 228 Ti[i] = ce.Ti[i]; 229 jac[i] = ce.jac[i]; 230 node[i] = ce.node[i]; 231 concentratedmass[i] = ce.concentratedmass[i]; 232 } 233 234 rho = ce.rho; //DR 2013-02-04 deleted rho from class element 235 } 236 Initialize()237 virtual void Initialize() 238 { 239 Body3D::Initialize(); 240 } 241 virtual void LinkToElements(); 242 virtual void BuildDSMatrices(); 243 SOS()244 virtual int SOS() const {return Dim()*NS();}; //size of K and M SOSowned()245 virtual int SOSowned() const {return sos2;}; //len(u) ES()246 virtual int ES() const {return 0;}; //size of first order explicit equations IS()247 virtual int IS() const {return 0;}; //implicit (algebraic) size 248 IsRigid()249 virtual int IsRigid() const {return 0;} SetConcentratedMass(double cm,int i)250 virtual void SetConcentratedMass(double cm, int i) {concentratedmass[i-1] = cm;} 251 252 // (AD) changed () to .Get() XGP(int iloc)253 virtual const double& XGP(int iloc) const {return GetXact(ltg.Get(iloc+SOS()));} XGPD(int iloc)254 virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg.Get(iloc+SOS()));} 255 // virtual const double& XGP(int iloc) const {return GetXact(ltg(iloc+SOS()));} 256 // virtual const double& XGPD(int iloc) const {return mbs->GetDrawValue(ltg(iloc+SOS()));} 257 258 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 259 //change for other shape functions and dimension: Dim()260 virtual int Dim() const {return 3;} //default value 261 virtual int NS() const; NNodes()262 virtual int NNodes() const {return 4;} NodeSize()263 virtual int NodeSize() const {return 12;} 264 virtual void GetS0(Vector& sf, const Vector3D& ploc) const; 265 virtual void GetS0x(Vector& sfx, const Vector3D& ploc) const; 266 virtual void GetDSMatrix0(const Vector3D& ploc, Matrix& sf) const; 267 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 268 269 //compute v = T*v 270 virtual void ApplyT(Vector& v) const; 271 272 virtual void ApplyTtp(Vector& v) const; 273 274 virtual void ApplyTtp(Matrix& m) const; 275 GetSMatrix(const Vector & sf,Matrix & sm)276 virtual void GetSMatrix(const Vector& sf, Matrix& sm) const 277 { 278 for (int j = 1; j <= Dim(); j++) 279 for (int i = 1; i<=sf.Length(); i++) 280 sm(j,i) = sf(i); 281 } 282 GetCoordinates(Vector & dc)283 virtual void GetCoordinates(Vector& dc) const 284 { 285 for (int i = 1; i <= SOS(); i++) 286 dc(i) = XG(i); 287 } 288 GetCoordinatesP(Vector & dc)289 virtual void GetCoordinatesP(Vector& dc) const 290 { 291 for (int i = 1; i <= SOS(); i++) 292 dc(i) = XGP(i); 293 } 294 GetDrawCoordinates(Vector & dc)295 virtual void GetDrawCoordinates(Vector& dc) const 296 { 297 for (int i = 1; i <= SOS(); i++) 298 dc(i) = XGD(i); 299 } 300 GetDrawCoordinatesP(Vector & dc)301 virtual void GetDrawCoordinatesP(Vector& dc) const 302 { 303 for (int i = 1; i <= SOS(); i++) 304 dc(i) = XGPD(i); 305 } 306 307 virtual Vector3D GetPos(const Vector3D& p_loc) const; 308 virtual Vector3D GetPos(const Vector3D& p_loc, const Vector& xg) const; 309 310 virtual Vector3D GetVel(const Vector3D& p_loc) const; 311 312 //from -1 to +1 313 virtual Vector3D GetPosx(const Vector3D& p_loc, const Vector& xg) const; 314 315 virtual void GetRot(Matrix3D& rot, const Vector& xg) const; 316 317 virtual Vector3D GetPosD(const Vector3D& p_loc) const; 318 319 //in reference element coordinates (-1..1) 320 virtual Vector3D GetPos0D(const Vector3D& p_loc) const; 321 322 virtual Vector3D GetVelD(const Vector3D& p_loc) const; 323 SetComputeCoordinates()324 virtual void SetComputeCoordinates() 325 { 326 for (int i = 1; i <= SOS(); i++) 327 xg(i) = XG(i); 328 } 329 330 virtual void GetH(Matrix& H); 331 332 virtual void EvalM(Matrix& m, double t); 333 GetJacobi(Matrix3D & jac,const Vector3D & p,const Matrix & DS,const Vector & x0)334 virtual void GetJacobi(Matrix3D& jac, const Vector3D& p, const Matrix& DS, const Vector& x0) const 335 { 336 int ns = NS(); 337 for (int j = 1; j <= Dim(); j++) 338 { 339 for (int i = 1; i <= Dim(); i++) 340 { 341 jac(i,j) = 0; 342 for (int k=1; k <= ns; k++) 343 { 344 jac(i,j) += DS(j,k)*x0((k-1)*Dim()+i); 345 } 346 } 347 } 348 //global_uo << "jac=" << jac << "\n"; 349 } 350 351 virtual void EvalF2(Vector& f, double t); 352 GetRotMatrix(const Vector3D & ploc)353 virtual Matrix3D GetRotMatrix(const Vector3D& ploc) const 354 { 355 //Compute Gradient ... 356 static Vector u; 357 u.SetLen(SOS()); 358 GetCoordinates(u); 359 360 ApplyT(u); 361 u -= e0; 362 Vector3D p0(ploc); 363 p0.Scale(0.5*lx,0.5*ly,0.5*lz); 364 Matrix3D G; 365 Gradu(p0, u, G); 366 return G+Matrix3D(1); 367 } GetRotMatrixP(const Vector3D & ploc)368 virtual Matrix3D GetRotMatrixP(const Vector3D& ploc) const 369 { 370 //Compute Gradient ... 371 static Vector u; 372 u.SetLen(SOS()); 373 GetCoordinatesP(u); 374 375 ApplyT(u); 376 Vector3D p0(ploc); 377 p0.Scale(0.5*lx,0.5*ly,0.5*lz); 378 Matrix3D G; 379 Gradu(p0, u, G); 380 return G; 381 } GetRotMatrixD(const Vector3D & ploc)382 virtual Matrix3D GetRotMatrixD(const Vector3D& ploc) const 383 { 384 //Compute Gradient ... 385 static Vector u; 386 u.SetLen(SOS()); 387 GetDrawCoordinates(u); 388 389 u = T*u; 390 u -= e0; 391 Vector3D p0(ploc); 392 p0.Scale(0.5*lx,0.5*ly,0.5*lz); 393 Matrix3D G; 394 GraduD(p0, u, G); 395 return G+Matrix3D(1); 396 } 397 398 //ploc -1 ... +1 Gradu(const Vector3D & ploc,const Vector & u,Matrix3D & gradu)399 virtual void Gradu(const Vector3D& ploc, const Vector& u, Matrix3D& gradu) const 400 { 401 static Matrix DS; 402 GetDSMatrix0(ploc,DS); 403 404 Matrix3D jac, jacinv; 405 GetJacobi(jac,ploc,DS,e0); 406 407 jac.GetInverse(jacinv); 408 jacinv = jacinv.GetTp(); 409 410 static Matrix grad; 411 grad.SetSize(Dim(),NS()); 412 Mult(jacinv, DS, grad); 413 414 gradu.SetAll(0); 415 int dim = Dim(); 416 int l; 417 for (int j = 1; j <= dim; j++) 418 { 419 for (int i = 1; i <= NS(); i++) 420 { 421 l = (i-1)*dim+j; 422 for (int k = 1; k <= dim; k++) 423 { 424 gradu(j,k) += grad(k,i)*u(l); 425 } 426 } 427 } 428 } 429 430 //for body loads: 431 //Computes f = d p_ref/d q * x ApplyDprefdq(Vector & f,const Vector3D & x)432 virtual void ApplyDprefdq(Vector& f, const Vector3D& x) 433 { 434 //fill in, f.Length is already set 435 UO() << "Not yet implemented\n"; 436 437 } 438 //Computes f = d rot_ref/d q * x, rot bedeutet rotation um x, y, und z-Achse ApplyDrotrefdq(Vector & f,const Vector3D & x)439 virtual void ApplyDrotrefdq(Vector& f, const Vector3D& x) 440 { 441 //fill in, f.Length is already set 442 UO() << "Not yet implemented\n"; 443 } 444 //only displacements, rotations makes no sense, even in rigid body 445 //->only for volumeloads (gravity ...) GetIntDuDq(Matrix & dudq)446 virtual void GetIntDuDq(Matrix& dudq) 447 { 448 GetH(dudq); 449 } GetdRotvdqT(const Vector3D & vloc,const Vector3D & ploc,Matrix & d)450 virtual void GetdRotvdqT(const Vector3D& vloc, const Vector3D& ploc, Matrix& d) 451 { 452 d.SetSize(NS()*Dim(),Dim()); 453 static Matrix d2; 454 d2.SetSize(NS()*Dim(),Dim()); 455 double diffpar = 1e-8; 456 GetdPosdqT(ploc+diffpar*vloc, d2); 457 d = d2; 458 GetdPosdqT(ploc,d2); 459 d -= d2; 460 d *= 1./diffpar; 461 462 } 463 GetdPosdqT(const Vector3D & ploc,Matrix & d)464 virtual void GetdPosdqT(const Vector3D& ploc, Matrix& d) 465 { 466 //p = S(p.x,p.y,p.z)*q; d(p)/dq 467 Vector3D p0=ploc; 468 p0.Scale(0.5*lx,0.5*ly,0.5*lz); 469 static Vector SV; 470 GetS0(SV, p0); 471 d.SetSize(NS()*Dim(),Dim()); 472 d.FillWithZeros(); 473 for (int i = 1; i <= NS(); i++) 474 { 475 d((i-1)*3+1,1)=SV(i); 476 d((i-1)*3+2,2)=SV(i); 477 d((i-1)*3+3,3)=SV(i); 478 } 479 ApplyTtp(d); 480 //UO() << "Not yet implemented\n"; 481 } 482 GetdRotdqT(const Vector3D & ploc,Matrix & d)483 void GetdRotdqT(const Vector3D& ploc, Matrix& d) 484 { 485 //*********************************** 486 //vector rz for y-rotation 487 //vector rz for x-rotation 488 //vector rx/ry for z-rotation 489 490 //damit gleich wie 491 //Compute Gradient ... 492 493 494 static Vector xg; 495 xg.SetLen(SOS()); 496 GetCoordinates(xg); 497 ApplyT(xg); 498 499 Vector3D p0(ploc); 500 p0.Scale(0.5*lx,0.5*ly,0.5*lz); 501 502 static Matrix DS; 503 GetDSMatrix0(p0,DS); 504 505 Matrix3D jac, jacinv; 506 GetJacobi(jac,p0,DS,e0); 507 508 jac.GetInverse(jacinv); 509 jacinv = jacinv.GetTp(); 510 511 static Matrix grad; 512 Matrix3D F; 513 grad.SetSize(Dim(),NS()); 514 Mult(jacinv, DS, grad); 515 516 F.SetAll(0); 517 int dim = Dim(); 518 int l; 519 for (int j = 1; j <= dim; j++) 520 { 521 for (int i = 1; i <= NS(); i++) 522 { 523 l = (i-1)*dim+j; 524 for (int k = 1; k <= dim; k++) 525 { 526 F(j,k) += grad(k,i)*xg(l); 527 } 528 } 529 } 530 531 static Vector Sx, Sy, Sz; 532 Sx.SetLen(NS()); 533 Sy.SetLen(NS()); 534 Sz.SetLen(NS()); 535 for (int i = 1; i <= NS(); i++) 536 { 537 Sx(i) = grad(1,i); 538 Sy(i) = grad(2,i); 539 Sz(i) = grad(3,i); 540 } 541 542 543 Vector3D rx(F(1,1),F(2,1),F(3,1)); 544 Vector3D ry(F(1,2),F(2,2),F(3,2)); 545 Vector3D rz(F(1,3),F(2,3),F(3,3)); 546 547 double dyx = Sqr(ry.X())+Sqr(ry.Y()); 548 double dxy = Sqr(rx.Y())+Sqr(rx.X()); 549 550 double dzx = Sqr(rz.X())+Sqr(rz.Z()); 551 //double dxz = Sqr(rx.Z())+Sqr(rx.X()); 552 553 double dzy = Sqr(rz.Y())+Sqr(rz.Z()); 554 //double dyz = Sqr(ry.Y())+Sqr(ry.Z()); 555 556 for (int i = 1; i <= NS(); i++) 557 { 558 //moment x/y is applied to surface normal vector r,z 559 //from delta gamma_x 560 d((i-1)*3+1,1) = 0; //only delta r.X() terms! 561 d((i-1)*3+2,1) =-(rz.Z()*Sz(i))/dzy; //only delta r.Y() terms! 562 d((i-1)*3+3,1) = (rz.Y()*Sz(i))/dzy; //only delta r.Z() terms! 563 564 //from delta gamma_y 565 d((i-1)*3+1,2) = (rz.Z()*Sz(i))/dzx; //only delta r.X() terms! 566 d((i-1)*3+2,2) = 0; //only delta r.Y() terms! 567 d((i-1)*3+3,2) =-(rz.X()*Sz(i))/dzx; //only delta r.Z() terms! 568 569 //from delta gamma_z 570 d((i-1)*3+1,3) =-(ry.Y()*Sy(i))/dyx - (rx.Y()*Sx(i))/dyx; //only delta r.X() terms! 571 d((i-1)*3+2,3) = (ry.X()*Sy(i))/dyx + (rx.X()*Sx(i))/dyx; //only delta r.Y() terms! 572 d((i-1)*3+3,3) = 0; //only delta r.Z() terms! 573 574 } 575 //ApplyTtp(d);??? 576 577 } 578 579 GetdPosdx(const Vector3D & ploc,Vector3D & dpdx)580 virtual void GetdPosdx(const Vector3D& ploc, Vector3D& dpdx) 581 { 582 Vector3D p0(ploc.X()/(0.5*lx),0,0); 583 584 SV.SetLen(NS()); 585 GetS0x(SV,p0); 586 587 static Vector xg; 588 xg.SetLen(SOS()); 589 GetCoordinates(xg); 590 ApplyT(xg); 591 dpdx = 0; 592 for (int i=1; i <= 3; i++) 593 { 594 for (int j=1; j <= NS(); j++) 595 { 596 dpdx(i) += SV(j) * xg((j-1)*3+i); 597 } 598 } 599 }; 600 601 //ploc -1 ... +1 GraduD(const Vector3D & ploc,const Vector & u,Matrix3D & gradu)602 virtual void GraduD(const Vector3D& ploc, const Vector& u, Matrix3D& gradu) const 603 { 604 static Matrix DS; 605 GetDSMatrix0(ploc,DS); 606 607 Matrix3D jac, jacinv; 608 GetJacobi(jac,ploc,DS,e0); 609 610 jac.GetInverse(jacinv); 611 jacinv = jacinv.GetTp(); 612 613 static Matrix grad; 614 grad.SetSize(Dim(),NS()); 615 Mult(jacinv, DS, grad); 616 617 gradu.SetAll(0); 618 int dim = Dim(); 619 int l; 620 for (int j = 1; j <= dim; j++) 621 { 622 for (int i = 1; i <= NS(); i++) 623 { 624 l = (i-1)*dim+j; 625 for (int k = 1; k <= dim; k++) 626 { 627 gradu(j,k) += grad(k,i)*u(l); 628 } 629 } 630 } 631 } 632 633 // variables, available for post-processing and sensing 634 virtual void GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables); 635 // computation of the variables 636 virtual double GetFieldVariableValue(const FieldVariableDescriptor & fvd, const Vector3D & local_position, bool flagD); 637 638 virtual void DrawElement(); 639 640 protected: 641 //mechanical: 642 double lx, ly, lz, Em, nu; 643 644 int node[4]; 645 int sos2; 646 double concentratedmass[4]; 647 648 //temporary storage for acceleration: 649 Vector xg, xgd; 650 Matrix massmatrix; 651 Matrix Hmatrix; 652 Matrix DS; //3*8 653 Vector SV; //8 654 Vector temp; 655 656 //integration points 657 Vector x1,x2,x3,w1,w2,w3; 658 int orderxy, orderz; 659 Matrix grad[ANCFPlateMaxIP]; 660 double jacdet[ANCFPlateMaxIP]; 661 662 Matrix Ti[4]; 663 Matrix T; //slope discontinuities 664 Matrix3D jac[4]; //slope discontinuities 665 Vector e0; //initial vector in e, not in p 666 667 double rho; //DR 2013-02-04 deleted rho from class element 668 }; 669 670 671 672 673 #endif 674 675 /* 676 B.SetSize(ns*dim,9); 677 B.SetAll(0); 678 Vector piola1v; 679 Matrix B; 680 */