1 //#**************************************************************
2 //# filename: ANCFBeamShear2D.cpp
3 //#
4 //# author: Karin & Astrid
5 //#
6 //# generated: 2009 / 2010
7 //# description: 2D ANCF beam element with shear deformation
8 //# comments: formulation is based on the following paper:
9 //# Gerstmayr, Matikainen, Mikkola: A geometrically exact beam element based on the ANCF
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 #include "rigid2d.h"
25 #include "femathhelperfunctions.h"
26 #include "material.h"
27 #include "ANCFBeamShear2D.h"
28 #include "Node.h"
29 //#include "graphicsconstants.h"
30 //#include "elementdataaccess.h"
31 //#include "solversettings_auto.h"
32
33 //SetANCFBeamShear-Fkt. for the linear element:
SetANCFBeamShear2D(const Vector & xc1,const Vector & xc2,int n1i,int n2i,int materialnumi,const Vector3D & si,const Vector3D & coli)34 void ANCFBeamShear2D::SetANCFBeamShear2D(const Vector& xc1, const Vector& xc2, int n1i, int n2i, int materialnumi,
35 const Vector3D& si, const Vector3D& coli)
36 {
37 nnodes=2;
38 n1=n1i; n2=n2i; sos2=0;
39 size = si;
40 n3=-1;
41
42 //x_init = xc1.Append(xc2);
43 //xg = xc1.Append(xc2);
44 xg.SetLen(SOS());
45 xg.SetAll(0.);
46 x_init.SetLen(SOS());
47 x_init.SetAll(0.);
48 x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!!
49
50 lx = size.X(); ly = size.Y(); lz = size.Z(); //lx is parameter to identify reference configuration
51 materialnum = materialnumi;
52 mass = lx*ly*lz*GetMaterial().Density(); //use lx, because this is the parameter of the real geometry ...
53 //lx=1;ly=1;lz=1;
54 col = coli;
55 BuildDSMatrices();
56
57 q0 = xc1.Append(xc2);
58 };
59
60 //SetANCFBeamShear-Fkt. for quadratic element:
SetANCFBeamShear2D(const Vector & xc1,const Vector & xc2,const Vector & xc3,int n1i,int n2i,int n3i,int materialnumi,const Vector3D & si,const Vector3D & coli)61 void ANCFBeamShear2D::SetANCFBeamShear2D(const Vector& xc1, const Vector& xc2, const Vector& xc3, int n1i, int n2i, int n3i, int materialnumi,
62 const Vector3D& si, const Vector3D& coli)
63 {
64 nnodes=3;
65 n1=n1i; n2=n2i; n3=n3i;
66 sos2=0;
67 size = si;
68
69 //x_init = xc1.Append(xc2);
70 //xg = xc1.Append(xc2);
71 xg.SetLen(SOS());
72 xg.SetAll(0.);
73 x_init.SetLen(SOS());
74 x_init.SetAll(0.);
75 x_init = x_init.Append(Vector(SOS())); //Velocity initial conditions can also be transformed by Tinv!!!
76
77 lx = size.X(); ly = size.Y(); lz = size.Z(); //lx is parameter to identify reference configuration
78
79 //size.X() = ComputeCurvedLength(x_init); //size.X() is true curve length, but other curve length is better for sliding joint!
80 //UO() << "linit=" << lx << ", ltrue=" << size.X() << "\n";
81 materialnum = materialnumi;
82
83 mass = lx*ly*lz*GetMaterial().Density(); //use lx, because this is the parameter of the real geometry ...
84 //lx=1;ly=1;lz=1;
85
86 col = coli;
87 //UO() << "Cable: rho=" << rho << ", E=" << Em << ", size=" << size << ", n1=" << n1 << ", n2=" << n2 << "\n";
88
89 BuildDSMatrices();
90
91 //q0 is generated from vector xc1(4),... :q0=xc1.Append(xc2), q0.Append(xc3)
92 q0.SetLen(SOS());
93 q0.SetAll(0.);
94 int i;
95
96 //should be realized with append. ...
97 for(i=1;i<=4;i++)
98 {
99 q0(i) = xc1(i);
100 }
101 for(i=5;i<=8;i++)
102 {
103 q0(i) = xc2(i-4);
104 }
105 for(i=9;i<=12;i++)
106 {
107 q0(i) = xc3(i-8);
108 }
109 //if (GetBeamEIy() == 0) // if Beam parameters were not initialized in Material..
110 //{
111 // GetMaterial().BeamEIy() = -1;
112 // GetMaterial().BeamEA() = 0;
113 // GetMaterial().BeamRhoA() = 0;
114 // GetMaterial().BeamA() = ly*lz;
115 //}
116 };
117
118
LinkToElements()119 void ANCFBeamShear2D::LinkToElements()
120 {
121 if (SOSowned() == 0)
122 {
123 //UO() << "Link nodes to elements in Cable\n";
124 const Node& node1 = GetMBS()->GetNode(n1);
125 const Node& node2 = GetMBS()->GetNode(n2);
126 for (int i=1; i <= node1.SOS(); i++)
127 {
128 AddLTG(node1.Get(i)); //LocalToGlobal-Element-list
129 }
130 for (int i=1; i <= node2.SOS(); i++)
131 {
132 AddLTG(node2.Get(i));
133 }
134 if(nnodes == 3)
135 {
136 const Node& node3 = GetMBS()->GetNode(n3);
137 for (int i=1; i <= node3.SOS(); i++)
138 {
139 AddLTG(node3.Get(i));
140 }
141 }
142 for (int i=1; i <= node1.SOS(); i++)
143 {
144 AddLTG(node1.Get(i+node1.SOS()));
145 }
146 for (int i=1; i <= node2.SOS(); i++)
147 {
148 AddLTG(node2.Get(i+node2.SOS()));
149 }
150 if(nnodes == 3)
151 {
152 const Node& node3 = GetMBS()->GetNode(n3);
153 for (int i=1; i <= node3.SOS(); i++)
154 {
155 AddLTG(node3.Get(i+node3.SOS()));
156 }
157 }
158 }
159 }
160
BuildDSMatrices()161 void ANCFBeamShear2D::BuildDSMatrices()
162 {
163 //orderx = 9; //max 9x5, otherwise array grad too small!!!!
164 //GetIntegrationRule(x1,w1,orderx);
165 };
166
167
168 //----------------
169 //Shapefunctions
170 //
171 //defined on scaled rectangular element: ploc=(\xi,\eta)
172 //-L/2<\xi<+L/2, -H/2<\eta<+H/2
173 //----------------
174
175 //GetShapes(sf,ploc) computes Shapefunctions at position ploc and saves values in vector sf
GetShapes(Vector & sf,const Vector2D & ploc) const176 void ANCFBeamShear2D::GetShapes(Vector& sf, const Vector2D& ploc) const //ploc=(x,y)
177 {
178 if(nnodes == 2)
179 {
180 sf.SetLen(NS()); //NS=2*nnodes, vector sf has 4 entries
181 sf(1)=1./lx*(lx/2.-ploc(1));
182 sf(2)=ploc(2)/lx*(lx/2.-ploc(1));
183 sf(3)=1./lx*(lx/2.+ploc(1));
184 sf(4)=ploc(2)/lx*(lx/2.+ploc(1));
185 }
186 else if(nnodes == 3)
187 {
188 sf.SetLen(NS()); //NS=2*nnodes, vector sf has 6 entries
189 sf(1)=-2./(lx*lx)*ploc(1)*(lx/2.-ploc(1));
190 sf(2)=-2./(lx*lx)*ploc(1)*ploc(2)*(lx/2.-ploc(1));
191 sf(3)=2./(lx*lx)*ploc(1)*(lx/2.+ploc(1));
192 sf(4)=2./(lx*lx)*ploc(1)*ploc(2)*(lx/2.+ploc(1));
193 sf(5)=-4./(lx*lx)*(ploc(1)-lx/2.)*(ploc(1)+lx/2.);
194 sf(6)=-4./(lx*lx)*ploc(2)*(ploc(1)-lx/2.)*(ploc(1)+lx/2.);
195 }
196 }
197
198 //GetSF(i,ploc) computes i-th Shapefunctions at ploc=(x,y)
GetSF(int i,const Vector2D & ploc) const199 double ANCFBeamShear2D::GetSF(int i, const Vector2D& ploc) const
200 {
201 if(nnodes == 2)
202 {
203 switch(i)
204 {
205 case 1: return 1./lx*(lx/2.-ploc(1));
206 case 2: return ploc(2)/lx*(lx/2.-ploc(1));
207 case 3: return 1./lx*(lx/2.+ploc(1));
208 case 4: return ploc(2)/lx*(lx/2.+ploc(1));
209 default: mbs->UO()<<"only 4 Shapefunctions\n"; return 0.; //in case i>4
210 }
211 }
212 else if(nnodes == 3)
213 {
214 switch(i)
215 {
216 case 1: return -2./(lx*lx)*ploc(1)*(lx/2.-ploc(1));
217 case 2: return -2./(lx*lx)*ploc(1)*ploc(2)*(lx/2.-ploc(1));
218 case 3: return 2./(lx*lx)*ploc(1)*(lx/2.+ploc(1));
219 case 4: return 2./(lx*lx)*ploc(1)*ploc(2)*(lx/2.+ploc(1));
220 case 5: return -4./(lx*lx)*(ploc(1)-lx/2.)*(ploc(1)+lx/2.);
221 case 6: return -4./(lx*lx)*ploc(2)*(ploc(1)-lx/2.)*(ploc(1)+lx/2.);
222 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.; //fin case i>6
223 }
224 }
225 return 0.;
226 }
227 //computes derivatives dS/dx at ploc
GetShapesx(Vector & sfx,const Vector2D & ploc) const228 void ANCFBeamShear2D::GetShapesx(Vector& sfx, const Vector2D& ploc) const
229 {
230 if(nnodes == 2)
231 {
232 sfx.SetLen(NS());
233 sfx(1)=-1./lx;
234 sfx(2)=-ploc(2)/lx;
235 sfx(3)=1./lx;
236 sfx(4)=ploc(2)/lx;
237 }
238 else if(nnodes == 3)
239 {
240 sfx.SetLen(NS());
241 sfx(1)=-1./lx+4.*ploc(1)/(lx*lx);
242 sfx(2)=-ploc(2)/lx+4.*ploc(1)*ploc(2)/(lx*lx);
243 sfx(3)=1./lx+4.*ploc(1)/(lx*lx);
244 sfx(4)=ploc(2)/lx+4.*ploc(1)*ploc(2)/(lx*lx);
245 sfx(5)=-8.*ploc(1)/(lx*lx);
246 sfx(6)=-8.*ploc(1)*ploc(2)/(lx*lx);
247 }
248 }
249
GetSFx(int i,const Vector2D & ploc) const250 double ANCFBeamShear2D::GetSFx(int i, const Vector2D& ploc) const
251 {
252 if(nnodes == 2)
253 {
254 switch(i)
255 {
256 case 1: return -1./lx;
257 case 2: return -ploc(2)/lx;
258 case 3: return 1./lx;
259 case 4: return ploc(2)/lx;
260 default: mbs->UO()<<"only 4 Shapefunctions\n"; return 0.;
261 }
262 }
263 else if(nnodes == 3)
264 {
265 switch(i)
266 {
267 case 1: return -1./lx+4.*ploc(1)/(lx*lx);
268 case 2: return -ploc(2)/lx+4.*ploc(1)*ploc(2)/(lx*lx);
269 case 3: return 1./lx+4.*ploc(1)/(lx*lx);
270 case 4: return ploc(2)/lx+4.*ploc(1)*ploc(2)/(lx*lx);
271 case 5: return -8.*ploc(1)/(lx*lx);
272 case 6: return -8.*ploc(1)*ploc(2)/(lx*lx);
273 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
274 }
275 }
276 return 0.;
277 }
278
279 //2nd derivatives d^2S/dx^2 are 0, GetShapesxx creates vector sfxx with zeros
GetShapesxx(Vector & sfxx,const Vector2D & ploc) const280 void ANCFBeamShear2D::GetShapesxx(Vector& sfxx, const Vector2D& ploc) const
281 {
282 if(nnodes == 2)
283 {
284 sfxx.SetLen(NS()); //NS=4
285 sfxx.SetAll(0); //all derivatives are 0
286 }
287 else if(nnodes == 3)
288 {
289 sfxx.SetLen(NS()); //NS=6
290 sfxx(1)=4./(lx*lx);
291 sfxx(2)=4.*ploc(2)/(lx*lx);
292 sfxx(3)=4./(lx*lx);
293 sfxx(4)=4.*ploc(2)/(lx*lx);
294 sfxx(5)=-8./(lx*lx);
295 sfxx(6)=-8.*ploc(2)/(lx*lx);
296 }
297 }
298
GetSFxx(int i,const Vector2D & ploc) const299 double ANCFBeamShear2D::GetSFxx(int i, const Vector2D& ploc) const
300 {
301 if(nnodes == 2)
302 {
303 return 0;
304 }
305 else if(nnodes == 3)
306 {
307 switch(i)
308 {
309 case 1: return 4./(lx*lx);
310 case 2: return 4.*ploc(2)/(lx*lx);
311 case 3: return 4./(lx*lx);
312 case 4: return 4.*ploc(2)/(lx*lx);
313 case 5: return -8./(lx*lx);
314 case 6: return -8.*ploc(2)/(lx*lx);
315 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
316 }
317 }
318 return 0.;
319 }
320
321 //computes derivatives dS/dy at ploc
GetShapesy(Vector & sfy,const Vector2D & ploc) const322 void ANCFBeamShear2D::GetShapesy(Vector& sfy, const Vector2D& ploc) const
323 {
324 if(nnodes == 2)
325 {
326 sfy.SetLen(NS());
327 sfy(1)=0;
328 sfy(2)=1./lx*(lx/2.-ploc(1));
329 sfy(3)=0;
330 sfy(4)=1./lx*(lx/2.+ploc(1));
331 }
332 else if(nnodes == 3)
333 {
334 sfy.SetLen(NS());
335 sfy(1)=0;
336 sfy(2)=-ploc(1)/lx+2*ploc(1)*ploc(1)/(lx*lx);
337 sfy(3)=0;
338 sfy(4)=ploc(1)/lx+2*ploc(1)*ploc(1)/(lx*lx);
339 sfy(5)=0;
340 sfy(6)=1.-4.*ploc(1)*ploc(1)/(lx*lx);
341 }
342 }
343
GetSFy(int i,const Vector2D & ploc) const344 double ANCFBeamShear2D::GetSFy(int i, const Vector2D& ploc) const
345 {
346 if(nnodes == 2)
347 {
348 switch(i)
349 {
350 case 1: return 0;
351 case 2: return 1./lx*(lx/2.-ploc(1));
352 case 3: return 0;
353 case 4: return 1./lx*(lx/2.+ploc(1));
354 default: mbs->UO()<<"only 4 Shapefunctions\n"; return 0.;
355 }
356 }
357 else if(nnodes == 3)
358 {
359 switch(i)
360 {
361 case 1: return 0;
362 case 2: return -ploc(1)/lx+2*ploc(1)*ploc(1)/(lx*lx);
363 case 3: return 0;
364 case 4: return ploc(1)/lx+2*ploc(1)*ploc(1)/(lx*lx);
365 case 5: return 0;
366 case 6: return 1.-4.*ploc(1)*ploc(1)/(lx*lx);
367 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
368 }
369 }
370 return 0.;
371 }
372
373 //computes mixed derivatives d^2S/dxdy at ploc
GetShapesxy(Vector & sfy,const Vector2D & ploc) const374 void ANCFBeamShear2D::GetShapesxy(Vector& sfy, const Vector2D& ploc) const
375 {
376 if(nnodes == 2)
377 {
378 sfy.SetLen(NS());
379 sfy(1)=0;
380 sfy(2)=-1./lx;
381 sfy(3)=0;
382 sfy(4)=1./lx;
383 }
384 else if(nnodes == 3)
385 {
386 sfy.SetLen(NS());
387 sfy(1)=0;
388 sfy(2)=-1./lx+4.*ploc(1)/(lx*lx);
389 sfy(3)=0;
390 sfy(4)=1./lx+4.*ploc(1)/(lx*lx);
391 sfy(5)=0;
392 sfy(6)=-8.*ploc(1)/(lx*lx);
393 }
394 }
395
GetSFxy(int i,const Vector2D & ploc) const396 double ANCFBeamShear2D::GetSFxy(int i, const Vector2D& ploc) const
397 {
398 if(nnodes == 2)
399 {
400 switch(i)
401 {
402 case 1: return 0;
403 case 2: return -1./lx;
404 case 3: return 0;
405 case 4: return 1./lx;
406 default: mbs->UO()<<"only 4 Shapefunctions\n"; return 0.;
407 }
408 }
409 else if(nnodes == 3)
410 {
411 switch(i)
412 {
413 case 1: return 0;
414 case 2: return -1./lx+4.*ploc(1)/(lx*lx);
415 case 3: return 0;
416 case 4: return 1./lx+4.*ploc(1)/(lx*lx);
417 case 5: return 0;
418 case 6: return -8.*ploc(1)/(lx*lx);
419 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
420 }
421 }
422 return 0.;
423 }
424
Gradu(const Vector2D & ploc,const Vector & u,Matrix3D & gradu) const425 void ANCFBeamShear2D::Gradu(const Vector2D& ploc, const Vector& u, Matrix3D& gradu) const
426 {
427
428 Matrix3D jac, jacinv;
429 GetJacobi(jac,ploc,q0);
430
431 jac.GetInverse(jacinv);
432 jacinv = jacinv.GetTp();
433
434 gradu.SetSize(2,2);
435 gradu.SetAll(0);
436
437 int dim = Dim();
438 int l;
439 for (int j = 1; j <= dim; j++)
440 {
441 for (int i = 1; i <= NS(); i++)
442 {
443 l = (i-1)*dim+j;
444 gradu(j,1) += GetSFx(i,ploc)*u(l);
445 gradu(j,2) += GetSFy(i,ploc)*u(l);
446 }
447 }
448 }
449
450 //Kappa
GetDeltaKappa(const double & x,const Vector & xg,Vector & dkappa,double & kappa) const451 void ANCFBeamShear2D::GetDeltaKappa(const double& x, const Vector& xg, Vector& dkappa, double& kappa) const
452 {
453
454 int dim = Dim();
455 int ns = NS();
456
457 static Vector SVx;
458 static Vector SVxx;
459 SVx.SetLen(NS());
460 SVxx.SetLen(NS());
461
462 Vector2D ploc(x,0.);
463 GetShapesx(SVx,ploc);
464 GetShapesxx(SVxx,ploc);
465
466 Vector2D rx(0.,0.);
467 rx = GetPosx2D(ploc);
468
469 Vector2D rxx(0.,0.);
470 rxx = GetPos2D(ploc);
471
472 double rxn = rx.Norm();
473
474 if (rxn == 0)
475 {
476 dkappa *= 0; kappa=0; return;
477 }
478 double rxcrxx = rx.Cross(rxx);
479 double f = rxcrxx;
480
481 double g = Sqr(rxn);
482
483 if (rxn == 0)
484 {
485 kappa=0; //f/0
486 }
487 else
488 {
489 kappa = f/g;
490 }
491
492 double g2inv = 1./Sqr(g);
493 double fn = f*g2inv;
494 fn *= 2; //==cable element kappa mode = 2
495
496 double gn = g*g2inv;
497
498 double t1;
499 double df, dg;
500
501 for (int i=1; i <= dim; i++)
502 {
503 for (int j=1; j <= ns; j++)
504 {
505 //dr,x/de x r,xx + r,x x dr,xx/de:
506 //delta r,x x d^2(r)/dx^2 + ... =>paper
507 switch (i) {
508 case 1: { t1 = SVx(j)*rxx.Y()-SVxx(j)*rx.Y(); break; }
509 case 2: { t1 = -SVx(j)*rxx.X()+SVxx(j)*rx.X(); break; }
510 default: ;
511 }
512
513 dg = (rx(i)*SVx(j)); //normed
514 df = t1;
515 dkappa((j-1)*dim+i) = df*gn-fn*dg; //normed
516 }
517 }
518 }
519
520 //----------------
521 //Velocity
522 //----------------
523
524 //flagD = 0 f�r Berechnungen, flagD = 1 for Visualization
525
GetVel2D(const Vector2D & p_loc,int flagD) const526 Vector2D ANCFBeamShear2D::GetVel2D(const Vector2D& p_loc, int flagD) const
527 {
528 Vector2D p(0.,0.);
529 for (int i = 1; i <= Dim(); i++) //i=1,2
530 {
531 for (int j = 1; j <= NS(); j++)//j=1,...,4
532 {
533 if(flagD==0)
534 p(i) += GetSF(j,p_loc)*XGP((j-1)*Dim()+i); //XGP=actual solution vector
535
536 else
537 p(i) += GetSF(j,p_loc)*XGPD((j-1)*Dim()+i);
538 }
539 }
540 return p;
541 };
542
543 //----------------
544 //Displacement
545 //----------------
546 //only for visualization
GetDisplacement2DD(const Vector2D & p_loc) const547 Vector2D ANCFBeamShear2D::GetDisplacement2DD(const Vector2D& p_loc) const
548 {
549 Vector2D p(0.,0.);
550 for (int i = 1; i <= Dim(); i++)
551 {
552 for (int j = 1; j <= NS(); j++)
553 {
554 p(i) += GetSF(j,p_loc)*(XGD((j-1)*Dim()+i));
555 }
556 }
557 return p;
558 };
559
560 //----------------
561 //GetPos
562 //----------------
563
564 //GetPos2D berechnet Vektor Positionsvektor r im deformierten Element
565 //GetPos2D computes r in deformed element
566 //xg = generalized coordinates (= our unknowns q)
GetPos2D(const Vector2D & p_loc,const Vector & xg) const567 Vector2D ANCFBeamShear2D::GetPos2D(const Vector2D& p_loc, const Vector& xg) const
568 {
569 Vector2D p(0.,0.);
570 for (int i = 1; i <= Dim(); i++)
571 {
572 for (int j = 1; j <= NS(); j++)
573 {
574 p(i) += GetSF(j,p_loc)*(xg((j-1)*Dim()+i)+q0((j-1)*Dim()+i)); //r=u+r_0=Sq+Sq_0 (q_0=initial values, x_init=0-Vector -> no displacement in the beginning))
575 }
576 }
577 return p;
578 };
579
580 //GetPos2D von oben, aber falls im Aufruf als zweite Eingabe kein Vektor, sondern ein Integer steht,
581 //wird diese Fkt f�r Grafik aufgerufen; globale Variable XG wird f�r Rechnung verwendet anstatt Eingabe xg
582
583 //same GetPos2D function as above, but with the second parameter flag D the function for the visualization is activated
GetPos2D(const Vector2D & p_loc,int flagD) const584 Vector2D ANCFBeamShear2D::GetPos2D(const Vector2D& p_loc, int flagD) const
585 {
586 Vector2D p(0.,0.);
587 for (int i = 1; i <= Dim(); i++)
588 {
589 for (int j = 1; j <= NS(); j++)
590 {
591 if(flagD==0)
592 p(i) += GetSF(j,p_loc)*(XG((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
593 else
594 p(i) += GetSF(j,p_loc)*(XGD((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
595 }
596 }
597 return p;
598 };
599
600 //Ableitungen von GetPos entsprechen unseren Richtungsvektoren
601 //GetPosx2D (Ableitung in Richtung der Balkenachse: dr/dx=r,x)
602 //GetPosx2D is the derivative with respect to the direction of the beam centerline
GetPosx2D(const Vector2D & p_loc,const Vector & xg) const603 Vector2D ANCFBeamShear2D::GetPosx2D(const Vector2D& p_loc, const Vector& xg) const
604 {
605 Vector2D p(0.,0.);
606 for (int i = 1; i <= Dim(); i++)
607 {
608 for (int j = 1; j <= NS(); j++)
609 {
610 p(i) += GetSFx(j,p_loc)*(xg((j-1)*Dim()+i)+q0((j-1)*Dim()+i)); //r,x=S,x*u+S,x*q_0
611 }
612 }
613 return p;
614 }
615
GetPosx2D(const Vector2D & p_loc,int flagD) const616 Vector2D ANCFBeamShear2D::GetPosx2D(const Vector2D& p_loc, int flagD) const
617 {
618 Vector2D p(0.,0.);
619 for (int i = 1; i <= Dim(); i++)
620 {
621 for (int j = 1; j <= NS(); j++)
622 {
623 if(flagD==0)
624 p(i) += GetSFx(j,p_loc)*(XG((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
625 else
626 p(i) += GetSFx(j,p_loc)*(XGD((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
627 }
628 }
629 return p;
630 }
631
632 //GetPosy2D (Ableitung in Richtung des Querschnitts: dr/dy=r,y)
633 //GetPosy2D is the derivative with respect to the direction of the cross section
GetPosy2D(const Vector2D & p_loc,const Vector & xg) const634 Vector2D ANCFBeamShear2D::GetPosy2D(const Vector2D& p_loc, const Vector& xg) const
635 {
636 Vector2D p(0.,0.);
637 for (int i = 1; i <= Dim(); i++)
638 {
639 for (int j = 1; j <= NS(); j++)
640 {
641 p(i) += GetSFy(j,p_loc)*(xg((j-1)*Dim()+i)+q0((j-1)*Dim()+i)); //r,y=S,y*u+S,y*q_0
642 }
643 }
644 return p;
645 }
646
GetPosy2D(const Vector2D & p_loc,int flagD) const647 Vector2D ANCFBeamShear2D::GetPosy2D(const Vector2D& p_loc, int flagD) const
648 {
649 Vector2D p(0.,0.);
650 for (int i = 1; i <= Dim(); i++)
651 {
652 for (int j = 1; j <= NS(); j++)
653 {
654 if(flagD==0)
655 p(i) += GetSFy(j,p_loc)*(XG((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
656 else
657 p(i) += GetSFy(j,p_loc)*(XGD((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
658 }
659 }
660 return p;
661 }
662
663 //not used at the moment:
664 //second derivatives of position vector r with respect to x
665
666 //Vector2D ANCFBeamShear2D::GetPosxx2D(const Vector2D& p_loc, const Vector& xg) const
667 // {
668 // Vector2D p(0.,0.);
669 // for (int i = 1; i <= Dim(); i++)
670 // {
671 // for (int j = 1; j <= NS(); j++)
672 // {
673 // p(i) += GetSFxx(j,p_loc)*(xg((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
674 // }
675 // }
676 // return p;
677 // }
678 //
679 //Vector2D ANCFBeamShear2D::GetPosxx2D(const Vector2D& p_loc, int flagD) const
680 // {
681 // Vector2D p(0.,0.);
682 // for (int i = 1; i <= Dim(); i++)
683 // {
684 // for (int j = 1; j <= NS(); j++)
685 // {
686 // if(flagD==0)
687 // p(i) += GetSFxx(j,p_loc)*(XG((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
688 // else
689 // p(i) += GetSFxx(j,p_loc)*(XGD((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
690 // }
691 // }
692 // return p;
693 // }
694
695
696 //GetPosxy2D: d2r/dydx=r,xy
697 //is used in DeltaThetax
GetPosxy2D(const Vector2D & p_loc,const Vector & xg) const698 Vector2D ANCFBeamShear2D::GetPosxy2D(const Vector2D& p_loc, const Vector& xg) const
699 {
700 Vector2D p(0.,0.);
701 for (int i = 1; i <= Dim(); i++)
702 {
703 for (int j = 1; j <= NS(); j++)
704 {
705 p(i) += GetSFxy(j,p_loc)*(xg((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
706 }
707 }
708 return p;
709 }
710
GetPosxy2D(const Vector2D & p_loc,int flagD) const711 Vector2D ANCFBeamShear2D::GetPosxy2D(const Vector2D& p_loc, int flagD) const
712 {
713 Vector2D p(0.,0.);
714 for (int i = 1; i <= Dim(); i++)
715 {
716 for (int j = 1; j <= NS(); j++)
717 {
718 if(flagD==0)
719 p(i) += GetSFxy(j,p_loc)*(XG((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
720 else
721 p(i) += GetSFxy(j,p_loc)*(XGD((j-1)*Dim()+i)+q0((j-1)*Dim()+i));
722 }
723 }
724 return p;
725 }
726
727
728 //computes position vector in the reference element: r_0
GetInitPos2D(const Vector2D & p_loc) const729 Vector2D ANCFBeamShear2D::GetInitPos2D(const Vector2D& p_loc) const
730 {
731 Vector2D p(0.,0.);
732 for (int i = 1; i <= Dim(); i++)
733 {
734 for (int j = 1; j <= NS(); j++)
735 {
736 p(i) += GetSF(j,p_loc)*(q0((j-1)*Dim()+i)); //r_0=Sq_0; (q_0 includes initial values)
737 }
738 }
739 return p;
740 };
741
742 //for visualization
GetPos2D0D(const Vector2D & p_loc) const743 Vector2D ANCFBeamShear2D::GetPos2D0D(const Vector2D& p_loc) const
744 {
745 Vector2D plocscaled;
746 plocscaled(1)=p_loc(1)*lx/2.;
747 plocscaled(2)=p_loc(2)*ly/2.;
748 return GetPos2D(plocscaled,1);
749 };
750
GetPos2D0D(const Vector2D & p_loc,double defscale) const751 Vector2D ANCFBeamShear2D::GetPos2D0D(const Vector2D& p_loc, double defscale) const
752 {
753 Vector2D plocscaled;
754 plocscaled(1)=p_loc(1)*lx/2.;
755 plocscaled(2)=p_loc(2)*ly/2.;
756 return GetInitPos2D(plocscaled) + defscale*GetDisplacement2DD(plocscaled);
757 };
758
759
760 //----------------
761 //directions T1, T2
762 //----------------
763 //T1 steht normal auf Querschnittsrichtung T2
764 //T1 is perpendicular to the direction of the cross section T2
GetT12D(const double & xbar) const765 Vector2D ANCFBeamShear2D::GetT12D(const double& xbar) const
766 {
767 Vector2D T1(0.,0.), T2(0.,0.);
768 Vector2D ploc(xbar,0.);
769 T2=GetPosy2D(ploc); //r,y
770 T2.Normalize(); //T1 has length 1
771 T1(1)=T2(2); //Rechts-Kippregel: vector T1 is perpendicular to T2
772 T1(2)=-T2(1);
773 return T1;
774 }
775
776 //direction of cross section: T2=r,y/Norm(r,y)
GetT22D(const double & xbar) const777 Vector2D ANCFBeamShear2D::GetT22D(const double& xbar) const
778 {
779 Vector2D T2(0.,0.);
780 Vector2D ploc(xbar,0.);
781 T2=GetPosy2D(ploc); //r,y
782 T2.Normalize(); //T2 has length 1
783 return T2;
784 }
785
786
787 //----------------
788 //Gamma, Theta'
789 //----------------
790 //Gamma1,Gamma2 (and corresponding Delta-Variations)
GetGamma12D(const double & xbar) const791 double ANCFBeamShear2D::GetGamma12D(const double& xbar) const
792 {
793 Vector2D ploc(xbar,0.);
794 return GetT12D(xbar)*GetPosx2D(ploc)-1; //Gamma1 = T1.r,x-1 (r,x=direction of beam axis, in Paper denoted by r_0
795 }
796
GetGamma22D(const double & xbar) const797 double ANCFBeamShear2D::GetGamma22D(const double& xbar) const
798 {
799 Vector2D ploc(xbar,0.);
800 return GetT22D(xbar)*GetPosx2D(ploc); //Gamma2 = T2.r,x
801 }
802
GetDeltaGamma22D(const double & xbar,Vector & DeltaGamma2) const803 void ANCFBeamShear2D::GetDeltaGamma22D(const double& xbar, Vector& DeltaGamma2) const
804 {
805 //Paper-equation (56)
806 //DeltaGamma2 = r,y.S,x/Norm + r,x.S,y/Norm - (r,x.r,y).r,y.S,y/Norm^3
807
808 DeltaGamma2.SetLen(SOS());
809 DeltaGamma2.SetAll(0);
810
811 Vector2D ry, rx, ploc(xbar,0.);
812 rx=GetPosx2D(ploc);
813 ry=GetPosy2D(ploc);
814
815 double ryNorm, ry3;
816 ryNorm=ry.Norm();
817 ry3=ryNorm*ryNorm*ryNorm;
818
819 if(ryNorm<1e-10)
820 return; //if denominator gets 0
821
822 for (int i = 1; i <= Dim(); i++) //i=1,2
823 {
824 for (int j = 1; j <= NS(); j++)//j=1,...,4
825 {
826 int k=(j-1)*Dim()+i;
827 DeltaGamma2(k)=ry(i)*GetSFx(j,ploc)/ryNorm + rx(i)*GetSFy(j,ploc)/ryNorm - (rx*ry)*ry(i)*GetSFy(j,ploc)/ry3;
828 }
829 }
830 }
831
GetDeltaGamma12D(const double & xbar,Vector & DeltaGamma1) const832 void ANCFBeamShear2D::GetDeltaGamma12D(const double& xbar, Vector& DeltaGamma1) const
833 {
834 //Paper-equation (57)
835 //DeltaGamma1 = \hat{r,y}.S,x/Norm - \hat{r,x}.S,y/Norm + (\hat{r,x}.r,y).r,y.S,y/Norm^3
836
837 DeltaGamma1.SetLen(SOS());
838 DeltaGamma1.SetAll(0);
839
840 Vector2D ry, rx, ploc(xbar,0.), rxhat, ryhat;
841 rx=GetPosx2D(ploc);
842 ry=GetPosy2D(ploc);
843 rxhat(1)=rx(2);
844 rxhat(2)=-rx(1);
845 ryhat(1)=ry(2);
846 ryhat(2)=-ry(1);
847
848 double ryNorm, ry3;
849 ryNorm=ry.Norm();
850 ry3=ryNorm*ryNorm*ryNorm;
851
852 if(ryNorm<1e-10)
853 return; //if denominator gets 0
854
855 for (int i = 1; i <= Dim(); i++) //i=1,2
856 {
857 for (int j = 1; j <= NS(); j++)//j=1,...,4
858 {
859 int k=(j-1)*Dim()+i;
860 DeltaGamma1(k)=ryhat(i)*GetSFx(j,ploc)/ryNorm - rxhat(i)*GetSFy(j,ploc)/ryNorm + (rxhat*ry)*ry(i)*GetSFy(j,ploc)/ry3;
861 }
862 }
863 }
864
865 //angle (for curvature): Theta'
GetThetax2D(const double & xbar) const866 double ANCFBeamShear2D::GetThetax2D(const double& xbar) const
867 {
868 //Paper-equation (50)
869 //Theta'= (r,y Cross r,xy)/Norm^2
870
871 Vector2D ploc(xbar,0.);
872 double f,g;
873 f=GetPosy2D(ploc).Cross(GetPosxy2D(ploc));
874 g=sqr(GetPosy2D(ploc).Norm());
875 if(fabs(g)<1e-10)
876 return 0; //if denominator gets 0
877 else
878 return f/g;
879 }
880
GetDeltaThetax2D(const double & xbar,Vector & DeltaThetax) const881 void ANCFBeamShear2D::GetDeltaThetax2D(const double& xbar, Vector& DeltaThetax) const
882 {
883 //Paper-equations (54) und (55)
884
885 DeltaThetax.SetLen(SOS());
886 DeltaThetax.SetAll(0);
887
888 double f,g;
889 Vector2D ploc(xbar,0);
890 f=GetPosy2D(ploc).Cross(GetPosxy2D(ploc));
891 g=sqr(GetPosy2D(ploc).Norm());
892 if(fabs(g)<1e-10)
893 return; //if denominator gets 0
894
895 Vector Deltaf, Deltag;
896 Deltaf.SetLen(SOS());
897 Deltaf.SetAll(0);
898 Deltag.SetLen(SOS());
899 Deltag.SetAll(0);
900
901 Vector2D ry, rxy;
902 ry=GetPosy2D(ploc);
903 rxy=GetPosxy2D(ploc);
904
905 //Deltaf(k)=(r,y Cross S,xy) - (r,xy Cross S,y)
906 //Deltag(k)=2*r,y*S,y
907
908 for (int i = 1; i <= Dim(); i++) //i=1,2
909 {
910 for (int j = 1; j <= NS(); j++)//j=1,...,4
911 {
912 int k=(j-1)*Dim()+i; //k = coordinate in Deltaf-vector
913 if(i==1)
914 Deltaf(k)=-ry(2)*GetSFxy(j,ploc)+rxy(2)*GetSFy(j,ploc);
915 else
916 Deltaf(k)=ry(1)*GetSFxy(j,ploc)-rxy(1)*GetSFy(j,ploc);
917
918 Deltag(k)=2.*ry(i)*GetSFy(j,ploc);
919 }
920 }
921
922 DeltaThetax=(g*Deltaf-f*Deltag);
923 DeltaThetax*=1./sqr(g);
924 }
925
926 //is used for \Pi^{thickness}:
927 //Eyy = transverse strain
GetEyy2D(const double & xbar) const928 double ANCFBeamShear2D::GetEyy2D(const double& xbar) const
929 {
930 Vector2D ploc(xbar,0);
931 Vector2D ry;
932 ry=GetPosy2D(ploc);
933 double Eyy;
934 Eyy = 0.5*(ry*ry)-0.5; //equation (51), Fehler: -1/2 statt -1, (bei Ableitung DeltaEyy wirkt sich dieser Faktor nicht aus)
935 return Eyy;
936 }
937
938 //DeltaEyy
GetDeltaEyy2D(const double & xbar,Vector & DeltaEyy) const939 void ANCFBeamShear2D::GetDeltaEyy2D(const double& xbar, Vector& DeltaEyy) const
940 {
941 Vector2D ploc(xbar,0);
942 Vector2D ry;
943 ry=GetPosy2D(ploc);
944
945 for (int i = 1; i <= Dim(); i++) //i=1,2
946 {
947 for (int j = 1; j <= NS(); j++)//j=1,...,4
948 {
949 int k=(j-1)*Dim()+i;
950 DeltaEyy(k)=ry(i)*GetSFy(j,ploc); //DeltaEyy = 0.5*2*r,y*Delta(r,y)=r,y.S,y
951 }
952 }
953 }
954
GetAngle2D(const Vector2D & ploc) const955 double ANCFBeamShear2D::GetAngle2D(const Vector2D& ploc) const
956 {
957 mbs->UO() << "called GetAngle2D";
958 //static Vector xg(SOS());
959 //xg.SetLen(SOS());
960
961 //for (int i = 1; i <= SOS(); i++)
962 // xg(i) = XG(i); //e=e1..e8
963
964 //Vector2D rx = GetPosx2D(ploc.X(),xg);
965 //double ang = atan2(rx.Y(),rx.X());
966
967 return 0.;//ang;
968 }
969 //
970 //double ANCFBeamShear2D::GetAngle2DP(const Vector2D& ploc) const
971 //{
972 // static Vector xg(SOS());
973 // xg.SetLen(SOS());
974 //
975 // for (int i = 1; i <= SOS(); i++)
976 // xg(i) = XG(i); //e=e1..e8
977 //
978 // double p0=ploc.X()/(0.5*lx);
979 // static Vector SVx;
980 //
981 // GetShapesx(SVx, p0);
982 // Vector2D rx(0.,0.);
983 // for (int i = 1; i <= Dim(); i++)
984 // {
985 // for (int j = 1; j <= NS(); j++)
986 // {
987 // rx(i) += SVx(j)*xg((j-1)*Dim()+i);
988 // }
989 // }
990 //
991 // double rxn = Sqr(rx.Norm());
992 // if (rxn == 0) rxn = 1;
993 //
994 // double angp = 0;
995 // for (int j = 1; j <= NS(); j++)
996 // {
997 // angp += XGP((j-1)*Dim()+1) * (-rx.Y()*SVx(j))/rxn;
998 // angp += XGP((j-1)*Dim()+2) * (SVx(j)*rx.X())/rxn;
999 // }
1000 //
1001 // return angp;
1002 //}
1003
GetdPosdqT(const Vector2D & ploc,Matrix & d)1004 void ANCFBeamShear2D::GetdPosdqT(const Vector2D& ploc, Matrix& d)
1005 {
1006 d.SetSize(NS()*Dim(),Dim());//12 x 2 (eigentlich: 6*(2x2))
1007 d.FillWithZeros();
1008
1009 for (int i = 1; i <= NS(); i++)
1010 {
1011 d((i-1)*Dim()+1,1)=GetSF(i,ploc);
1012 d((i-1)*Dim()+2,2)=GetSF(i,ploc);
1013 }
1014 }
1015
1016 //H = int(S,dV,V)
GetH(Matrix & H)1017 void ANCFBeamShear2D::GetH(Matrix& H)
1018 {
1019 if (Hmatrix.Getrows() == SOS())
1020 {
1021 H = Hmatrix;
1022 return;
1023 }
1024 else
1025 {
1026 int dim = Dim();
1027 int ns = NS();
1028
1029 H.SetSize(ns*dim,dim);
1030 H.SetAll(0);
1031 Matrix3D jac;
1032
1033 ConstVector<1> x2, w2;
1034 GetIntegrationRule(x1,w1,3); //3x1x1 !!!!!
1035 GetIntegrationRule(x2,w2,1);
1036
1037 for (int i1=1; i1<=x1.GetLen(); i1++)
1038 {
1039 for (int i2=1; i2<=x2.GetLen(); i2++)
1040 {
1041 Vector2D ploc(x1(i1)*lx*0.5,x2(i2)*ly*0.5);
1042 GetJacobi(jac,ploc,q0);
1043 double jacdet = jac.Det() *0.25*lx*ly;
1044 //mbs->UO() << "jacobian = " << jac << "\njacdet = " << jacdet << "\n";
1045 double fact = fabs (jacdet) * w1(i1)*w2(i2) * lz;
1046
1047 for (int i=0; i<ns; i++)
1048 {
1049 for (int j=1; j<=dim; j++)
1050 {
1051 H(i*dim+j,j)+=fact*GetSF(i+1, ploc);
1052 }
1053 }
1054 }
1055 }
1056 Hmatrix = H;
1057 }
1058 }
1059
1060
1061 //----------------
1062 //Massmatrix M
1063 //----------------
1064 //M = int(rho*((S)^T).S, dV,V)
EvalM(Matrix & m,double t)1065 void ANCFBeamShear2D::EvalM(Matrix& m, double t)
1066 {
1067 if (massmatrix.Getcols() == SOS())
1068 {
1069 m = massmatrix;
1070 return;
1071 }
1072 else
1073 {
1074 int dim = Dim();
1075 int ns = NS();
1076
1077 Vector x1,w1;
1078 int order=4;
1079
1080 GetIntegrationRule(x1,w1,order);
1081
1082 //we have a general cross section; instead of Em() we use Rho()
1083 double rhoI_w = Rho()*lz*Cub(ly)/12.;
1084 double A = ly*lz;
1085 double rhoA = Rho()*A;
1086
1087 double Shape_i, Shape_j;
1088 double Shapey_i, Shapey_j;
1089
1090 for (int i1=1; i1<=x1.GetLen(); i1++)//sum over integration points
1091 {
1092 double xi = x1(i1); //unit element
1093 double x = xi*lx*0.5; //scaled rect. reference element: \xi
1094 Vector2D ploc(x,0);
1095
1096 Vector2D rx0 = GetPosx2D(ploc);
1097 double rxn = rx0.Norm(); //curvature in reference element
1098 double det = 0.5*lx*rxn; //element transformation
1099
1100 for (int i_dim=1; i_dim<=dim; i_dim++)
1101 {
1102 for(int i_ns=1; i_ns<=ns; i_ns++)
1103 {
1104 int i=(i_ns-1)*dim+i_dim;
1105 Shape_i=GetSF(i_ns,ploc);
1106 Shapey_i=GetSFy(i_ns,ploc);
1107
1108 for (int j_dim=1; j_dim<=dim; j_dim++)
1109 {
1110 for(int j_ns=1; j_ns<=ns; j_ns++)
1111 {
1112 int j=(j_ns-1)*dim+j_dim;
1113 Shape_j=GetSF(j_ns,ploc);
1114 Shapey_j=GetSFy(j_ns,ploc);
1115 if(j_dim==i_dim)//i not equal j->Si*Sj=0
1116 {
1117 m(i,j)+=w1(i1)*det*rhoA*Shape_i*Shape_j+w1(i1)*det*rhoI_w*Shapey_i*Shapey_j;
1118 }
1119 }
1120 }
1121 }
1122 }
1123
1124 }
1125 massmatrix = m;
1126 //UO() << "m-invert=" << m.Invert2() << "\n";
1127 }
1128 };
1129 //last entry: number of rows, number of columns -> 0
1130
1131 //variational formulation: Delta\Pi (equation (52))
1132 //we calculate residual
EvalF2(Vector & f,double t)1133 void ANCFBeamShear2D::EvalF2(Vector& f, double t)
1134 {
1135 Body2D::EvalF2(f,t);
1136 TMStartTimer(22); //CPU timing
1137 int sos = SOS();
1138 ConstVector<12> fadd;//element residual
1139
1140 int ns = NS();
1141 int dim = Dim();
1142
1143 fadd.SetLen(SOS());
1144 fadd.SetAll(0);
1145
1146 ConstVector<12> temp;
1147 temp.SetLen(SOS());
1148 temp.SetAll(0);
1149
1150 SV.SetLen(NS());
1151
1152 //choose formulation:
1153 //Paper: A geom. exact beam element based on the ANCF (Gerstmayr, Matikainen, Mikkola)
1154 //int BeamFormulation = 0;//Finnland Paper: Simo&Vu-Quoc
1155 //int BeamFormulation = 1;//continuum mechanics based formulation (locking)
1156 //int BeamFormulation = 2;//continuum mechanics based formulation (locking compensation)
1157
1158 double EI_w = Em()*lz*Cub(ly)/12.; //Fl�chentr�gheitsmoment: moment of inertia of area
1159 double A = ly*lz; //Querschnittsfl�che: cross section area
1160 double EA = Em()*A; //E-Modul*area
1161 //Shear-Modul: G = E/(2+2nu)
1162 double GAks = A*ks*Em()/(2.+2.*Nu());
1163 //double GAks = A*Em()/(2.+2.*Nu());//GA (without ks)
1164 SetComputeCoordinates();
1165
1166 if(BeamFormulation == 0)//Finnland Paper
1167 {
1168
1169 //if (IsBeamParameters())
1170 //{
1171 // EI_w = GetBeamEIy();
1172 // EA = GetBeamEA();
1173 // GAks = GetMaterial().BeamGAky();
1174 //}
1175 //mbs->UO()<<"EI="<<EI_w<<", EA="<<EA<<", GAks="<<GAks<<"\n";
1176 //mbs->UO()<<"Nu="<<Nu()<<"\n";
1177
1178 //Integration order:
1179 //for Gau� Integration:
1180 int order_Gamma, order_Theta;
1181 //for Lobatto Integration:
1182 int order_Thickness;
1183
1184 //int_order
1185 if(nnodes == 2)
1186 {
1187 if(RI == 0)
1188 {
1189 order_Gamma = 4;
1190 order_Theta = 4;
1191 }
1192 else if(RI ==1)
1193 {
1194 order_Gamma = 1;
1195 order_Theta = 1;
1196 }
1197 order_Thickness =2;
1198
1199 }
1200 else if(nnodes == 3)
1201 {
1202 if(RI == 0)
1203 {
1204 order_Gamma = 6;
1205 order_Theta = 6;
1206 }
1207 else if(RI == 1)
1208 {
1209 order_Gamma = 3;
1210 order_Theta = 3;
1211 }
1212 order_Thickness = 3; //this gives order 4 in Lobatto
1213 }
1214 //Lobatto:
1215 // order=1 and order=2 -> both cases 2 Integration points on the borders
1216 // order=3 -> 3 integration points, 2 on the borders, 1 on the mid point
1217
1218 static Vector xG,xT,wG,wT,xE,wE; //xE,wE for Lobatto Integration
1219
1220 //Integration points and weights:
1221 GetIntegrationRule(xG,wG,order_Gamma);
1222 GetIntegrationRule(xT,wT,order_Theta);
1223 GetIntegrationRuleLobatto(xE,wE, order_Thickness);
1224 //mbs->UO()<<"Integrationspunkte Lobatto:"<<xE<<"\n";//How many Lobatto-points do we have?
1225 SetComputeCoordinates(); //values from global XG -> xg
1226
1227 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1228 //Integration
1229 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1230 static Vector x0;
1231 x0.SetLen(SOS());
1232 x0.SetAll(0.);
1233
1234 //Simo-Vu-Quoc-Strain:
1235
1236 //Theta*DeltaTheta (bending stiffness):
1237 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
1238 {
1239 double xi = xT(i1); //unit element
1240 double x = xi*lx*0.5; //scaled rect. reference element: xi
1241 Vector2D ploc(x,0);
1242
1243 Vector2D rx0 = GetPosx2D(ploc, x0);
1244 double rxn = rx0.Norm(); //Vorkr�mmung im Referenzelement: curvature in reference element
1245 double det = 0.5*lx*rxn; //element transformation
1246
1247 double fact_EI= EI_w*wT(i1)*det; //Det from element transformation
1248
1249 GetDeltaThetax2D(x,temp); //get DeltaTheta'-Werte and write in temp-vector
1250
1251 temp *= fact_EI*GetThetax2D(x); //temp=EI*Theta'*DeltaTheta'
1252 //mbs->UO()<<"temp="<<temp<<"\n";
1253 fadd += temp;
1254 }
1255
1256 //Gamma1*DeltaGamma1 (axial stiffness):
1257 for (int i1=1; i1<=xG.GetLen(); i1++)
1258 {
1259 double xi = xG(i1); //unit element
1260 double x = xi*lx*0.5; //scaled rect. reference element: xi
1261 Vector2D ploc(x,0);
1262
1263 Vector2D rx0 = GetPosx2D(ploc, x0);
1264 double rxn = rx0.Norm(); //Vorkr�mmung im Referenzelement: curvature in reference element
1265 double det = 0.5*lx*rxn; //element transformation
1266
1267 double fact_EA= EA*wG(i1)*det; //Det from element transformation
1268
1269 GetDeltaGamma12D(x,temp); //get DelatGamma1-Werte and write in temp-vector
1270
1271 temp *= fact_EA*GetGamma12D(x); //temp=EA*Gamma1*DeltaGamma1
1272 fadd += temp;
1273 }
1274
1275 //Gamma2*DeltaGamma2 (shear stiffness):
1276 for (int i1=1; i1<=xG.GetLen(); i1++)
1277 {
1278 double xi = xG(i1); //unit element
1279 double x = xi*lx*0.5; //scaled rect. reference element: xi
1280 Vector2D ploc(x,0);
1281
1282 Vector2D rx0 = GetPosx2D(ploc, x0);
1283 double rxn = rx0.Norm(); //Vorkr�mmung im Referenzelement: curvature in reference element
1284 double det = 0.5*lx*rxn; //element transformation
1285
1286 double fact_GAks= GAks*wG(i1)*det; //Det from element transformation
1287 //double fact_GAks= 1./(1./GAks+(lx*lx)/(12*EI_w))*wG(i1)*det; //Substitution from Felippa-document
1288
1289 GetDeltaGamma22D(x,temp); //get DeltaGamma2-Werte and write in temp-vector
1290
1291 temp *= fact_GAks*GetGamma22D(x); //temp=ks*GA*Gamma2*DeltaGamma2
1292 fadd += temp;//fadd=fadd+temp
1293 }
1294
1295 //thickness-Strain:
1296
1297 //Eyy*DeltaEyy:
1298 for (int i1=1; i1<=xE.GetLen(); i1++)
1299 {
1300 double xi = xE(i1); //unit element
1301 double x = xi*lx*0.5; //scaled rect. reference element: xi
1302 Vector2D ploc(x,0);
1303
1304 Vector2D rx0 = GetPosx2D(ploc, x0);
1305 double rxn = rx0.Norm();
1306 double det = 0.5*lx*rxn;
1307 //double update, Eyy, DeltaEyy;
1308 double fact_Eyy= EA*wE(i1)*det;
1309
1310 GetDeltaEyy2D(x,temp); //get DeltaEyy-Werte and write in temp-vector
1311
1312 temp *= fact_Eyy*GetEyy2D(x); //temp=EA*Eyy*DeltaEyy
1313 fadd += temp;
1314
1315 //linearized thickness-strain:
1316 //linearisiertes Eyy:
1317
1318 //if(i1==1)
1319 //{
1320 //Eyy=GetSFy(2,ploc)*XG(4); //XG = Vektor der FG
1321 //DeltaEyy=GetSFy(2,ploc);
1322 //update = fact_Eyy*Eyy*DeltaEyy; //temp=EA*Eyy(lin)*DeltaEyy(lin)
1323 //fadd(4) += update;
1324 //}
1325 //else
1326 //{
1327 //Eyy=GetSFy(4,ploc)*XG(8);
1328 //DeltaEyy=GetSFy(4,ploc);
1329 //update = fact_Eyy*Eyy*DeltaEyy; //temp=EA*Eyy(lin)*DeltaEyy(lin)
1330 //fadd(8) += update;
1331 //}
1332
1333 }
1334 }
1335 //= 1;//continuum mechanics based formulation (locking)
1336 //= 2;//continuum mechanics based formulation (locking compensation)
1337 else if(BeamFormulation == 1 || BeamFormulation == 2)
1338 {
1339
1340 //Integration Order
1341 int order_Poisson;
1342 int order_Thickness;
1343
1344 if(nnodes == 2)
1345 {
1346 if(RI == 0)
1347 {
1348 order_Poisson = 4;
1349 }
1350 else if(RI ==1)
1351 {
1352 order_Poisson = 1;
1353 }
1354 order_Thickness =2;
1355
1356 }
1357 else if(nnodes == 3)
1358 {
1359 if(RI == 0)
1360 {
1361 order_Poisson = 6;
1362 }
1363 else if(RI == 1)
1364 {
1365 order_Poisson = 3;
1366 }
1367 order_Thickness = 3;
1368 }
1369 //mbs->UO()<<"order="<<order_Poisson<<"\n";
1370 //mbs->UO()<<"order="<<order_Thickness<<"\n";
1371
1372 //UO() << "standard ANCForig\n";
1373
1374 Matrix3D Id;
1375 Id.Set22(1.,0.,0.,1.);
1376
1377 Matrix3D strain, piola1, F;
1378 F.SetSize(2,2);
1379
1380 temp.SetLen(SOS());
1381 fadd.SetLen(SOS());
1382 fadd.SetAll(0);
1383
1384 int poissoncorrection = 0; //1==reduced integration of poisson part
1385 if (BeamFormulation == 2) poissoncorrection = 1;
1386
1387 static Vector x1,x2,w1,w2;
1388
1389 for (int kk=1; kk <= 1+poissoncorrection; kk++)
1390 {
1391 Matrix3D Dm;
1392
1393 double ks = 10.*(1.+Nu())/(12.+11.*Nu()); //for rectangular cross-section only!
1394
1395 if (poissoncorrection)
1396 {
1397 if (BeamFormulation == 2)
1398 {
1399 Dm.SetAll(0);
1400 if (kk == 1)//D^0 (equ. (25))
1401 {
1402 //poisson ratio zero:
1403 Dm(1,1) = Em();
1404 Dm(2,2) = Em();
1405 Dm(3,3) = ks*Em()/(2.*(1.+Nu()));//=ks*G
1406
1407 //GetDMatrix(Dm,nu,Em);
1408 GetIntegrationRule(x1,w1,order_Poisson);//x_i ... integration points, w_i ... integration weights
1409 GetIntegrationRule(x2,w2,order_Thickness);
1410 }
1411 else if (kk == 2)//D^v (equ. (26))
1412 {
1413 Dm(1,1) = Em()/(1.-Sqr(Nu()))-Em();//=Em*nu^2/(1-nu^2)
1414 Dm(1,2) = Nu()*Em()/(1.-Sqr(Nu()));
1415 Dm(2,1) = Nu()*Em()/(1.-Sqr(Nu()));
1416 Dm(2,2) = Em()/(1.-Sqr(Nu()))-Em();
1417
1418 GetIntegrationRule(x1,w1,order_Poisson);//x_i ... integration points, w_i ... integration weights
1419 GetIntegrationRule(x2,w2,1);//only 1 int. point along cross section
1420 }
1421 }
1422 }
1423 else
1424 {
1425 GetIntegrationRule(x1,w1,order_Poisson);//x_i ... integration points, w_i ... integration weights
1426 GetIntegrationRule(x2,w2,order_Thickness);
1427 double nu2 = Nu();
1428
1429 double f = Em()/(1.-Sqr(nu2));
1430
1431 Dm(1,1)=f;
1432 Dm(1,2)=nu2*f;
1433 Dm(1,3)=0;
1434
1435 Dm(2,1)=nu2*f;
1436 Dm(2,2)=f;
1437 Dm(2,3)=0;
1438
1439 Dm(3,1)=0;
1440 Dm(3,2)=0;
1441 Dm(3,3)=.5*f*(1.-nu2);
1442 }
1443
1444 int kx1 = x2.Length();//number of integration points
1445
1446 for (int i1=1; i1<=x1.GetLen(); i1++)
1447 {
1448 for (int i2=1; i2<=x2.GetLen(); i2++)
1449 {
1450 //TMStartTimer(20);
1451 int i,j,k;
1452 int ind = (i1-1)*kx1+(i2-1);
1453 Vector2D p(x1(i1)*0.5*lx,x2(i2)*0.5*ly);
1454
1455 // compute F
1456 //F.SetAll(0);
1457 //F(1,1)=GetPosx2D(p).X();
1458 //F(2,1)=GetPosx2D(p).Y();
1459 //F(1,2)=GetPosy2D(p).X();
1460 //F(2,2)=GetPosy2D(p).Y();
1461
1462 Matrix3D jac;
1463 GetJacobi(jac,p,q0);
1464 double jacdet = jac.Det();
1465
1466 Matrix3D jacinv;
1467 jac.GetInverse(jacinv);
1468 jacinv = jacinv.GetTp();
1469
1470 static Matrix grad0, grad;//am Einheitselement
1471 grad0.SetSize(Dim(),NS());
1472 for(i=1; i<=NS(); i++)
1473 {
1474 grad0(1,i)=GetSFx(i,p);
1475 grad0(2,i)=GetSFy(i,p);
1476 }
1477 Mult(jacinv, grad0, grad);
1478 Gradu(p,xg,F);//Verschiebungsgrad
1479 F(1,1) += 1;//Positionsgrad
1480 F(2,2) += 1;
1481
1482 int linear=0;
1483 if (linear)
1484 {
1485 F -= Id;
1486 strain = 0.5*(F+F.GetTp());
1487 Vector3D s3(strain(1,1), strain(2,2), 2*strain(1,2));
1488 Vector3D stress = Dm*s3;
1489 piola1.Set22(stress.X(),stress.Z(),stress.Z(),stress.Y()); //linearized!
1490 }
1491 else if (BeamFormulation == 1 || BeamFormulation == 2)
1492 {
1493 //Green-Lagrange strain tensor
1494 //strain = 0.5 * (F.GetTp() * F - I);
1495
1496 //F.GetATA2(strain); //does not exist in 2D
1497 strain = 0.5 * (F.GetTp() * F);
1498 strain(1,1) -= 0.5; strain(2,2) -= 0.5;
1499
1500 //this formula does not work for plane stress/strain:
1501 //piola1 = F * ((2*mu) * strain + (la * strain.Trace())*Id );
1502 Vector3D s3(strain(1,1), strain(2,2), 2.*strain(1,2));
1503 Vector3D stress;
1504 stress = Dm * s3;
1505
1506 piola1.Set22(stress.X(),stress.Z(),stress.Z(),stress.Y()); //linearized!
1507
1508 piola1 = F * piola1;//2.PK
1509 }
1510
1511 for (int j=1; j <= dim; j++)
1512 {
1513 for (int i = 0; i < ns; i++)
1514 {
1515 //temp(dim*i+j) = grad[ind](1, i+1)*piola1(j,1) + grad[ind](2, i+1)*piola1(j,2);
1516 temp(dim*i+j) = grad(1, i+1)*piola1(j,1) + grad(2, i+1)*piola1(j,2);
1517 }
1518 }
1519
1520 //fadd.MultAdd(fabs (jacdet[ind]) * lz * w1(i1)*w2(i2),temp);
1521 fadd.MultAdd(fabs(jacdet)*0.5*lx*0.5*ly * lz * w1(i1)*w2(i2),temp);
1522 }
1523 }
1524 }
1525
1526 }
1527 //mbs->UO()<<"fadd="<<fadd<<"\n";
1528 //Residual:
1529 f -= fadd; //f=f-fadd, (Ku=f -> Residual=f-Ku, Ku=fadd)
1530
1531 //mbs->UO()<<"general.coord.: "<<XG(1)<<","<<XG(4)<<"\n";
1532
1533 TMStopTimer(22);
1534
1535 //if (GetMassDamping() != 0)
1536 //{
1537 // // +++++++ damping: +++++++
1538 // for (int i = 1; i <= SOS(); i++)
1539 // xg(i) = XGP(i);
1540
1541 // if (massmatrix.Getcols() == SOS())
1542 // {
1543 // Mult(massmatrix,xg,temp);
1544 // }
1545 // else
1546 // {
1547 // Matrix dmat;
1548 // dmat.SetSize(SOS(),SOS());
1549 // EvalM(dmat,t);
1550 // massmatrix = dmat;
1551 // Mult(dmat,xg,temp);
1552 // }
1553 // //double k=1;
1554 // //if (GetMBS()->GetTime() < 0.4) k=500;
1555 // temp *= GetMassDamping();
1556 // f -= temp;
1557 //}
1558
1559 };
1560
1561
1562
1563 //++++++++++++++++++++++++ Plasticity ++++++++++++++++++++++
1564
PostNewtonStep(double t)1565 double ANCFBeamShear2D::PostNewtonStep(double t) // changes plastic strains, returns yieldfunction(sigma)/Em
1566 {
1567 return 0;
1568 }
1569
PostprocessingStep()1570 void ANCFBeamShear2D::PostprocessingStep()
1571 {
1572 }
1573
1574
1575 //--------------------------------------------------------
1576 //for Displacement/Stress/Strain/Stress Resultants-Visualization:
1577 //---------------------------------------------------------
1578
GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables)1579 void ANCFBeamShear2D::GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables)
1580 {
1581 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_displacement,
1582 FieldVariableDescriptor::FVCI_z);
1583 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_stress,
1584 FieldVariableDescriptor::FVCI_z, true);
1585 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_stress_mises);
1586 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_total_strain,
1587 FieldVariableDescriptor::FVCI_y, true);
1588 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_force_axial);
1589 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_force_transversal);
1590 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_moment_bending);
1591 }
1592
GetFieldVariableValue(const FieldVariableDescriptor & fvd,const Vector2D & local_position,bool flagD)1593 double ANCFBeamShear2D::GetFieldVariableValue(const FieldVariableDescriptor & fvd, const Vector2D & local_position, bool flagD)
1594 {
1595 if(!flagD)
1596 return FIELD_VARIABLE_NO_VALUE;
1597
1598 //Displacements:
1599 //ploc_xi_0 is from -1 .. +1
1600 Vector2D ploc;
1601 ploc(1)=local_position(1)*lx*0.5; //ploc_xi_0 -> ploc_\xi
1602 ploc(2)=local_position(2)*ly*0.5; //ploc_eta_0 -> ploc_\eta
1603
1604 xgd.SetLen(SOS());
1605 GetDrawCoordinates(xgd); //xgd=XGD
1606
1607 if(fvd.VariableType() == FieldVariableDescriptor::FVT_displacement)
1608 return fvd.GetComponent(GetDisplacement2DD(ploc)); // displacements
1609
1610 //Stress and Strain:
1611
1612 double Thetax = GetThetax2D(ploc(1)); //sp�ter: event. hier Vorkr�mmung abziehen: -Vorkr�mmung;
1613 double Gamma1 = GetGamma12D(ploc(1));
1614 double Gamma2 = GetGamma22D(ploc(1));
1615
1616 Matrix3D strain(0);
1617 Vector2D rx, ry;
1618 rx=GetPosx2D(ploc, xgd);
1619 ry=GetPosy2D(ploc, xgd);
1620 strain(1,1) = 0.5*(rx(1)*rx(1)+rx(2)*rx(2)-1);
1621 strain(2,1) = 0.5*(rx(1)*ry(1)+rx(2)*ry(2));
1622 strain(1,2) = strain(2,1);
1623 strain(2,2) = 0.5*(ry(1)*ry(1)+ry(2)*ry(2)-1);
1624
1625 Matrix3D stress(0);
1626 double fact=Em()/(1.+Nu())/(1.-2.*Nu());
1627 stress(1,1)=fact*((1-Nu())*strain(1,1)+Nu()*strain(2,2));
1628 stress(2,2)=fact*(Nu()*strain(1,1)+(1-Nu())*strain(2,2));
1629 stress(1,2)=fact*(1-2*Nu())*strain(1,2);
1630 stress(2,1)=stress(1,2);
1631 stress(3,3)=fact*(Nu()*strain(1,1)+Nu()*strain(2,2));
1632
1633 double A = ly*lz;
1634 double EI_w = Em()*lz*Cub(ly)/12.;
1635 double EA = Em()*A;
1636 // for Gamma2:
1637 // G=E/(2+2nu)
1638 double GAks = ks*Em()/(2+2*Nu());
1639
1640 switch(fvd.VariableType())
1641 {
1642 case FieldVariableDescriptor::FVT_stress_mises: return stress.Mises();
1643 case FieldVariableDescriptor::FVT_stress: return fvd.GetComponent(stress);
1644 case FieldVariableDescriptor::FVT_total_strain: return fvd.GetComponent(strain);
1645 case FieldVariableDescriptor::FVT_beam_force_axial: return EA * Gamma1;
1646 case FieldVariableDescriptor::FVT_beam_force_transversal: return GAks * Gamma2;
1647 case FieldVariableDescriptor::FVT_beam_moment_bending: return EI_w * Thetax;
1648 }
1649
1650 return FIELD_VARIABLE_NO_VALUE;
1651 }
1652
1653
1654 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1655 //++++++++++++++++++++++++++++++++�nderungen bis hier++++++++++++++++++++++++++++++++++++++++
1656 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1657
DrawElement()1658 void ANCFBeamShear2D::DrawElement()
1659 {
1660 mbs->SetColor(col);
1661
1662 // draw elements at z = 0.1
1663 Vector3D offset(0,0,0.);
1664
1665 double lx1 = lx; double ly1 = ly*GetMBS()->GetMagnifyYZ();
1666
1667 int linemode = 1; //0=no lines, 1=outline+color, 2=outline, 3=elementline+color, 4=elementline
1668 if (GetMBS()->GetIOption(110) && !GetMBS()->GetIOption(111))
1669 {
1670 linemode = 2;
1671 }
1672 if (!GetMBS()->GetIOption(110) && GetMBS()->GetIOption(111))
1673 {
1674 linemode = 0;
1675 }
1676 if (!GetMBS()->GetIOption(110) && !GetMBS()->GetIOption(111))
1677 {
1678 linemode = 4;
1679 }
1680
1681 int colormode = 0;
1682 if (GetMBS()->GetActualPostProcessingFieldVariable() != NULL)
1683 colormode = 1;
1684
1685 double def_scale = GetMBS()->GetDOption(105); // deformation scaling
1686
1687 int drawlinemode = GetMBS()->GetIOption(117); // Draw Flat Elements, in this case, function is plotted normal to element axis
1688
1689 if ( (1||colormode) && !drawlinemode) // plot coloured elements
1690 {
1691 lx1 /= lx;
1692 ly1 /= ly;
1693
1694 double modeval = 0;
1695 int xgset = 0;
1696
1697 double tilex = GetMBS()->GetIOption(137);
1698 double tiley = GetMBS()->GetIOption(138);
1699
1700 TArray<Vector3D> points((int)(tilex+1)*(int)(tiley+1));
1701 TArray<double> vals((int)(tilex+1)*(int)(tiley+1));
1702 double v=0;
1703
1704 points.SetLen(0); vals.SetLen(0);
1705 Vector2D p0, p0_mag, vx, vy, vy_mag;
1706 int tileyn = (int)tiley;
1707 int tilexn = (int)tilex;
1708
1709 p0 = Vector2D(-lx1,-1.);
1710 vx = Vector2D(2.*lx1/tilexn,0);
1711 vy = Vector2D(0,2./tileyn);
1712
1713 p0_mag = Vector2D(-lx1,-ly1);
1714 vy_mag = Vector2D(0,2.*ly1/tileyn);
1715
1716 for (double iy = 0; iy <= tileyn+1e-10; iy++)
1717 {
1718 for (double ix = 0; ix <= tilexn+1e-10; ix++)
1719 {
1720 Vector2D ploc = (p0+ix*vx+iy*vy);
1721 Vector2D ploc_mag = (p0_mag+ix*vx+iy*vy_mag);
1722 Vector2D pg = GetPos2D0D(ploc_mag, def_scale);
1723
1724 Vector3D p(ToP3D(pg)+offset);
1725 points.Add(p);
1726 if (colormode)
1727 v = GetFieldVariableValue(*GetMBS()->GetActualPostProcessingFieldVariable(), ploc, true);
1728 vals.Add(v);
1729 }
1730 }
1731 mbs->DrawColorQuads(points,vals,(int)tilexn+1,(int)tileyn+1,colormode,linemode);
1732 }
1733 else if (drawlinemode) // plot function on center line over beam axis
1734 {
1735 double scalex = lx1/lx;
1736 double scaley = ly1/ly;
1737
1738 int tilex = GetMBS()->GetIOption(137);
1739
1740 TArray<Vector3D> points(tilex+1), normals(tilex+1), drawpoints(tilex+1);
1741 TArray<double> vals(tilex+1);
1742 double v=0;
1743
1744 points.SetLen(0); vals.SetLen(0); normals.SetLen(0); drawpoints.SetLen(0);
1745 Vector2D p0, vx, vy;
1746
1747 p0 = Vector2D(-scalex,0);
1748 vx = Vector2D(2.*scalex/tilex,0);
1749 vy = Vector2D(0.,0.);
1750
1751 int cnt = 0;
1752 for (double ix = 0.; ix <= tilex+1e-10; ix++)
1753 {
1754 // point on the middle line
1755 Vector2D ploc = (p0+ix*vx);
1756 Vector2D pg = GetPos2D0D(ploc, def_scale);
1757
1758 Vector3D p(ToP3D(pg)+offset);
1759 points.Add(p);
1760
1761 // normal direction (numerical differentiation for the moment)
1762 Vector2D ploc_eps = ploc + 1e-5*vx;
1763 Vector2D pg_eps = GetPos2D0D(ploc_eps, def_scale);
1764
1765 Vector3D normal = Vector3D(-pg_eps.Y()+pg.Y(), pg_eps.X()-pg.X(), 0.);
1766 double scale = ly1/normal.Norm();
1767 normal *= scale;
1768 normals.Add(normal);
1769
1770 // computation of value
1771 if (colormode)
1772 v = GetFieldVariableValue(*GetMBS()->GetActualPostProcessingFieldVariable(), ploc, true);
1773 vals.Add(v);
1774
1775 MBS* constmbs = const_cast<MBS*>(mbs);
1776 constmbs->UpdateFEMinMaxCol(v);
1777 }
1778
1779 for (int cnt=1; cnt<=tilex+1; cnt++)
1780 {
1781 // scale normal according to minimum/maximum
1782 double scale = Maximum(fabs(GetMBS()->GetTImincol()), fabs(GetMBS()->GetTImaxcol()));
1783 if (scale==0) scale = 1.;
1784 // point on the function line
1785 drawpoints.Add(points(cnt)+(vals(cnt)/scale)*normals(cnt));
1786
1787 // draw arrow for value of function, in normal direction to middle line, colored matching legend
1788 double diff = mbs->GetFEmaxcol()-mbs->GetFEmincol();
1789 Vector3D fecolor = mbs->FEColor( (vals(cnt)-mbs->GetFEmincol())/diff );
1790 mbs->MyDrawArrow(points(cnt), drawpoints(cnt), fecolor);
1791 if (cnt>1)
1792 {
1793 // red line for plotting function value
1794 mbs->MyDrawLine(drawpoints(cnt-1), drawpoints(cnt), ly1*1e2, col);
1795 // middle line of element
1796 mbs->MyDrawLine(points(cnt-1), points(cnt), ly1*5e1);
1797 // upper and lower line of element, unscaled
1798 mbs->MyDrawLine(points(cnt-1)+(0.5/scaley)*normals(cnt-1), points(cnt)+(0.5/scaley)*normals(cnt), ly1*5e1);
1799 mbs->MyDrawLine(points(cnt-1)-(0.5/scaley)*normals(cnt-1), points(cnt)-(0.5/scaley)*normals(cnt), ly1*5e1);
1800 }
1801
1802 }
1803
1804 TArray<Vector3D> nodepoints(4);
1805 nodepoints(1) = ToP3D(GetPos2D0D(Vector2D(-scalex, 1.), def_scale)) + offset;
1806 nodepoints(2) = ToP3D(GetPos2D0D(Vector2D(-scalex,-1.), def_scale)) + offset;
1807 nodepoints(3) = ToP3D(GetPos2D0D(Vector2D( scalex,-1.), def_scale)) + offset;
1808 nodepoints(4) = ToP3D(GetPos2D0D(Vector2D( scalex, 1.), def_scale)) + offset;
1809 //mbs->DrawPolygon(nodepoints, 1, ly*ly1*1e-1);
1810 // left and right line of element
1811 mbs->MyDrawLine(nodepoints(1), nodepoints(2), ly1*4e1);
1812 mbs->MyDrawLine(nodepoints(3), nodepoints(4), ly1*4e1);
1813 }
1814 else
1815 {
1816 double tiling = 24;
1817
1818 Vector2D p1,p2,p3,p4;
1819 mbs->SetColor(col);
1820
1821 for (double i = 0; i < tiling; i++)
1822 {
1823 double l1 = -lx1*0.5+lx1*i/tiling;
1824 double l2 = -lx1*0.5+lx1*(i+1)/tiling;
1825 if (i == 0)
1826 {
1827 p1 = GetPos2D0D(Vector2D(l1/(0.5*lx), 1.), def_scale);
1828 p2 = GetPos2D0D(Vector2D(l1/(0.5*lx),-1.), def_scale);
1829
1830 }
1831 else
1832 {
1833 p1 = p4;
1834 p2 = p3;
1835 }
1836 p3 = GetPos2DD(Vector2D(l2,-ly*0.5));
1837 p4 = GetPos2DD(Vector2D(l2, ly*0.5));
1838 mbs->DrawQuad(ToP3D(p4),ToP3D(p3),ToP3D(p2),ToP3D(p1));
1839 }
1840 }
1841 };
1842
1843
1844